A unified model for transient subglacial water pressure and basal sliding

Abstract

Changes in water pressure at the beds of glaciers greatly modify their sliding rate, affecting rates of ice mass loss and sea level change. However, there is still no agreement about the physics of subglacial sliding or how water affects it. Here, we present a new simplified physical model for the effect of transient subglacial hydrology on basal ice velocity. This model assumes that a fraction of the glacier bed is connected by an active hydrologic system that, when averaged over an appropriate scale, is governed by two parameters with limited spatial variability. The sliding model is reminiscent of Budd's empirical sliding law but with fundamental differences including a dependence on the fractional area of the active hydrologic system. With periodic surface meltwater forcing, the model displays classic diffusion-wave behavior, with a downstream time lag and decay of subglacial water pressure perturbations. Testing the model against Greenland observations suggests that, despite its simplicity, it captures key features of observed proglacial discharges and ice velocities with reasonable physical parameter values. Given these encouraging findings, including this sliding model in predictive ice-sheet models may improve their ability to predict time-evolving velocities and associated sea level change and reduce the related uncertainties.

1. Introduction

Debate over what physics governs the sliding of glaciers over their bed and how this depends on subglacial water has existed since the first glacier sliding theories of Weertman (1957) and Lliboutry (1958). The elegant Weertman sliding law states that $u_{\rm b} = C\tau _{\rm b}^m$, where u _b is the basal velocity, C is a constant related to the roughness of the bed, τ _b is the basal shear stress and m is an exponent in the range of 2–4, related to Glen's flow law exponent n (Weertman, 1957). Weertman's theory predicts that sliding velocity is independent of the amount or pressure of subglacial water. In contrast, Lliboutry's more complex, but perhaps more observationally motivated theory, relies heavily on the role of water pressure. Lliboutry's theory is also more intricate, and testing it requires independent field observations of subglacial water pressure, which were exceedingly rare in Lliboutry's lifetime (Fowler, 2010) and are still limited today (Andrews and others, 2014).

About 60 years later, there is still no widespread agreement about the physics of subglacial sliding or how water affects it. That said, understanding of such processes has progressed significantly, particularly with the recognition that till sediments under high water pressure can fail in a standard solid friction Coulomb regime (Iverson and others, 1998; Tulaczyk and others, 2000; Tsai and others, 2015; Zoet and Iverson, 2020), and that sliding with subglacial cavities under high pressures may lead to higher velocities in a qualitatively similar manner (Iken, 1981; Fowler, 1987; Schoof, 2005; Gagliardini and others, 2007). However, the use of empirical sliding laws that incorporate water pressure in an ad hoc manner remains widespread (Budd and others, 1979; Schoof, 2010; Gladstone and others, 2017; Stearns and Van der Veen, 2018).

Observations have demonstrated that transient water input can drive significant increases in glacier surface velocity (Iken and others, 1983; Iken and Bindschadler, 1986; Zwally and others, 2002; Bartholomew and others, 2010; Andrews and others, 2014; Cowton and others, 2016; Smith and others, 2021), but also that sustained water input can lead to slower velocities, particularly later in the melt season (Tedstone and others, 2015). The way in which water inputs affect subglacial water pressure, and hence sliding velocities, is thought to be dependent on the hydraulic conductivity at glacier bed. Modeling of subglacial hydraulic conductivity has progressed significantly (e.g. Rothlisberger, 1972; Kamb, 1987; Clarke, 1996; Flowers and Clarke, 2002; Schoof, 2010; Schoof and others, 2012; Hewitt, 2013; Flowers, 2015), but it remains challenging to reconcile observations of basal water pressure with modeled values.

To help address the above shortcomings, particularly for short (e.g., diurnal or monthly) timescales, we introduce a simplified physical model for transient subglacial hydrology and its expected effect on subglacial sliding. Our model is similar in some ways to previous studies (Schoof, 2010; Schoof and others, 2012; Hewitt, 2013; Flowers, 2015) in its explicit inclusion and importance of spatial and temporal heterogeneity of the hydrologic system, and thus its effect on average sliding speeds, yet is simple enough that its few parameters are constrainable through sparse field observations. Application of the model to field observations of meltwater input, output and associated ice velocity changes from the Greenland ice sheet results in good fits when just a few physical parameters are adjusted within an expected range of reasonable values. The relative success of the two-parameter hydrologic model and simplified sliding law provides an encouraging path forward for use with large-scale predictive modeling efforts (e.g. Goelzer and others, 2020; Seroussi and others, 2020).

2. Model description

2.1. Sliding law with heterogeneous subglacial water pressure

