Sea Ice Thermodynamics

Sea ice thermodynamics governs the freezing and melting of ice through heat exchange with the atmosphere above and the ocean below. ClimaSeaIce implements a slab sea ice model, which represents the ice as a single, vertically-integrated layer with uniform properties.

The slab model

In the slab approximation, the ice is characterized by its thickness $h$, concentration $\aleph$, and a single temperature at the top surface $T_u$. The bottom of the ice sits at the ice-ocean interface, where temperature equals the melting point.

The evolution of ice thickness follows the Stefan condition at both interfaces:

\[\frac{\partial h}{\partial t} = w_u + w_b\]

where $w_u$ and $w_b$ are the interface velocities at the top (atmosphere-ice) and bottom (ice-ocean) interfaces, respectively.

Phase transitions and the liquidus

The melting temperature of sea ice depends on its salinity through the liquidus relationship. ClimaSeaIce uses a linear liquidus:

\[T_m(S) = T_0 - m S\]

where $T_m$ is the melting temperature, $T_0$ is the freshwater melting temperature (0°C by default), $S$ is the salinity, and $m$ is the liquidus slope (0.054 °C/psu by default).

The PhaseTransitions struct encapsulates the thermodynamic properties governing ice-water transitions:

using ClimaSeaIce.SeaIceThermodynamics: PhaseTransitions

phase_transitions = PhaseTransitions()

PhaseTransitions{Float64}
├── density: 917.0
├── heat_capacity: 2000.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)

Temperature-dependent latent heat

The latent heat of fusion $\mathscr{L}$ is the energy per unit mass of pure ice required to transform ice into liquid water (or released during freezing). It varies with temperature through a Stefan correction:

\[\mathscr{L}(T) = \mathscr{L}_0 + \left(\frac{\rho_\ell c_\ell}{\rho_i} - c_i\right)(T - T_0)\]

where $\rho_i, c_i$ are the microscopic pure-ice density and heat capacity, $\rho_\ell, c_\ell$ are the liquid water density and heat capacity, $\mathscr{L}_0$ is the reference latent heat at temperature $T_0$, and $T$ is the current temperature. All quantities here are material properties of pure H₂O and identical for sea ice and snow (the ice crystals are the same substance in both media).

To obtain energy per unit volume of a porous medium (snow or sea ice), multiply by the bulk density of that medium. This separation of the mass-budget multiplier from the per-mass Stefan correction matches the CICE / CLM convention.

using ClimaSeaIce.SeaIceThermodynamics: latent_heat

# Per-mass latent heat at different temperatures
L_at_0C = latent_heat(phase_transitions, 0.0)
L_at_minus10C = latent_heat(phase_transitions, -10.0)
println("Latent heat at 0°C: ", L_at_0C, " J/kg")
println("Latent heat at -10°C: ", L_at_minus10C, " J/kg")

Latent heat at 0°C: 334000.0 J/kg
Latent heat at -10°C: 308360.2748091603 J/kg

Heat fluxes and the Stefan condition

The slab model computes thickness changes from the balance of heat fluxes at each interface.

Internal conductive flux

Within the ice slab, heat is conducted from the warmer bottom to the colder top surface. The conductive flux is:

\[Q_i = -k \frac{T_u - T_b}{h}\]

where $k$ is the thermal conductivity (default 2 W/(m·K) for freshwater ice), $T_u$ is the top surface temperature, $T_b$ is the bottom temperature (at the melting point), and $h$ is the ice thickness.

The ConductiveFlux struct implements this:

using ClimaSeaIce.SeaIceThermodynamics: ConductiveFlux

# Thermal conductivity of 2 W/(m·K)
conductive_flux = ConductiveFlux(Float64; conductivity = 2.0)

ConductiveFlux{Float64}(2.0)

Interface velocities

At each interface, the difference between incoming and outgoing heat fluxes drives melting or freezing:

\[w_u = \frac{Q_x - Q_i}{\mathscr{L}(T_u)}, \quad w_b = \frac{Q_i - Q_b}{\mathscr{L}(T_b)}\]

where:

• $Q_x$ is the external (atmospheric) heat flux into the top surface
• $Q_i$ is the internal conductive flux
• $Q_b$ is the oceanic heat flux at the bottom

