Public Documentation

ClimaSeaIce.jl

ClimaSeaIce.jl is Julia software for simulating the freezing, melting, and horizontal motion of sea ice on CPUs and GPUs. It is designed for climate-scale sea-ice modeling, supports standalone simulations, and can be coupled to ocean models built with Oceananigans.jl.

What ClimaSeaIce provides

• Sea-ice thermodynamics for freezing, melting, conductive heat transfer, and surface and basal heat-flux parameterizations.
• Sea-ice dynamics with explicit and split-explicit momentum solvers.
• Multiple rheological closures, including viscous and elasto-visco-plastic options.
• Support for standalone sea-ice models and models coupled to Oceananigans-based ocean simulations.

Installation

Install the latest registered release with

using Pkg
Pkg.add("ClimaSeaIce")

To use the development version from GitHub, use

using Pkg
Pkg.add(url = "https://github.com/CliMA/ClimaSeaIce.jl")

Documentation

Documentation is available at:

• Stable documentation
• Development documentation

The documentation includes model setup, physics descriptions, examples, and internal implementation notes.

Quick start

The main entry point is SeaIceModel. Thermodynamic parameterizations live under SeaIceThermodynamics, and momentum closures and rheologies live under SeaIceDynamics and Rheologies.

Example scripts are available in examples/.

Citing

If you use ClimaSeaIce.jl in research, teaching, or derived software, please cite the Zenodo record:

Silvestri, S. et al. (2026). CliMA/ClimaSeaIce.jl: v0.5.6 (v0.5.6). Zenodo. https://doi.org/10.5281/zenodo.16143708

SeaIceModel(grid;
            clock                       = Clock{eltype(grid)}(time = 0),
            ice_consolidation_thickness = 0.05, # m
            ice_salinity                = 0,    # psu
            sea_ice_density             = 900,  # kg m⁻³, bulk sea-ice
            snow_density                = 330,  # kg m⁻³, bulk snow
            phase_transitions           = PhaseTransitions(eltype(grid)),
            top_heat_flux               = nothing,
            bottom_heat_flux            = 0,    # W m⁻²
            velocities                  = nothing,
            advection                   = nothing,
            tracers                     = (),
            timestepper                 = :SplitRungeKutta3,
            boundary_conditions         = NamedTuple(),
            ice_thermodynamics          = sea_ice_slab_thermodynamics(grid),
            snow_thermodynamics         = nothing,
            snowfall                    = 0,
            dynamics                    = nothing,
            forcing                     = NamedTuple())

Construct a sea-ice model on grid with slab thermodynamics, optional momentum equations, and optional snow thermodynamics.

The model evolves sea-ice thickness h, concentration ℵ, optional salinity S, optional snow thickness hs, and, when dynamics is provided, the horizontal ice velocity components u and v.

Arguments

• grid: The computational grid.

Keyword Arguments

• clock: Model clock. Defaults to Clock{eltype(grid)}(time = 0).
• ice_consolidation_thickness: Threshold thickness above which the slab is treated as consolidated ice. Default: 0.05 meters.
• ice_salinity: Sea-ice salinity field or constant. When non-constant it is included in the prognostic state as S. Default: 0 psu.
• sea_ice_density: Bulk sea-ice density. Default: 900 kg m⁻³.
• snow_density: Bulk snow density. Default: 330 kg m⁻³.
• phase_transitions: Thermodynamic phase-transition parameters shared by the ice and snow slabs. Default: PhaseTransitions(eltype(grid)).
• top_heat_flux: External top heat flux, or tuple of fluxes, passed to the ice upper boundary condition. Default: nothing.
• bottom_heat_flux: External bottom heat flux passed to the ice lower boundary condition. Default: 0 W m⁻².
• velocities: Pre-existing velocity fields. If omitted, they are allocated by the constructor, possibly on an extended grid required by the selected dynamics solver. Default: nothing.
• advection: Advection scheme for prognostic tracers. Default: nothing.
• tracers: Additional prognostic tracers.
• timestepper: Time stepper specification passed to TimeStepper. Default: :SplitRungeKutta3.
• boundary_conditions: Boundary conditions for allocated prognostic fields.
• ice_thermodynamics: Ice thermodynamics model. Default: sea_ice_slab_thermodynamics(grid).
• snow_thermodynamics: Optional snow thermodynamics model. When provided, the model allocates prognostic snow thickness hs.
• snowfall: Snowfall forcing applied only when snow_thermodynamics is present. May be a constant, Field, or FieldTimeSeries.
• dynamics: Optional sea-ice dynamics model, for example, SeaIceMomentumEquation(grid).
• forcing: Additional model forcing passed through model_forcing.

