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 Simulations and the like.
using Oceananigans
using Oceananigans.Units
using ClimaSeaIce
using CairoMakieGenerate 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ₐ)
endsensible_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 entriesThe 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.615 seconds)
[ Info: Executing initial time step...
[ Info: ... initial time step complete (1.472 seconds).
[ Info: Simulation is stopping after running for 5.374 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.