Negative velocities indicate melting; positive velocities indicate freezing.

Top surface boundary conditions

The top surface temperature $T_u$ can be determined in two ways.

Melting-constrained flux balance

The default approach, MeltingConstrainedFluxBalance, solves for the temperature that equilibrates the external and internal heat fluxes:

\[\sum_n Q_{x,n}(T_u) - Q_i(T_u) = 0\]

subject to the constraint that the surface cannot exceed the melting temperature:

\[T_u \leq T_m(S)\]

When the surface would otherwise be warmer than the melting point, it is clamped to $T_m(S)$ and the excess heat drives surface melting.

using ClimaSeaIce.SeaIceThermodynamics: MeltingConstrainedFluxBalance

top_bc = MeltingConstrainedFluxBalance()

MeltingConstrainedFluxBalance{ClimaSeaIce.SeaIceThermodynamics.HeatBoundaryConditions.NonlinearSurfaceTemperatureSolver}(ClimaSeaIce.SeaIceThermodynamics.HeatBoundaryConditions.NonlinearSurfaceTemperatureSolver())

Prescribed temperature

Alternatively, the surface temperature can be prescribed directly using PrescribedTemperature:

using ClimaSeaIce.SeaIceThermodynamics: PrescribedTemperature

# Fix the surface at -5°C
top_bc = PrescribedTemperature(-5.0)

PrescribedTemperature{Float64}(-5.0)

Bottom boundary condition

The default bottom boundary condition, IceWaterThermalEquilibrium, assumes the ice-ocean interface is at the salinity-dependent melting temperature:

\[T_b = T_m(S)\]

This is appropriate when the ocean mixed layer is well-mixed and maintains thermal equilibrium with the ice bottom.

External heat fluxes

External heat fluxes drive the thermodynamic evolution of sea ice.

Radiative emission

The RadiativeEmission flux represents longwave radiation emitted by the ice surface:

\[Q_r = \varepsilon \sigma (T + T_{\text{ref}})^4\]

where $\varepsilon$ is the emissivity, $\sigma$ is the Stefan-Boltzmann constant, and $T_{\text{ref}}$ converts from Celsius to Kelvin.

Custom flux functions

The FluxFunction wrapper allows arbitrary user-defined heat fluxes that can depend on position, time, surface temperature, and model fields:

using ClimaSeaIce.SeaIceThermodynamics: FluxFunction

# A simple sensible heat flux depending on surface temperature
sensible_heat(Tu, clock, fields, parameters) = parameters.coefficient * (parameters.T_air - Tu)
flux = FluxFunction(sensible_heat; parameters = (coefficient = 15.0, T_air = -10.0))

FluxFunction of sensible_heat (generic function with 1 method) with parameters @NamedTuple{coefficient::Float64, T_air::Float64}

Putting it together: SlabThermodynamics

The SlabThermodynamics struct describes a single slab layer (sea ice or snow). It carries the top surface temperature field, the top/bottom heat boundary conditions, the internal conductive flux coefficient, and the concentration-evolution rule. Phase-transition parameters (densities, latent heat, liquidus) are not stored here — they live on the SeaIceModel as phase_transitions, shared between the ice and the (optional) snow layer.

using Oceananigans
using ClimaSeaIce
using ClimaSeaIce.SeaIceThermodynamics: MeltingConstrainedFluxBalance

grid = RectilinearGrid(size=(10, 10), x=(0, 1e5), y=(0, 1e5), topology=(Periodic, Periodic, Flat))

ice_thermodynamics = SlabThermodynamics(grid;
    top_heat_boundary_condition = MeltingConstrainedFluxBalance(),
    internal_heat_flux = ConductiveFlux(Float64; conductivity = 2.0))

ice_thermodynamics

SlabThermodynamics
└── top_surface_temperature: 10×10×1 Field{Oceananigans.Grids.Center, Oceananigans.Grids.Center, Nothing} reduced over dims = (3,) on Oceananigans.Grids.RectilinearGrid on CPU

The bulk sea-ice density is passed to SeaIceModel separately via the sea_ice_density keyword (default 900 kg/m³), and is stored as a Field{Center, Center, Nothing} on the model. It is distinct from the microscopic pure-ice density inside PhaseTransitions, which refers to the pure ice crystal and not to the porous sea-ice matrix.

