Extending the lumped subglacial–englacial hydrology model of Bartholomaus and others (2011)

Bartholomaus and others (2011) model the subglacial and englacial hydrology of Kennicott Glacier, Alaska, USA, and the glacier’s response to a lateral lake outburst flood that occurs each summer, i.e. a jökulhlaup. The model uses observed input flux and highly simplified subglacial and englacial morphology to reproduce the hydrograph of the Kennicott River, other hydrographs and the glacier motion (sliding). The model is ‘lumped’, i.e. the entire hydrological system is represented by one cell, with the advantage that parameter identification is possible given the observed input and output fluxes.

The goal of this correspondence is to extend their model in three specific ways: (1) to show a (lumped) pressure equation which follows from their equations, (2) to show how to extend their model to the distributed flowline case, with attention to the form of the distributed pressure equation, and (3) to state a water-amount-vs-pressure relation that applies in the steady-state case of their model. Thus this correspondence is largely deductive, showing what holds in the Bartholomaus and others (2011) theory, but with a distributed model as an extension.

Their model (the ‘Bartholomaus model’) would seem to be significantly different from more recent distributed models, specifically the theory of Hewitt (2011), Hewitt and others (2012) and Schoof and others (2012). Some similarities exist among them, however. All of these theories describe the evolution of subglacial linked-cavity systems and include physical opening and closing processes for these cavities. On the other hand, the just-cited distributed theories are entirely subglacial, while the Bartholomaus model has both subglacial and englacial water storage.

The distributed version of the lumped pressure equation derivable in the Bartholomaus model (Eqn (8) below) is a parabolic regularization of the elliptic pressure equation in Schoof and others (2012). In other words, we show that the pressure equation in Schoof and others (2012) is the zero englacial porosity limit of the distributed version of the pressure equation implicit in the Bartholomaus model. A connection between the Bartholomaus model and distributed models is observed by Hewitt (2013), but this is limited to the addition of an englacial storage term in the mass conservation equation. Werder and others (2013) include a parabolic pressure equation, with englacial void ratio in the coefficient of the pressure rate term, similar to what is implicit in the Bartholomaus model. Exposing these connections of the Bartholomaus model to more recent literature is a motivation for this correspondence.

In the Bartholomaus model, the total volume of liquid water stored in the glacier is S(t), and this is split into englacial S _en(t) and subglacial S _sub(t) portions:

(1)

Mass conservation in the model is the simple statement (Bartholomaus and others, 2008)

(2)

In the Kennicott Glacier application, fluxes Q _in and Q _out are observed.

The subglacial cavities have geometry determined by bedrock bumps that have horizontal spacing λ_x, λ_y , height h and width w _c. These combine to give a dimensionless capacity parameter f = hw _c/(λ _x λ _y); the value f = 0.05 is used for Kennicott Glacier. Each cavity has cross-sectional area A _c(t) and thus volume w _c A _c. The glacier occupies a rectangle of dimensions L × W in the map-plane so that the number of cavities is ν = (LW) /(λ_x λ_y ). It follows that the subglacial storage volume is S _sub = (w _c A _c)v = fLWA _c/h, proportional to A _c in this lumped case.

Englacial water is assumed to fill a well-connected system of crevasses and moulins up to a level z _w(t) above the bedrock, in a system that has macroporosity ϕ, so the englacial storage is S _en = LWϕz_w. Denote the subglacial water pressure by P(t). In the Bartholomaus model, P is hydrostatic, so knowledge of P is equivalent to knowledge of englacial storage, because of the assumed efficient connection to the subglacial system; the englacial system acts as a piezometer. In fact, noting the relation S _en = LWϕz _w above, and denoting the density of fresh water by ρ _w and gravity by g, we have

(3)

Now, in the Bartholomaus model the cavity cross-sectional area A _c evolves by physical opening and closure processes. The rate of production of subglacial water through wall melt is denoted here simply by Z, with sign Z > 0 when cavities are enlarging. The melt rate Z is further parameterized in the Bartholomaus model, but the details are unimportant to the derivation here. Let C _c = (2A)/nⁿ , where A and n are the ice softness parameter and power in the Glen ice flow law (Glen, 1955), respectively. Denote the sliding speed by u _b and let P _o = ρ _i gH be the overburden pressure, where ρ _i is the density of ice. In these terms the cavity evolution equation is

(4)

The three terms on the right are opening by cavitation, opening by wall melt, and closure by creep, respectively (cf. equation (4) in Bartholomaus and others, 2011).

Equations (1–4) combine to give an evolution equation for the pressure, as follows. From Eqn (3) we can write the pressure rate of change dP/dt in terms of the englacial storage rate of change dS _en/dt. Then using the time-derivative of Eqn (1), and both Eqn (2) and the proportionality between S _sub and A _c, we can rewrite in terms of fluxes and cavity area:

(5)

Finally, Eqn (4) eliminates dA _c/dt:

(6)

Though it is not stated there, Eqn (6) follows from the equations in Bartholomaus and others (2011). However, equation (12) in Bartholomaus and others (2011) is a restriction of Eqn (6) to the case of no creep closure and no input/output fluxes.

