Skip to content

Module enthalpy

Brief summary

The enthalpy module determines the thermal profile in the ice. It uses an enthalpy-based approach that solves for both cold and temperate ice regions simultaneously (Aschwanden et al., 2012)1. The enthalpy module requires the iceflow module as it depends on ice velocity. In turn, the enthalpy determines the basal friction coefficient and the Arrhenius factor, which are used for ice-flow computation. The parameters of the module are described here.

New variable names

Some enthalpy variables and parameters have been renamed for clarity. The following state variables have been changed:

Previous name New name
basalMeltRate basal_melt_rate
bheatflx basal_heat_flux
tillwat h_water_till
Tpmp T_pmp
Epmp E_pmp
surftemp T_s
surfenth E_s

See the parameters section for parameter name changes.

The implementation is largely inspired by PISM (Aschwanden et al., 2012)1. Other references have also helped verify the implementation through benchmarks (Kleiner et al., 2015; Wang et al., 2020)2 3. The enthalpy module will be described in further detail in the in-preparation paper for IGM 3 (Jouvet et al., 2026)4.

Physical model

The thermal state of glaciers can be complex, with some regions at the pressure-melting point while others are not. Two ice regimes are typically defined:

  • cold ice: \(T < T_\mathrm{pmp}\) and \(\omega = 0\);
  • temperate ice: \(T = T_\mathrm{pmp}\) and \(0 < \omega \le 1\).

Here, \(T\) denotes ice temperature and \(\omega\) is the water content. The pressure-melting-point temperature \(T_\mathrm{pmp}\) depends on pressure through the Clausius-Clapeyron relation, so we may write \(T_\mathrm{pmp}=T_\mathrm{pmp}(p)\) with \(p\) the ice pressure. To avoid dealing with two variables \((T, \omega)\), it is convenient to introduce the enthalpy as follows:

\[ E = \left\{ \begin{aligned} &E_\mathrm{pmp} + c_\mathrm{i}\,(T - T_\mathrm{pmp}), & \quad \text{(cold ice)}\\ &E_\mathrm{pmp} + L\,\omega,& \quad \text{(temperate ice)} \end{aligned} \right. \]

where \(c_\mathrm{i}\) and \(L\) are the specific heat capacity and latent heat of ice, respectively. Importantly, each value of \(E\) corresponds to a unique state \((T, \omega)\), so the values of these variables can be inferred once \(E\) is known.

Since enthalpy can be defined up to an additive constant, we can choose the value of \(E_\mathrm{pmp}\); here, we choose it in a such a way that the enthalpy is zero at a reference temperature \(T_\mathrm{ref}\): $$ E_\mathrm{pmp} = c_\mathrm{i}(T_\mathrm{pmp} - T_\mathrm{ref}). $$

Based on energy conservation, the governing equation for enthalpy takes the form of the following partial differential equation:

\[ \rho_\mathrm{i} \left(\dfrac{\partial E}{\partial t} + u\,\dfrac{\partial E}{\partial x} + v\,\dfrac{\partial E}{\partial y} + w\,\dfrac{\partial E}{\partial z}\right) = \dfrac{\partial}{\partial z} \left(K \,\dfrac{\partial E}{\partial z}\right) + \Phi - \rho_\mathrm{w} L D_\mathrm{w} (\omega), \]

together with suitable boundary conditions (Aschwanden et al., 2012)1. Here, \(\rho_\mathrm{i}\) is the ice density, \(\rho_\mathrm{w}\) is the liquid water density, \(\mathbf{v}=(u, v, w)\) are the velocity components in the ice along each direction, \(K\) is the effective enthalpy diffusivity, \(\Phi=\Phi(\mathbf{v})\) is the strain heating, and \(D_\mathrm{w}\) is a drainage function. The enthalpy diffusivity depends on the thermal state of the ice; here we follow the usual approach in ice-sheet modeling by writing

\[ K = \left\{ \begin{aligned} &k_\mathrm{i}/c_\mathrm{i}, & \quad \text{(cold ice)}\\ &\epsilon\,k_\mathrm{i}/c_\mathrm{i},& \quad \text{(temperate ice)} \end{aligned} \right. \]

where \(k_\mathrm{i}\) is the thermal conductivity of ice and \(\epsilon \ll 1\) denotes the ratio of temperate to cold ice diffusivity.

Numerical solution

At each time step, the enthalpy module performs the following operations:

  1. Computation of surface conditions \((T_\mathrm{s}, E_\mathrm{s})\).
  2. Computation of pressure-melting point conditions \((T_\mathrm{pmp}, E_\mathrm{pmp})\).
  3. Solution of the equation for \(E\).
  4. Computation of the thermal state \((T, \omega)\).
  5. Computation of the Arrhenius factor \(A\).
  6. Computation of till hydrology conditions.
  7. Computation of till friction conditions.

The computationally intensive step is solving the equation for \(E\). To achieve this efficiently, we use an operator splitting method:

  • first, solve the horizontal advection equation with an explicit backward Euler method;
  • then, solve the vertical advection-diffusion equation with an implicit solver.

The implicit solver consists of the Thomas algorithm for tridiagonal linear systems, applied simultaneously to each vertical ice column.

Coupling with ice flow

The enthalpy module builds upon the iceflow module. To activate the coupling, set the following option:

processes.iceflow.physics.sliding.law = weertman

To ensure proper functionality, also activate vertical_iceflow, use a relatively fine vertical discretization, and ensure sufficient retraining.

Parameters

The complete default configuration file can be found here: enthalpy.yaml.

Structure of the parameters:

enthalpy
├── numerics
│   └── ...
├── solver
│   └── ...
├── thermal
│   └── ...
├── surface
│   └── ...
├── arrhenius
│   └── ...
├── till
│   └── ...
└── drainage
    └── ...

Description of the parameters:

numerics
Name Description Default value Units
numerics.Nz Number of vertical layers for the enthalpy discretization. 10 —
numerics.vert_spacing Parameter controlling the vertical spacing; 1.0 means uniform, larger values concentrate points near the bed. 4.0 —
solver
Name Description Default value Units
solver.allow_basal_refreezing Allow negative basal melt rates (refreezing). True —
solver.correct_w_for_melt Correct vertical velocity for basal melting. True —
solver.override_basal_at_pmp Force basal enthalpy to be at pressure-melting point. False —
thermal
Name Description Default value Units
thermal.T_min Minimum temperature bound for enthalpy solver. 243.15 K
thermal.T_pmp_ref Reference pressure melting point temperature. 273.15 K
thermal.T_ref Reference temperature for enthalpy calculations. 223.15 K
thermal.c_ice Specific heat capacity of ice. 2009.0 J kg\( ^{-1} \) K\( ^{-1} \)
thermal.L_ice Latent heat of fusion for ice. 334000.0 J kg\( ^{-1} \)
thermal.k_ice Thermal conductivity of ice. 2.1 W m\( ^{-1} \) K\( ^{-1} \)
thermal.K_ratio Ratio of temperate to cold ice thermal diffusivity. 1e-1 —
thermal.basal_heat_flux_ref Reference geothermal heat flux (used when not specified in state). 0.065 W m\( ^{-2} \)
thermal.beta Clausius-Clapeyron constant. 7.9e-08 K Pa\( ^{-1} \)
surface
Name Description Default value Units
surface.T_offset Temperature offset for surface conditions. 0.0 °C
arrhenius
Name Description Default value Units
arrhenius.A_cold Pre-exponential factor for cold ice. 3.985e-13 Pa\( ^{-n} \) s\( ^{-1} \)
arrhenius.A_warm Pre-exponential factor for warm ice. 1916.0 Pa\( ^{-n} \) s\( ^{-1} \)
arrhenius.Q_cold Activation energy for cold ice. 60000.0 J mol\( ^{-1} \)
arrhenius.Q_warm Activation energy for warm ice. 139000.0 J mol\( ^{-1} \)
arrhenius.T_threshold Temperature threshold between cold and warm regimes. 263.15 K
arrhenius.omega_coef Water content enhancement coefficient. 181.25 —
arrhenius.omega_max Maximum water content for enhancement. 0.01 —
arrhenius.R Universal gas constant. 8.314 J mol\( ^{-1} \) K\( ^{-1} \)
till
Name Description Default value Units
till.friction.phi Till friction angle (uniform case). 30.0 °
till.friction.phi_min Friction angle at bed_min. 15.0 °
till.friction.phi_max Friction angle at bed_max. 45.0 °
till.friction.bed_min Minimum bed elevation for interpolation. None m
till.friction.bed_max Maximum bed elevation for interpolation. None m
till.friction.tauc_min Minimum yield stress. 100000.0 Pa
till.friction.tauc_max Maximum yield stress. 10000000000.0 Pa
till.friction.tauc_ice_free Yield stress for ice-free areas. 1000000.0 Pa
till.friction.u_ref Reference sliding velocity. 100.0 m yr\( ^{-1} \)
till.hydro.h_water_till_max Maximum till water layer thickness. 2.0 m
till.hydro.N_ref Reference effective pressure. 1000.0 Pa
till.hydro.e_ref Reference void ratio. 0.69 —
till.hydro.C_c Till compressibility coefficient. 0.12 —
till.hydro.delta Minimum effective pressure fraction. 0.02 —
till.hydro.drainage_rate Till drainage rate. 0.001 m yr\( ^{-1} \)
drainage
Name Description Default value Units
drainage.water_density Water density constant. 1000.0 kg m\( ^{-3} \)
drainage.drain_ice_column Enable water drainage through ice column. True —
drainage.omega_target Target water content after drainage. 0.01 —
drainage.omega_threshold_1 First water content threshold. 0.01 —
drainage.omega_threshold_2 Second water content threshold. 0.02 —
drainage.omega_threshold_3 Third water content threshold. 0.03 —

Contributors: G. Jouvet, T. Gregov, L. Bacchin.

  1. Aschwanden, A., Bueler, E., Khroulev, C., & Blatter, H. (2012). An enthalpy formulation for glaciers and ice sheets. Journal of Glaciology, 58(209), 441--457. https://doi.org/10.3189/2012jog11j088 ↩↩↩

  2. Kleiner, T., Rückamp, M., Bondzio, J. H., & Humbert, A. (2015). Enthalpy benchmark experiments for numerical ice sheet models. The Cryosphere, 9(1), 217--228. https://doi.org/10.5194/tc-9-217-2015 ↩

  3. Wang, Y., Zhang, T., Xiao, C., Ren, J., & Wang, Y. (2020). A two-dimensional, higher-order, enthalpy-based thermomechanical ice flow model for mountain glaciers and its benchmark experiments. Computers & Geosciences, 141, 104526. https://doi.org/10.1016/j.cageo.2020.104526 ↩

  4. Jouvet, G., Cook, S., Cordonnier, G., Finley, B., Henz, A., Herrmann, O., Maussion, F., Mey, J., Scherler, D., & Welty, E. (2026). Concepts and capabilities of the instructed glacier model. https://doi.org/10.31223/x5t99c ↩