Muninn.jl API

ElevationLUT{F <: AbstractFloat}

Lookup table of positive degree days and solid precipitation as a function of the elevation offset ΔS = S - ref_hgt, one column per mass balance window.

Both quantities are exactly piecewise linear in ΔS:

PDD(ΔS)  = Σ_d max(0, T_d + g_d·ΔS)
snow(ΔS) = Σ_d p_d·clamp((2 - T_d - g_d·ΔS)/2, 0, 1)

so linear interpolation on a uniform table is exact except within one cell of a breakpoint. Monthly mass balance error from the table is 7.1e-5 m at dΔS = 1 m; the lookup costs about an eighth of one SIA2D! call, versus ten times that for the daily loop it replaces.

Fields

• ΔS_min::F, ΔS_max::F: Elevation offset range covered, in m.
• dΔS::F, inv_dΔS::F: Table resolution in m, and its reciprocal.
• PDD::Matrix{F}: Positive degree days, (n_elevations, n_windows).
• snow::Matrix{F}: Solid precipitation in mm, same shape.

MBcache{F <: AbstractFloat}

Everything a mass balance model needs to be evaluated as a source term inside the ice flow RHS, precomputed once per glacier at init_cache time.

The RHS is called many times per mass balance window and must not touch Rasters, so the climate is sliced up front into ClimateWindows and, for models whose rate depends on the surface only through ΔS, collapsed further into an ElevationLUT.

An empty cache (no windows, empty table) is built when mass balance is off, keeping ModelCache concretely typed either way — a run without mass balance pays neither the memory nor the build time. Use mb_cache_active to tell them apart.

Fields

• windows::Vector{ClimateWindow{F}}: Climate per mass balance window, in time order.
• lut::ElevationLUT{F}: Elevation lookup table, empty for :general models.
• ṁ::Matrix{F}: Buffer holding the most recently computed rate, in m/yr.
• t₀::F, step_MB::F: Start of the run and the mass balance step, in decimal years.
• ref_hgt::F: Reference elevation of the climate data, in m.
• temp_bias::F: Temperature bias baked into lut, in °C.
• H_acc::F, H_acc_w::F: Ice thickness at which accumulation switches on, and the width of that switch, in m. The ramp is centred on H_acc: it is zero below H_acc - H_acc_w/2 and full above H_acc + H_acc_w/2.
• H_abl::F: Ramp width for ablation, in m, rising from zero at a bare bed.

MBmodel <: AbstractModel

An abstract type representing a mass balance model in the Muninn package. This serves as a base type for all specific mass balance models, ensuring they adhere to a common interface and can be used interchangeably within the ODINN framework.

TImodel <: MBmodel

An abstract type representing a temperature index mass balance models within the ODINN framework. This type serves as a parent type for more specialized mass balance models, ensuring they adhere to a common interface defined by the MBmodel abstract type.

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.

TImodel2{F <: AbstractFloat}

A type representing a temperature-index model with parameters for snow and ice degree-day factors.

Keyword arguments

• DDF_snow::F: Degree-day factor for snow in m w.e. °C⁻¹ d⁻¹, which determines the melt rate of snow per degree above the melting point.
• DDF_ice::F: Degree-day factor for ice in m w.e. °C⁻¹ d⁻¹, which determines the melt rate of ice per degree above the melting point.

Type Parameters

• F: A subtype of AbstractFloat, representing the numeric type used for the model parameters.

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.

TImodel2(params::Sleipnir.Parameters; DDF_snow::F = 3.0/1000.0, DDF_ice::F = 6.0/1000.0) where {F <: AbstractFloat}

Create a temperature-index model with two degree-day factors (TImodel2) for mass balance calculations.

Arguments

• params::Sleipnir.Parameters: The parameters object containing simulation settings.
• DDF_snow::F: Degree-day factor for snow in m w.e. °C⁻¹ d⁻¹ (default: 3.0/1000.0 = 0.003).
• DDF_ice::F: Degree-day factor for ice in m w.e. °C⁻¹ d⁻¹ (default: 6.0/1000.0 = 0.006).

Returns

• TI2_model: An instance of the TImodel2 with the specified parameters.

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

MB_rate!(ṁ, H, cache::MBcache, mb_model::MBmodel, glacier, t)

Fallback for mass balance models that have no RHS form yet.

MB_rate!(ṁ, H, cache::MBcache, mb_model::TImodel1, glacier, t)

