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)

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 γ² = ζ * π² * (Δt / mᵢ) / Az is a stability parameter with $Az$ is the area of the grid cell, $mᵢ$ the ice mass, and $Δt$ the time step.

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 value is 1e-10.
• min_relaxation_parameter: Minimum value for the relaxation parameter α. Default value is 30.
• max_relaxation_parameter: Maximum value for the relaxation parameter α. Default value is 500.

a simple explicit solver

SeaIceMomentumEquation(grid; 
                       coriolis=nothing,
                       rheology=ElastoViscoPlasticRheology(eltype(grid)),
                       auxiliary_fields=NamedTuple(),
                       ocean_velocities=nothing,
                       solver=ExplicitSolver(),
                       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().
• ocean_velocities: The ocean surface velocities used to limit the 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.

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.