The process by which subglacial water affects sliding is complicated and not agreed upon, but there are a few physical principles that are incontrovertible. First, whenever and wherever subglacial water pressure is locally equal to or above the ice overburden pressure, the glacier bed (whether hard bedded or till bedded) is expected to be fully supported by the water pressure and therefore unable to withstand any local shear stress. The simplest physical model that incorporates this idea is the Coulomb friction model that asserts τ _b ≤ f(σ − p), where τ _b is the basal shear stress, f is the friction coefficient, σ is the ice overburden pressure and p is the water pressure. Notably, when p = σ the local shear stress must be zero. Second, a number of observations (e.g. Andrews and others, 2014) strongly suggest that there is spatial heterogeneity in hydraulic conductivity and subglacial water pressures. It is thus useful to separately consider regions that have sufficient water pressure to be locally decoupled from the bed, for example in a Coulomb manner, and other regions of the bed where either low water pressures or significant bed roughness (which turns local normal stresses into larger-scale regional shear stress) prevents failure so that enhanced ice deformation or regelation govern ice motion, as in Weertman sliding. The simplest subglacial sliding law that incorporates the idea that basal shear stresses are limited either by Coulomb friction or by Weertman (enhanced ice deformation and regelation) mechanisms is

(1)$$\tau _{\rm b} = {\rm min} \big [ {{( {u_{\rm b}/C} ) }^{1/m}, \;f( {\sigma -p} ) } \big ] $$

as suggested by Tsai and others (2015) and Zoet and Iverson (2020), and which is similar to other suggestions (e.g. Brondex and others, 2017). Equation (1) specifically accounts for basal shear stresses being limited by the Coulomb limit in till and by the Weertman mechanisms immediately above that. This law, however, does not account for the limiting of Coulomb failure by bed roughness, which would imply an additional condition for switching from the Coulomb to the Weertman type of shear stress.

If one would like to be agnostic regarding the reasons for switching behavior, one could simply assume a given area of the bed to potentially have significant lowering of shear strength due to the Coulomb mechanism and that the shear stress is otherwise given by the Weertman law, regardless of the pressure at the bed. With this assumption, the law would be