Customising the internal flux

The value you pass via internal_heat_flux = ... to SlabThermodynamics (and to its sea_ice_slab_thermodynamics and snow_slab_thermodynamics helpers) is the raw coefficient of the layer's internal flux — not a FluxFunction. The tendency kernel wraps it in a FluxFunction on the fly, threading in model.phase_transitions.liquidus and the slab's bottom boundary condition as parameters. This avoids baking the shared liquidus into each slab at construction time.

Four shapes are supported out of the box:

A ConductiveFlux (default for both ice and snow)

using ClimaSeaIce.SeaIceThermodynamics: ConductiveFlux

ice_thermodynamics = SlabThermodynamics(grid;
    internal_heat_flux = ConductiveFlux(Float64; conductivity = 2.0))

SlabThermodynamics
└── top_surface_temperature: 10×10×1 Field{Oceananigans.Grids.Center, Oceananigans.Grids.Center, Nothing} reduced over dims = (3,) on Oceananigans.Grids.RectilinearGrid on CPU

Computes Q_i = -k (T_u - T_b) / h_i. The most common case.

A bare Function

Pass a kernel function directly. Its signature must match the standard FluxFunction calling convention — (i, j, grid, Tu, clock, fields, parameters) -> Q:

@inline function radiative_internal_flux(i, j, grid, Tu, clock, fields, parameters)
    hi = @inbounds fields.h[i, j, 1]
    # parameters.flux is the function itself; parameters.liquidus and
    # parameters.bottom_heat_boundary_condition are injected by the slab.
    return ifelse(hi ≤ 0, zero(hi), -5.0 * Tu)  # toy radiative model
end

ice_thermodynamics_rad = SlabThermodynamics(grid;
    internal_heat_flux = radiative_internal_flux)

SlabThermodynamics
└── top_surface_temperature: 10×10×1 Field{Oceananigans.Grids.Center, Oceananigans.Grids.Center, Nothing} reduced over dims = (3,) on Oceananigans.Grids.RectilinearGrid on CPU

A fully-assembled FluxFunction

If you have already packaged the kernel, its parameters, and the temperature-dependence flag in a FluxFunction, you can pass it directly. The slab stores it unchanged; the tendency kernel does not inject its own liquidus / bottom_heat_boundary_condition:

using ClimaSeaIce.SeaIceThermodynamics: FluxFunction

user_flux = FluxFunction(radiative_internal_flux;
                         parameters = (flux = nothing, absorption = 0.8),
                         top_temperature_dependent = true)

ice_thermodynamics_fn = SlabThermodynamics(grid;
    internal_heat_flux = user_flux)

SlabThermodynamics
└── top_surface_temperature: 10×10×1 Field{Oceananigans.Grids.Center, Oceananigans.Grids.Center, Nothing} reduced over dims = (3,) on Oceananigans.Grids.RectilinearGrid on CPU

This is the right choice when your flux needs parameters that the default liquidus/bottom-BC injection cannot supply.

A custom flux struct (extend flux_kernel)

For more structured flux descriptions, define a concrete type and provide one new method:

struct NonLinearConductiveFlux{T}
    k0 :: T
    α  :: T   # temperature sensitivity
end

@inline function nonlinear_conductive_flux(i, j, grid, Tu, clock, fields, parameters)
    flux = parameters.flux           # ::NonLinearConductiveFlux
    bottom_bc = parameters.bottom_heat_boundary_condition
    liquidus = parameters.liquidus
    Tb = bottom_temperature(i, j, grid, bottom_bc, liquidus)
    hi = @inbounds fields.h[i, j, 1]
    k_eff = flux.k0 * (1 + flux.α * Tu)
    return ifelse(hi ≤ 0, zero(hi), -k_eff * (Tu - Tb) / hi)
end

# Extend the dispatch so the slab knows how to wrap NonLinearConductiveFlux
ClimaSeaIce.SeaIceThermodynamics.flux_kernel(::NonLinearConductiveFlux) = nonlinear_conductive_flux

SlabThermodynamics(grid;
    internal_heat_flux = NonLinearConductiveFlux(2.0, 0.01))

