A freezing bucket
A common laboratory experiment freezes an insultated bucket of water from the top down, using a metal lid to keep the top of the bucket at some constant, very cold temperature. In this example, we simulate such a scenario using the SeaIceModel. Here, the bucket is perfectly insulated and infinitely deep, like many buckets are: if the Simulation is run for longer, the ice will keep freezing, and freezing, and will never run out of water. Also, the water in the infinite bucket is (somehow) all at the same temperature, in equilibrium with the ice-water interface (and therefore fixed at the melting temperature). Yep, this kind of thing happens all the time.
We start by using Oceananigans to bring in functions for building grids and Simulations and the like.
using Oceananigans
using Oceananigans.UnitsNext we using ClimaSeaIce to get some ice-specific names.
using ClimaSeaIceAn infinitely deep bucket with a single grid point
Perhaps surprisingly, we need just one grid point to model an possibly infinitely thick slab of ice with SlabSeaIceModel. We would only need more than 1 grid point if our boundary conditions vary in the horizontal direction.
grid = RectilinearGrid(size=(), topology=(Flat, Flat, Flat))1×1×1 RectilinearGrid{Float64, Oceananigans.Grids.Flat, Oceananigans.Grids.Flat, Oceananigans.Grids.Flat} on Oceananigans.Architectures.CPU with 0×0×0 halo
├── Flat x
├── Flat y
└── Flat z Next, we build our model of an ice slab freezing into a bucket. We start by defining a constant internal ConductiveFlux with ice_conductivity
conductivity = 2 # kg m s⁻³ K⁻¹
internal_heat_flux = ConductiveFlux(; conductivity)ConductiveFlux{Float64}(2.0)Note that other units besides Celsius can be used, but that requires setting model.phase_transitions` with appropriate parameters. We set the ice heat capacity and density as well,
ice_heat_capacity = 2100 # J kg⁻¹ K⁻¹
ice_density = 900 # kg m⁻³
phase_transitions = PhaseTransitions(; ice_heat_capacity, ice_density)PhaseTransitions{Float64, ClimaSeaIce.SeaIceThermodynamics.LinearLiquidus{Float64}}(900.0, 2100.0, 999.8, 4186.0, 334000.0, 0.0, ClimaSeaIce.SeaIceThermodynamics.LinearLiquidus{Float64}(0.0, 0.054))We set the top ice temperature,
top_temperature = -10 # ᵒC
top_heat_boundary_condition = PrescribedTemperature(-10)PrescribedTemperature{Int64}(-10)Construct the icethermodynamics of sea ice, for this we use a simple slab sea ice representation of icethermodynamics
ice_thermodynamics = SlabSeaIceThermodynamics(grid;
internal_heat_flux,
phase_transitions,
top_heat_boundary_condition)SlabSeaIceThermodynamics
└── top_surface_temperature: ConstantField(-10)We also prescribe a frazil ice heat flux that stops when the ice has reached a concentration of 1. This heat flux represents the initial ice formation from a liquid bucket.
@inline frazil_ice_formation(i, j, grid, Tuᵢ, clock, fields) = - (1 - fields.ℵ[i, j, 1]) # W m⁻²
bottom_heat_flux = FluxFunction(frazil_ice_formation)FluxFunction of frazil_ice_formationThen we assemble it all into a model.
model = SeaIceModel(grid; ice_thermodynamics, bottom_heat_flux)SeaIceModel{Oceananigans.Architectures.CPU, Oceananigans.Grids.RectilinearGrid}(time = 0 seconds, iteration = 0)
├── grid: 1×1×1 RectilinearGrid{Float64, Oceananigans.Grids.Flat, Oceananigans.Grids.Flat, Oceananigans.Grids.Flat} on Oceananigans.Architectures.CPU with 0×0×0 halo
├── ice_thermodynamics: SlabThermodynamics
├── advection: Nothing
└── external_heat_fluxes:
├── top: FluxFunction of slab_internal_heat_flux with parameters (flux=ConductiveFlux{Float64}, liquidus=ClimaSeaIce.SeaIceThermodynamics.LinearLiquidus{Float64}, bottom_heat_boundary_condition=ClimaSeaIce.SeaIceThermodynamics.HeatBoundaryConditions.IceWaterThermalEquilibrium{Int64})
└── bottom: FluxFunction of frazil_ice_formationNote that the default bottom heat boundary condition for SlabSeaIceThermodynamics is IceWaterThermalEquilibrium with freshwater. That's what we want!
model.ice_thermodynamics.heat_boundary_conditions.bottomClimaSeaIce.SeaIceThermodynamics.HeatBoundaryConditions.IceWaterThermalEquilibrium{Int64}(0)Ok, we're ready to freeze the bucket for 10 straight days. The ice will start forming suddenly due to the frazil ice heat flux and then eventually grow more slowly.
simulation = Simulation(model, Δt=10minute, 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 entriesCollecting data and running the simulation
Before simulating the freezing bucket, we set up a Callback to create a timeseries of the ice thickness saved at every time step.
# Container to hold the data
timeseries = []
# Callback function to collect the data from the `sim`ulation
function accumulate_timeseries(sim)
h = sim.model.ice_thickness
ℵ = sim.model.ice_concentration
push!(timeseries, (time(sim), first(h), first(ℵ)))
end
# Add the callback to `simulation`
simulation.callbacks[:save] = Callback(accumulate_timeseries)Callback of accumulate_timeseries on IterationInterval(1)Now we're ready to hit the Big Red Button (it should run pretty quick):
run!(simulation)[ Info: Initializing simulation...
[ Info: ... simulation initialization complete (2.610 seconds)
[ Info: Executing initial time step...
[ Info: ... initial time step complete (887.578 ms).
[ Info: Simulation is stopping after running for 0 seconds.
[ Info: Simulation time 10 days equals or exceeds stop time 10 days.Visualizing the result
It'd be a shame to run such a "cool" simulation without looking at the results. We'll visualize it with Makie.
using CairoMakietimeseries is a Vector of Tuple. So we have to do a bit of processing to build Vectors of time t and thickness h. It's not much work though:
t = [datum[1] for datum in timeseries]
h = [datum[2] for datum in timeseries]
ℵ = [datum[3] for datum in timeseries]
V = h .* ℵ1441-element Vector{Float64}:
0.0
1.7967665389365299e-6
0.0007223000840850853
0.0014427517483228223
0.002158006379959062
0.002868134797340176
0.003573218484895081
0.004273339343467654
0.004968578034121583
0.005659013552915546
⋮
0.30803096270311536
0.3081473858342732
0.3082637655349162
0.3083801018527354
0.3084963948353363
0.30861264453023896
0.30872885098487823
0.30884501424660404
0.3089611343626815Just for fun, we also compute the velocity of the ice-water interface:
dVdt = @. (h[2:end] .* ℵ[2:end] - h[1:end-1] .* ℵ[1:end-1]) / simulation.Δt1440-element Vector{Float64}:
2.9946108982275498e-9
1.2008388625769145e-6
1.2007527737295617e-6
1.1920910527270663e-6
1.1835473623018568e-6
1.175139479258175e-6
1.1668680976209543e-6
1.1587311510898817e-6
1.1507258646566054e-6
1.1428492656136997e-6
⋮
1.9411101574930608e-7
1.9403855192973163e-7
1.939661677383131e-7
1.9389386303201677e-7
1.9382163766817896e-7
1.9374949150441355e-7
1.9367742439879708e-7
1.9360543620968357e-7
1.9353352679579716e-7All that's left, really, is to put those lines! in an Axis:
set_theme!(Theme(fontsize=24, linewidth=4))
fig = Figure(size=(1600, 700))
axh = Axis(fig[1, 1], xlabel="Time (days)", ylabel="Ice thickness (cm)")
axℵ = Axis(fig[1, 2], xlabel="Time (days)", ylabel="Ice concentration (-)")
axV = Axis(fig[1, 3], xlabel="Ice Volume (cm)", ylabel="Freezing rate (μm s⁻¹)")
lines!(axh, t ./ day, 1e2 .* h)
lines!(axℵ, t ./ day, ℵ)
lines!(axV, 1e2 .* V[1:end-1], 1e6 .* dVdt)
save("freezing_bucket.png", fig)
If you want more ice, you can increase simulation.stop_time and run!(simulation) again (or just re-run the whole script).
This page was generated using Literate.jl.