PhaseTransitions(FT=Oceananigans.defaults.FloatType;
                 density               = 917,    # kg m⁻³
                 heat_capacity         = 2000,   # J / (kg ᵒC)
                 liquid_density        = 999.8,  # kg m⁻³
                 liquid_heat_capacity  = 4186,   # J / (kg ᵒC)
                 reference_latent_heat = 334e3,  # J kg⁻³
                 reference_temperature = 0,      # ᵒC
                 liquidus = LinearLiquidus(FT))  # default assumes psu, ᵒC

Return a representation of transitions between the solid and liquid phases of salty water: in other words, the freezing and melting of sea ice.

PhaseTransitions stores the thermodynamic parameters shared by the slab sea-ice and snow parameterizations in SeaIceModel, including densities, heat capacities, a reference latent heat, and the liquidus relation.

The latent heat of fusion $ℒ(T)$ (more simply just "latent heat") is a function of temperature $T$ via

\[ρ ℒ(T) = ρ ℒ₀ + (ρ_ℓ c_ℓ - ρ c) (T - T₀)\]

where $ρ$ is the solid density, $ρ_ℓ$ is the liquid density, $c$ is the solid heat_capacity, $c_ℓ$ is the liquid heat capacity, and $T₀$ is a reference temperature, all of which are assumed constant.

The default liquidus assumes that salinity has practical salinity units (psu) and that temperature is degrees Celsius.

SlabThermodynamics(grid; kw...)

A minimal slab representation of a single sea-ice or snow layer.

The object stores:

• the prognostic top surface temperature,
• top and bottom heat boundary conditions,
• an internal heat-flux model (for example ConductiveFlux), and
• the concentration-evolution rule used when thermodynamic growth or melt changes ice volume.

Shared thermodynamic material properties such as densities, latent heat, and the liquidus relation are stored at the SeaIceModel level and threaded into the tendency kernels via model.phase_transitions.

sea_ice_slab_thermodynamics(grid; kw...)

Construct a SlabThermodynamics with default parameters appropriate for sea ice: conductivity = 2 W/(m K). Bulk density and all phase-transition parameters live on SeaIceModel (as sea_ice_density and phase_transitions respectively).

snow_slab_thermodynamics(grid;
                         conductivity = 0.31,
                         kw...)

Construct a SlabThermodynamics with default parameters appropriate for snow: conductivity = 0.31 W/(m K). Bulk density and all phase-transition parameters live on SeaIceModel (as snow_density and phase_transitions respectively).

FluxFunction(func; parameters=nothing, top_temperature_dependent=false)

Return a FluxFunction representing a flux across an air-ice, air-snow, or ice-water interface.

The wrapped callable func must have one of the following signatures:

flux = func(i, j, grid, clock, top_temperature, model_fields)

when isnothing(parameters), or

flux = func(i, j, grid, clock, top_temperature, model_fields, parameters)

when !isnothing(parameters).

Set top_temperature_dependent = true when the flux must be recomputed during the diagnostic solve for the surface temperature.

Example

@inline sensible_heat_flux(i, j, grid, Tₛ, clock, fields, coefficient) =
    coefficient * (fields.T_air[i, j, 1] - Tₛ)

Q = FluxFunction(sensible_heat_flux; parameters = 15.0,
                 top_temperature_dependent = true)

MeltingConstrainedFluxBalance(top_surface_temperature_solver=NonlinearSurfaceTemperatureSolver())

Return a boundary condition that determines the top (or "upper surface") temperature $Tᵤ$ to equilibrate the top heat fluxes,

\[Qₓ₁(Tᵤ) + Qₓ₂(Tᵤ) + ⋯ - Qᵢ = Σᴺ Qₓₙ(Tᵤ) - Qᵢ = δQ ,\]

where $Qᵢ$ is the intrinsic flux into the top from within the ice (typically, a conductive flux), $Qₓₙ$ represent external fluxes into the air above the ice, $δQ$ is the residual flux, and $Tᵤ$ is the (upper) top temperature. $Tᵤ$ is evaluated under the constraint that

\[Tᵤ ≤ Tₘ(S),\]

where $Tₘ(S)$ is the melting temperature, which is a function of the ice salinity at the top, $S$. When $Tᵤ < Tₘ(S)$, the top is frozen and $δQ = 0$. When the constraint operates, such that $Tᵤ = Tₘ(S)$, the top is melting and the residual flux is non-zero.

\[δQ ≡ Σᴺ Qₙ(Tₘ) - Qᵢ(Tₘ).\]

The residual flux is consumed by the cost of transforming ice into liquid water, and is related to the rate of change of ice thickness, $h$, by

\[\frac{dhₛ}{dt} = δQ / ℒ(Tᵤ)\]

