Public Documentation

Ocean 🌊 Sea ice component of CliMa's Earth system model.

PhaseTransitions(FT=Float64,
                 ice_density           = 917,   # kg m⁻³
                 ice_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⁻³
                 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.

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 ice_density, $ρ_ℓ$ is the liquid density, $cᵢ$ is the heat capacity of ice, and $c_ℓ$ is the heat capacity of liquid, 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.

SlabSeaIceThermodynamics(grid; kw...)

Pretty simple model for sea ice.

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

Return FluxFunction representing a flux across an air-ice, air-snow, or ice-water interface. The flux is computed by func with the signature

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

if isnothing(parameters), or

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

if !isnothing(parameters). If func is top_temperature_dependent, then it will be recomputed during a diagnostic solve for the top temperature.

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=Float64; kw...)

Returns a flux representing radiative emission from a surface.

ElastoViscoPlasticRheology(FT::DataType = Float64; 
                           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())

Constructs an ElastoViscoPlasticRheology object representing a "modified" elasto-visco-plastic rheology for slab sea ice dynamics that follows the implementation of Kimmritz et al (2016). 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★ h exp( - C ⋅ ( 1 - ℵ ))$ where $P★$ 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 in Kimmritz et al (2016) In particular: α = sqrt(γ²) where γ² = ζ * cα * (Δt / mᵢ) / Az is a stability parameter with $Az$ is the area of the grid cell, $mᵢ$ the ice mass, $Δt$ the time step, and $cα$ a numerical stability parameter which controls the stregth 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 thougth of as a $pseudo substep number'' or a$relaxation parameter''. If we are using a subcycling solver, if α ≪ number of substeps, the convergence will be faster.

Arguments

• grid: the SlabSeaIceModel grid

Keyword Arguments

• ice_compressive_strength: parameter expressing compressive strength (in Nm²). Default: 27500.
• ice_compaction_hardening: exponent coefficient for compaction hardening. Default: 20.
• yield_curve_eccentricity: eccentricity of the elliptic yield curve. Default: 2.
• Δ_min: Minimum value for the visco-plastic parameter. Limits the maximum viscosity of the ice, transitioning the ice from a plastic to a viscous behaviour. Default: 1e-10.
• min_relaxation_parameter: Minimum value for the relaxation parameter α. Default: 30.
• max_relaxation_parameter: Maximum value for the relaxation parameter α. Default: 500.
• relaxation_strength: parameter controlling the strength of the relaxation parameter. The maximum value is π², see Kimmritz et al (2016). Default: π² / 2.
• pressure_formulation: can use ReplacementPressure or IceStrength. The replacement pressure formulation avoids ice motion in the absence of forcing. Default: ReplacementPressure.

a simple explicit solver

SeaIceMomentumEquation(grid; 
                       coriolis = nothing,
                       rheology = ElastoViscoPlasticRheology(eltype(grid)),
                       auxiliary_fields = NamedTuple(),
                       top_momentum_stress    = nothing,
                       bottom_momentum_stress = nothing,
                       free_drift = nothing,
                       solver = SplitExplicitSolver(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:

\[ ∂u τₒ τₐ -- + f x u = ∇ ⋅ σ + -- + -- ∂t mᵢ mᵢ\]

where the terms (left to right) represent (1) the time derivative of the ice velocity, (2) the coriolis force. (3) the divergence of internal stresses, (4) the ice-ocean boundary stress, and (5) the ice-atmosphere boundary stress.

Arguments

• grid: The computational grid.

Keyword Arguments

• coriolis: Parameters for the background rotation rate of the model.
• rheology: The sea ice rheology model, default is ElastoViscoPlasticRheology(eltype(grid)).
• auxiliary_fields: A named tuple of auxiliary fields, default is an empty NamedTuple().
• free_drift: The free drift velocities used to limit sea ice momentum when the mass or the concentration are below a certain threshold. Default is nothing (indicating that the free drift velocities are zero).
• solver: The momentum solver to be used.
• minimum_concentration: The minimum sea ice concentration above which the sea ice velocity is dynamically calculated, default is 1e-3.
• minimum_mass: The minimum sea ice mass per area above which the sea ice velocity is dynamically calculated, default is 1.0 kg/m².

SemiImplicitStress(FT = Float64; 
                   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

\[τᵤ = ρₑ Cᴰ sqrt((uₑ - uᵢⁿ)² + (vₑ - vᵢⁿ)²) (uₑ - uᵢⁿ⁺¹)\]

\[τᵥ = ρₑ Cᴰ sqrt((uₑ - uᵢⁿ)² + (vₑ - vᵢⁿ)²) (vₑ - vᵢⁿ⁺¹)\]

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: Float64).

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(; 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{T, B}

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

The only supported configuration is when either the top_momentum_stess 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, rheology::AbstractExplicitRheology, Δ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.