Private types and functions

ForwardEulerTimeStepper(grid, prognostic_fields;
                        implicit_solver = nothing,
                        Gⁿ = map(similar, prognostic_fields))

Return a first-order Forward-Euler timestepper (ForwardEulerTimeStepper) on grid, with tracers. The tendency fields Gⁿ, usually equal to the prognostic_fields passed as positional argument, can be specified via optional kwargs.

The first-order Forward-Euler timestepper steps forward the state Uⁿ by Δt via

Uⁿ⁺¹ = Uⁿ + Δt * Gⁿ

where Uⁿ is the state at the $n$-th timestep and Gⁿ is the tendency at the $n$-th timestep.

dynamic_time_step!(model::RKSeaIceModel, Δt)

Update ice concentration ℵ and content 𝓋 = ℵ·h using the advective tendencies stored in model.timestepper.Gⁿ for a Runge–Kutta substep, then recover the thickness h = 𝓋/ℵ.

Unlike the Forward Euler version, this method updates the state from the cached previous state Ψ⁻ (stored by cache_current_fields!) rather than from the already-updated fields:

\[\begin{align*} 𝓋ⁿ⁺¹ = ℵⁿ hⁿ + Δt \, Gⁿ_𝓋 \\ ℵⁿ⁺¹ = ℵⁿ + Δt \, Gⁿ_ℵ \end{align*}\]

where $hⁿ$ and $ℵⁿ$ are retrieved from model.timestepper.Ψ⁻.

The kernel _dynamic_step_tracers! also handles:

• Clipping negative content and concentration values
• Recovering the thickness h = 𝓋/ℵ
• Ridging: when ℵ > 1, concentration is capped at 1 and thickness is adjusted to conserve ice content

Arguments

• model: A SeaIceModel using SplitRungeKuttaTimeStepper, following the low-storage framework of Silvestri et al. 2026.
• Δt: The time increment for this substep.

See also: rk_substep!, cache_current_fields!.

References

• Silvestri, S., Campin, J.-M., Wagner, G. L., Constantinou, N. C., Lee, X. K., and Ferrari, R. (2026). A low-storage Runge-Kutta framework for nonlinear free-surface ocean models. J. Adv. Model. Earth Sy. Submitted in April 2026. doi:10.22541/essoar.15002225/v1.

cache_current_fields!(model::RKSeaIceModel)

Cache the current prognostic fields (ice thickness h, ice concentration ℵ, and any additional tracers) into the timestepper's Ψ⁻ storage before performing a Runge–Kutta substep.

This function is called by Oceananigans' SplitRungeKuttaTimeStepper at the beginning of each full time step to store the state Uⁿ that is needed for computing the RK3 weighted average in the low-storage framework of Silvestri et al. 2026.

See also: rk_substep!, dynamic_time_step!.

References

• Silvestri, S., Campin, J.-M., Wagner, G. L., Constantinou, N. C., Lee, X. K., and Ferrari, R. (2026). A low-storage Runge-Kutta framework for nonlinear free-surface ocean models. J. Adv. Model. Earth Sy. Submitted in April 2026. doi:10.22541/essoar.15002225/v1.

rk_substep!(model::RKSeaIceModel, Δτ, callbacks)

Perform a single Runge–Kutta substep for the sea-ice model, advancing the state by Δτ.

The substep consists of three sequential operations:

1. Dynamic step: Compute advective tendencies and update ice thickness h and concentration ℵ via dynamic_time_step!.

2. Thermodynamic step: Apply column physics (melting/freezing) via thermodynamic_time_step!. This step is performed all at once since thermodynamics is local column physics.

3. Momentum step: Advance ice velocities using either an implicit or split-explicit scheme via time_step_momentum!.

This function is called by Oceananigans' SplitRungeKuttaTimeStepper for each of the three RK3 substeps within a full time step, following the low-storage framework of Silvestri et al. 2026.

Arguments

• model: A SeaIceModel using SplitRungeKuttaTimeStepper.
• Δτ: The substep time increment (a fraction of the full time step Δt).
• callbacks: Callbacks to execute during the substep (currently unused).

See also: cache_current_fields!, dynamic_time_step!.

References

• Silvestri, S., Campin, J.-M., Wagner, G. L., Constantinou, N. C., Lee, X. K., and Ferrari, R. (2026). A low-storage Runge-Kutta framework for nonlinear free-surface ocean models. J. Adv. Model. Earth Sy. Submitted in April 2026. doi:10.22541/essoar.15002225/v1.

LinearLiquidus(FT = Oceananigans.defaults.FloatType;
               salinity_slope = 0.0542,            # ᵒC / (g kg⁻¹)
               depth_slope = 7.89e-4,              # ᵒC / m
               freshwater_melting_temperature = 0) # ᵒC

Return a liquidus that depends linearly on salinity $S$ and on the height $z$ relative to the sea surface,

\[Tₘ(S, z) = T₀ - m S + λ z ,\]

where $T₀$ is the freshwater_melting_temperature, $m$ is the salinity_slope, and $λ$ is the depth_slope. With $m > 0$ and $λ > 0$ the melting temperature decreases as salinity increases and as depth increases ($z < 0$ below the sea surface). The height $z$ stands in for the water pressure $p ≈ - ρ g z$ and is therefore measured from the sea surface: at the base of floating ice it is set by the weight of the overlying ice and snow (see IceWaterThermalEquilibrium).

melting_temperature(liquidus, S) returns the surface value, melting_temperature(liquidus, S, z) the value at z.

The defaults fit the TEOS-10 freezing point of air-saturated seawater, in Conservative Temperature, within 0.02 K over S = 28-36 g kg⁻¹ and the top 1000 m, and assume that salinity is in g kg⁻¹ and temperature is in degrees Celsius.

