Module `subglacial_hydrology`

This IGM module is the single source of truth for subglacial hydrology. It computes and maintains state.effective_pressure (MPa) — the basal effective pressure \(N = p_i - p_w\) consumed by the Budd, regularized Coulomb, and Mohr-Coulomb sliding laws — and, when mode: till_storage is selected, also evolves the subglacial till water layer \(W\) (state.h_water_till, m water) forward in time.

State variables

Reads: thk

Writes: effective_pressure

The module is opt-in. If the only sliding law in use is Weertman, \(N\) cancels out of the cost and the module should not be activated. Conversely, when Budd, regularized Coulomb, or Mohr-Coulomb is active, subglacial_hydrology is required — the iceflow energy validator refuses to run without effective_pressure in its input list.

Closure modes

The mode parameter selects how \(N\) is computed at each update:

constant_one — \(N = 1\) MPa everywhere. Useful as a sanity check, or to recover Weertman-like behaviour with the Budd cost when paired with N_ref = 1.
percentage — \(N = (1 - \texttt{percentage}) \cdot \rho_i\, g\, h\). A fixed fraction of the ice-overburden pressure.
ocean_connected — \(N = \rho_i g h - \rho_w g \max(w - b,\, 0)\), where \(w\) is the local water level (state.water_level) and \(b\) is the bedrock elevation (state.topg). Assumes the basal hydraulic system is connected to the ocean / proglacial water body. Recommended for marine / tidewater settings. Requires state.water_level — see the local and load_ncdf input modules.
from_input — leave state.effective_pressure as it was loaded from the input data (e.g. an inversion product). The module acts as a pass-through.
till_storage — full Tulaczyk till hydrology: evolves state.h_water_till each time step and derives \(N\) from the resulting saturation. See the physical model below. Requires enthalpy to run before subglacial_hydrology.

After the chosen formula, \(N\) is clamped from below by N_min (default 1 kPa) to avoid the singularity in the sliding-law cost when \(N \to 0\).

Physical model (`mode: till_storage`)

At every time step the module first integrates the till water ODE:

\[W^{n+1} = \text{clip}\!\left(W^n + \Delta t \left(\frac{\rho_i}{\rho_w}\,\dot{m}_b - d_r\right),\; 0,\; W_{\max}\right)\]

where \(\dot{m}_b\) is the basal melt rate (m ice yr⁻¹) from state.basal_melt_rate, \(\rho_i / \rho_w\) converts ice melt to water-equivalent volume, \(d_r\) is the constant drainage rate (till_storage.drainage_rate), and \(W_{\max}\) is the till capacity (till_storage.h_water_till_max). The result is zeroed where ice thickness is zero.

The saturation \(s = W / W_{\max}\) then drives the Tulaczyk (2000) effective-pressure parameterisation:

\[N = \min\!\left(p_i,\; N_\text{ref} \left(\frac{\delta\, p_i}{N_\text{ref}}\right)^s \cdot 10^{e_\text{ref}(1-s)/C_c}\right)\]

Unit convention

\(N\) is stored in MPa to match the rest of the iceflow stress quantities (slidingco, viscosity costs, surface stress). Literature reference values for \(N_\text{ref}\) are typically \(\mathcal{O}(1\,\text{MPa})\) — see (Brondex et al., 2019)¹, (Pollard & DeConto, 2012)², (Pattyn, 2017)³.

Module ordering

subglacial_hydrology must run before iceflow so that state.effective_pressure is up to date when the sliding cost is assembled. When mode: ocean_connected is used it must also run after the input module that created state.water_level. When mode: till_storage is used it must run after enthalpy so that state.basal_melt_rate is fresh:

override /processes:
  - time
  - iceflow
  - thk
  - enthalpy           # provides state.basal_melt_rate (needed for till_storage)
  - arrhenius
  - subglacial_hydrology   # after enthalpy, before iceflow

Parameters

Default configuration file (subglacial_hydrology.yaml):

