Ice sheets, glacial erosion and meltwater#

goSPL represents glaciers with a cheap, robust diagnostic glacial-erosion model. Each step the ELA surface mass balance is routed downhill into an ice discharge, from which an ice thickness, a basal sliding velocity, glacial abrasion, till/moraine deposition, meltwater and an ice load are derived — one linear solve, no time integration. The model is enabled by adding an ice section to the input file.

Note

Why a diagnostic, not full ice dynamics. A true Shallow-Ice-Approximation (SIA) thickness solve was implemented and then removed. The ice margin \(H \ge 0\) is a free-boundary (obstacle) problem, on which an implicit F(H)=0 thickness solve diverged on real continental runs; and an ice-dynamics solve is stiff and over-thickens km-scale continental ice at coarse resolution over goSPL’s long (\(10^2\)–\(10^4\) yr) timesteps. The diagnostic below is physical and robust at any resolution and over those long steps, which is why it is the only model goSPL ships.

Surface mass balance#

The accumulation/ablation source \(\dot{m}\) (m ice/yr) at each cell is a fraction of the local precipitation set by the surface elevation relative to two time-varying altitudes — the equilibrium-line altitude (hela) and the ice-cap altitude (hice):

\[\mathrm{\dot{m}_i} = \mathrm{P_i} \cdot \min\!\left(1, \frac{\eta_i - h_{ela}}{h_{ice} - h_{ela}}\right)\]

Above the ice-cap altitude the cell accumulates the full local precipitation as ice; between the ELA and the ice cap the fraction ramps linearly from 0 to 1; below the ELA the ramp goes negative and the source becomes ablation (melt). The accumulation part is scaled by accum_factor (a precipitation→ice conversion fraction) and optionally capped at accum_max (m ice/yr); ablation is amplified by melt.

A degenerate-configuration guard zeroes the ice and returns immediately when hice <= hela or when the maximum surface elevation lies below the ELA.

Note

For global models a single ELA is unphysical — it ranges from ~5000–6000 m in the tropics to near sea level at the poles. hela, hice and hterm may therefore each be a per-vertex map rather than a scalar, and may vary in time through a glaciers time series (mirroring the precipitation climate block). The mass-balance ramp is then evaluated per node with its local ELA / ice-cap altitude, so the same run can grow tropical summit glaciers and polar ice sheets simultaneously. See the user guide for the input syntax.

Ice routing and discharge#

The accumulated ice is routed downhill on the epsilon-filled bed using a multiple-flow-direction (MFD) algorithm — the same flow-matrix / KSP machinery that builds the river flow accumulation, with the number of directions set by icedir. A single linear solve turns the per-cell accumulation into an ice discharge \(Q\) (m3/yr) at each cell: the volume of ice passing through it per year. There is no time integration of an ice-thickness PDE.

Ice thickness (Bahr scaling)#

The ice thickness \(H\) follows from a Bahr-type width–area (here discharge) scaling of the discharge:

\[H = \mathrm{eheight}\,\cdot\,\mathrm{fwidth}\,\cdot\,Q^{\,0.3}\]

with the Bahr thickness factor eheight and width factor fwidth. Thicker ice forms where the routed discharge is larger (the trunk of the ice stream), tapering to zero where there is no ice.

Basal sliding velocity (Glen sliding law)#

From that thickness and the bed-surface slope, goSPL derives the basal sliding velocity \(u_b\) from Glen’s sliding law (the ice_velocity Fortran kernel):

\[u_b \;\propto\; H^{\,n-1}\,|\nabla s|^{\,n-1}\,\nabla s, \qquad s = z_\mathrm{bed} + H\]

with sliding coefficient slide and Glen exponent \(n\) = glen (usually 3). The velocity is physically bounded. It is written to the output as iceUb and drives glacial abrasion.

Glacial abrasion (vertical and lateral)#

When abrasion.Kg > 0, sliding ice abrades the bed (deepening valleys) at the rate

\[E_g = K_g\,|u_b|^{\,l}\]

