Freezing of a lake
In this example we simulate the freezing of a lake in winter. The lake is represented by four points that start at 1ᵒC and are cooled down by an atmosphere with temperatures of -20, -10, -1, and -0.1ᵒC. The lake is 10 m deep and not subjected to radiative transfer (this lake is covered in a wind tunnel under which we blow some cold air).
We start by using Oceananigans
to bring in functions for building grids and Simulation
s and the like.
using Oceananigans
using Oceananigans.Units
using ClimaSeaIce
using CairoMakie
Generate a 1D grid for difference ice columns subject to different solar insolation
grid = RectilinearGrid(size=4, x=(0, 1), topology=(Periodic, Flat, Flat))
4×1×1 RectilinearGrid{Float64, Oceananigans.Grids.Periodic, Oceananigans.Grids.Flat, Oceananigans.Grids.Flat} on Oceananigans.Architectures.CPU with 3×0×0 halo
├── Periodic x ∈ [0.0, 1.0) regularly spaced with Δx=0.25
├── Flat y
└── Flat z
The sensible heat flux from the atmosphere is represented by a FluxFunction
.
parameters = (
transfer_coefficient = 1e-3, # Unitless
atmosphere_density = 1.225, # kg m⁻³
atmosphere_heat_capacity = 1004, #
atmosphere_temperature = [-20, -10, -5, 0], # ᵒC
atmosphere_wind_speed = 5 # m s⁻¹
)
(transfer_coefficient = 0.001, atmosphere_density = 1.225, atmosphere_heat_capacity = 1004, atmosphere_temperature = [-20, -10, -5, 0], atmosphere_wind_speed = 5)
Flux is positive (cooling by fluxing heat up away from upper surface) when Tₐ < Tᵤ:
@inline function sensible_heat_flux(i, j, grid, Tᵤ, clock, fields, parameters)
Cₛ = parameters.transfer_coefficient
ρₐ = parameters.atmosphere_density
cₐ = parameters.atmosphere_heat_capacity
Tₐ = parameters.atmosphere_temperature[i]
uₐ = parameters.atmosphere_wind_speed
return Cₛ * ρₐ * cₐ * uₐ * (Tᵤ - Tₐ)
end
sensible_heat_flux (generic function with 1 method)
We also evolve a bucket freshwater lake that cools down and freezes from below generating fresh sea-ice (or lake-ice in this case?). We set the Δt of the lake to 10 minutes. This time step will be used to also for the sea-ice model.
lake = (
lake_density = 1000, # kg m⁻³
lake_heat_capacity = 4000, #
lake_temperature = [1.0, 1.0, 1.0, 1.0], # ᵒC
lake_depth = 10, # m
atmosphere = parameters,
Δt = 10minutes
)
(lake_density = 1000, lake_heat_capacity = 4000, lake_temperature = [1.0, 1.0, 1.0, 1.0], lake_depth = 10, atmosphere = (transfer_coefficient = 0.001, atmosphere_density = 1.225, atmosphere_heat_capacity = 1004, atmosphere_temperature = [-20, -10, -5, 0], atmosphere_wind_speed = 5), Δt = 600.0)
We build a flux function that serves three purposes:
- computing the cooling of the lake caused by the atmosphere
- If the lake temperature is low enough, freezing the lake from above
- and adding the heat flux to the bottom of the ice
@inline function advance_lake_and_frazil_flux(i, j, grid, Tuᵢ, clock, fields, parameters)
atmos = parameters.atmosphere
lake = parameters
Cₛ = atmos.transfer_coefficient
ρₐ = atmos.atmosphere_density
cₐ = atmos.atmosphere_heat_capacity
Tₐ = atmos.atmosphere_temperature[i]
uₐ = atmos.atmosphere_wind_speed
Tₒ = lake.lake_temperature
cₒ = lake.lake_heat_capacity
ρₒ = lake.lake_density
Δ = lake.lake_depth
ℵ = fields.ℵ[i, j, 1]
Δt = lake.Δt
Qₐ = Cₛ * ρₐ * cₐ * uₐ * (Tₐ - Tₒ[i]) * (1 - ℵ)
Tₒ[i] = Tₒ[i] + Qₐ / (ρₒ * cₒ) * Δt
Qᵢ = ρₒ * cₒ * (Tₒ[i] - 0) / Δt * Δ # W m⁻²
Qᵢ = min(Qᵢ, zero(Qᵢ)) # We only freeze, not melt
Tₒ[i] = ifelse(Qᵢ == 0, Tₒ[i], zero(Qᵢ))
return Qᵢ
end
aerodynamic_flux = FluxFunction(sensible_heat_flux; parameters)
top_heat_flux = (aerodynamic_flux)
bottom_heat_flux = FluxFunction(advance_lake_and_frazil_flux; parameters=lake)
model = SeaIceModel(grid;
ice_consolidation_thickness = 0.05, # m
top_heat_flux,
bottom_heat_flux)
SeaIceModel{Oceananigans.Architectures.CPU, Oceananigans.Grids.RectilinearGrid}(time = 0 seconds, iteration = 0)
├── grid: 4×1×1 RectilinearGrid{Float64, Oceananigans.Grids.Periodic, Oceananigans.Grids.Flat, Oceananigans.Grids.Flat} on Oceananigans.Architectures.CPU with 3×0×0 halo
├── ice_thermodynamics: SlabThermodynamics
├── advection: Nothing
└── external_heat_fluxes:
├── top: FluxFunction of sensible_heat_flux with parameters (transfer_coefficient=0.001, atmosphere_density=1.225, atmosphere_heat_capacity=1004, atmosphere_temperature=4-element Vector{Int64}, atmosphere_wind_speed=5)
└── bottom: FluxFunction of advance_lake_and_frazil_flux with parameters (lake_density=1000, lake_heat_capacity=4000, lake_temperature=4-element Vector{Float64}, lake_depth=10, atmosphere=(transfer_coefficient=0.001, atmosphere_density=1.225, atmosphere_heat_capacity=1004, atmosphere_temperature=4-element Vector{Int64}, atmosphere_wind_speed=5), Δt=600.0)
We initialize all the columns with open water (0 thickness and 0 concentration)
set!(model, h=0, ℵ=0)
simulation = Simulation(model, Δt=lake.Δt, stop_time=10days)
Simulation of SeaIceModel
├── Next time step: 10 minutes
├── Elapsed wall time: 0 seconds
├── Wall time per iteration: NaN days
├── Stop time: 10 days
├── Stop iteration: Inf
├── Wall time limit: Inf
├── Minimum relative step: 0.0
├── Callbacks: OrderedDict with 3 entries:
│ ├── stop_time_exceeded => Callback of stop_time_exceeded on IterationInterval(1)
│ ├── stop_iteration_exceeded => Callback of stop_iteration_exceeded on IterationInterval(1)
│ └── wall_time_limit_exceeded => Callback of wall_time_limit_exceeded on IterationInterval(1)
├── Output writers: OrderedDict with no entries
└── Diagnostics: OrderedDict with no entries
The data is accumulated in a timeseries for visualization.
timeseries = []
function accumulate_timeseries(sim)
T = model.ice_thermodynamics.top_surface_temperature
h = model.ice_thickness
ℵ = model.ice_concentration
To = lake.lake_temperature
push!(timeseries, (time(sim),
h[1, 1, 1], ℵ[1, 1, 1], T[1, 1, 1],
h[2, 1, 1], ℵ[2, 1, 1], T[2, 1, 1],
h[3, 1, 1], ℵ[3, 1, 1], T[3, 1, 1],
h[4, 1, 1], ℵ[4, 1, 1], T[4, 1, 1],
To[1], To[2], To[3], To[4]))
end
simulation.callbacks[:save] = Callback(accumulate_timeseries)
run!(simulation)
[ Info: Initializing simulation...
[ Info: ... simulation initialization complete (3.460 seconds)
[ Info: Executing initial time step...
[ Info: ... initial time step complete (1.450 seconds).
[ Info: Simulation is stopping after running for 0 seconds.
[ Info: Simulation time 10 days equals or exceeds stop time 10 days.
Extract and visualize data
t = [datum[1] for datum in timeseries]
h1 = [datum[2] for datum in timeseries]
ℵ1 = [datum[3] for datum in timeseries]
T1 = [datum[4] for datum in timeseries]
h2 = [datum[5] for datum in timeseries]
ℵ2 = [datum[6] for datum in timeseries]
T2 = [datum[7] for datum in timeseries]
h3 = [datum[8] for datum in timeseries]
ℵ3 = [datum[9] for datum in timeseries]
T3 = [datum[10] for datum in timeseries]
h4 = [datum[11] for datum in timeseries]
ℵ4 = [datum[12] for datum in timeseries]
T4 = [datum[13] for datum in timeseries]
L1 = [datum[14] for datum in timeseries]
L2 = [datum[15] for datum in timeseries]
L3 = [datum[16] for datum in timeseries]
L4 = [datum[17] for datum in timeseries]
set_theme!(Theme(fontsize=18, linewidth=3))
fig = Figure(size=(1000, 900))
axT = Axis(fig[1, 1], xlabel="Time (days)", ylabel="Top temperature (ᵒC)")
axℵ = Axis(fig[2, 1], xlabel="Time (days)", ylabel="Ice concentration (-)")
axh = Axis(fig[3, 1], xlabel="Time (days)", ylabel="Ice thickness (m)")
axL = Axis(fig[4, 1], xlabel="Time (days)", ylabel="Lake temperature (ᵒC)")
lines!(axT, t / day, T1)
lines!(axℵ, t / day, ℵ1)
lines!(axh, t / day, h1)
lines!(axL, t / day, L1)
lines!(axT, t / day, T2)
lines!(axℵ, t / day, ℵ2)
lines!(axh, t / day, h2)
lines!(axL, t / day, L2)
lines!(axT, t / day, T3)
lines!(axℵ, t / day, ℵ3)
lines!(axh, t / day, h3)
lines!(axL, t / day, L3)
lines!(axT, t / day, T4)
lines!(axℵ, t / day, ℵ4)
lines!(axh, t / day, h4)
lines!(axL, t / day, L4)
save("freezing_in_winter.png", fig)
This page was generated using Literate.jl.