Furthermore, Eqn (6) suggests how to extend the Bartholomaus model from ‘lumped’ to ‘distributed’. Consider a one-dimensional glacier under which the water is flowing in the positive x direction. Recalling W is the transverse width, replace L by ∆x, the along-flow length of one cell in a finite-difference or finite-volume scheme. Then Eqn (6) becomes

(7)

which still describes the pressure in one cell. Because Q _in is the upstream input flux into the cell while Q _out is the downstream output, the continuum limit of Eqn (7) is therefore clear, with (Q _out – Q _in)/ ∆x →∂Q/∂x as ∆x → 0. Thus a distributed flowline form of the Bartholomaus model includes the partial differential equation

(8)

A distributed extension of the Bartholomaus model is not viable without a Darcy or other expression for Q. Such a flux parameterization was not needed in the Kennicott Glacier application because the input and output flux data form complete boundary conditions for a one-cell model. However, using any reasonable Darcy-type formulation for the flux Q, such as the power laws (2.10) of Schoof and others (2012), Eqn (8) becomes a nonlinear parabolic equation for the pressure. Furthermore, when using a Darcy flux expression, the ϕ → 0 (singular) limit of Eqn (8) is an elliptic equation for the water pressure, namely equation (2.12) in Schoof and others (2012).

Equation (8) should be solved in a time-dependent, coupled system with Eqn (4). The state space of the resulting model, i.e. the values that must be given as initial conditions for these evolution equations, is the pair of variables (A _c, P), where both A _c (equivalently, the subglacial water layer thickness η) and P depend on time t and space x.

Note that Eqn (8) should be solved subject to inequalities 0 ≤P ≤ (ρ _w/ρ _i)P _o. This is because the emergence of water at the surface of the glacier, from the efficiently connected englacial system, bounds the subglacial pressure. (Though Bartholomaus and others (2011) report a supraglacial geyser lifting water a few meters above the glacier surface, which violates the given upper bound, the use of the bound in a model supposes, reasonably, that such over-pressure geysers are short-lived and small in magnitude relative to over-burden.) Compare the physical justification of related bounds 0 ≤P ≤P _o in Schoof and others (2012), in a model which includes the ϕ → 0 limit of Eqn (8).

An additional implication of the Bartholomaus model, again not stated in Bartholomaus and others (2011), concerns the pressure in steady state. Specifically, d/dt = 0 in Eqn (4) gives this relationship between pressure P and cavity area A _c:

(9)

In the Bartholomaus model, the cavities described by A _c are assumed to always be water-filled (cf. Schoof and others, 2012). Thus Eqn (9) is equivalent to saying that in steady state the pressure is a function P(η) of the area-averaged thickness η of the subglacial water layer. Equation (9) also suggests the behavior of coupled Eqns (4) and (6) (or (4) and (8) in the distri: buted case) at large englacial porosity values.

The above reasons might support the a priori inclusion of a water-amount-vs-pressure functional relationship like Eqn (9) in a model, for use only on longer-timescale questions, but this is certainly not recommended for numerical models which are intended to follow the timescales that drove the construction of the Bartholomaus model. However, a power law P = P _o(η/η _crit)^7/2 with η _crit a positive constant, is used by Flowers and Clarke (2002) in non-steady circumstances also. Equation (9) is not a power law, however, though it is at least increasing with η like the Flowers and Clarke (2002) relation.

It is interesting to observe from Eqn (9) that the steady water pressure does not depend on the englacial macro-porosity ϕ. Though the englacial pressure is parameterized by P = ρ _w gz _w in Bartholomaus and others (2011), its steady value is entirely determined by the balance between sliding, wall melt and creep closure in the subglacial system. Thus the englacial system is passive in determining subglacial steady state, while the value of is critical in determining timescales of the hydrological, and thus glacier sliding, response (Van Pelt, 2013).

A two-horizontal-dimension version of coupled equations (4) and (8), using a full-cavity assumption as in the Bartholomaus model, and using the general Darcy flux relation (2.10) from Schoof and others (2012), is numerically solved in the ‘distributed’ hydrology model code which is part of the Parallel Ice Sheet Model (PISM) in its 0.6 release (February 2014) (see Bueler and Van Pelt, 2014). The steady-state result (Eqn (9)) is not used in the PISM implementation, other than when a user supplies to PISM the initial values for water amount (i.e. A _c or η in the above equations) without supplying initial values for pressure; the result of Eqn (9) is as good an initial guess for pressure as any other in this data-deficient situation.

We have observed that implementing a parabolic pressure equation like Eqn (8), subject to either of the above bounds on pressure, is substantially simpler than numerically solving the variational inequality form of the corresponding elliptic pressure equation (Schoof and others, 2012). In particular, our work shows that physical pressure bounds can be preserved, a part of the linked-cavities model dismissed as ‘prohibitively expensive’ by Werder and others (2013), by explicit time-stepping of such a parabolic pressure equation (Bueler and Van Pelt, 2014).