Class IceMesh#

class flow.iceplex.IceMesh(*args, **kwargs)[source]#

Diagnostic glacial-erosion model.

Rather than time-integrating an ice-dynamics PDE, this class derives a cheap, stable diagnostic of the glacial state each step and uses it to drive glacial erosion, sediment transport and ice loading. For the equilibrium-line-altitude (ELA) surface mass balance \(\dot{m}\):

The accumulation \(\dot{m}^+\) is routed downhill on the (epsilon-filled) bed using a multiple-flow-direction algorithm — the same machinery as the river flow accumulation — into an ice discharge \(Q\) (m³/yr); one linear solve, no time integration.
The ice thickness follows a Bahr width–area scaling of the discharge, \(H = e_h\, f_w\, Q^{0.3}\).
The basal sliding velocity comes from Glen’s sliding law on that thickness and the bed-surface slope (the ice_velocity kernel, \(u_b \propto H^{n-1}|\nabla s|^{n-1}\nabla s\)) — physically bounded, so erosion does not spike at flow-convergence cells.

From these the model computes velocity-based glacial abrasion (\(E_g = K_g |u_b|^l\), with an optional lateral valley-wall term for U-shaping), ice-transported till deposited as moraine where the ice melts out (conserving rock volume), the discharge-conserving glacial meltwater re-injected into the river flow accumulation, and the ice load feeding the flexural isostasy.

The diagnostic is robust and physical at any resolution and over goSPL’s long (\(10^2\)–\(10^4\) yr) timesteps, where a true ice-dynamics solve is stiff and over-thickens km-scale continental ice (see docs/DESIGN_ICE_SHEET.md for the history).

Useful links: - Braun et al. (1999) - Herman & Braun (2008) - Egholm (2012) - Deal & Prasicek (2020) - Hergarten (2021) - Liebl et al. (2023)

Methods

`glacialTill`()	Glacial till — production, transport and deposition.
`iceAccumulation`()	Main glacial update.

Initialise

__init__(*args, **kwargs)

Initialisation of the IceMesh class.

Public Methods

`iceAccumulation`()	Main glacial update.
`glacialTill`()	Glacial till — production, transport and deposition.

Private Methods

`_iceMassBalance`(elaH, iceH)	Bed elevation zbed and the ELA surface mass balance mdot (m ice/yr, no meltfac): accumulation above the ELA, ablation below.
`_matrixIceFlow`([dir_ice])	Build the multiple-flow-direction (MFD) routing matrix for ice on the (epsilon-filled) bed surface — used by the diagnostic `mfd` flow model to accumulate the ELA mass balance into an ice discharge.
`_iceFlowMFD`(elaH, iceH, iceT)	Diagnostic glacial driver (no ice-dynamics solve).
`_glacialLateralErosion`()	Lateral glacial erosion rate (m/yr) — valley-wall abrasion by adjacent fast ice, the explicit term that widens glaciated valleys toward a U-profile (vertical abrasion alone only deepens the trough).
`_glacialTillStrata`(dz_ero, dep_w, owned)	Stratigraphic / dual-lithology coupling for glacial till.
`_glacialMeltwater`(zbed, mdot)	Glacial meltwater volume (m^3/yr) delivered to the rivers.
`_routeTill`(Vero, Vtot, owned)	Catchment-aware till routing on the ice-surface flow network (`till.route: True`).

Public functions#

IceMesh.iceAccumulation()[source]#

Main glacial update.

Computes the diagnostic glacial state (_iceFlowMFD()) from the ice-cap altitude, the equilibrium-line altitude (ELA) and the glacier terminus. Each of these may be a uniform scalar (time-varying) or a per-vertex map — the latter for global models where the ELA varies strongly with latitude. It also snapshots the ice load for the flexural isostasy when enabled.

IceMesh.glacialTill()[source]#

Glacial till — production, transport and deposition.

Glacial abrasion (\(E_g = K_g |u_b|^l\)) erodes rock under sliding ice; the abraded material becomes till carried by the ice and deposited where the ice melts out — the ablation zone / terminal moraine. The bed is therefore lowered under fast ice and raised in the ablation zone: a transport, not a source/sink.

Two modes, selected by whether stratigraphy is on:

Bulk bed (no stratigraphy): the abraded volume Vtot is redistributed to the ablation cells weighted by the meltwater rate (iceMeltL) as a bed-to-bed transport, so Σ deposited == Vtot and the net bed-volume change is zero by construction.
Stratigraphic (_glacialTillStrata()): the till is removed from the layers it was abraded from and re-deposited as a fresh moraine layer in the ablation zone — split into the coarse/fine lithology fractions under dual lithology. Conservation is on the solid phase (per fraction), so the bed bulks up by the porosity contrast (uncompacted till vs compacted source rock).

Both modes are guarded by a dedicated till-conservation test (the dual-lithology lesson: a volume the total budget can’t see needs its own guard).

Active only with till.on and Kg > 0; a no-op otherwise.

Deposition distribution. By default the till is spread across the whole ablation zone weighted by the meltwater rate. With till.route: True it is instead routed down the ice-surface flow network and melts out toward each catchment’s terminus (_routeTill()) — appropriate for high-resolution regional runs where individual glacier catchments and termini are resolved. Both produce a per-cell deposition weight (summing to one) consumed identically by the bulk and stratigraphic paths, so conservation is unchanged.

Private functions#

