Models

Model{IFM <: AbstractEmptyModel, MBM <: Union{<:AbstractEmptyModel, Vector{<:AbstractEmptyModel}}, TC <: AbstractEmptyModel}

A mutable struct that represents a model with three components: iceflow, mass balance, and machine learning.

Model(
    iceflow::IFM,
    mass_balance::MBM,
    trainable_components::TC,
) where {IFM <: AbstractEmptyModel, MBM <: Union{<:AbstractEmptyModel, Vector{<:AbstractEmptyModel}}, TC <: AbstractEmptyModel}

Model(;iceflow, mass_balance) = Model(iceflow, mass_balance, nothing)

Initialize Model (no machine learning model).

Keyword arguments

• iceflow::IFM: Represents the iceflow component, which is an instance of IFM.
• mass_balance::MBM: Represents the mass balance component. Either a single MBM instance or a Vector{MBM} of per-glacier models.
• trainable_components::TC: Represents the trainable components, which is an instance of TC.

Type Parameters

• IFM: A subtype of AbstractEmptyModel representing the type of the iceflow model.
• MBM: Either a subtype of AbstractEmptyModel (single shared model) or a Vector of such subtypes (per-glacier models). The field stores the exact runtime type — no Union in the field — which is required for AD compatibility.
• TC: A subtype of AbstractEmptyModel representing the type of the trainable components.

Model(;
    iceflow::Union{<: AbstractModel, Nothing} = nothing,
    mass_balance::Union{<: AbstractModel, Nothing} = nothing,
    regressors::Union{NamedTuple, Nothing} = nothing
)

Creates a new model instance using the provided iceflow, mass balance, and machine learning components.

Arguments

• iceflow::Union{IFM, Nothing}: The iceflow model to be used. Can be a single model or nothing.
• mass_balance::Union{MBM, Nothing}: The mass balance model to be used. Can be a single model or nothing.
• regressors::Union{NamedTuple, Nothing}: The regressors to be used in the laws.

Returns

• model: A new instance of Sleipnir.Model initialized with the provided components.

Note

Since Model is a keyword argument function, it is defined in Sleipnir. Otherwise we would have to define it multiple times, making the ODINN precompilation impossible because of method overwriting.

SIA2Dmodel(A, C, n, Y, U, n_H, n_∇S)
SIA2Dmodel(params; A, C, n, Y, U, n_H, n_∇S)

Create a SIA2Dmodel, representing a two-dimensional Shallow Ice Approximation (SIA) model.

The SIA model describes glacier flow under the assumption that deformation and basal sliding dominate the ice dynamics. It relies on:

• Glen's flow law for internal deformation, with flow rate factor A and exponent n,

• A sliding law governed by coefficient C,

• Optionally the user can provide either:

  • A specific diffusive velocity U such that D = U * H
  • A modified creep coefficient Y that takes into account the ice thickness such that D = (C + Y * 2/(n+2)) * (ρ*g)^n * H^(n_H+1) * |∇S|^(n_∇S-1) where n_H and n_∇S are optional parameters that control if the SIA should use the n law or not. This formulation is denoted as the hybrid diffusivity in the code.

This struct stores the laws used to compute these three parameters during a simulation. If not provided, default constant laws are used based on glacier-specific values.

Arguments

