Snow mass balance

compute_meltwater_outflow(
    hydraulics::ConstantSnowHydraulics,
    θ_liq::NF,
    W_snow::NF
) where NF

Meltwater outflow M_r (m/s, SWE) from a snowpack with snow water equivalent W_snow (m). Liquid water in excess of the capillary retention L_c drains from the snowpack with a cubic conductivity [5] (in excess-saturation form), K_sat · S*³ with S* = max(θ_liq − L_c, 0) / (1 − L_c), where θ_liq is the liquid fraction of the water substance. This Darcy rate does not depend on the amount of water in the pack, so it is additionally limited by the excess liquid water depth max(θ_liq − L_c, 0) · W_snow drained over the timescale τ_m:

M_r = min(K_sat · S*³, max(θ_liq − L_c, 0) · max(W_snow, 0) / τ_m)

Outflow vanishes smoothly as θ_liq → L_c or W_snow → 0, so a thin pack cannot drain more water than it stores; explicit time steps should not exceed τ_m.

References

• [5] Tarboton, Chowdhury and Jackson (1994)

compute_meltwater_outflow_derivative(
    hydraulics::ConstantSnowHydraulics,
    θ_liq::NF,
    W_ice::NF
) where NF

Derivative dM_r/dθ_liq (m/s) of the meltwater outflow (compute_meltwater_outflow) with respect to the liquid fraction θ_liq at a fixed ice content W_ice (m SWE), i.e. with the snow water equivalent W_snow = W_ice / (1 − θ_liq). It is the derivative of the active branch of the outflow: 3 K_sat S*² / (1 − L_c) where the Darcy rate is the smaller one, and W_ice (1 − L_c) / (τ_m (1 − θ_liq)²) where the storage limit is; zero below the capillary retention L_c. Requires θ_liq < 1.

snow_meltwater_flux(
    i,
    j,
    grid,
    fields,
    snow::SingleLayerSnow
)

Snow meltwater outflow M_r (m/s SWE) at grid cell i, j, the meltwater_flux auxiliary: the outflow at the current liquid fraction under explicit integration (see compute_meltwater_flux!), and the drainage of the last implicit step otherwise.

compute_meltwater_flux!(
    out,
    i,
    j,
    grid,
    fields,
    snow::SingleLayerSnow,
    ::Explicit
)

Store the storage-limited Darcy-type meltwater outflow M_r (m/s SWE) at grid cell i, j in meltwater_flux, evaluated at the diagnosed liquid water fraction and the snow water equivalent of the pack (see compute_meltwater_outflow). A no-op under implicit integration of snow_water_equivalent, where the implicit step sets it.

implicit_operator(
    snow::SingleLayerSnow,
    state,
    grid,
    ::Val{:snow_water_equivalent}
) -> DiagonalOperator

Allocate the DiagonalOperator of the implicit meltwater drainage of snow_water_equivalent, solved in the liquid fraction snow_liquid_fraction (see implicit_variable).

Drainage removes liquid water at the zero-enthalpy reference (liquid water at 0 °C), so it changes neither the snow energy Ū nor the ice content W_i = (1 − θ) W = −Ū / (ρ_w L_f), clamped to [0, W]. With Ū⁺ and W⁺ = Wⁿ + Δt G_exp from the explicit stage (snowfall, rain, and sublimation), the backward Euler step of the drainage M_r(θ, W) (compute_meltwater_outflow) is a scalar equation per column in the liquid fraction θ, via W = W_i / (1 − θ):

G(θ) = W_i / (1 − θ) + Δt M_r(θ, W_i / (1 − θ)) − W⁺ = 0.

Each Picard iteration takes one Newton step: linearizing G about the iterate θᵐ, G(θ) ≈ G(θᵐ) + G′(θᵐ) (θ − θᵐ), and setting this to zero gives the diagonal system G′(θᵐ) θᵐ⁺¹ = G′(θᵐ) θᵐ − G(θᵐ) (see the snow mass balance docs for the convergence argument). Each iteration calls the hooks:

1. linearize!: the diagonal G′(θᵐ) and right-hand side G′(θᵐ) θᵐ − G(θᵐ), at the liquid fraction θᵐ of the current iterate Wᵐ, bounded to [0, θ⁺];
2. solve!: θᵐ⁺¹ in snow_liquid_fraction;
3. implicit_update!: Wᵐ⁺¹ = W_i / (1 − θᵐ⁺¹) with min(L_c, θ⁺) ≤ θᵐ⁺¹ ≤ θ⁺, since the iteration starts from Wⁿ, which may lie below the root, and meltwater_flux ← (W⁺ − Wᵐ⁺¹) / Δt, the drainage passed on to the surface hydrology in the next step;
4. implicit_closure!: the snow closure at the new iterate.