IceMesh._iceMassBalance(elaH, iceH)[source]#: Bed elevation zbed and the ELA surface mass balance mdot (m ice/yr, no meltfac): accumulation above the ELA, ablation below.

IceMesh._matrixIceFlow(dir_ice=1)[source]#: Build the multiple-flow-direction (MFD) routing matrix for ice on the (epsilon-filled) bed surface — used by the diagnostic mfd flow model to accumulate the ELA mass balance into an ice discharge. Mirrors the water flow-direction matrix; deterministic slope tie-break via gid.

IceMesh._iceFlowMFD(elaH, iceH, iceT)[source]#

Diagnostic glacial driver (no ice-dynamics solve).

ELA surface mass balance (_iceMassBalance: accumulation above the ELA, ablation below — including the accum_factor/accum_max controls).
Route the accumulation downhill through the MFD matrix into an ice discharge iceFA (one linear solve; no time integration).
Thickness iceHL from a Bahr width–area scaling of the discharge.
Basal sliding velocity iceUbL from Glen’s sliding law on that thickness + bed-surface slope (ice_velocity), the abrasion driver E_g = K_g |u_b|^l. Physically bounded (∝ H^{n-1}|∇s|^{n-1}∇s), unlike a raw balance velocity Q/(H·W) which blows up with catchment size and spikes at flow-convergence cells.

No stiffness, no CFL — one routing solve per step.

IceMesh._glacialLateralErosion()[source]#

Lateral glacial erosion rate (m/yr) — valley-wall abrasion by adjacent fast ice, the explicit term that widens glaciated valleys toward a U-profile (vertical abrasion alone only deepens the trough).

Returns a per-cell rate (≥0), nonzero only on subaerial ‘wall’ cells (little ice of their own) standing within the ice column of a faster-sliding neighbour — see the ice_lateral_erosion kernel. No-op (zeros) unless ice.abrasion.Kl > 0.

IceMesh._glacialTillStrata(dz_ero, dep_w, owned)[source]#

Stratigraphic / dual-lithology coupling for glacial till.

The abraded bed lowering dz_ero (m, ≤0) is removed from the stratigraphic pile with erodeStrat(), which returns the uncompacted solid removed per fraction (thCoarse / thFine, m/yr). That solid is the till: its total uncompacted volume is redistributed to the ablation cells weighted by the meltwater rate and laid down as a fresh moraine layer with deposeStrat(), carrying the abraded fine fraction (a single, ice-mixed composition).

Solid mass is conserved per fraction by construction: the fine volume deposited equals the fine volume eroded (so the dual-lithology _fineEroded / _fineDeposited budget stays balanced), and the coarse volume likewise. The bed bulks up by the porosity contrast between the uncompacted till and the compacted source rock.

erodeStrat / deposeStrat overwrite the transient routing arrays (thCoarse / thFine / depoFineFrac) that the fluvial sedChange consumes later this step, so they are saved and restored around the till calls.

IceMesh._glacialMeltwater(zbed, mdot)[source]#

Glacial meltwater volume (m^3/yr) delivered to the rivers.

Two models, selected by ice.melt_conserve:

Discharge-conserving (default): the accumulation (the water that fell as ice above the ELA, \(\dot{m}^+\,A\)) is routed down the ice-surface flow network and released as meltwater where the ice melts out (fraction \(f\), forced to 1 at the margin/terminus), so \(\sum \mathrm{melt} = \sum \mathrm{accumulation}\). Over goSPL’s long timesteps a land-terminating glacier is ~steady, so all the ice that accumulates leaves as meltwater at the snout — this closes the glacial water budget (the precipitation removed from runoff above the ELA returns downstream). Transport-with-loss on the ice network, same MPI-correct flow-matrix / KSP machinery as _routeTill.
Precip-scaled (melt_conserve: False): the local ablation rate \(\max(-\dot{m},0)\,A\) where ice is present — cheaper, but it loses water (the ablation generally under-returns the accumulation).

IceMesh._routeTill(Vero, Vtot, owned)[source]#

Catchment-aware till routing on the ice-surface flow network (till.route: True).

The abraded till is transported down the ice surface \(s = z_\mathrm{bed} + H\) along steepest descent (the same direction the ice flux follows) and melts out progressively toward each catchment’s terminus, building moraine where the ice actually ends rather than smearing it across the whole ablation zone. This matters at high (sub-km) resolution where individual glacier catchments and termini are resolved.

It is a transport-with-loss on the ice network, solved with the same MPI-correct flow-matrix / KSP machinery as the river accumulation:

\[L_i = A_i + \sum_{u \to i} w_{ui}\,(1 - f_u)\,L_u, \qquad D_i = f_i\,L_i\]

where \(A_i\) is the local abraded volume, \(L_i\) the till load passing through cell \(i\), and \(f_i \in [0,1]\) the melt-out fraction — the share of the local ice column lost to ablation this step, \(f_i = \min(1, \dot{a}_i\,\Delta t / H_i)\). At the ice margin (a cell whose steepest-descent receiver is ice-free) \(f\) is forced to 1 so no till leaks onto bare ground, which makes the transport exactly conservative: \(\sum_i D_i = \sum_i A_i\).

Parameters:

Vero – per-cell abraded volume (m^3, local array).
Vtot – total abraded volume over owned nodes (m^3).
owned – boolean mask of owned (non-ghost) nodes.

Returns:

per-cell deposition weight (local array, Σ over owned = 1).