Write the instantaneous mass balance rate ṁ(H, t), in m/yr, into ṁ.

Rate form of compute_MB: that returns the balance accumulated over a window, for calibration and diagnostics; this returns the rate at which it accrues, to add to ∂H/∂t. Integrating the rate over a window with H frozen reproduces compute_MB exactly; letting H evolve within the window is the operator-splitting error this avoids.

ṁ depends on H through S = B + H, so ∂ṁ/∂H (the elevation feedback) reaches the gradient with respect to ice flow parameters like A even when no mass balance parameter is trained.

DDF and prcp_fac come from mb_model on every call, so they can change without invalidating the lookup table; temp_bias is baked into the table and checked against the cache instead.

ṁ is overwritten, not accumulated into — the caller combines it with the rest of the RHS.

MB_rate_∂H!(∂ṁ, H, cache::MBcache, mb_model, glacier, t)

Fill ∂ṁ with the derivative of the mass balance rate with respect to the ice thickness, in-place.

ṁ at a cell depends only on that cell's H, through S = B + H and through the ramp, so the Jacobian is diagonal and a vector-Jacobian product is just an elementwise multiply.

From ṁ = rate(ΔS) · ramp(H) and ∂ΔS/∂H = 1, the product rule gives

\[∂ṁ/∂H = \frac{∂rate}{∂ΔS} ramp + rate \frac{∂ramp}{∂H}\]

Both terms vanish below the ice margin, so the derivative is continuous there. Only manual adjoints need this; the automatic sensitivity path differentiates MB_rate! directly.

MB_timestep!(cache, model::Model, glacier::G, step::F, t) where {F <: AbstractFloat, G <: AbstractGlacier}

Simulates a mass balance timestep for a given glacier model.

Arguments

• cache: The model cache to update.
• model::Model: The glacier model containing iceflow and mass balance information.
• glacier::G: The glacier object containing climate and DEM data.
• step::F: The timestep duration.
• t: The current time.

Description

This function performs the following steps:

1. Computes the period from the previous timestep to the current time.
2. Retrieves cumulative climate data for the specified period.
3. Downscales the climate dataset to a 2D grid based on the glacier's DEM.
4. Computes the mass balance for the glacier and updates the model's iceflow mass balance.

MB_timestep(model::Model, glacier::G, step::F, t::F) where {F <: AbstractFloat, G <: AbstractGlacier}

Calculate the mass balance (MB) for a glacier over a given timestep.

Keyword arguments

• model::Model: The model containing mass balance parameters.
• glacier::G: The glacier object containing climate and DEM data.
• step::F: The timestep duration.
• t::F: The current time.

Returns

• MB::Matrix{F}: The computed mass balance matrix for the given timestep.

Details

1. Computes the period between the current time t and the previous step t - step.
2. Retrieves cumulative climate data for the specified period.
3. Downscales the climate data to a 2D grid based on the glacier's DEM.
4. Computes the mass balance using the downscaled climate data.

build_elevation_lut(windows, ΔS_min, ΔS_max, temp_bias; dΔS = LUT_dΔS_DEFAULT)

Build an ElevationLUT covering [ΔS_min, ΔS_max] for the given climate windows.

temp_bias is baked into the table, so a table is only valid for the mass balance model it was built from. DDF and prcp_fac are deliberately not baked in: they scale the melt and accumulation terms linearly and are applied at lookup time, which keeps the table valid if they change.

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.

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

compute_MB(
    mb_model::TImodel1,
    climate_2D_period::Climate2Dstep,
    step::AbstractFloat,
)

Compute the mass balance (MB) for a given mass balance model and climate period.

Arguments

• mb_model::TImodel1: The mass balance model containing parameters such as accumulation factor (acc_factor) and degree-day factor (DDF).
• climate_2D_period::Climate2Dstep: The climate data for a specific period, including snow accumulation (snow) and positive degree days (PDD).
• step::AbstractFloat: The step used to update MB. This scales the MB so that the accumulation and degree-day factors are scaled monthly.

Returns

• A numerical array representing the computed mass balance, calculated as the difference between the product of the accumulation factor and snow, and the product of the degree-day factor and positive degree days.

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.

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.

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.

lut_lookup(lut::ElevationLUT, ΔS, k::Int)

Positive degree days and solid precipitation at elevation offset ΔS in window k, by linear interpolation in the table. Calls lut_lookup_∂ and discards the derivatives — the interpolation already computes them as intermediates, so this costs nothing extra.