where $ℒ(Tᵤ)$ is the latent heat, equal to the different between the higher internal energy of liquid water and the lower internal energy of solid ice, at the temperature $Tᵤ$.

PrescribedTemperature()

heat boundary condition indicating that temperature is prescribed on the boundary.

RadiativeEmission(FT=Oceananigans.defaults.FloatType;
                  emissivity = 1,
                  stefan_boltzmann_constant = 5.67e-8,
                  reference_temperature = 273.15)

Return a flux representing radiative emission from a surface.

Auxiliaries(rheology, grid)

A struct holding any auxiliary fields and kernels needed for the computation of sea ice stresses.

ElastoViscoPlasticRheology(FT = Oceananigans.defaults.FloatType;
                           ice_compressive_strength = 27500,
                           ice_compaction_hardening = 20,
                           yield_curve_eccentricity = 2,
                           minimum_plastic_stress = 2e-9,
                           min_relaxation_parameter = 50,
                           max_relaxation_parameter = 300,
                           relaxation_strength = π^2,
                           pressure_formulation = ReplacementPressure())

Construct an ElastoViscoPlasticRheology object representing a modified elasto-visco-plastic rheology for sea-ice dynamics following Kimmritz et al. 2017. The stress tensor is computed following the constitutive relation:

\[σᵢⱼ = 2η ϵ̇ᵢⱼ + [(ζ - η) (ϵ̇₁₁ + ϵ̇₂₂) - P / 2] δᵢⱼ\]

where $ϵ̇ᵢⱼ$ are the strain rates, $η$ is the shear viscosity, $ζ$ is the bulk viscosity, and $P$ is the ice strength (acting as the isotropic part of the stress tensor) parameterized as $P_\star h \exp[ - C ( 1 - ℵ )]$ where $P_\star$ is the ice_compressive_strength, $C$ is the ice_compaction_hardening, $h$ is the ice thickness, and $ℵ$ is the ice concentration.

The stresses are substepped using a dynamic substepping coefficient $α$ that is spatially varying and computed dynamically as done by Kimmritz et al. 2017. In particular: $α = \sqrt{γ²}$, where $γ² = ζ c_α (Δt / mᵢ) / A_z$ is a stability parameter with $A_z$ is the area of the grid cell, $mᵢ$ the ice mass, $Δt$ the time step, and $c_α$ a numerical stability parameter which controls the strength of $γ²$.

The stresses are substepped with:

\[σᵢⱼᵖ⁺¹ = σᵢⱼᵖ + (σᵢⱼᵖ⁺¹ - σᵢⱼᵖ) / α\]

This formulation allows fast convergence in regions where $α$ is small. Regions where $α$ is large correspond to regions where the ice is more solid and the convergence is slower. $α$ can be thought of as a "pseudo substep number" or a "relaxation parameter". If we are using a subcycling solver, then if $α$ ≪ number of substeps, the convergence is faster.

Arguments

• grid: The computational grid.

Keyword Arguments

• ice_compressive_strength: Parameter expressing compressive strength (in N m⁻²). Default: 27500.
• ice_compaction_hardening: Exponent coefficient for compaction hardening. Default: 20.
• yield_curve_eccentricity: Eccentricity of the elliptic yield curve. Default: 2.
• minimum_plastic_stress: Minimum value for the visco-plastic parameter. Limits the maximum viscosity of the ice, transitioning the rheology from plastic to viscous behavior. Default: 2e-9.
• min_relaxation_parameter: Minimum value for the relaxation parameter α. Default: 50.
• max_relaxation_parameter: Maximum value for the relaxation parameter α. Default: 300.
• relaxation_strength: Parameter controlling the strength of the relaxation parameter. The maximum value is π²; see Kimmritz et al. 2017. Default: π².
• pressure_formulation: Either ReplacementPressure() or IceStrength(). The replacement-pressure formulation avoids ice motion in the absence of forcing. Default: ReplacementPressure().

References

• Kimmritz, M., Losch, M., and Danilov, S. (2017). A comparison of viscous-plastic sea ice solvers with and without replacement pressure. Ocean Modelling, 115, 59-69. doi:10.1016/j.ocemod.2017.05.006.

a simple explicit solver

SeaIceMomentumEquation(grid;
                       coriolis = nothing,
                       rheology = ElastoViscoPlasticRheology(eltype(grid)),
                       top_momentum_stress    = nothing,
                       bottom_momentum_stress = nothing,
                       free_drift = nothing,
                       solver = SplitExplicitSolver(grid; substeps=150),
                       minimum_concentration = 1e-3,
                       minimum_mass = 1.0)

Constructs a SeaIceMomentumEquation object that controls the dynamical evolution of sea-ice momentum. The sea-ice momentum obey the following evolution equation:

