GitHub

Freezing bucket example

A common laboratory experiment freezes an insulated 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: 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). This example demonstrates how to:

  • use SlabSeaIceThermodynamics with prescribed boundary conditions,
  • configure internal heat conduction,
  • use FluxFunction for frazil ice formation,
  • collect and visualize time series data.

Install dependencies

First let's make sure we have all required packages installed.

using Pkg
pkg"add Oceananigans, ClimaSeaIce, CairoMakie"

Using ClimaSeaIce.jl

Write

using Oceananigans
using Oceananigans.Units
using ClimaSeaIce
using CairoMakie

to load ClimaSeaIce functions and objects into our script.

An infinitely deep bucket with a single grid point

Perhaps surprisingly, we need just one grid point to model a 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, Flat, Flat, Flat} on CPU with 0×0×0 halo
├── Flat x 
├── Flat y 
└── Flat z

Model configuration

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 # W m⁻¹ 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}
├── ice_density: 900.0
├── ice_heat_capacity: 2100.0
├── liquid_density: 999.8
├── liquid_heat_capacity: 4186.0
├── reference_latent_heat: 334000.0
├── reference_temperature: 0.0
└── liquidus: LinearLiquidus(freshwater_melting_temperature = 0.0, slope = 0.054)

We set the top ice temperature:

top_temperature = -10 # °C
top_heat_boundary_condition = PrescribedTemperature(top_temperature)
PrescribedTemperature{Int64}(-10)

Constructing the thermodynamics

We construct the ice thermodynamics using a simple slab sea ice representation:

ice_thermodynamics = SlabSeaIceThermodynamics(grid;
                                              internal_heat_flux,
                                              phase_transitions,
                                              top_heat_boundary_condition)
SlabSeaIceThermodynamics
└── top_surface_temperature: ConstantField(-10)

Frazil ice formation

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_formation (generic function with 1 method)

Building the model

Then 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, Flat, Flat, Flat} on CPU with 0×0×0 halo
├── timestepper: SplitRungeKuttaTimeStepper(3)
├── ice_thermodynamics: SlabThermodynamics
├── advection: Nothing
└── external_heat_fluxes: 
    ├── top: FluxFunction of slab_internal_heat_flux (generic function with 2 methods) with parameters @NamedTuple{flux::ConductiveFlux{Float64}, liquidus::ClimaSeaIce.SeaIceThermodynamics.LinearLiquidus{Float64}, bottom_heat_boundary_condition::ClimaSeaIce.SeaIceThermodynamics.HeatBoundaryConditions.IceWaterThermalEquilibrium{Int64}}
    └── bottom: FluxFunction of frazil_ice_formation (generic function with 1 method)

Note that the default bottom heat boundary condition for SlabSeaIceThermodynamics is IceWaterThermalEquilibrium with freshwater. That's what we want!

model.ice_thermodynamics.heat_boundary_conditions.bottom
ClimaSeaIce.SeaIceThermodynamics.HeatBoundaryConditions.IceWaterThermalEquilibrium{Int64}(0)

Running a simulation

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
├── run_wall_time: 0 seconds
├── run_wall_time / 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

Collecting data

Before simulating the freezing bucket, we set up a Callback to create a time series of the ice thickness saved at every time step:

timeseries = []

function accumulate_timeseries(sim)
    h = sim.model.ice_thickness
    ℵ = sim.model.ice_concentration
    push!(timeseries, (time(sim), first(h), first(ℵ)))
end

simulation.callbacks[:save] = Callback(accumulate_timeseries)
Callback of accumulate_timeseries on IterationInterval(1)

Now we're ready to run the simulation:

run!(simulation)
[ Info: Initializing simulation...
[ Info:     ... simulation initialization complete (1.146 seconds)
[ Info: Executing initial time step...
[ Info:     ... initial time step complete (719.980 ms).
[ Info: Simulation is stopping after running for 2.136 seconds.
[ Info: Simulation time 10 days equals or exceeds stop time 10 days.

Visualizing the results

It'd be a shame to run such a "cool" simulation without looking at the results. We'll visualize it with Makie:

timeseries 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.0007753433778587996
 0.0014957931316414758
 0.002210862057632628
 0.002920744476221077
 0.00362555428902386
 0.004325386595955387
 0.005020328691172857
 0.005710463347073555
 ⋮
 0.3080398908749131
 0.3081563106764367
 0.3082726870510988
 0.308389020046584
 0.3085053097104911
 0.30862155609033376
 0.30873775923354035
 0.30885391918745425
 0.30897003599933415

Just 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.Δt
1440-element Vector{Float64}:
 2.9946108982275498e-9
 1.2892443521997717e-6
 1.2007495896377936e-6
 1.191781543318587e-6
 1.183137364314082e-6
 1.1746830213379713e-6
 1.166387178219212e-6
 1.1582368253624497e-6
 1.1502244265011628e-6
 1.1423446039537195e-6
 ⋮
 1.9410546025888652e-7
 1.9403300253924744e-7
 1.9396062443684712e-7
 1.9388832580865175e-7
 1.9381610651190522e-7
 1.9374396640440643e-7
 1.936719053443244e-7
 1.9359992318982813e-7
 1.9352801979982682e-7

All 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.