lut_lookup_∂(lut::ElevationLUT, ΔS, k)

Positive degree days and solid precipitation at elevation offset ΔS in window k, by linear interpolation in the table, together with their derivatives with respect to ΔS. The table is piecewise linear, so the slope of the bracketing segment is the exact derivative, not an approximation of it.

Throws a DomainError, rather than clamping, if ΔS is outside the tabulated range or not finite — clamping would quietly return a wrong mass balance.

lut_∂ṁ_bound(lut::ElevationLUT, mb_model, step_MB, H_abl, H_acc_w)

Upper bound on |∂ṁ/∂H| over the whole table, in yr⁻¹.

This is the mass balance contribution to the spectral radius of the ice flow right hand side. It's a state-independent bound, not the exact maximum at the current state: a stabilised solver only needs an upper bound to size its stages, and the adjoint hands the RHS an augmented [H; θ] vector that a state-dependent version couldn't read anyway.

From ṁ = rate(ΔS)·ramp(H), with ramp ≤ 1 and ramp' ≤ 1.5/w for a ramp of width w, taking the narrower of the accumulation and ablation ramps:

\[|∂ṁ/∂H| ≤ \max|∂rate/∂ΔS| + \max|rate| · 1.5 / \min(H_{abl}, H_{acc,w})\]

Both maxima come straight from the table, whose rows are the tabulated elevations and whose slopes between them are exact.

mb_S_dependence(mb_model::MBmodel)

How the mass balance rate of mb_model depends on the ice surface S, which decides how it can be evaluated inside the ice flow RHS.

• :elevation_only: the rate at a cell depends on S only through the scalar offset ΔS = S - ref_hgt. The whole model then collapses to a one-dimensional function of ΔS per time window, which is what makes a lookup table possible.
• :general: the rate depends on S in a way that does not reduce to ΔS (a neural network reading several fields, for instance) and must be evaluated directly.

:general is the conservative default; models opt in to :elevation_only.

mb_cache_active(cache)

Whether cache actually carries precomputed climate, i.e. whether mass balance is being evaluated inside the ice flow RHS for this simulation. nothing, the cache of a mass balance model that has no RHS form, is never active.

pdd_snow_exact(window::ClimateWindow, ΔS, temp_bias)

Exact daily sums of positive degree days and solid precipitation at elevation offset ΔS.

This is the reference the lookup table approximates, and the routine that fills it. It is the same daily temperature-index computation downscale_2D_climate! performs, evaluated at a single elevation instead of over a grid.

smoothstep(x)

Cubic Hermite ramp x²(3 - 2x), clamped to [0, 1].

Fades the mass balance rate in and out with ice thickness. A hard mask (as apply_MB_mask! applies) is a step discontinuity in H, which as a source term is unresolvable by an adaptive controller and undifferentiable for an adjoint. This ramp is C¹: its derivative 6x(1 - x) is zero at both ends, so it joins the constant parts smoothly. It is not C², since the second derivative jumps at the ends. Its derivative is at most 1.5 over its width.

smoothstep_∂(x)

Derivative of smoothstep, zero outside [0, 1].

window_index(cache::MBcache, t)

Index of the mass balance window containing time t.

Window k spans (t₀ + (k-1)·step, t₀ + k·step]. Times at or before t₀ map to the first window and times past the end map to the last.

init_mb_cache(mb_model::TImodel1, simulation, glacier_idx, θ)

Build the MBcache for one glacier. Returns an empty cache when mass balance is off.

The lookup table spans the elevations the glacier can reach — lowest bed to highest initial surface — padded either side, to absorb thickening beyond the initial state or a bed emerging below the current minimum. Going outside that range throws rather than extrapolating.

plot_cumulative_mb(
    mb_model::MBmodel,
    glacier::AbstractGlacier,
    t_start::F,
    t_end::F;
    step::F = F(1.0 / 12.0),
    kwargs...,
) where {F <: AbstractFloat}

Plot the mean annual mass balance map (m w.e. yr⁻¹) for a single glacier over [t_start, t_end] using a calibrated mass balance model and static geometry.

Extends Sleipnir.plot_cumulative_mb with a dispatch that accepts an MBmodel and a glacier directly. The cumulative MB is computed via compute_cumulative_MB and divided by (t_end - t_start) to obtain an annual rate before plotting. All kwargs are forwarded (e.g. title, colormap).