(Kg and l in the abrasion block). The abrasion is masked to subaerial, ice-covered cells (no marine abrasion; \(u_b\) is already zero where there is no ice).

An optional lateral (valley-wall) erosion term widens glaciated valleys toward a U-profile. With abrasion.Kl > 0, each wall cell flanking fast ice is eroded at \(K_l\,|u_{b,\mathrm{neighbour}}|^{\,\mathrm{lat\_l}}\), tapered by how much of the wall is in contact with the neighbouring ice column. The eroded wall rock joins the same conserved till → moraine budget as the vertical abrasion. lat_l (default = l) is its velocity exponent. Both Kg and Kl default to 0 (abrasion off).

When till handling is off, the abraded material is added directly to the erosion–deposition rate as an incision and flows into the fluvial sediment system through the standard erosion path in all SPL flavours (SPL, nlSPL, soilSPL).

Glacial till and moraine deposition#

When till.on is set (the default), the abraded rock is not released straight into the rivers but carried as till by the ice and deposited as moraine where the ice melts out — the ablation zone / terminus. Rock volume is conserved (rock is moved, not created or destroyed). Two modes are used depending on whether stratigraphy is active:

  • Bulk bed (no stratigraphy): the abraded volume is redistributed as a bed-to-bed transport, so the net bed-volume change is exactly zero.

  • Stratigraphic: the till is removed from the stratigraphic layers it was abraded from and re-deposited as a fresh moraine layer in the ablation zone. Under dual lithology the moraine is split into coarse and fine fractions carrying the abraded (ice-mixed) fine fraction, so the per-fraction solid budget stays balanced. The bed bulks up by the porosity contrast between the uncompacted till and the compacted source rock.

Independently of the bulk/stratigraphic split, the deposition distribution has two options. With till.route: True (the default) the till is routed down the ice-surface flow network: it is transported along steepest descent of the surface \(s = z_\mathrm{bed} + H\) and melts out progressively toward each catchment’s terminus, building moraine at the actual ice margins. This is a transport-with-loss 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 = \min(1, \dot a_i \Delta t / H_i)\) the melt-out fraction (forced to 1 at the ice margin so no till leaks onto bare ground). With till.route: False the global abraded volume is instead spread across the whole ablation zone weighted by the local meltwater rate — cheaper (no extra solve) but it decouples erosion and deposition across separate ice masses. Both options conserve mass and are protected by dedicated tests.

Meltwater re-injection into the river network#

By default (melt_conserve: True) the glacial meltwater delivered to the rivers is discharge-conserving. The precipitation that fell as ice above the ELA is routed down-glacier and released as liquid meltwater where the ice melts out, so the total meltwater equals the total accumulation — closing the glacial water budget so downstream basins do not under-predict discharge. The river accumulation step in flowAccumulation removes the precipitation that was captured as ice from the river source and adds this meltwater back, so cells downstream of glacier termini see the corresponding meltwater discharge instead of losing it.

With melt_conserve: False goSPL reverts to the local precipitation-scaled ablation rate,

\[\mathrm{m_i} = \max(-\dot{m}_i, 0)\;\Omega_i\;[H_i > 0]\]

which is cheaper but generally returns less water than fell as ice.

Ice loading and flexure#

When flexure is enabled, the diagnostic ice thickness is used as the ice load in the flexural-isostasy computation: the change in ice thickness between steps is converted to an equivalent load (scaled by \(\rho_i / \rho_c\)) and applied through the existing flexure path, so a growing ice sheet drives isostatic subsidence and deglaciation drives rebound.

Output#

Four ice fields are written to the HDF5 / XDMF output whenever the ice model is on:

  • iceH — ice thickness (m); restored on restart;

  • iceUb — basal sliding speed (m/yr);

  • iceMelt — ablation meltwater (m3/yr) re-injected into the rivers (the glacial contribution to downstream discharge);

  • iceAbr — glacial abrasion rate \(E_g = K_g|u_b|^{l}\) (m/yr); zero where abrasion is off (Kg = 0).