• A: Law for the flow rate factor. Defaults to a constant value from the glacier.
• C: Law for the sliding coefficient. Defaults similarly.
• n: Law for the flow law exponent. Defaults similarly.
• p: Law for the sliding law exponent (basal drag). Defaults similarly.
• q: Law for the sliding law exponent (normal stress). Defaults similarly.
• Y: Law for the hybrid diffusivity. Providing a law for Y discards the laws of A, C and n.
• U: Law for the diffusive velocity. Defaults behavior is to disable it and in such a case it is computed from A, C and n. Providing a law for U discards the laws of A, C, n and Y.
• n_H::F: The exponent to use for H in the SIA equation when using the Y law (hybrid diffusivity). It should be nothing when this law is not used.
• n_∇S::F: The exponent to use for ∇S in the SIA equation when using the Y law (hybrid diffusivity). It should be nothing when this law is not used.
• Y_is_provided::Bool: Whether the diffusivity is provided by the user through the hybrid diffusivity Y or it has to be computed from the SIA formula from A, C and n.
• U_is_provided::Bool: Whether the diffusivity is provided by the user through the diffusive velocity U or it has to be computed from the SIA formula from A, C and n.
• n_H_is_provided::Bool: Whether the H exponent is prescribed by the user, or the one of the n law has to be used. This flag is used only when a law for Y is used.
• n_∇S_is_provided::Bool: Whether the ∇S exponent is prescribed by the user, or the one of the n law has to be used. This flag is used only when a law for Y is used.
• apply_A_in_SIA::Bool: Whether the value of the A law should be computed each time the SIA is evaluated.
• apply_C_in_SIA::Bool: Whether the value of the C law should be computed each time the SIA is evaluated.
• apply_n_in_SIA::Bool: Whether the value of the n law should be computed each time the SIA is evaluated.
• apply_p_in_SIA::Bool: Whether the value of the p law should be computed each time the SIA is evaluated.
• apply_q_in_SIA::Bool: Whether the value of the q law should be computed each time the SIA is evaluated.
• apply_Y_in_SIA::Bool: Whether the value of the Y law should be computed each time the SIA is evaluated.
• apply_U_in_SIA::Bool: Whether the value of the U law should be computed each time the SIA is evaluated.

SIA2D!(
    dH::Matrix{R},
    H::Matrix{R},
    simulation::SIM,
    t::R,
    θ,
) where {R <:Real, SIM <: Simulation}

Simulates the evolution of ice thickness in a 2D shallow ice approximation (SIA) model. Works in-place.

Arguments

• dH::Matrix{R}: Matrix to store the rate of change of ice thickness.
• H::Matrix{R}: Matrix representing the ice thickness.
• simulation::SIM: Simulation object containing model parameters and state.
• t::R: Current simulation time.
• θ: Parameters of the laws to be used in the SIA. Can be nothing when no learnable laws are used.

Details

This function updates the ice thickness H and computes the rate of change dH using the shallow ice approximation in 2D. It retrieves necessary parameters from the simulation object, enforces positive ice thickness values, updates glacier surface altimetry and computes surface gradients. It then applies the necessary laws that are not updated via callbacks (A, C, n or U depending on the use-case) and computes the flux components, and flux divergence.

Notes

• The function operates on a staggered grid for computing gradients and fluxes.
• Surface elevation differences are capped using upstream ice thickness to impose boundary conditions.
• The function modifies the input matrices dH and H in-place.

See also SIA2D

TImodel1{F <: AbstractFloat}

A structure representing a temperature index model with degree-day factor and precipitation correction factor.

Keyword arguments

• DDF::F: Degree-day factor (m w.e. °C⁻¹ d⁻¹), which converts positive degree days to melt.
• prcp_fac::F: Dimensionless precipitation correction factor applied as a multiplier to the snowfall field before computing accumulation. A value of 1.0 (default) leaves precipitation unchanged; values > 1 increase accumulation (useful when the climate input underestimates solid precipitation), values < 1 reduce it.

Type Parameters

• temp_bias::F: Uniform temperature bias (°C) added to the glacier's climate before computing melt and snow/rain partitioning. 0.0 (default) leaves the climate unchanged. Used as the third calibration lever when DDF and prcp_fac alone cannot bracket the geodetic mass-balance target.

• F: A subtype of AbstractFloat representing the type of the factors.

Note: The unit conversion from mm to m w.e. is handled internally via the constant PRECIP_UNIT_CONVERSION (1/1000) and is not a tunable parameter.

TImodel1(params::Sleipnir.Parameters; DDF::F = 7.0/1000.0, prcp_fac::F = 1.0, temp_bias::F = 0.0) where {F <: AbstractFloat}

Create a temperature index model with one degree-day factor (DDF) with the given parameters.

Arguments

• params::Sleipnir.Parameters: The simulation parameters.
• DDF::F: Degree-day factor in m w.e. °C⁻¹ d⁻¹ (default is 7.0/1000.0 = 0.007).
• prcp_fac::F: Dimensionless precipitation correction factor (default is 1.0).
• temp_bias::F: Uniform temperature bias in °C (default is 0.0).

Returns

• TI1_model: An instance of TImodel1 with the specified parameters.

Note: Precipitation unit conversion (mm → m w.e.) is handled internally via PRECIP_UNIT_CONVERSION.

