Snow energy balance

struct SnowEnergyTemperatureClosure{NF} <: Terrarium.AbstractEnergyClosure

Energy–temperature closure for snow volumes. For SingleLayerSnow, the depth-averaged snow temperature T_snow (°C) and liquid water fraction θ_liq are recovered from the depth-integrated internal energy Ū_snow (J/m²) using the medium-agnostic FreeWater enthalpy relations, treating the bulk snowpack as an ice-water-air mixture. The internal energy is defined as:

\[U(T) = T_{ ext{snow}} \times C(T) - \rho_{snow} L_{sl} (1 - F(T))\]

with C(T) the temperature-dependent volumetric heat capacity of the snowpack (J/m³/K), ρ_snow L_sl = ρ_w L_sl θ the volumetric latent heat of fusion (J/m³), and F(T) = θ_liq/θ the fraction of the total (liquid water + ice) volumetric water content that is liquid.

compute_snow_volumetric_energy(Ū_snow, d_snow, d_min) -> Any

Volumetric snow internal energy U_snow = Ū_snow/max(d_snow, d_min) [J/m³] from the depth-integrated energy Ū_snow [J/m²] and the snow depth d_snow [m]. The depth is floored at the minimum thermal thickness d_min (see min_snow_conduction_thickness) rather than a machine-eps offset. This bounds the snow temperature recovered downstream: with only an eps offset, any residual Ū_snow over a vanishing d_snow gives a huge U_snow and hence a snow temperature far below physical bounds, which then corrupts the (cover-fraction-blended) skin temperature and basal heat flux. Flooring at d_min gives a thin snowpack a bounded effective heat capacity C_snow·d_min, so it stores negligible energy and is thermally transient — the ground heat flux passes essentially unmediated to the soil (f_snow → 0 in the blend) while the pack still stores mass and modifies the surface albedo.

compute_snow_volumetric_heat_capacity(
    _::Terrarium.AbstractSnow{NF},
    constants::PhysicalConstants,
    ρ_snow,
    liq
) -> Any

Bulk volumetric heat capacity C_snow [J/m³/K] of the snowpack treated as an ice–liquid–air mixture, given the bulk snow density ρ_snow [kg/m³] and the liquid water fraction liq ∈ [0,1] of the water substance. The water-substance mass per unit snow volume is ρ_snow, of which (1 − liq) is ice and liq is liquid; the remaining void space is dry air. The corresponding constituent volume fractions are

\[\begin{aligned} θ_{ice} &= ρ_{snow}·(1 − liq)/ρ_i,\ θ_{liq} &= ρ_{snow}·liq/ρ_w,\ θ_{air} &= 1 − θ_{ice} − θ_{liq} \end{aligned}\]

and the heat capacity is the volume-weighted sum over the ice and liquid constituents, C_snow = cp_i·ρ_i·θ_ice + cp_w·ρ_w·θ_liq. The air's own sensible-heat storage (≈0.1%) is neglected; note the air is still reflected in the bulk density ρ_snow < ρ_ice.

compute_snow_precip_heat_flux(
    _::Terrarium.AbstractSnow,
    constants::PhysicalConstants,
    P_s,
    R_on_snow,
    T_air
) -> Any

Advected heat flux [W/m²] carried into the snowpack by precipitation, relative to liquid water at 0 °C (the U = 0 reference of the FreeWater enthalpy closure). Fresh snow P_s arrives as ice, which sits L_sl below the liquid reference, plus sensible heat for T_air < 0; rain-on-snow R_on_snow arrives as liquid carrying only its sensible heat for T_air > 0. The latent heat released when rain refreezes in a cold pack is captured implicitly by the enthalpy closure, so it is not added here (adding L_sl would double-count relative to the liquid-water reference).

compute_snow_energy_tendency(
    i,
    j,
    grid,
    fields,
    snow::SingleLayerSnow,
    atmos::Terrarium.AbstractAtmosphere,
    constants::PhysicalConstants
) -> Any

Depth-integrated snow energy tendency [W/m²] at grid cell i, j (all fluxes positive upward):

dŪ_snow/dt = Q_base − Q_top + Q_precip + Q_subl

where Q_top/Q_base are the surface/basal heat fluxes, Q_precip the advected precipitation heat (see compute_snow_precip_heat_flux), and Q_subl an advective correction for sublimation.

The sublimation correction Q_subl = ρ_w·L_sl·E_subl is required because the latent heat flux carries the full sublimation enthalpy ρ_w·L_sg·E_subl, whereas the mass leaving the snowpack departs as ice, whose specific enthalpy relative to the liquid-water reference is −L_sl. Adding back ρ_w·L_sl·E_subl leaves the snowpack with a net loss of ρ_w·(L_sg − L_sl)·E_subl = ρ_w·L_lg·E_subl, the vaporization enthalpy carried by the departing vapor.