\[\frac{∂\boldsymbol{u}}{∂t} + \boldsymbol{f} × \boldsymbol{u} = \frac{\boldsymbol{\nabla} \cdot \boldsymbol{\sigma}}{mᵢ} + \frac{\boldsymbol{\tau}ₒ}{mᵢ} + \frac{\boldsymbol{\tau}ₐ}{mᵢ}\]

where $∂\boldsymbol{u}/∂t$ is the time derivative of the ice velocity, $\boldsymbol{f}$ is the Coriolis parameter, $\boldsymbol{\nabla} \cdot \boldsymbol{\sigma} / mᵢ$ is the divergence of internal stresses, $\boldsymbol{\tau}ₒ/mᵢ$ is the ice-ocean boundary stress, $\boldsymbol{\tau}ₐ/mᵢ$ is the ice-atmosphere boundary stress, and $mᵢ = ρᵢ h ℵ$ is the ice mass per unit area.

Arguments

• grid: The computational grid.

Keyword Arguments

• coriolis: Parameters for the background rotation rate of the model.
• rheology: The sea-ice rheology model. Default: ElastoViscoPlasticRheology(eltype(grid)).
• top_momentum_stress: Atmosphere-to-ice momentum stress, or an object that can be materialized into one. Default: nothing.
• bottom_momentum_stress: Ocean-to-ice momentum stress, or an object that can be materialized into one. Default: nothing.
• free_drift: The free drift velocities used when nonzero sea ice mass or concentration are below the dynamical momentum thresholds. Default is nothing.
• solver: Momentum solver used to advance the velocity field. Default: SplitExplicitSolver(grid; substeps = 150).
• minimum_concentration: Minimum sea-ice concentration above which the velocity is evolved dynamically. Below this threshold, nonzero ice moves with free drift and roundoff-level concentration cells are set to zero. Default: 1e-3.
• minimum_mass: Minimum sea-ice mass per area above which the velocity is evolved dynamically. Below this threshold, nonzero ice moves with free drift and roundoff-level mass cells are set to zero. Default: 1.0 kg/m².

SemiImplicitStress(FT = Oceananigans.defaults.FloatType;
                   uₑ = ZeroField(FT),
                   vₑ = ZeroField(FT),
                   ρₑ = 1026.0,
                   Cᴰ = 5.5e-3)

A structure representing the semi-implicit stress between the sea ice and an external fluid (either the ocean or the atmosphere), calculated as:

\[\begin{align*} τᵤ & = ρₑ Cᴰ \sqrt{(uₑ - uᵢⁿ)² + (vₑ - vᵢⁿ)²} (uₑ - uᵢⁿ⁺¹) \\ τᵥ & = ρₑ Cᴰ \sqrt{(uₑ - uᵢⁿ)² + (vₑ - vᵢⁿ)²} (vₑ - vᵢⁿ⁺¹) \end{align*}\]

where $uₑ$ and $vₑ$ are the external velocities, $uᵢⁿ$ and $vᵢⁿ$ are the sea ice velocities at the current time step, and $uᵢⁿ⁺¹$ and $vᵢⁿ⁺¹$ are the sea ice velocities at the next time step.

Arguments

• FT: The field type of the velocities (optional, default: Oceananigans.defaults.FloatType).

Keyword Arguments

• uₑ: The external x-velocity field.
• vₑ: The external y-velocity field.
• ρₑ: The density of the external fluid.
• Cᴰ: The drag coefficient.

SplitExplicitSolver(grid::AbstractGrid; substeps=120)

Creates a SplitExplicitSolver that controls the dynamical evolution of sea-ice momentum by subcycling substeps times in between each ice_thermodynamics / tracer advection time step.

The default number of substeps is 120.

StressBalanceFreeDrift(; top_momentum_stress = nothing,
                         bottom_momentum_stress = nothing)

A free drift parameterization that computes the free drift velocities as a balance between top and bottom stresses $τₐ ≈ τₒ$.

The only supported configuration is when either the top_momentum_stress or the bottom_momentum_stress are a SemiImplicitStress. The model will compute the free drift velocity exactly assuming that the other stress does not depend on the sea ice velocity.

Can be used to limit the sea ice velocity when the mass or the concentration are below a certain threshold, or as a dynamics model itself that substitutes the sea ice momentum equation calculation everywhere.

time_step_momentum!(model, dynamics::SplitExplicitMomentumEquation, Δt)

function for stepping u and v in the case of explicit solvers. The sea-ice momentum equations are characterized by smaller time-scale than sea-ice ice_thermodynamics and sea-ice tracer advection, therefore explicit rheologies require substepping over a set number of substeps.