The benefit over the bare-function approach is that the flux coefficient is introspectable (show, GPU adapt, user copying, etc.) and carries its own type parameters.

When coupled to a SeaIceModel, these thermodynamics automatically compute thickness tendencies based on the configured heat fluxes and boundary conditions.

Snow thermodynamics

ClimaSeaIce supports a layered snow model where a snow layer sits on top of the ice slab. Snow insulates the ice by adding thermal resistance in series, reducing the conductive flux and slowing ice growth.

Snow layer physics

When snow is present, the combined conductive flux through the snow+ice slab uses resistors in series:

\[F_c = \frac{T_b - T_u}{h_s / k_s + h_i / k_i}\]

where $h_s, k_s$ are the snow thickness and conductivity, and $h_i, k_i$ are the ice thickness and conductivity. Snow typically has $k_s \approx 0.31$ W/(m·K) compared to $k_i \approx 2$ W/(m·K) for ice, making it an effective insulating layer.

The snow surface temperature $T_u$ is determined by balancing the atmospheric heat fluxes against the combined conductive flux, using the same MeltingConstrainedFluxBalance approach as bare ice. The melting point for snow is 0°C (rather than the salinity-dependent liquidus for ice).

Interface temperature

The snow-ice interface temperature $T_{si}$ is computed analytically from the surface temperature using the thermal resistance ratio:

\[T_{si} = T_b + (T_u - T_b) \frac{R_i}{R_s + R_i}\]

where $R_i = h_i / k_i$ and $R_s = h_s / k_s$. This temperature is passed to the ice thermodynamics as a prescribed boundary condition, so the ice layer is completely unaware of the snow layer above it.

Snow processes

The snow model includes three processes that modify the snow thickness each time step:

1. Surface melt: When the surface temperature reaches the melting point, any excess conductive flux first melts snow before reaching the ice. The snow absorbs energy $\rho_s L_s \Delta h_s$ before excess energy melts ice from the top.

2. Snowfall accumulation: Snow precipitates at a prescribed rate (kg/m²/s) and accumulates only where ice is present ($\aleph > 0$).

3. Snow-ice formation (flooding): When the snow load pushes the ice surface below the waterline (negative freeboard), seawater floods the snow layer and converts it to ice. The freeboard criterion is:

   \[h_f = h_i (1 - \rho_i / \rho_w) - h_s \rho_s / \rho_w < 0\]

   The flooding step formulated like this is energy-conserving

Snow volume conservation

Snow thickness $h_s$ represents the thickness over the ice-covered fraction $\aleph$. The total snow energy per unit area is $\rho_s L_s h_s \aleph$, consistent with how ice energy is $\mathscr{L} h \aleph$. When ice concentration changes (e.g., during consolidation), the snow thickness is adjusted to conserve the snow volume $h_s \times \aleph$.

Setting up a snow-covered model

Pass a second SlabThermodynamics via snow_thermodynamics (or use the snow_slab_thermodynamics helper, which defaults the snow conductivity to 0.31 W m⁻¹ K⁻¹). The model reads the bulk snow density from the snow_density keyword (default 330 kg/m³).

using ClimaSeaIce

grid = RectilinearGrid(size=(), topology=(Flat, Flat, Flat))

snow_thermodynamics = snow_slab_thermodynamics(grid)

model = SeaIceModel(grid;
    snow_thermodynamics,
    snow_density = 330,
    snowfall = 3e-6) # kg/m²/s

set!(model, h=1, ℵ=1, hs=0.1) # 10 cm snow on 1 m ice

The combined snow+ice conductive flux is constructed inline inside the layered thermodynamic kernel from each layer's conductivity, and the snow-ice interface temperature $T_{si}$ is computed analytically at each step from the resistance ratio.

Dynamics-only configurations

When running a pure-dynamics simulation (no phase physics), pass ice_thermodynamics = nothing. SeaIceModel still carries the bulk sea_ice_density field that the momentum equations need, so the dynamics path works without a thermodynamic step:

model_dyn = SeaIceModel(grid;
    ice_thermodynamics = nothing,
    sea_ice_density = 900)

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: Nothing
├── snow_thermodynamics: Nothing
├── advection: Nothing
└── external_heat_fluxes: 
    ├── top: Nothing
    └── bottom: Int64