A fully liquid pack (W_i = 0, θ = 1) drains at min(K_sat, (1 − L_c) W / τ_m), so its step has the exact solution max(W⁺ − Δt K_sat, W⁺ / (1 + Δt (1 − L_c) / τ_m)), independent of the iterate.

linearize!(
    operator::DiagonalOperator,
    rhs,
    state,
    grid,
    snow::SingleLayerSnow,
    constants::PhysicalConstants,
    explicit_water,
    Δt,
    ::Terrarium.Iteration{:snow_water_equivalent}
)

Fill the diagonal and right-hand side of the Newton step for the liquid fraction of the meltwater drainage (compute_implicit_meltwater_coefficients!), given the explicit stage W⁺ = explicit_water.

implicit_update!(
    state,
    grid,
    snow::SingleLayerSnow,
    constants::PhysicalConstants,
    operator::DiagonalOperator,
    explicit_water,
    Δt,
    ::Terrarium.Iteration{:snow_water_equivalent}
)

Update snow_water_equivalent from the solved liquid fraction and store the drainage of the step in meltwater_flux (compute_implicit_snow_water_update!).

compute_snow_ice_content(
    i,
    j,
    grid,
    fields,
    explicit_water,
    constants::PhysicalConstants
)

Ice content W_i = −Ū⁺ / (ρ_w L_f) (m SWE) of the snowpack at grid cell i, j, clamped to [0, W⁺], from the explicit stage W⁺ = explicit_water and the snow energy Ū⁺; unchanged by the meltwater drainage.

compute_implicit_meltwater_coefficients!(
    out,
    i,
    j,
    grid,
    fields,
    explicit_water,
    snow::SingleLayerSnow{NF},
    constants::PhysicalConstants,
    Δt
) where NF

Diagonal G′(θᵐ) and right-hand side G′(θᵐ) θᵐ − G(θᵐ) at grid cell i, j of the Newton step for the liquid fraction, with G(θ) = W_i / (1 − θ) + Δt M_r(θ, W_i / (1 − θ)) − W⁺, the storage-limited outflow M_r (compute_meltwater_outflow), the ice content W_i (compute_snow_ice_content), and W⁺ = explicit_water. The linearization point is the liquid fraction 1 − W_i / Wᵐ of the current iterate, bounded to [0, θ⁺]. The coefficients are the identity for a pack without ice or snow, which compute_implicit_snow_water_update! treats separately.

compute_implicit_snow_water_update!(
    out,
    i,
    j,
    grid,
    fields,
    explicit_water,
    snow::SingleLayerSnow{NF},
    constants::PhysicalConstants,
    Δt
) where NF

Snow water equivalent W = W_i / (1 − θ) at grid cell i, j from the solved liquid fraction θ = snow_liquid_fraction, bounded to [min(L_c, θ⁺), θ⁺] with the capillary retention L_c, so that W ≤ W⁺ = explicit_water and a Newton step that undershoots the root at the switch between the branches of the outflow drains neither ice nor capillary water, and the drainage meltwater_flux = (W⁺ − W) / Δt. A pack without ice (W_i = 0) takes the exact solution W = max(W⁺ − Δt K_sat, W⁺ / (1 + Δt (1 − L_c) / τ_m)) of its step, and a step without snow (W⁺ ≤ 0) leaves W = W⁺.

compute_snow_water_tendency(
    i,
    j,
    grid,
    fields,
    snow::SingleLayerSnow{NF},
    atmos::Terrarium.AbstractAtmosphere
) where NF
compute_snow_water_tendency(
    i,
    j,
    grid,
    fields,
    snow::SingleLayerSnow{NF},
    atmos::Terrarium.AbstractAtmosphere,
    water_timestepping::Timestepping
) where NF

Snow water equivalent (SWE) tendency (m/s) at grid cell i, j:

dW_snow/dt = S + R_snow − M − E_subl

where S is snowfall, R_snow = f_snow · rainfall the rain intercepted by the snow-covered fraction, M the Darcy meltwater outflow (see snow_meltwater_flux), and E_subl the sublimation rate. M is omitted under implicit integration of snow_water_equivalent, where the implicit step drains the meltwater.

compute_snow_tendencies!(
    tendencies,
    i,
    j,
    grid,
    fields,
    snow::SingleLayerSnow,
    atmos::Terrarium.AbstractAtmosphere,
    constants::PhysicalConstants
)
compute_snow_tendencies!(
    tendencies,
    i,
    j,
    grid,
    fields,
    snow::SingleLayerSnow,
    atmos::Terrarium.AbstractAtmosphere,
    constants::PhysicalConstants,
    water_timestepping::Timestepping
)

Accumulate the snow water-equivalent and depth-integrated energy tendencies at grid cell i, j from the mass and energy balances (see compute_snow_water_tendency and compute_snow_energy_tendency).