subglacial_hydrology:
  # How to compute state.effective_pressure (MPa). Modes:
  #   constant_one    — N = 1 MPa everywhere (use with Budd N_ref=1 for Weertman-like behaviour)
  #   percentage      — N = (1 - percentage) * rho_ice * g * h
  #   ocean_connected — N = rho_ice*g*h - rho_water*g*max(water_level - topg, 0)
  #   from_input      — leave state.effective_pressure as already loaded
  #   till_storage    — Tulaczyk till hydrology (Van Pelt & Bueler 2015);
  #                     evolves state.h_water_till each time step, then derives
  #                     N from saturation. Requires `enthalpy` to run before
  #                     `subglacial_hydrology` in the process order.
  #
  # Unit is MPa to match the rest of the iceflow stresses (slidingco, σ, τ).
  # Literature N_ref values are typically O(1 MPa) (Brondex et al. 2019,
  # Pollard & DeConto 2012, Pattyn 2017).
  mode: ocean_connected
  percentage: 0.0      # only used if mode == percentage
  N_min: 1.0e-3        # MPa: floor (1 kPa) applied after the closed-form modes

  till_storage:
    # Till water layer evolution (Tulaczyk 2000 drainage ODE).
    # Only active when mode == till_storage.
    h_water_till_max: 2.0     # m water: maximum till water thickness (capacity)
    drainage_rate:   0.001    # m water / yr: basal drainage rate to the groundwater system
    water_density:   1000.0   # kg m^-3: density of liquid water
    # Tulaczyk effective-pressure parameters.
    # h_water_till_max above drives saturation s = h_water_till / h_water_till_max.
    N_ref:  1.0e-3   # MPa (legacy: 1000 Pa = 1e-3 MPa)
    e_ref:  0.69
    C_c:    0.12
    delta:  0.02

Structure of the parameters:

subglacial_hydrology
├── mode
├── percentage
├── N_min
└── till_storage
    └── ...

Description of the parameters:

Name	Description	Default value	Units
`mode`	Closure used to compute `state.effective_pressure` (MPa). One of: `constant_one` (N = 1 MPa everywhere), `percentage` (N = (1 - percentage) * rho_ice * g * h), `ocean_connected` (N = rho_icegh - rho_watergmax(water_level - topg, 0); requires `state.water_level`), `from_input` (leave `state.effective_pressure` as already loaded from the input data), or `till_storage` (evolves `state.h_water_till` each time step via the Tulaczyk (2000) ODE, then derives N from saturation; requires `enthalpy` to run before `subglacial_hydrology`).	ocean_connected	—
`percentage`	Fraction of the ice-overburden pressure subtracted in `mode: percentage`. 0 means N equals overburden; 1 means N is zero. Ignored for other modes.	0.0	—
`N_min`	Lower floor applied to N after the chosen formula (default 1.0e-3 MPa = 1 kPa) to avoid the singularity in the Budd / Coulomb sliding cost.	0.001	MPa

till_storage

Name	Description	Default value	Units
`till_storage.h_water_till_max`		2.0
`till_storage.drainage_rate`		0.001
`till_storage.water_density`		1000.0
`till_storage.N_ref`		0.001
`till_storage.e_ref`		0.69
`till_storage.C_c`		0.12
`till_storage.delta`		0.02

Example usage

Ocean-connected basal pressure with Budd sliding

# @package _global_

inputs:
  local:
    filename: input.nc
    water_level:
      include: True
      value: 0.0          # uniform sea level at 0 m

processes:
  subglacial_hydrology:
    mode: ocean_connected
    N_min: 1.0e-3
  iceflow:
    physics:
      sliding:
        law: budd
        N_ref: 1.0      # MPa
        q_exponent: 1.0
  thk: {}
  time:
    start: 2000.0
    end: 2100.0
    save: 5.0

Till hydrology with Mohr-Coulomb sliding

# @package _global_

processes:
  enthalpy: {}
  arrhenius: {}
  subglacial_hydrology:
    mode: till_storage
    till_storage:
      h_water_till_max: 2.0      # m — till capacity
      drainage_rate: 0.001       # m/yr — basal drainage
      N_ref: 1.0e-3              # MPa
      e_ref: 0.69
      C_c: 0.12
      delta: 0.02
  iceflow:
    physics:
      sliding:
        law: mohr_coulomb

Or override on the command line:

igm_run +experiment/params processes.subglacial_hydrology.mode=percentage processes.subglacial_hydrology.percentage=0.96

Contributors: Guillaume Jouvet, Thomas Gregov.

Brondex, J., Gillet-Chaulet, F., & Gagliardini, O. (2019). Sensitivity of centennial mass loss projections of the amundsen basin to the friction law. The Cryosphere, 13(1), 177--195. https://doi.org/10.5194/tc-13-177-2019 ↩
Pollard, D., & DeConto, R. M. (2012). Description of a hybrid ice sheet-shelf model, and application to antarctica. Geoscientific Model Development, 5(5), 1273--1295. https://doi.org/10.5194/gmd-5-1273-2012 ↩
Pattyn, F. (2017). Sea-level response to melting of antarctic ice shelves on multi-centennial timescales with the fast elementary thermomechanical ice sheet model (f.ETISh v1.0). The Cryosphere, 11(4), 1851--1878. https://doi.org/10.5194/tc-11-1851-2017 ↩

Module subglacial_hydrology