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 (Reference Weertman1957) and Lliboutry (Reference Lliboutry1958). 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, Reference Weertman1957). 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, Reference Fowler2010) and are still limited today (Andrews and others, Reference Andrews2014).
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, Reference Iverson, Hooyer and Baker1998; Tulaczyk and others, Reference Tulaczyk, Kamb and Engelhardt2000; Tsai and others, Reference Tsai, Stewart and Thompson2015; Zoet and Iverson, Reference Zoet and Iverson2020), and that sliding with subglacial cavities under high pressures may lead to higher velocities in a qualitatively similar manner (Iken, Reference Iken1981; Fowler, Reference Fowler1987; Schoof, Reference Schoof2005; Gagliardini and others, Reference Gagliardini, Cohen, Raback and Zwinger2007). However, the use of empirical sliding laws that incorporate water pressure in an ad hoc manner remains widespread (Budd and others, Reference Budd, Keage and Blundy1979; Schoof, Reference Schoof2010; Gladstone and others, Reference Gladstone2017; Stearns and Van der Veen, Reference Stearns and van der Veen2018).
Observations have demonstrated that transient water input can drive significant increases in glacier surface velocity (Iken and others, Reference Iken, Rothlisberger, Flotron and Haeberli1983; Iken and Bindschadler, Reference Iken and Bindschadler1986; Zwally and others, Reference Zwally2002; Bartholomew and others, Reference Bartholomew2010; Andrews and others, Reference Andrews2014; Cowton and others, Reference Cowton, Nienow, Sole, Bartholomew and Mair2016; Smith and others, Reference Smith2021), but also that sustained water input can lead to slower velocities, particularly later in the melt season (Tedstone and others, Reference Tedstone2015). 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, Reference Rothlisberger1972; Kamb, Reference Kamb1987; Clarke, Reference Clarke1996; Flowers and Clarke, Reference Flowers and Clarke2002; Schoof, Reference Schoof2010; Schoof and others, Reference Schoof, Hewitt and Werder2012; Hewitt, Reference Hewitt2013; Flowers, Reference Flowers2015), 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, Reference Schoof2010; Schoof and others, Reference Schoof, Hewitt and Werder2012; Hewitt, Reference Hewitt2013; Flowers, Reference Flowers2015) 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, Reference Goelzer2020; Seroussi and others, Reference Seroussi2020).
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, Reference Andrews2014) 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
as suggested by Tsai and others (Reference Tsai, Stewart and Thompson2015) and Zoet and Iverson (Reference Zoet and Iverson2020), and which is similar to other suggestions (e.g. Brondex and others, Reference Brondex, Gagliardini, Gillet-Chaulet and Durand2017). 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
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, Reference Nye1953; Rothlisberger, Reference Rothlisberger1972) 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, Reference Iken1981; Tsai and Rice, Reference Tsai and Rice2010). Thus, a reasonable approximation is to simply set p ≈ σ over the region R hydro or
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 0 and the area of R hydro is A, then a simple regional force balance over A 0 (e.g. Cuffey and Paterson, Reference Cuffey and Paterson2010) 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
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, Reference Cuffey and Paterson2010), 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
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
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, Reference Budd, Keage and Blundy1979, Reference Budd, Jenssen and Smith1984) $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 km2) 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, Reference Malvern1969).
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, Reference Lliboutry1958; Boulton and Jones, Reference Boulton and Jones1979; Kamb, Reference Kamb1987; Flowers, Reference Flowers2015). 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, Reference Rothlisberger1972; Tsai and Rice, Reference Tsai and Rice2010). 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, Reference Flowers and Clarke2002) or whether it is primarily determined by melting, viscous creep or sliding (e.g. Hewitt, Reference Hewitt2013). 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, Reference Rothlisberger1972; Tsai and Rice, Reference Tsai and Rice2010) is a good approximation, where
or
where U is the average fluid velocity over the height of the conduit, f 0 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, Reference Tsai and Rice2010). Other turbulent flow laws could also be assumed (e.g. Schoof and others, Reference Schoof, Hewitt and Werder2012), 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
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, Reference Fitts2013), 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, Reference Flowers2015). 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
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, Reference Fitts2013), η 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, Reference Cuffey and Paterson2010; Tsai and Rice, Reference Tsai and Rice2010). 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, Reference Hewitt2013). 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, Reference Hewitt2013; Flowers, Reference Flowers2015), 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
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
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, Reference Fitts2013) 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, Reference Rothlisberger1972).
Substituting Eqn (11) into Eqn (10) produces a single equation governing p(x,t)
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, Reference Fitts2013), and used in some classic subglacial hydrology models (e.g. Boulton and Jones, Reference Boulton and Jones1979). 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 (Reference Iken, Rothlisberger, Flotron and Haeberli1983) or the many other models in which creep contributes to closure or opening (e.g., Hewitt, Reference Hewitt2013). A balance between the last two terms produces the standard Rothlisberger (Reference Rothlisberger1972) 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, Reference Schoof, Hewitt and Werder2012; Flowers, Reference Flowers2015). 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, Reference Fitts2013). 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, Reference Mandelis2001). This is most easily seen by assuming the ansatz
where k is the wavenumber and ω is the frequency. Upon substitution into Eqn (12), this implies that k 2 = iω/κ or $k = ( {1 + i} ) \sqrt {\omega /2\kappa }$
. Substituting this back into Eqn (13) implies
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, Reference Clarke1996), 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, Reference Cuffey and Paterson2010).
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
Solving Eqn (15) (see ‘Code availability’ section) with a diurnal sinusoidal meltwater input forcing Q in = Q ss + Q 0sin(ωt) with ɛ = 0 and κ = 400 km2 d−1 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 km2 d−1, 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).
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 km2 d−1, ɛ = 0, diurnal forcing. (c) Flux and (d) pressure using κ = 4 km2 d−1, ɛ = 0, diurnal forcing. (e) Flux and (f) pressure using κ = 400 km2 d−1, ɛ = 0 for a 7-day period forcing. (g) Flux and (h) pressure using κ = 400 km2 d−1, ɛ = 10 d−1, 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 m3 s−1/(kPa km−1), L = 42 km, Q ss = 18 m3 s−1, and Q 0 = 12 m3 s−1.
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, Reference Clarke1996; Flowers and Clarke, Reference Flowers and Clarke2002; Hewitt, Reference Hewitt2013).
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 (Reference Smith2021) 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, Reference Rennermalm2013, Reference Rennermalm2017; van As and others, Reference van As2019) (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.
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 km2 d−1 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 (Reference Smith2021) 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.
Results modeled using observed Q in, κ = 600 km2 d−1, ɛ = 0. (a) Fluxes at moulin input (blue) and proglacially (red), (b) observed (squares, with uncertainties bounded in gray; see Smith and others (Reference Smith2021)) 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 0 − 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−3.
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, Reference Maier, Humphrey, Harper and Meierbachtol2019) 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, Reference Hewitt2013; Flowers, Reference Flowers2015). 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 km2 d−1 and ɛ = 4 d−1, along with sliding law parameters k A/[ρ igH(A 0 − 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−1 compared to 6.5 m d−1 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.
Same as Figure 3 but for results modeled using observed Q in, κ = 1400km2 d−1, ɛ = 4 d−1. (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 0 − 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.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, Reference Weertman1957); 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−1, respectively), suggesting that tradeoffs between m and kA 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, Reference Smith2021). 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−1 range and the hydraulic gradient follows the hydrostatic gradient (~1 km/40 km), then K can be estimated to be ~40 m s−1 near the moulin. Since κ = K/S, a value of κ in the range of 400–1400 km2 d−1 then implies that S is in the range of 0.002–0.009 m−1. 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, Reference Smith2021). This estimated S along with the constraint ɛ ≤ 4 d−1 suggests that η ≥ 1010 Pa s. This is consistent with effective viscosities from Glen's flow law of ~2 × 1011 Pa s for a perturbation stress of 10% hydrostatic (~1 MPa) (Cuffey and Paterson, Reference Cuffey and Paterson2010), 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 κ = 600km2 d−1 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 (Reference van As2017) 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 (Reference van As2017) depend on fluctuation timescales in a predictable way, and are not directly comparable to tracer velocity measurements (e.g. Chandler and others, Reference Chandler2013).
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, Reference Robel, Tsai, Minchew and Simons2017) and with better constraints on the subglacial hydrologic system (e.g. Rada and Schoof, Reference Rada and Schoof2018), 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, Reference Goelzer2020; Seroussi and others, Reference Seroussi2020) 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 (Reference Smith2021), available from https://doi.org/10.1029/2020GL091418.
Acknowledgements
This research was partially funded by the NASA Cryospheric Science Program (grant 80NSSC19K0942) managed by Dr Thorsten Markus.
Author contributions
VC Tsai designed the study, developed the model and led the writing. LC Smith, AS Gardner and H Seroussi contributed to the data analysis and to editing of the paper.
Appendix A: Budd-type sliding law
Equation (6) is reminiscent of Budd's classic empirical sliding ‘law’ (Budd and others, Reference Budd, Keage and Blundy1979)
(where q is a fitting parameter) with a number of important differences. Although both have a −p in the denominator and thus have velocities that increase with increasing water pressure, in Eqn (6) the sensitivity of the pressure dependence is not in proportion to the degree away from hydrostatic, but instead is proportional to k A/(A 0 − A ss). Furthermore, whereas Budd's law has the (local) basal stress in the numerator, Eqn (6) has the regionally averaged basal stress (or driving stress, when the shallow ice approximation is used) in the numerator. We note that if the relationship between A and p was nonlinear (unlike Eqn (5)) then Eqn (6) would also be nonlinear, and if A was not algebraically related to p then one would not be able write out an algebraic formula for the sliding law as in Eqn (6). Finally, we reiterate that the pressure dependence in Eqn (6) is indirect, caused only by the expected change in area of the active hydrologic system by relatively modest changes in water pressure, and the pressure everywhere is itself assumed close enough to overburden for its changes to be irrelevant.
The assumed linear proportionality ΔA = k AΔp in Eqn (5) is an approximation that is expected to breakdown for large changes in pressure, e.g. when a different regime of bed roughness is encountered. If instead it were assumed that ΔA = f(Δp) where f is an arbitrary function that satisfies f(0) = 0, then one could Taylor expand f to obtain
so that the modified Eqn (6) would have nonlinear contributions from Δp in its denominator and the sliding law would then be more complex than the Budd-type of law. The Budd-type contribution, however, would still represent the first-order contribution for small perturbations in Δp.
The physical interpretation of k A is perhaps best understood as being related to the proportionality of storage with water pressure. For a model in which water storage is only proportional to water pressure (i.e. ɛ = 0 in Eqn (12)), the storage change can be separated into a portion due to changes in average height and a portion due to changes in average width of the original conduits. Changes in the width would naturally result in changes to A, with the simplest picture being that changes in A are equal to the average length of conduits multiplied by the change in width. Although it is beyond the scope of the current study to do so, one could thus attempt to predict the value of kA prognostically from a detailed subglacial hydrologic model. Since storage changes can potentially include a term that is not proportional to water pressure, one could envision a more complex version of Eqn (5) that is not algebraic but is expressed as a differential equation like in Eqn (12). A sliding law constructed with such a non-algebraic version of Eqn (5) would also not be algebraic and would instead have a time dependence that depends on the history of water pressure rather than the instantaneous water pressure. Although this is possible, we believe the current evidence does not require such a complex assumption, and we assume that A depends algebraically on p in the active hydrologic system.
Appendix B: Hydraulic conductivity
In standard hydrology notation, the conductivity K is the proportionality constant that links changes in hydraulic head gradients to changes in water flux per unit area. In our notation, and for a 1-D gradient, then
or k Q = (KA s/ρg). If Q is governed by Eqn (8), then one would not expect a linear scaling of Q with pressure gradient and thus the value of k Q may depend on the pressure gradient. Taking a partial derivative of Eqn (8) with respect to ∂p/∂x yields
so that
Unfortunately, most of the physical parameters in this equation are poorly constrained, and we therefore opt to determine k Q from observations of flux at a known pressure gradient, rather than attempting to determine it from first principles. To do this, we rewrite Eqn (8) in steady state as Q ss = 2k Q(−∂p/∂x)ss. We note that changes in flux also depend on changes to h and W in addition to changes in ∂p/∂x but that such flux changes are of second order given the expected magnitude of S (see Section 3.3) compared with the order-one fluctuations in ∂p/∂x (changes comparable to the steady state). Thus, we ignore these additional contributions to Q to concentrate on the most important contribution.
To simplify Eqn (12), we have approximated
This approximation is valid as long as the length scale of pressure fluctuations is significantly shorter than the length scale of fluctuations in KA s. For cases modeled in this study, this can be verified post facto; in more general cases, the un-approximated equation can be used at the expense of a more complex model.
Appendix C: Time lag of sinusoidally forced diffusion equation
If p(x, t) = ei(kx−ωt) with $k = ( {1 + i} ) \sqrt {\omega /2\kappa }$
, then
Comparing this with p(x,t), the 1 − i is at −45° in the complex plane so that the flux at x = 0 is phase advanced by one-eighth of a period with respect to p(0,t), the pressure at x = 0.
Since it is the average pressure over some area that is expected to be the direct cause of basal velocity, it is of interest to determine the phase lag of the average pressure perturbation. Given the calculation above, if only pressures near x = 0 are important for the measured velocity then one would expect a phase lag of one-eighth of a period. On the other extreme, if measured velocities are affected by all downstream pressures, then one can calculate the average phase delay as the phase of
The 1 + i is at 45° in the complex plane so that the average phase delay is one-eighth of a period lagged with respect to p(0,t). Thus, in total, the phase delay of the average pressure must be between 1/8 and ¼ of a period relative to the periodic flux forcing, depending on what fraction of the downstream area is responsible for the signal of interest. If ɛ > 0, then the phase delay can be reduced below 1/8 of a period.
Open access