ProportionalEvolution()

Parameterize sea-ice concentration changes by partitioning thermodynamic volume change between lateral and vertical growth using the proportional-evolution rule introduced by Hibler 1979.

References

• Hibler, W. D. III (1979). A Dynamic Thermodynamic Sea Ice Model. Journal of Physical Oceanography, 9(4), 815-846. doi:10.1175/1520-0485(1979)009<0815:ADTSIM>2.0.CO;2.

flux_kernel(flux)

Return the kernel function that computes the slab's internal heat flux given a raw flux-coefficient flux. This is a public extension point: users who define a custom flux struct should add a method

ClimaSeaIce.SeaIceThermodynamics.flux_kernel(::MyFlux) = my_flux_kernel

where my_flux_kernel(i, j, grid, Tu, clock, fields, parameters) returns a heat flux with parameters.flux::MyFlux.

Minimal example:

struct MyFlux{T}
    coefficient :: T
end

@inline function my_flux_kernel(i, j, grid, Tu, clock, fields, parameters)
    flux = parameters.flux
    return flux.coefficient * (fields.h[i, j, 1] - Tu)
end

ClimaSeaIce.SeaIceThermodynamics.flux_kernel(::MyFlux) = my_flux_kernel

Built-in dispatches:

• ConductiveFlux → slab_internal_heat_flux (single-layer Fourier)
• IceSnowConductiveFlux → ice_snow_conductive_flux (resistors in series)
• Function → returned directly (the function is its own kernel)

internal_flux_function(flux, phase_transitions, bottom_heat_boundary_condition)

Wrap a raw internal-flux coefficient (ConductiveFlux, IceSnowConductiveFlux, user Function, or user struct) in the FluxFunction shape expected by the surface-temperature solver and by getflux. The wrapper is built at the tendency-kernel level so that phase_transitions and bottom_heat_boundary_condition can be read from model.phase_transitions and the slab's heat BCs without threading those values through SlabThermodynamics at construction time.

If flux is already a FluxFunction, it is returned unchanged — the user is assumed to have fully assembled the wrapper themselves, and the kernel does not inject its own parameters.

To plug a custom flux-coefficient struct into the bare-ice slab, extend the flux_kernel dispatch:

struct MyInternalFlux{T}
    some_parameter :: T
end

# Signature must match the standard FluxFunction kernel:
#   (i, j, grid, Tu, clock, fields, parameters) -> Q
@inline function my_internal_flux(i, j, grid, Tu, clock, fields, parameters)
    flux = parameters.flux         # ::MyInternalFlux
    # ...use `flux.some_parameter`, `fields.h`, etc...
end

ClimaSeaIce.SeaIceThermodynamics.flux_kernel(::MyInternalFlux) = my_internal_flux

Once flux_kernel dispatches, MyInternalFlux can be passed to any SlabThermodynamics/sea_ice_slab_thermodynamics/snow_slab_thermodynamics constructor as internal_heat_flux = MyInternalFlux(...).

The layered snow+ice path additionally assumes that each layer's flux coefficient carries a .conductivity field and that the coupling is resistors-in-series. Custom flux types for a layered column are out of scope for this refactor.

latent_heat(phase_transitions::PhaseTransitions, T)

Return the per-mass latent heat of fusion of pure ice at temperature T,

\[ℒ(T) = ℒ₀ + \left(\frac{ρ_ℓ c_ℓ}{ρ} - c\right)(T - T₀) ,\]

where $ρ$, $c$ are the microscopic pure-ice density and heat capacity, $ρ_ℓ$, $c_ℓ$ are the liquid density and heat capacity, and $T₀$ is the reference temperature at which the reference latent heat $ℒ₀$ is defined.

This is the per-mass form of the volumetric expression $ρ ℒ(T) = ρ ℒ₀ + (ρ_ℓ c_ℓ - ρ c)(T - T₀)$ (divided through by $ρ$).

The returned quantity is per unit mass of pure ice. To obtain energy per unit volume of a porous medium (snow or sea ice), multiply by the bulk density of that medium.

IceWaterThermalEquilibrium(; salinity = 0)

Represents an ice-water interface in heat equilibrium, such that the bottom temperature $T_b$ is equal to the melting temperature at the base of the ice,

\[T_b = Tₘ(S, z_b) ,\]

where $S$ is the salinity at the ice-water boundary and $z_b$ is the height of the base of floating ice relative to the sea surface (see ice_base_height), set by the ice and snow load and by the liquid_density of the PhaseTransitions.

Both freezing and melting may occur at an ice-water boundary.

ice_base_height(hi, hs, ρi, ρs, ρw)

Return the height of the base of floating ice of thickness hi and density ρi, covered by snow of thickness hs and density ρs, relative to the sea surface,

\[z_b = - \frac{ρ_i h_i + ρ_s h_s}{ρ_w} ,\]

where $ρ_w$ is the density of the water. The water pressure at $z_b$ is the weight of the overlying ice and snow, independently from the elevation of the free surface.

EnthalpyMethodSeaIceModel(; grid,
                            closure = default_closure(grid),
                            ice_heat_capacity = 2090.0 / reference_density,
                            water_heat_capacity = 3991.0 / reference_density,
                            fusion_enthalpy = 3.3e5 / reference_density,
                            boundary_conditions = NamedTuple())

Return a thermodynamic model for ice sandwiched between an atmosphere and ocean on an Eulerian grid.

Calculate the right-hand-side of the free surface displacement (η) equation.

initialize_rheology!(model, rheology::ElastoViscoPlasticRheology)

Initialize the elasto-visco-plastic rheology. In this step we calculate the ice strength given the ice mass (thickness and concentration).

Compute the sea ice $u$-velocity tendencies.

Compute the sea ice $v$-velocity tendencies.