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 [86].
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
horizontallytransporting 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 massconserving. 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 mapplane but much faster than 

A version of the 

A massconserving 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 Temperatureindex scheme computes runoff
using computed melt and the assumption that a fraction of this melt refreezes, all other
models assume that all melt turns into runoff.)
Option 
Description 


Specifies a NetCDF file which contains a timedependent field 

Maximum effective thickness for water stored in till. 

Water accumulates in the till at the basal melt rate 

Replace the conservationofenergy 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 [88]; see also
[10].
In particular, denoting tillwat
by \(W_{till}\), the tillstored 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
MohrCoulomb 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 griddependent.
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 slowlyvarying 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 steadystate subglacial water flux for technical details.
The massconserving model: hydrology routing
¶
In this model the water is conserved in the mapplane. Water does get put into the till,
with the same maximum value hydrology
.tillwat_max
, but excess water is
horizontallytransported. 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 [89]. Ultimately the
transportable water will reach the ice sheet grounding line or icefreeland 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],
[90]
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 CFLtype timestep 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 timestepping 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 spatiallydistributed timeseries 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 