calibrate_MB_model(
    model::Sleipnir.Model,
    glaciers::Vector{<:AbstractGlacier},
    params::Parameters,
) -> Sleipnir.Model

High-level entry point that calibrates the mass balance model of model per glacier against geodetic observations, returning a new model.

Calibration cannot happen in place: for TImodel1 it replaces a single model by one model per glacier, which changes the type of the mass_balance field. Callers must use the return value; discarding it silently keeps the uncalibrated model.

The behaviour dispatches on the mass balance model type:

• TImodel1: returns a new Model whose mass_balance is a per-glacier vector of TImodel1s, each fitted against its glacier's geodetic mass balance with calibrate_ti_model. Glaciers without dhdtData keep the original (uncalibrated) model and a warning is emitted. If no glacier carries geodetic data, the model is returned unchanged.
• any other MBmodel: no-op — no calibration routine is defined, so the model is returned untouched. Add a _calibrate_MB_model(model, ::YourType, …) method to support a new model type.

An already-vectorized (per-glacier) mass balance model is left unchanged. Any keyword arguments are forwarded to the type-specific calibrator (e.g. calibrate_ti_model for TImodel1). This is the function the Prediction and Inversion constructors call when params.simulation.calibrate_MB is true.

calibrate_ti_model(
    glacier::AbstractGlacier,
    params::Parameters;
    DDF_bounds::Tuple{Sleipnir.Float, Sleipnir.Float} = (
        Sleipnir.Float(params.physical.DDF_min),
        Sleipnir.Float(params.physical.DDF_max)),
    prcp_fac_bounds::Tuple{Sleipnir.Float, Sleipnir.Float} = (
        Sleipnir.Float(params.physical.prcp_fac_min),
        Sleipnir.Float(params.physical.prcp_fac_max)),
    temp_bias_bounds::Tuple{Sleipnir.Float, Sleipnir.Float} = (
        Sleipnir.Float(params.physical.temp_bias_min),
        Sleipnir.Float(params.physical.temp_bias_max)),
    density_ratio::Sleipnir.Float = Sleipnir.Float(1.0),
    calibration_period::Union{Nothing, Tuple{Sleipnir.Float, Sleipnir.Float}} = nothing,
    prcp_fac::Union{Symbol, Real} = :from_winter_prcp,
    step::Sleipnir.Float = Sleipnir.Float(1.0 / 12.0),
) -> TImodel1

Calibrate TImodel1 for a single glacier against the geodetic mass balance stored in glacier.dhdtData (e.g. the 2000-2020 observations from Hugonnet et al. 2021).

The calibration follows a 3-step cascade analogous to the OGGM v1.6 approach, using Brent's method for root-finding at each step:

1. DDF step: With prcp_fac fixed (glacier-specific from winter precipitation by default, see the prcp_fac keyword) and temp_bias = 0.0, find the degree-day factor that matches the geodetic MB target. If the target can be bracketed within DDF_bounds, this step alone produces the calibrated model.

2. prcp_fac step (fallback): If the DDF search cannot bracket the target, DDF is fixed at its boundary value and prcp_fac is varied within prcp_fac_bounds. A warning is emitted when this fallback is used.

3. temp_bias step (fallback): If both previous steps fail, DDF and prcp_fac are fixed at their best boundary values and a uniform temperature bias (°C) is varied within temp_bias_bounds. This handles glaciers where the climate forcing is systematically biased for the glacier's hypsometry (e.g. high-elevation accumulation overestimation in coarse reanalysis data).

A static glacier geometry is assumed throughout (no ice-flow dynamics).

Arguments

• glacier::AbstractGlacier: Glacier whose dhdtData field contains the observed geodetic mass balance (see DhdtData).
• params::Parameters: Simulation parameters (provides physical bounds).

Keyword arguments