(2)$$\tau _{\rm b} = \left\{{\matrix{ {\min [ {{( {u_{\rm b}/C} ) }^{1/m}, \;f( {\sigma -p} ) } ] } {{\rm \;\; if\;}x\in R_{{\rm hydro}}} \cr {{( {u_b/C} ) }^{1/m}{\rm \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;}} {{\rm otherwise}} \cr } } \right., \;$$

where x is the spatial location of the piece of bed of interest, and R _hydro is the region where the hydrological system is expected to be active. Here, we define the active hydrological system as any region where water pressure affects the local shear stress, including any such areas of the bed affected by hydrostatic or tidal ocean water pressure.

We note that Eqn (2) can be further simplified by considering that the hydrologic system is likely to persist only where the water pressure is close to the overburden pressure. If water pressures are significantly lower than overburden, then ice deformation will tend to close the hydrologic system (e.g. Nye, 1953; Rothlisberger, 1972) and make shear stress roughness (i.e. Weertman) dominated; if water pressure is significantly larger than overburden, then the ice would either lift up or deform until the water pressure reaches overburden (Iken, 1981; Tsai and Rice, 2010). Thus, a reasonable approximation is to simply set p ≈ σ over the region R _hydro or

(3)$$\tau _{\rm b} = \left\{{\matrix{ {0{\rm \;\;\;\;\;\;\;\;\;\;\;\;\;\;if\;}x\in R_{{\rm hydro}}, \;} \cr {{( {u_{\rm b}/C} ) }^{1/m}{\rm \;\;\;\;\;\;otherwise} .} \cr } } \right.$$

As defined in this way, the area represented by R _hydro should include all subglacial cavities, subglacial channels and also any till-covered areas of the bed that have reached failure due to high water pressure, but should not include ice–bed coupled areas between cavities. If the total regional bed area of interest is A ₀ and the area of R _hydro is A, then a simple regional force balance over A ₀ (e.g. Cuffey and Paterson, 2010) implies that the basal stress over the area unaffected by hydrology is increased relative to what it would have been otherwise (equal to the driving stress) by the ratio of the hydrologically affected area to the total area, or

(4)$$\tau _{{\rm avg}}A_0 = \tau _{\rm b}( {A_0-A} ) , \;$$

where τ _avg is the regionally averaged basal shear stress and τ _b in Eqn (4) refers to the basal stress over the area unaffected by hydrology. In the commonly applied shallow ice approximation (e.g. Cuffey and Paterson, 2010), this average basal shear stress must balance the driving stress, i.e. τ _avg = τ _d where τ _d = ρ _igHα is the driving stress, ρ _i is the ice density, g is the gravity, H is the ice thickness and α is the slope at the ice surface.

Finally, to account for a temporally variable hydrologic system, we would expect that increases in pressure would lead to an increased area for R _hydro, and the simplest approximation would be to assume a linear relationship between changes in area and pressure to accommodate instantaneous area changes as

(5)$$\Delta A = A-A_{{\rm ss}} = k_A\Delta p = k_A( {\,p-p_{{\rm ss}}} ) , \;$$

where A _ss is the steady-state area affected by hydrology, p _ss is the steady-state pressure (potentially equal to the local overburden σ = ρ _igH so that p _ss = σ), and k _A = ∂A/∂p is the proportionality constant linking changes in A to changes in p (see Appendix A). We recognize that there could potentially be a more complex relationship between A and p, particularly over long timescales or for large variations or with thresholds related to the bed roughness (see Appendix A). Continuing with Eqn (5) and substituting it along with Eqn (4) into Eqn (3) implies

(6)$$u_{\rm b} = C\left[{\displaystyle{{A_0\tau_{{\rm avg}}} \over {A_0-A_{{\rm ss}} + k_Ap_{{\rm ss}}-k_Ap}}} \right]^m.$$

Interestingly, this sliding law which results from assuming Weertman sliding on some area of the bed and a Coulomb limited shear stress of τ _b = 0 on the remaining area is reminiscent of Budd's classic empirical sliding ‘law’ (Budd and others, 1979, 1984) $u_{\rm b} = C\tau_{\rm b}^m /{( {\sigma -p} ) }^q$ with a number of important differences (see Appendix A). Equation (6) therefore provides the first physical justification of this type of sliding law that we are aware of; that is, Eqn (6) expresses a regional force balance rather than a local stress balance since it accounts for regions with two distinct types of basal properties. Nevertheless, Eqn (6) may be thought of as a ‘local’ sliding law in the sense that it may be appropriate once used as an average over a large enough area (e.g. 1 km²) to account for a statistically representative portion of the subglacial hydrologic system, yet a small enough area to be used in a spatially variable way as a sliding law in an ice-sheet model, much in the same way that the continuum hypothesis is utilized in classical continuum mechanics (Malvern, 1969).

2.2. Subglacial hydrology model

Equation (6) provides a prediction for the basal sliding velocity so long as the subglacial water pressure is known. We thus next turn to the question of how to estimate this pressure. Subglacial hydrologic modeling has a long history, and has gone through phases in which various authors have advocated for physics related to cavities, channels and porous media flow (Lliboutry, 1958; Boulton and Jones, 1979; Kamb, 1987; Flowers, 2015). Today, there appears to be reasonable agreement that mass conservation should be satisfied, and that the most active part of the hydrologic system is likely flowing turbulently (Rothlisberger, 1972; Tsai and Rice, 2010). Unfortunately, there is less agreement about equally important aspects of the physics, most notably whether storage of water is directly related to water pressure (e.g. Flowers and Clarke, 2002) or whether it is primarily determined by melting, viscous creep or sliding (e.g. Hewitt, 2013). For reasons detailed below, we believe all of these mechanisms to be physically reasonable, particularly when discussing hydrological changes on seasonal (including diurnal) timescales, and we describe a framework in which all of their impacts on basal sliding are accounted for and their relative magnitudes evaluated.

For turbulent, hydraulically rough flows, the Manning–Strickler approximation (Rothlisberger, 1972; Tsai and Rice, 2010) is a good approximation, where

(7a)$$-\displaystyle{{\partial p} \over {\partial x}} = \displaystyle{{\,f_0\rho a^{1/3}} \over {4h^{4/3}}}U^2$$

or

(7b)$$U = \displaystyle{2 \over {\,f_0^{1/2} \rho ^{1/2}a^{1/6}}}h^{2/3}\left({-\displaystyle{{\partial p} \over {\partial x}}} \right)^{1/2}, \;$$

where U is the average fluid velocity over the height of the conduit, f ₀ is a constant, ρ is the water density, a is the channel roughness length, h is the height of the conduit and ∂p/∂x is the pressure gradient. p is assumed to decrease as x increases; the sign of the left-hand side of Eqn (7a) becomes positive if the pressure gradient is reversed. We note that the Darcy–Weisbach equation is equivalent to the Manning–Strickler approximation when the expected h dependence of the friction factor in rough turbulent flow is accounted for (Tsai and Rice, 2010). Other turbulent flow laws could also be assumed (e.g. Schoof and others, 2012), and would lead to qualitatively similar model results. Assuming h is the characteristic (average) height of these hydrological features, and the total width of these flow features (e.g. within a given catchment) is W, the total flux is then given by

(8)$$Q = hWU = \displaystyle{2 \over {\,f_0^{1/2} \rho ^{1/2}a^{1/6}}}h^{5/3}W\left({-\displaystyle{{\partial p} \over {\partial x}}} \right)^{1/2}.$$

Equation (8) assumes an approximately rectangular cross section for subglacial flow, but a geometrical correction can be made for scenarios in which the cross section is more circular. Using a single cross-sectional average value for U(x), h(x) and W(x) is also an approximation that may result in a factor of 2 or more difference in total flux, but Eqn (8) should represent the correct physical scalings for the active hydrologic system.

Many subglacial hydrologic models assume that the steady-state radius of conduits, and hence water storage, depends primarily on the hydraulic pressure gradient due to the turbulent melting caused, and that this melting and opening of conduits is balanced by creep closure. In addition to these changes in water storage caused by local melting, we argue that transient instantaneous changes in water storage are also caused by changes in water pressure itself. Given likely spatial heterogeneity of the active subglacial hydrological system, we expect the area changes described by Eqn (5) (i.e. widening of basal channels proportional to water pressure) to result in a component of water storage that is proportional to instantaneous water pressure. This type of storage change is standard in groundwater hydrology (e.g. Fitts, 2013), where it sometimes has a very different (e.g. poroelastic or unconfined porous medium) interpretation, but can nonetheless be described by the same proportionality and has previously been used in a number of subglacial hydrology studies (Flowers, 2015). Additional instantaneous separation of ice from its bed would cause additional contribution to instantaneous storage from water pressure. Accounting for all three effects, and referencing pressure and melting to their steady-state values, then

(9)$$\displaystyle{{\partial S_{{\rm Tot}}} \over {\partial t}} = \displaystyle{{S\Delta V} \over {\rho g}}\displaystyle{{\partial p} \over {\partial t}} + \displaystyle{{\Delta V} \over \eta }( {\,p-p_{{\rm ss}}} ) + ( {M-M_{{\rm ss}}} ) \Delta V, \;$$

where S _Tot is the total storage volume, ΔV is the unit volume of this storage, S is called specific storage in the groundwater literature (fractional additional water storage per unit change in hydraulic head) (Fitts, 2013), η is an effective viscosity that relates the pressure in excess of the steady-state pressure p _ss to the rate of change in total storage and M is the local melt rate of ice per unit volume, with steady-state value M _ss. Equation (9) depends on variability in p and M from their steady-state values, where these values are the mean states around which perturbations occur, and could refer to monthly or annual average values when considering shorter timescale variability. Transient changes expressed by Eqn (9) can be evaluated without knowing the steady-state values themselves, as long as the variability is known. When M is caused by melting from turbulent frictional heating, it is given by (ΔV/ρL _i)U( − (∂p/∂x)) where L _i is the latent heat of ice (e.g. Cuffey and Paterson, 2010; Tsai and Rice, 2010). We note that η could be nonlinear and have a pressure dependence, for example if the change in storage is due to nonlinear deformation of ice (e.g. Glen's flow law) (Hewitt, 2013). For simplicity, we are agnostic about pressure dependence of η and later assume it is independent of pressure. Surprisingly, despite ample discussion of the relative merits of including the three terms on the right-hand side of Eqn (9) (Hewitt, 2013; Flowers, 2015), to our knowledge these three terms have never been included together with an attempt to determine their relative importance.

Setting the local changes in water storage volume (Eqn (9)) equal to the amount of water change from flux in minus the flux out (Eqn (8)), from moulin delivery, and from melt gives

(10)$$\displaystyle{{\partial p( {x, \;t} ) } \over {\partial t}} + \varepsilon ( {\,p( {x, \;t} ) -p_{{\rm ss}}( x ) } ) + \displaystyle{{\Delta \rho g} \over {SA_{\rm s}}}( {M( {x, \;t} ) -M_{{\rm ss}}( x ) } ) = {-}\displaystyle{{\rho g} \over {SA_{\rm s}}}\displaystyle{{\partial Q( {x, \;t} ) } \over {\partial x}} + \displaystyle{{\rho g} \over {SA_{\rm s}}}q( {x, \;t} ) , \;$$

where A _s is the cross-sectional area through which flow occurs (A _s = hW for a perfectly rectangular flow opening), ɛ = ρg/(Sη) is a scaled inverse viscosity, Δρ = ρ − ρ _i, q is local moulin water delivery per unit length and Q is given by Eqn (8).

If there are time-dependent changes in ∂p/∂x away from a mean value, the flux will also fluctuate from its mean. For small perturbations, we can linearize Eqn (8) about its mean state so that

(11)$$Q-Q_{{\rm ss}}\approx \displaystyle{{\partial Q} \over {\partial ( {-( {\partial p/\partial x} ) } ) }}\left[{-\displaystyle{{\partial p} \over {\partial x}} + {\left({\displaystyle{{\partial p} \over {\partial x}}} \right)}_{{\rm ss}}} \right].$$

Although this approximation may become inaccurate as the perturbation grows, at least it should have the right physics encapsulated. Defining the coefficient k _Q = ∂Q/∂(−∂p/∂x), we can link this to terminology from the hydrological literature (e.g. Fitts, 2013) as k _Q = KA _s/(ρg) where K is the effective hydraulic conductivity (see Appendix B). For most of the ice sheet (except for very near the terminus), the ice thickness gradient is nearly linear, with surface slopes in the range of 0.1–3° (0.001–0.05 radians), resulting in a nearly constant steady-state pressure head gradient(−∂p/∂x)_ss/(ρg)in the range of 0.001–0.05. The coefficient k _Q for perturbations is 1/2 of the coefficient relating the steady-state flux Q _ss to the steady-state pressure gradient (∂p/∂x)_ss (see Appendix B), so that the observed Q _ss and (−∂p/∂x)_ss ≈ ρ _igH/L constrains k _Q, where H is the thickness of ice at the location of meltwater input (e.g. a moulin) and L is the horizontal distance to the terminus. We note that the pressure must be atmospheric at the glacier terminus (for a land-terminating glacier), and the steady-state pressure at the moulin must be close to the hydrostatic ice pressure if the moulin is to exist in steady state (e.g., Rothlisberger, 1972).

Substituting Eqn (11) into Eqn (10) produces a single equation governing p(x,t)

(12)$$\displaystyle{{\partial p} \over {\partial t}} = \kappa \displaystyle{{\partial ^2p} \over {\partial x^2}}-\varepsilon ( {\,p-p_{{\rm ss}}} ) + \hat{q}( x ) , \;$$

where κ = K/S is the hydraulic diffusivity, $\hat{q} = \kappa [ {( {q-q_{{\rm ss}}} ) - ( {\Delta \rho /\rho } ) ( {M-M_{{\rm ss}}} ) } ] /k_Q$ is the normalized net input of water (per unit length) and spatial variability in KA _s is assumed to be over a longer length scale than that of p (see Appendix B). When ɛ = 0, Eqn (12) expresses the standard diffusion equation that is commonly appropriate in groundwater hydrology (Fitts, 2013), and used in some classic subglacial hydrology models (e.g. Boulton and Jones, 1979). When ɛ > 0, Eqn (12) expresses the idea that the rate of water storage increase can be additionally caused by water pressure, as for example in Iken and others (1983) or the many other models in which creep contributes to closure or opening (e.g., Hewitt, 2013). A balance between the last two terms produces the standard Rothlisberger (1972) steady-state channel theory. If hydraulic properties are isotropic then the 2-D generalization of Eqn (12) is simply obtained by replacing the second derivative in x by the 2-D Laplacian.

Importantly, κ and ɛ have limited spatial variability (see Appendix B) and we assume them to be constant and unknown (to be determined). k _Q = LQ _ss/(2ρ _igH) and $\hat{q}$ are assumed known through observations. Boundary conditions (such as a known influx of water at a known moulin location, and atmospheric pressure at the terminus) are also assumed to be known. When boundary fluxes or $\hat{q}$ are unknown, proxies for such inputs such as surface runoff data could be used as long as they are appropriately corrected. Thus, unlike many other subglacial hydrology models which have numerous under-constrained parameters, there are only two free parameters in this hydrologic model.

Finally, we note that we have chosen to express Eqn (12) as an equation for water pressure rather than for average channel height or area as is more commonly done (Schoof and others, 2012; Flowers, 2015). This choice offers a number of advantages. First, most sliding laws, including Eqn (6) proposed here, depend explicitly on pressure rather than subglacial channel height, making it more useful as a direct output of the model. Second, it is unlikely that we will ever have direct measurements of the spatial (vertical and horizontal) extent of the subglacial hydrologic system, so it is advantageous for a model to be as agnostic about these unobservable features as possible. Although the present model still depends implicitly on some of these unobserved variables, those variables are encapsulated in just a few unknown parameters (κ, ɛ, k _Q, $\hat{q}$) unlike more traditional formulations. Relatedly, we have chosen to express all variability as departures from steady-state conditions. This eliminates the need for knowledge of steady-state values when making transient predictions of the model.

3. Results

3.1. Idealized synthetic examples

Before applying the model more generally, we first present a few idealized model results that enable a better understanding of the hydrologic part of the model (Eqn (12)). When ɛ = 0, Eqn (12) simplifies to the pure diffusion equation that is well known in groundwater hydrology (Fitts, 2013). If this is either forced with a periodic boundary condition (i.e. periodic meltwater input) or with periodic local input through $\hat{q}( x )$, the solution for p is known to be given by a diffusion wave (Mandelis, 2001). This is most easily seen by assuming the ansatz

(13)$$p( {x, \;t} ) = p( {k, \;\omega } ) {\rm e}^{i( {kx-\omega t} ) }, \;$$

where k is the wavenumber and ω is the frequency. Upon substitution into Eqn (12), this implies that k ² = iω/κ or $k = ( {1 + i} ) \sqrt {\omega /2\kappa }$. Substituting this back into Eqn (13) implies

(14)$$p( {x, \;t} ) = p( {k, \;\omega } ) {\rm e}^{-\sqrt {( {\omega /2\kappa } ) } x}{\rm e}^{i\left({\sqrt {( {\omega /2\kappa } ) } x-\omega t} \right)}.$$

In other words, p has a wave-like structure that propagates with speed $\sqrt {2\omega \kappa }$ but also decays with distance with decay length $\sqrt {2\kappa /\omega }$. To some degree, this decay and time-lagged structure has previously been modeled and understood to potentially exist for subglacial channels (e.g. Clarke, 1996), but to our knowledge the specific diffusion wave solution expressed in Eqn (14) has never been recognized, and is qualitatively analogous to the diffusion of heat into near-surface ice layers (Cuffey and Paterson, 2010).

The specific boundary value problem that corresponds with a known arbitrary input flux of water Q _in(t), a constant proglacial atmospheric pressure p _atm, and assuming no (minimal) net creation or input of subglacial water between the input and output ($\hat{q} = 0$) is

(15a)$$\displaystyle{{\partial p} \over {\partial t}} = \kappa \displaystyle{{\partial ^2p} \over {\partial x^2}}-\varepsilon ( {\,p-p_{{\rm ss}}} ) , \;$$

(15b)$$\displaystyle{{\partial p( {0, \;t} ) } \over {\partial x}} = {-}Q_{{\rm in}}( t ) /k_Q, \;$$

(15c)$$p( {L, \;t} ) = p_{{\rm atm}}, \;$$

(15d)$$p( {x, \;0} ) = p_{{\rm ss}}( x ) .$$

Solving Eqn (15) (see ‘Code availability’ section) with a diurnal sinusoidal meltwater input forcing Q _in = Q _ss + Q ₀sin(ωt) with ɛ = 0 and κ = 400 km² d⁻¹ yields the results for flux and pressure shown in Figures 1a, b, and has exactly the anticipated diffusion wave structure superimposed on a linear gradient in pressure. Also as expected, decreasing the diffusivity (Figs 1c, d) results in a slower wave speed but also a much faster decay of the signal so that for κ = 4 km² d⁻¹, there is virtually no evidence of the sinusoidal diurnal signal 10 km downstream or farther of the diurnal input. Increasing the timescale of the forcing (decreasing the frequency) results in a slower wave speed and slower decay (Figs 1e, f). When ɛ = 0, the steady-state p(0,t) phase lags Q _in by exactly one eighth of a period (0.125 d) (Appendix C), with pressures downstream becoming progressively lagged according to the diffusion wave speed $\sqrt {2\omega \kappa }$. Averaged over infinite distance, the average phase lag is exactly an additional one-eighth of a period (Appendix C). Thus, with ɛ = 0, the minimum phase lag is one-eighth of a period (0.125 d = 3 h) and the maximum average phase lag is one-quarter of a period (0.25 d = 6 h).

Fig. 1. Result of solving Eqn (15) with sinusoidal meltwater input Q _in(t). Black dashed lines show the analytical solution of the wave speed for the diffusion wave model when ɛ = 0. (a) Flux (or scaled pressure gradient) and (b) pressure using κ = 400 km² d⁻¹, ɛ = 0, diurnal forcing. (c) Flux and (d) pressure using κ = 4 km² d⁻¹, ɛ = 0, diurnal forcing. (e) Flux and (f) pressure using κ = 400 km² d⁻¹, ɛ = 0 for a 7-day period forcing. (g) Flux and (h) pressure using κ = 400 km² d⁻¹, ɛ = 10 d⁻¹, diurnal forcing. In all panels, steady-state pressure p _ss(x) is assumed linear with downstream distance and proglacial pressure equal to p _atm, k _Q = 0.045 m³ s⁻¹/(kPa km⁻¹), L = 42 km, Q _ss = 18 m³ s⁻¹, and Q ₀ = 12 m³ s⁻¹.

Within reasonable parameter ranges for κ and ɛ, the general solution to Eqn (15) is qualitatively similar to the solution for ɛ = 0 except with a different phase and amplitude decay (see Figs 1g, h), and is evaluated numerically with finite differences (see ‘Code availability’ section). As expected, including the ɛ term advances the signal phase relative to ɛ = 0 by up to one-eighth of a period, with the maximum phase advance having p(0,t) perfectly in phase with Q _in. This extremum behavior can be understood by considering the limit where the ∂p/∂t term can be ignored compared to the other terms; in this case, the solution for p is exponentially decaying in space with a length scale $\sqrt {\kappa /\varepsilon }$ and is always perfectly in phase with Q _in (i.e. an infinite apparent wave speed). We also note that all time lags are deterministic outputs of the model once κ and ɛ are fixed, unlike other hydrological models where time lags are adjustable with additional storage parameters (Clarke, 1996; Flowers and Clarke, 2002; Hewitt, 2013).

3.2. Comparison of the model to field observations from Greenland

To test the combined hydrologic and sliding model against observations, we compare its predictions with field observations from Smith and others (2021) in southwest Greenland. This study collected hourly in situ discharge measurements in a large supraglacial river (‘Rio Behar’), together with local ice surface velocities from GPS, just upstream of a major moulin (see Fig. 2). Downstream proglacial discharge (~50 times larger than Rio Behar) was simultaneously monitored at two permanent gaging stations (Rennermalm and others, 2013, 2017; van As and others, 2019) (Fig. 2). Local ice surface velocities may be modulated by other moulin inputs in addition to Rio Behar (and thus contributions to $\hat{q}$), and the proglacial discharge also incorporates many moulin inputs from a much larger drainage area. Nonetheless, it is useful to compare the predicted output from a model constrained solely by observed input water flux for Q _in and otherwise as described in Eqn (15). Given the known fluxes over the 1-week timescale of the observations, we expect the contribution to $\hat{q}$ from temporal perturbations in melting of channels to be negligible. To spin-up the model, we assume a diurnal fluctuation in Q _in of similar magnitude for 4 d prior to the first observation.

Fig. 2. Validation site in southwest Greenland (left). ① and ② are locations of Watson River proglacial and AK4 proglacial discharge gaging stations, respectively. ③ is the location where GPS measurements of ice surface motion and measurements of moulin supraglacial input were made (along supraglacial river Rio Behar). Acoustic Doppler Current Profiler (ADCP: right) measuring moulin input at ③ (Rio Behar). Elevation contours shown with white lines and surface velocities shown with red shading.

Solving Eqn (15) initially assuming κ = 600 km² d⁻¹ and ɛ = 0 gives results shown in Figure 3. The predicted proglacial output is qualitatively similar to the observed proglacial discharge, with a predicted decline of the perturbation amplitude from ~70% at the inlet to ~7% proglacially and time lag of the discharge of ~12 h. Comparing this to the observations by Smith and others (2021) of the perturbation amplitude declining from ~70 to 9%, the model is reasonable, especially considering that the proglacial discharge is affected by numerous other moulins other than the one fed by Rio Behar and thus should not be expected to be quantitatively matched by the model, even if the model framework were correct.

Fig. 3. Results modeled using observed Q _in, κ = 600 km² d⁻¹, ɛ = 0. (a) Fluxes at moulin input (blue) and proglacially (red), (b) observed (squares, with uncertainties bounded in gray; see Smith and others (2021)) and predicted (line) ice velocities at the moulin location using the sliding law of Eqn (6) with τ _avg = τ _d and best-fit parameters k _A/[ρ _igH(A ₀ − A _ss)] = 0.07 and m = 4.1, (c) fluxes in the hydrologic model domain, (d) pressures in the hydrologic model domain. As above, p _ss is assumed linear with downstream distance, L = 42 km, H = 934 m, ρ _i = 920 kg m⁻³.

A more accurate quantitative comparison can be expected for the predicted basal sliding velocity. Although the observed ice surface velocities have a component due to ice deformation that is not modeled here, the deformation velocity is expected to be relatively small in this region (Maier and others, 2019) and temporal variability in deformation velocity is expected to be negligible in comparison with temporal variability in sliding velocity. In the predictions shown, we have applied the sliding model (Eqn (6)) at the observation point, using the predicted subglacial pressure at that location of the hydrologic model domain, and have not solved for a fully consistent global ice flow model. Comparing this approximate prediction with the observations, there is a small but significant time lag of ~1.5 h between the observed velocities and the local velocity predicted at the moulin location with ɛ = 0 (see Fig. 3b). When ɛ = 0, there is exactly a 3-h (1/8-period) delay in the subglacial pressure at the input location relative to the input flux, with pressures downstream being further delayed by the diffusion wave speed $\sqrt {2\omega \kappa }$ (where ω is the forcing frequency, see Fig. 1). Thus, the observed lag in peak velocity of ~1.5 h with respect to the peak input flux suggests that a non-zero viscous Iken-type of opening should be used (ɛ > 0). On the other extreme, the limit of (∂p/∂t) = 0 is also clearly not correct either, since the observed time lag is non-zero and implies that the instantaneous storage term is also important. These conclusions differ from those of previous authors who have variously suggested that only one or the other kind of subglacial storage change needs to be accounted for (Hewitt, 2013; Flowers, 2015). We recall that for the 10-d timescales modeled here, the expected departure from steady-state production of meltwater is estimated to be small but that it would be increasingly important over longer seasonal and multi-annual timescales.

Given the uncertainties on the values of κ and ɛ, we can also determine the values that are able to best predict the observed velocities assuming both the hydrological model (Eqn (15)) and the sliding law (Eqn (6)) are correct. This exercise yields an estimate of κ = 1400 km² d⁻¹ and ɛ = 4 d⁻¹, along with sliding law parameters k _A/[ρ _igH(A ₀ − A _ss)] = 0.05 and m = 4.0, and the results for this model are shown in Figure 4. These results are only modestly different from those of Figure 3 (RMSE of 5.5 m d⁻¹ compared to 6.5 m d⁻¹ for Fig. 3), with the main difference being the 1.5-h advance of the phase of the predicted velocity, which fits the observations better. The proglacial discharge is also significantly less delayed, with a lag of ~4.5 h (instead of 12 h), and has a slightly larger perturbation amplitude of ~10% (instead of 7%). The proglacial timing being significantly different than observed is likely due to local supraglacial contributions to the proglacial flux.

Fig. 4. Same as Figure 3 but for results modeled using observed Q _in, κ = 1400km² d⁻¹, ɛ = 4 d⁻¹. (a) Fluxes at moulin input (blue) and proglacially (red), (b) observed (squares, with uncertainties bounded in gray) and predicted (line) ice velocities at the moulin location using the sliding law of Eqn (6) with best-fit parameters k _A/[ρ _igH(A ₀ − A _ss)] = 0.05 and m = 4.0, (c) fluxes in model domain, (d) pressures in model domain. As above, p _ss is assumed linear, L = 42 km, H = 934 m, ρ _i = 920 kg m⁻³.

3.3. Inferences from model parameter value ranges

Both the sliding law parameters (Eqn (6)) and the parameters of the hydrologic model (Eqn (12)) have physical meaning, and it is useful to compare the values above with expectations from other physical considerations. Such comparisons may help us understand how applicable the parameter values may be to other regions or time periods that have not been analyzed.

The best-fitting sliding law exponent m = 4.0–4.1 found in Figures 3 and 4 is close to the expected m = 2–4 from a Weertman type of sliding law (Weertman, 1957); this agreement provides additional evidence that the sliding law of Eqn (6) may be physically realistic. Fixing m = 3 or m = 1 yield inverted models with good fits and similar RMSEs (5.5 and 5.6 m a⁻¹, respectively), suggesting that tradeoffs between m and k_A cannot be fully resolved. The best-fitting values of k _A suggest that perturbations in pressure that are 10% of hydrostatic result in changes to the active hydrological area on the order of 0.5% of the total area. Although there are no good observations of how much the active hydrological area changes with time, the order of magnitude size of this area change is consistent with the expected order of magnitude of the water storage changes (Smith and others, 2021). Thus, the value of k _A is physically plausible (see Appendix A).

The values of the hydrologic parameters κ and ɛ provide constraints on the hydraulic conductivity (K), specific storage (S) and effective viscosity (η). If fluid velocities are in the m s⁻¹ range and the hydraulic gradient follows the hydrostatic gradient (~1 km/40 km), then K can be estimated to be ~40 m s⁻¹ near the moulin. Since κ = K/S, a value of κ in the range of 400–1400 km² d⁻¹ then implies that S is in the range of 0.002–0.009 m⁻¹. In other words, for a 1 m change in hydraulic head, water storage is expected to change instantaneously by 0.2–0.9%. This is consistent with observations that inferred fluctuations of hundreds of meters of hydraulic head result in multi-centimeter-scale vertical displacements (Smith and others, 2021). This estimated S along with the constraint ɛ ≤ 4 d⁻¹ suggests that η ≥ 10¹⁰ Pa s. This is consistent with effective viscosities from Glen's flow law of ~2 × 10¹¹ Pa s for a perturbation stress of 10% hydrostatic (~1 MPa) (Cuffey and Paterson, 2010), and therefore also lends credence to the physical suitability of the model.

Since the model assumes that temporal variability is in the form of perturbations from a steady state, it is important to remember that there are explicit and implicit assumptions made about steady-state parameter values such as steady-state u _b, Q _ss, p _ss, A _ss and M _ss and their interdependences. For situations with different steady-state regimes, parameter values may be significantly different and potentially time variable over longer timescales (e.g. since the steady-state flux Q _ss and melt rate M _ss vary seasonally), and there is no guarantee that the best-fitting parameters for one time period and location will be the best when applied to a different time period or location. Thus, we caution against long-term application of the model until parameter values can be estimated during all appropriate regimes.

Finally, despite the uncertainties on parameter values, we note that the model and our estimates of κ and ɛ have implications for longer-term meltwater flux variability and downstream time lags and dampening that should be expected of subglacial water flows. Importantly, the model structure implies that time lags and dampening are not constants determined solely by hydrologic parameters but instead depend crucially on the timescale of flux variability. With κ = 600km² d⁻¹ and ɛ = 0, the time lag for diurnal fluctuations is 0.27 h per km downstream distance and the decay length is $\sqrt {2\kappa /\omega } = 14$ km, whereas for week-long fluctuations the predicted time lag increases to 0.73 h per km downstream distance and the decay length to 37 km. Thus, longer timescale flux variability, like that observed by van As and others (2017) in Greenland, is expected to be delayed longer than shorter timescale variability (by a ratio of the square root of the timescale, Fig. 1). Pressures and ice velocity fluctuations are predicted to have the same feature. These conclusions suggest that time lag measurements like those of van As and others (2017) depend on fluctuation timescales in a predictable way, and are not directly comparable to tracer velocity measurements (e.g. Chandler and others, 2013).

4. Conclusions

We have found that physical considerations of the expected spatial and temporal variability of the subglacial hydrologic system lead naturally to a simple Budd-like sliding law that is linked to a physically based, two-parameter subglacial hydrologic model. Predictions from the new model agree well with observations of diurnally modulated hydrologic fluxes and ice velocities at a field site in southwest Greenland, and suggest that despite its simplicity, the model captures important physics describing diurnal fluctuations in ice speed. To best fit observations, the model requires subglacial water storage changes to have both an instantaneous pressure dependence as well as a viscous component, unlike previously proposed hydrologic models. The small number of parameters needed to constrain the model, and the fact that the ranges of these parameters are within the expected physical range, indicate that the model may be useful as the physical basis for subglacial parameterizations in predictive ice-sheet models.

Although further research is needed to validate the new model in different settings, potentially including situations with grounding line migration (e.g. Robel and others, 2017) and with better constraints on the subglacial hydrologic system (e.g. Rada and Schoof, 2018), the model accomplishes the goal of producing a straightforward physical prediction of sliding velocities. Furthermore, it can be forced with quantities (moulin fluxes) that are directly observable or easily estimated through proxies (such as modeled runoff data). As such, the model reconciles the Weertman and Lliboutry schools of thought, and in a manner that is less dependent on unknown parameters than other hydrologic models or sliding laws. The success of the new simplified physics-based model provides an encouraging path forward for use with large-scale predictive ice-sheet modeling efforts (e.g. Goelzer and others, 2020; Seroussi and others, 2020) and may thus eventually help reduce uncertainties in predictions of future sea level rise that have been limited by models with poorly constrained parameters.

Code availability

Code for running models is available at https://doi.org/10.5281/zenodo.5231534.

Data

All datasets used here are available in the Supporting information of Smith and others (2021), available from https://doi.org/10.1029/2020GL091418.