Subglacial hydrology¶
At the present time, two simple subglacial hydrology models are implemented and
documented in PISM, namely -hydrology null (and its modification steady) and
-hydrology routing; see Table 15 and [91].
In both models, some of the water in the subglacial layer is stored locally in a layer of
subglacial till by the hydrology model. In the routing model water is conserved by
horizontally-transporting the excess water (namely bwat) according to the gradient of
the modeled hydraulic potential. In both hydrology models a state variable tillwat is
the effective thickness of the layer of liquid water in the till; it is used to compute
the effective pressure on the till (see Controlling basal strength). The pressure of the
transportable water bwat in the routing model does not relate directly to the
effective pressure on the till.
Note
Both null and routing described here provide all diagnostic quantities needed
for mass accounting, even though null is not mass-conserving. See
Mass accounting in subglacial hydrology models for details.
Option |
Description |
|---|---|
|
The default model with only a layer of water stored in till. Not mass conserving in
the map-plane but much faster than |
|
A version of the |
|
A mass-conserving horizontal transport model in which the pressure of transportable
water is equal to overburden pressure. The till layer remains in the model, so this
is a “drained and conserved plastic bed” model. The state variables are |
See Table 16 for options which apply to all hydrology models. Note
that the primary water source for these models is the energy conservation model which
computes the basal melt rate basal_melt_rate_grounded. There is, however, also option
-hydrology_input_to_bed_file which allows the user to add water directly into the
subglacial layer, in addition to the computed basal_melt_rate_grounded values. Thus
-hydrology_input_to_bed_file allows the user to model drainage directly to the bed
from surface runoff, for example. Also option -hydrology_bmelt_file allows the user
to replace the computed basal_melt_rate_grounded rate by values read from a file,
thereby effectively decoupling the hydrology model from the ice dynamics
(esp. conservation of energy).
To use water runoff computed by a PISM surface model, set
hydrology.surface_input_from_runoff. (The Temperature-index scheme computes runoff
using computed melt and the assumption that a fraction of this melt re-freezes, all other
models assume that all melt turns into runoff.)
Option |
Description |
|---|---|
|
Specifies a NetCDF file which contains a time-dependent field |
|
Maximum effective thickness for water stored in till. |
|
Water accumulates in the till at the basal melt rate |
|
Replace the conservation-of-energy basal melt rate |
The default model: -hydrology null¶
In this model the water is not conserved but it is stored locally in the till up to a
specified amount; the configuration parameter hydrology.tillwat_max sets that
amount. The water is not conserved in the sense that water above the
hydrology.tillwat_max level is lost permanently. This model is based on the
“undrained plastic bed” concept of [93]; see also
[10].
In particular, denoting tillwat by \(W_{till}\), the till-stored water layer effective
thickness evolves by the simple equation
where \(m=\) basal_melt_rate_grounded (kg \(\text{m}^{-2}\,\text{s}^{-1}\)), \(\rho_w\)
is the density of fresh water, and \(C\) hydrology_tillwat_decay_rate. At all times
bounds \(0 \le W_{till} \le W_{till}^{max}\) are satisfied.
This -hydrology null model has been extensively tested in combination with the
Mohr-Coulomb till (section Controlling basal strength) for modelling ice streaming (see
[2] and [10], among others).
“Steady state flow” model¶
This model (-hydrology steady) adds an approximation of the subglacial water flux to
the null model. It does not model the steady state distribution of the subglacial
water thickness.
Here we assume that the water input from the surface read from the file specified using
hydrology.surface_input.file instantaneously percolates to the base of the ice
and enters the subglacial water system.
We also assume that the subglacial drainage system instantaneously reaches its steady
state. Under this assumption the water flux can be approximated by using an auxiliary
problem that is computationally cheaper to solve, compared to using the mass conserving
routing model (at least at high grid resolutions).
We solve
on the grounded part of the domain with the initial state \(u_0 = \tau M\), where \(\tau\) is
the scaling of the input rate (hydrology.steady.input_rate_scaling) and \(M\) is
the input rate read from an input file.
The velocity field \(\V\) approximates the steady state flow direction. In this implementation \(\V\) is the normalized gradient of
Here \(H\) is the ice thickness, \(b\) is the bed elevation, \(g\) is the acceleration due to gravity and \(\rho_i\), \(\rho_w\) are densities of ice and water, respectively.
Note
Note that \(\psi\) ignores subglacial water thickness, producing flow patterns that are more concentrated than would be seen in a model that includes it. This effect is grid-dependent.
This implies that this code cannot model the distribution of the flow along the grounding line of a particular outlet glacier but it can tell us how surface runoff is distributed among the outlets.
The term \(\Delta \psi\) is the adjustment needed to remove internal minima from the “raw”
potential, filling any “lakes” it might have. This modification of \(\psi\) is performed
using an iterative process; set hydrology.steady.potential_n_iterations to
control the maximum number of iterations and hydrology.steady.potential_delta to
affect the amount it is adjusted by at each iteration.
The equation (28) is advanced forward in time until \(\int_{\Omega}u
< \epsilon\int_{\Omega} u_0\), where \(\epsilon\) (hydrology.steady.volume_ratio)
is a fraction controlling the accuracy of the approximation: lower values of \(\epsilon\)
would result in a better approximation and a higher computational cost (more iterations).
Set hydrology.steady.n_iterations to control the maximum number of these
iterations.
This model restricts the time step length in order to capture the temporal variability of the forcing: the flux is updated at least once for each time interval in the forcing file.
In simulations with a slowly-varying surface input rate the flux has to be updated every
once in a while to reflect changes in the flow pattern coming from changing geometry. Use
hydrology.steady.flux_update_interval (years) to set the update frequency.
See Computing steady-state subglacial water flux for technical details.
The mass-conserving model: -hydrology routing¶
In this model the water is conserved in the map-plane. Water does get put into the till,
with the same maximum value hydrology.tillwat_max, but excess water is
horizontally-transported. An additional state variable bwat, the effective thickness
of the layer of transportable water, is used by routing. This transportable water will
flow in the direction of the negative of the gradient of the modeled hydraulic potential.
In the routing model this potential is calculated by assuming that the transportable
subglacial water is at the overburden pressure [94]. Ultimately the
transportable water will reach the ice sheet grounding line or ice-free-land margin, at
which point it will be lost. The amount that is lost this way is reported to the user.
In this model tillwat also evolves by equation (27), but several
additional parameters are used in determining how the transportable water bwat flows
in the model; see Table 17. Specifically, the horizontal
subglacial water flux is determined by a generalized Darcy flux relation [18],
[95]
where \(\bq\) is the lateral water flux, \(W=\) bwat is the effective thickness of the
layer of transportable water, \(\psi\) is the hydraulic potential, and \(k\), \(\alpha\),
\(\beta\) are controllable parameters (Table 17).
In the routing model the hydraulic potential \(\psi\) is determined by
where \(P_o=\rho_i g H\) is the ice overburden pressure, \(g\) is gravity, \(\rho_i\) is ice density, \(\rho_w\) is fresh water density, \(H\) is ice thickness, and \(b\) is the bedrock elevation.
For most choices of the relevant parameters and most grid spacings, the routing model
is at least two orders of magnitude more expensive computationally than the null
model. This follows directly from the CFL-type time-step restriction on lateral flow of a
fluid with velocity on the order of centimeters to meters per second (i.e. the subglacial
liquid water bwat). (By comparison, much of PISM ice dynamics time-stepping is
controlled by the much slower velocity of the flowing ice.) Therefore the user should
start with short runs of order a few model years. We also recommend daily or even
hourly reporting for scalar and spatially-distributed time-series to see hydrology
model behavior, especially on fine grids (e.g. \(< 1\) km).
Option |
Description |
|---|---|
|
\(=k\) in formula (29). |
|
In the boundary strip water is removed and this is reported. This option specifies the width of this strip, which should typically be one or two grid cells. |
|
\(=\beta\) in formula (29). |
|
\(=\alpha\) in formula (29). |
| Previous | Up | Next |