• DDF_bounds: Search interval (DDF_min, DDF_max) in m w.e. °C⁻¹ d⁻¹. Default: from params.physical.
• prcp_fac_bounds: Search interval for the precipitation correction factor (dimensionless). Default: from params.physical.
• temp_bias_bounds: Search interval for the temperature bias (°C). Default: from params.physical.
• density_ratio: Conversion factor applied to glacier.dhdtData.dhdt. Use params.physical.ρ / params.physical.ρ_w ≈ 0.9 when the data are in m ice yr⁻¹, or 1.0 (default) for m w.e. yr⁻¹ (Hugonnet et al. 2021).
• calibration_period: Time window (t_start, t_end) in fractional years. Defaults to glacier.dhdtData.t.
• prcp_fac: Precipitation factor used in the DDF step. :from_winter_prcp (default) derives a glacier-specific factor from winter precipitation via Sleipnir.get_winter_prcp_factor; a Real value fixes it (e.g. 2.5 for OGGM's global W5E5 default). Only used in step 1; the fallback steps still search prcp_fac_bounds.
• step: Integration timestep in fractional years. Default: 1/12 (monthly).

Returns

• A TImodel1{Sleipnir.Float} with calibrated DDF, prcp_fac, and temp_bias.

This is the per-glacier building block used by calibrate_MB_model to build a per-glacier vector of calibrated models.

References

Hugonnet, R. et al. (2021). Accelerated global glacier mass loss in the early twenty-first century. Nature, 592, 726–731. https://doi.org/10.1038/s41586-021-03436-z

Maussion, F., Butenko, A., Eis, J., Fourteau, K., Jarosch, A. H., Landmann, J., Oesterle, F., Recinos, B., USias, S., Valsecchi, L., Marzeion, B., and Cogley, J. G. (2019). The Open Global Glacier Model (OGGM) v1.1. Geoscientific Model Development, 12, 909–941. https://doi.org/10.5194/gmd-12-909-2019

get_winter_prcp_factor(glacier::AbstractGlacier, params::Parameters; prcp_fac_bounds) -> Float

Compute a glacier-specific precipitation correction factor from mean winter precipitation, following OGGM's decide_winter_precip_factor.

Arguments

• glacier::AbstractGlacier: Glacier providing rgi_id, cenlat (hemisphere) and climate.climate_data_source (:W5E5 or :ERA5).
• params::Parameters: Simulation parameters, used to locate the climate file.
• prcp_fac_bounds: Lower/upper clip bounds (default (0.1, 10.0), matching OGGM).

Description

Reads the full daily climate record (1979–2019, the window OGGM's coefficients were fit on), averages winter daily precipitation (NH: Oct–Apr, SH: Apr–Oct) in kg/m²/day, and maps it to a precipitation factor: log fit for W5E5, linear fit for ERA5. If the daily file is missing, falls back to 2.5. Note: we average daily values directly (OGGM averages monthly daily-rates; equivalent within noise).

get_mb_model(mass_balance, glacier_idx::Integer = 1)

Return the mass balance model for a given glacier from a Model's mass_balance field.

When mass_balance is a per-glacier vector of models (e.g. one calibrated TImodel1 per glacier), the glacier_idx-th entry is returned. When it is a single shared model, that model is returned regardless of glacier_idx. This is the single accessor used throughout the ecosystem so that callers never need to branch on whether the mass balance model is vectorized.

compute_mean_annual_MB(
    mb_model::TImodel1,
    glacier::AbstractGlacier,
    t_start::F,
    t_end::F;
    step::F = F(1.0 / 12.0),
) where {F <: AbstractFloat} -> Sleipnir.Float

Glacier-wide mean annual mass balance (m w.e. yr⁻¹) over [t_start, t_end] with static geometry. Spatially averages the output of compute_cumulative_MB over the ice-covered area and normalises to an annual rate.

compute_cumulative_MB(
    mb_model::TImodel1,
    glacier::AbstractGlacier,
    t_start::F,
    t_end::F;
    step::F = F(1.0 / 12.0),
) where {F <: AbstractFloat} -> Matrix{Sleipnir.Float}

Accumulate the gridded mass balance (m w.e.) over [t_start, t_end] using a static glacier geometry. Returns the full 2D cumulative MB field, suitable for spatial visualisation via plot_cumulative_mb.

See compute_mean_annual_MB for the glacier-wide scalar summary.

NeuralNetwork{
    ChainType <: Lux.Chain,
    ComponentVectorType <: ComponentVector,
    NamedTupleType <: NamedTuple,
} <: FunctionalModel

Feed-forward neural network.

Fields

• architecture::ChainType: Flux.Chain neural network architecture
• θ::ComponentVectorType: Neural network parameters
• st::NamedTupleType: Neural network status