Note that no explicit meltwater energy term appears because meltwater drains as liquid water at 0 °C, which is the zero-enthalpy reference (U = 0) of the FreeWater closure, so it carries no enthalpy out of the snowpack.

The conductive/sublimation terms Q_base - Q_top + Q_subl are gated to zero when the pack is empty at the start of the step (W == 0), since they are otherwise evaluated against a degenerate zero-depth layer (Q_top/Q_base conduction targets floored at min_conduction_thickness, giving an arbitrarily large flux relative to a genuinely-thin pack's tiny heat capacity). Q_prcp is not gated: it is a pure source term from the precipitation rate alone, well-defined regardless of the existing pack depth. Gating it too — as a previous implementation did, multiplying the whole dUdt including Q_prcp by W > 0 — zeroed the energy input on the very first step any snow accumulates onto bare ground (W is read before that step's mass increment), leaving the freshly-fallen, physically cold snow with zero recorded enthalpy. Under the FreeWater closure (U = 0 ≡ liquid water at 0 °C) that reads as instantaneously fully melted, triggering spurious immediate Darcy drainage regardless of how cold the air actually was.

compute_snow_basal_heat_flux(
    i,
    j,
    grid,
    fields,
    snow::SingleLayerSnow,
    soil::Terrarium.AbstractSoil,
    constants::PhysicalConstants
) -> Any

Snow→soil basal conductive heat flux Q_base [W/m²] at grid cell i, j, positive upward (soil → snow), as the series-resistance conduction between the soil's top half-cell and the snow layer: Q_base = (T_soil − T_snow) / (Δz_soil/(2κ_soil) + d_snow/(2κ_snow)). Both the local soil top-layer thermal conductivity κ_soil and the snow conductivity κ_snow are recovered from their respective composition/density schemes rather than stored. The snow-side conduction thickness is floored at min_snow_conduction_thickness as d_snow → 0, matching the floor used elsewhere in the snow thermodynamics; without the soil-side term this reduces to the previous snow-resistance-only closure, which is only valid once the soil resistance is genuinely negligible next to the (floored) snow's.

compute_snow_soil_heat_flux(
    i,
    j,
    grid,
    fields,
    snow::SingleLayerSnow,
    constants::PhysicalConstants,
    soil::Terrarium.AbstractSoil
) -> Any

Blended soil-top heat flux [W/m²] at grid cell i, j: the snow-cover-fraction-weighted combination of the snow→soil basal conductive flux Q_base (see compute_snow_basal_heat_flux) and the explicit bare-ground conductive flux G (ground_heat_flux, already the unblended per-bare-ground-area quantity — see the ImplicitSkinTemperature-specific compute_ground_heat_flux in skin_temperature.jl), f_snow·Q_base + (1 − f_snow)·G.

compute_snow_sublimation_flux(
    i,
    j,
    grid,
    fields,
    _::Nothing,
    atmos,
    constants,
    skinT
) -> Any

Snow-surface sublimation rate [m/s SWE] at grid cell i, j. The snow surface is treated as saturated: a bulk-aerodynamic vapor flux Δq/rₐ evaluated at the skin temperature, with Δq taken over ice for a sub-freezing surface (the saturation humidity already dispatches over ice for T ≤ 0 — see saturation_specific_humidity_vapor). The water-vapor mass flux ρₐ·Δq/rₐ is converted to a snow-water-equivalent rate via ρ_w. Zero without snow (i.e. when snow === nothing).

compute_snow_interface_fluxes!(
    state,
    grid,
    _::Nothing,
    args...
)

Launch compute_snow_interface_fluxes! to diagnose the snow↔surface/soil coupling fluxes from the snow state and the surface energy balance outputs: the blended soil-top heat flux (soil_heat_flux), the snow-top conductive flux (surface_heat_flux, drives the snowpack's own energy tendency), and the snow surface sublimation rate (sublimation, the snow-fraction bulk-aerodynamic vapor flux at the converged skin temperature). Must run after the surface energy balance (which sets ground_heat_flux and the skin temperature). No-op when there is no snowpack (snow === nothing).

compute_snow_interface_fluxes!(
    out,
    i,
    j,
    grid,
    fields,
    snow::SingleLayerSnow,
    seb::Terrarium.AbstractSurfaceEnergyBalance,
    soil::Terrarium.AbstractSoil,
    constants::PhysicalConstants,
    atmos::Terrarium.AbstractAtmosphere
)

Diagnose the snow↔surface/soil coupling fluxes at grid cell i, j, run after the surface energy balance: the blended soil-top heat flux (see compute_snow_soil_heat_flux), the snow-top conductive flux (see compute_snow_surface_heat_flux), and the snow surface sublimation rate (see compute_snow_sublimation_flux, the same snow-fraction vapor flux the surface energy balance uses for the latent-flux partition).