A Consistent Framework for Coupling Basal Friction With Subglacial Hydrology on Hard-Bedded Glaciers.

Below hard-bedded glaciers, both basal friction and distributed subglacial drainage are thought to be controlled by a network of cavities. Previous coupled hydro-mechanical models, however, describe cavity-driven friction and hydraulic transmissivity independently, resulting in a physically inconsistent cavity evolution between the two components of the models. Here, we overcome this issue by describing the hydro-mechanical system using a common cavity-evolution description, that governs both transient friction and hydraulic transmissivity. We show that our coupling approach is superior to previous formulations in explaining a unique observation record of glacier sliding speed from the French Alps. We find that, at multi-day to multi-decadal timescales, sliding speed can be expressed as a direct function of basal shear stress and water discharge, without accounting for water pressure, which simply adjusts to maintain the cavitation ratio needed to accommodate the water supply.

Introduction

Flow variability of ice sheets and glaciers adds a large uncertainty to projections of their evolution and their future contribution to sea level rise (e.g., Mouginot et al., 2019; Ritz et al., 2015; Shepherd et al., 2019). Ice flow variability arises from the complex relationship between sliding speed, stress balance, water pressure, and temperature at the glacier base, all of which, in addition, depend on the properties of the substrate beneath the ice (Cuffey & Paterson, 2010). In particular, the difficulty in computing basal water pressure (e.g., Downs et al., 2018; Flowers, 2015) limits the predictive power of current ice sheet models (Ritz et al., 2015), and therefore the ability to project the future of the cryosphere under climate change.

Two‐way coupled models of ice flow and subglacial hydrology, in which sliding velocity has an effect on subglacial hydrology and vice‐versa (e.g., Hewitt, 2013; Hoffman & Price, 2014; Pimentel et al., 2010), provide useful tools to test the sensitivity of ice dynamics to melt water supply. These models are also needed to evaluate the subglacial hydrology and friction theories by confronting modeled with observed ice velocities (Brinkerhoff et al., 2021). Ice flow and subglacial hydrology models are usually linked by a friction law that relates water pressure, basal shear stress, and sliding velocity, and an equation linking the sliding speed to the efficiency of the distributed drainage system (e.g., Bueler & van Pelt, 2015; Gagliardini & Werder, 2018; Hoffman & Price, 2014). The distributed subglacial drainage system under hard‐bedded glaciers consists of a network of connected cavities (Iken & Bindschadler, 1986; Kamb, 1987; Walder, 1986) which, in current models, is represented by a macro‐porous sheet (Bueler & van Pelt, 2015; De Fleurian et al., 2014; Flowers, 2015; Hewitt, 2011; Schoof et al., 2012; Werder et al., 2013). It is commonly assumed that the thickness of this water sheet evolves through time according to a competition between the opening of cavities due to sliding over bedrock undulations and their closure by creep deformation (Schoof et al., 2012; Walder, 1986).

Although the water‐sheet thickness is linked to the cavitation process which, in turn, controls the basal friction (Fowler, 1986; Lliboutry, 1968), its evolution is typically considered independent of the friction law (Hewitt, 2013). Previous two‐way‐coupled modeling approaches use two distinct parameterizations to describe the influence of cavities on sliding speed and the water‐sheet thickness (e.g., Brinkerhoff et al., 2021; Bueler & van Pelt, 2015; Gagliardini & Werder, 2018; Hewitt, 2013; Hoffman & Price, 2014; Sommers et al., 2018). This leads to the simultaneous existence of two independent descriptions of the cavity size, and thus to inconsistent coupling. Moreover, sliding laws used in previous approaches neglect rate‐weakening friction (Fowler, 1986; Gagliardini et al., 2007; Helanow et al., 2020; Schoof, 2005) and assume steady‐state cavity behavior, a condition which may not be fulfilled given that the cavity evolution timescale (few days) is expected to be longer than the water pressure variations timescale (few hours).

In this study, we overcome the limitations listed above by developing a framework in which existing theories are coupled together such that the frictional state and the drainage efficiency are both controlled by a consistent transient behavior of cavities. We evaluate the model performance against a unique 28‐year record of sliding speed and water discharge from underneath the Argentière Glacier in the French Alps (Gimbert et al., 2021; Vincent & Moreau, 2016). This data set provides unprecedented constraints on subglacial parameters and allows characterizing the behavior of the hydro‐mechanical system from multi‐day to multi‐decadal time scales.

Materials and Methods

Unifying the Description of Drainage and Friction

Under a steady‐state situation, cavity geometry is at equilibrium with the sliding velocity u and the effective pressure (where and denote the ice and water pressure, respectively) such that basal shear stress is only a function of and N (Diego et al., 2022; Gagliardini et al., 2007; Schoof, 2005). This is however no longer the case in a transient situation, where cavity geometry does not necessarily have the time to fully adjust to changing sliding velocities and effective pressures. In this case, the friction law is expected to be of the form (Diego et al., 2022; Iken, 1981), where θ is a variable describing the cavity geometry. Although calculations of force balance at the sliding interface suggest that a transient sliding law should incorporate an instantaneous dependency on the effective pressure (Iken, 1981; Schoof, 2005); here, we neglect this aspect (as also done in Thøgersen et al. (2019) and Tsai et al. (2021)) and we use a friction law of the form . Following Thøgersen et al. (2019), we assume that θ is a dimensionless cavitation state that allows the transient friction law to be expressed as where A the Weertman friction coefficient (m yr−1 MPa−) and m an exponent. Note that we modified the original definition of by Thøgersen et al. (2019) to obtain tending to 0 when no cavitation occurs, and that the formulation in Equation (12), (13) is equivalent to equation 4 in Tsai et al. (2021), although the envisioned underlying physics and scales may differ. The cavitation state can be seen as a state variable in an analogy with the “rate and state” friction (Ruina, 1983), although here θ corresponds specifically to the cavitation process. θ monotonically increases with any variables describing cavity size (e.g., length or height), and is expected to be related to the water sheet thickness h considered in subglacial hydrology models, as described later in this section. As shown more specifically below, the particularity of the friction law in Equation (12), (13), in comparison to other “rate and state” approaches applied in glaciology (Goldberg et al., 2014; Lipovsky & Dunham, 2016; Minchew & Meyer, 2020; Zoet et al., 2021), is that its steady‐state form is equivalent to the friction law established by Gagliardini et al. (2007). In addition to including the current state knowledge on steady‐state, hard‐bed friction, the formulation in Equation (12), (13) is compatible with the experimental findings of Zoet et al. (2021) that a step increase in velocity leads first to a drag increase and then, is followed by a transient decrease due to cavities evolving toward a new steady‐state configuration.

Gagliardini et al. (2007) express steady friction through a single explicit dependency of bed shear stress on effective pressure as where C is a coefficient describing the maximum shear stress supported by the bedrock and q is an exponent. The steady‐state cavitation θ ∗ (Thøgersen et al., 2019) is obtained by combining Equations 1 and 2:

We note that at the transition to rate weakening (Gagliardini et al., 2007) and that the steady cavitation state at this transition can thus be expressed as

Adapting the formulation introduced by Schoof et al. (2012) for the evolution of the water sheet thickness to impose a steady cavitation state that equals θ ∗ as defined in Equation (1), (2), we obtain the following formulations for the opening and closing velocities of cavities: where l is the mean distance (m) between bedrock bumps responsible for the cavitation process and is the rate at which evolves (s−1). The temporal evolution of θ is thus described by

In this study, we use the subglacial hydrological model developed by Werder et al. (2013), which is mainly based on Schoof et al. (2012), with the modification that hydrological transmissivity is determined from the cavitation state defined by the friction law (Equation (1), (2)). The water sheet thickness h (m) and the sheet conductivity are thus expressed through their dependence on the common variable . We assume a direct relationship between and h of the form: with h the average bedrock bump height (m) and p an exponent. The evolution equation of h becomes: which is qualitatively similar to the one used in previous approaches (e.g., Schoof et al., 2012; Werder et al., 2013), although it involves additional exponents. We define as well the sheet conductivity k as a function of or h as where is the intrinsic sheet conductivity and p 2 is an exponent.

Deriving a Discharge‐Driven Sliding Law at Steady‐State

Assuming that the hydraulic potential gradient is constant with time (Shreve's [1972] approximation with constant surface slope and ice thickness gradient), it is possible to derive a steady‐state relationship between subglacial discharge and effective pressure (Delaney et al., 2019; Hewitt & Fowler, 2008; Schoof, 2010; Walder & Fowler, 1994) or cavitation state in our case. This can be used to evaluate sliding velocities without solving the hydro‐mechanical problem.

Following Schoof et al. (2012) combined with Equations 8 and 10 to express the sheet thickness and conductivity as a function of , the discharge in the cavity network (m2 s−1) can be expressed as a function of as: where is the hydraulic potential (MPa), (α , β ) are constant exponents and bold variables are vectors. Combining Equation (1), (2) with Equation (1), (2), we obtain the following sliding law as a function of water discharge: where . This sliding law can be viewed as a Weertman‐type law for which the friction coefficient is a function of water discharge. We note that this approximation is expected to hold only over time spans that are long enough for the cavities to adjust and that discharge should be averaged over that same time spans (one to few days).

Instead of using our modified Equations 1 and 11, a steady‐state expression can be derived from formulations by Schoof et al. (2012) and Gagliardini et al. (2007) (Equation (1), (11)) by combining the relationship given by a steady‐state water sheet thickness, the friction law and the sheet discharge (see Supporting Information). We obtain the following relationship: where and is a constant setting the characteristic time of cavity closure. Both Equations 12 and 13 are in fact equivalent to the steady‐state friction law (Equation 2) but use a different estimation of the effective pressure based on two different evolution equations of the sheet thickness. Equation (1), (11) differs from Equation (1), (11) in the way the power exponent on is linked to . Implications on the value of and the description of hydrologic transmissivity in models are discussed in the results section.

Numerical Model

We implement the coupled problem of ice flow and subglacial hydrology in the finite element tool Elmer/Ice (Gagliardini et al., 2013) which solves for both ice flow and hydrology models. The Stokes equation describing ice flow is solved together with the subglacial hydrology problem following the previous implementation of GlaDS (Werder et al., 2013) in Elmer/Ice (Gagliardini & Werder, 2018), in which we implemented the framework described in Section 2.1. The model is detailed in Supporting Information and all the variables and parameters are summarized in Tables S1 and S2 in Supporting Information S1.

Application to Argentière Glacier, French Alps

We take advantage of the unique simultaneous records of sliding speed, subglacial discharge, and ice thickness available from the Argentière Glacier (Mont Blanc Range, France) to evaluate the performance of our approach and quantify the unknown parameters. The topography of the Argentière Glacier has been continuously monitored since 1975 (Vincent et al., 2009) and the existence of subglacial infrastructure allows a continuous measurement of both subglacial discharge and glacier sliding speed (1990–2020) directly at the glacier base (Gimbert et al., 2021; Vincent & Moreau, 2016). The sliding speed is measured using an instrument called the “cavitometer,” which consists of a fixed wheel that rolls as the ice roof of a natural cavity slides above it. The sliding speed is measured with an accuracy of about 1 cm/day (Vincent & Moreau, 2016). Discharge and sliding speed have been recorded at a daily resolution until 2018, when the previous analog logging system was replaced by a digital solution, and the recording interval has been changed to 30 min.

To focus specifically on the interaction between sliding and hydrology and to limit contributions from other factors such as bed topography, we solve the numerical problem on a glacier with a simplified geometry, represented by a 0.6 × 5 km2 ice slab over a bedrock with uniform 6° slope representing the average situation of the Argentière Glacier. The surface topography is constructed using the “plastic approximation” from a given driving stress (Cuffey & Paterson, 2010) and evolves at a daily time scale following the reconstructed basal shear stress at the location of the sliding speed measurements (Gimbert et al., 2021) (see Text S1 and Figure S1 in Supporting Information S1). This methodology allows to impose basal shear stresses as expected in the real setup even though the modeled geometry is simplified.

Results

Steady‐State Cavitation

The novelty of our study is that the hydraulic transmissivity evolves according to the friction law through its link with the cavitation state (Equations 8 and 10). The sheet thickness and the hydraulic conductivity thus increase only when cavitation is affecting the sliding velocity (see Figure 1 for the steady‐state case ). In previous approaches (e.g., Bueler & van Pelt, 2015; Gagliardini & Werder, 2018; Hewitt, 2013; Hoffman & Price, 2014; Pimentel et al., 2010), independent parametrization of the friction law and the water sheet thickness evolution could lead to inconsistent configurations where the water‐sheet thickness is large while cavitation does not occur in the friction law (shown as the “Weertman range” in Figures 1 and 2). The additional constraints given by the coupling reduce the range of effective pressure that allows the water to drain. In particular, it links effective pressure to basal shear stress through imposing averaged effective pressure to be close to for the cavitation to occur and the water drainage to happen (Figures 1a and 1b). The distributed drainage system becomes able to accommodate the melt water supply only in a narrow range of effective pressure which is sensitive to the basal shear stress value.

Figure 1

Friction law and hydrological properties of the water‐sheet‐like drainage system at steady state as a function of sliding velocity and effective pressure determined from Equations 3, 8, and 10. (a) Normalized thickness of the water sheet underneath a glacier (contour lines and associated colored areas) and (b) hydraulic conductivity . The steady‐state friction law (Equation (3), (8), (10)) for basal shear stresses MPa is shown as black bold lines, with red stars showing the transition between rate‐strengthening and rate‐weakening. See Table S1 in Supporting Information S1 for parameter definitions and values.

Figure 2

Modeled and observed sliding velocity at Argentière Glacier, France. (a) Observed mean sliding velocity (normalized by ) for linearly spaced observed daily mean water discharge bins over the whole time series (blue line) and split by month (colored line). The black and yellow lines show predicted sliding velocity using Shreve's (1972) approximation and Equations 12 and 14, respectively. (b) Observed daily sliding velocities (orange line) compared to modeled sliding velocities using Shreve's (1972) approximation and Equation (3), (8), (10) (blue line) and using the newly introduced consistent hydro‐mechanical coupling (black line). (c) Same as (b) but compared to modeled sliding velocities using a standard coupling (Gagliardini & Werder, 2018) (green line). (d–f) Modeled as a function of measured daily velocities for the period 2000–2018 using Equation (3), (8), (10) (d), the consistent hydo‐mechanical coupling (e), and the standard coupling (f).

Evaluating the Discharge‐Driven Sliding Law Using Observations at the Argentière Glacier

Assuming that is proportional to the measured total discharge, Equation (3), (8), (10) provides a relationship between sliding speed and water discharge that can be compared to observations from the Argentière Glacier. Observed mean sliding velocity per given discharge intervals is well predicted by our framework providing that  = 7.2 (Figure 2a). The model also captures the intra‐monthly sliding variations during the melting period (Figures 2b and 2c), confirming that the equilibrium between discharge and the cavitation state is reached over a time scale of a few days (Bartholomaus et al., 2011). This shows that changes in cavitation state accommodate the variability in water supply at multi‐day and seasonal time scales, while changes in the hydraulic potential gradient, which are not accounted for here, likely only accommodate the variability in water supply at subdaily time scale. This supports that, on average, effective pressure and sliding velocity are driven by the cavitation state needed to accommodate water input and that the hydraulic potential gradient can be considered as constant in time. The relationship obtained from previous approaches (Equation (3), (8), (10)) shows similar agreement with data for  = 0.22 (Figure 2a). However, given that q ≥ 1, a value of ≥ 5 is needed to explain the observations. Such a high value is unlikely, given that ≈ 5/4 is appropriate for turbulent flow as described by the Darcy‐Weisbach law and ≈ 3 is appropriate for laminar sheet flow (Hewitt, 2013). We attribute the significant non‐linearity observed in the relationship between discharge and sliding velocity to the existence of a threshold cavitation state from which the cavities start to connect (large value of ) rather than a high value of . This result aligns with direct field observations of permeability acting as a binary switch, between connected and unconnected drainage networks (Andrews et al., 2014; Rada & Schoof, 2018).

Transient Coupled Numerical Model

Solving the full hydro‐mechanical coupled problem including channel dynamics and transient cavitation on an idealized setup allows to study daily to seasonal variations of sliding velocity that can be compared with the Argentière Glacier observations to calibrate unknown hydrological parameters and test the potential effect of channels. Parameter estimations are based on best fit between observed and modeled sliding velocity (see Supporting Information). In particular, this allows us to constrain separately the exponents p  = 0.6 and p  = 6.5 (Figure S2a in Supporting Information S1) and highlight the behavior of the cavity system regarding its hydrological properties. The low value of p indicates that the cavity volume increases as soon as the cavitation state increases while the high value of p indicates that the cavity network undergoes an abrupt transition in terms of connectivity as a certain threshold of cavitation state is reached (Figure S2b in Supporting Information S1).

Our coupled model captures remarkably well both the long‐term (multi‐decadal) and short‐term (multi‐day) variations in sliding velocity observed at the Argentière Glacier (Figures 2b and 2e and Figure S3a in Supporting Information S1). The model also accurately reproduces the average characteristics of the observed seasonal variability of sliding speed as a function of discharge (Figures S3b–S3d in Supporting Information S1). The fact that modeled sliding velocities remain close to those estimated using Equation (3), (8), (10) and Shreve's (1972) approximation during the melting period (Figure 2b) supports our suggestion that the system tends to adapt its hydraulic transmissivity through enhanced cavitation rather than by increasing pressure gradient to accommodate the water supply at multi‐days and seasonal timescales.

To compare our predictions with those using other approaches, we performed the same simulation using the standard coupling in GlaDS (Werder et al., 2013) as implemented in Elmer/Ice (Gagliardini & Werder, 2018) where the sheet thickness evolution and the resultant effective pressure are computed independently of the friction law (see Text S4 in Supporting Information S1). Results show that the best match with observations leads to a minimum RMSE of 5.2 cm day−1 whereas our approach gives a considerably lower RMSE of 3.1 cm day−1 (Figure 2f and Figure S4 in Supporting Information S1). In addition, best results using the standard coupling are obtained for unusual parameter values of and Pa−3 s−1 which are far from the more realistic values according to the Darcy‐Weisbach law and Pa−3 s−1 for temperate ice viscosity. In particular, the multi‐day variability is not well represented due to the water pressure variability being too high as a result of not being well regulated by the response of drainage efficiency to cavitation changes (Figure 2c). Even more importantly, the change in amplitude of the seasonal cycle in response to the multi‐decadal change in ice thickness is not well represented using GlaDS compared to using our present approach (Figure 2c). In both approaches, the cavitation state (our study) and the water sheet thickness (GlaDS) remain similar every summer because both are set by the water discharge. In our approach, this condition translates to the quantity being conserved, such that the ratio between winter and summer velocities varies proportionally to , consistent with observations (Gimbert et al., 2021). In contrast, the standard coupling fails to reproduce this behavior. With respect to the friction law only, using Equation (3), (8), (10) instead of a more conventional N‐dependent friction law (Gagliardini et al., 2007; Schoof, 2005) has the considerable advantage that the momentum balance can be solved regardless of the effective pressure value, since the law is rate strengthening at any given time. This allows to use rate‐weakening friction laws, and to solve the hydrological problem under high water input rates in which water pressure can significantly rise (even at overburden) during a short amount of time. However, at short timescale (subdaily), the data shows an instantaneous response of sliding speed to water pressure. This suggests that friction reacts to change in effective pressure at fixed cavitation ratio and that Equation (3), (8), (10) should include a dependency on effective pressure in order to properly capture subdaily sliding velocity variations.

Figure S3d in Supporting Information S1 highlights a seasonal hysteresis that is captured by our model due to seasonal glacier thickness variation and associated driving stress changes (Gimbert et al., 2021). The normalization of the sliding speed by removes this hysteresis (Figure 2a) and shows how the efficiency of the drainage at low water pressure in the Argentière Glacier does not increase during the melting season, as would be expected if channelization was occurring. This can be explained by a dominance of drainage through the cavity network which is not influenced by the development of the Röthlisberger‐channel (R‐channel) network. However, the simulation shows that, in winter, when the glacier slides in the Weertman range (no cavitation), the cavitation state becomes insignificant and the transmissivity of the water‐sheet‐like drainage vanishes, making the water to primarily drain through residual R‐channels. We find that these channels are in equilibrium with their discharge and impose a constant effective pressure independent of ice thickness in January/February/March (Figure S6 in Supporting Information S1). This result confirms the suggestion by Gimbert et al. (2021) that effective pressure is constant across winters.

Discussion and Conclusion

In this study, we introduce a missing link between modeling subglacial hydrology and glacier sliding by using the friction law to describe the basal hydraulic transmissivity. We demonstrate that, at the multi‐day time scale, sliding speed is set by the adjustment of the cavitation state needed to accommodate the water supply to the glacier bed. It follows that, at this time scale, water pressure is not a relevant variable to be estimated through complex hydrological models. Rather, melt‐driven subglacial discharge through the cavity system would be a more appropriate variable to predict changes in sliding speed. In the absence of an efficient drainage system, the analytical solution in Equation (3), (8), (10) enables estimating an upper bound for ice‐sheet and glacier acceleration associated with increase in melt water supply into the future. We show that the sensitivity of sliding speed to water discharge is a function of the maximum drainage capacity of the cavity network (m2 s−1) defined by and the exponent (Figure 3a). Assuming that the values of the terms and inferred in this study are representative of hard‐bedded glaciers in general and that the hydraulic potential gradient is dominated by the surface slope, we can determine this sensitivity as a function of surface slope only (Figure 3b). For example, given that surface slopes of the Greenland ice sheet are smaller than 2°, a surface runoff increase by a factor of 2–4 expected over the next century (Fettweis et al., 2013) would, according to our findings, lead to a maximum increase in sliding speed of about 10%–25% (Figure 3b).

Figure 3

Sliding speed sensitivity to discharge increase determined from Equation (3), (8), (10) using Shreve's (1972) approximation. (a) Sliding speed increase (%) in response to doubling water discharge as a function of the exponent and the maximum drainage capacity (see Table S2 in Supporting Information S1 for parameters definition). The black dot indicates the parameter set quantified for Argentière Glacier in this study. (b) Sliding speed increase (%) as a function of surface slope and discharge increase factor for the parameters determined at Argentière Glacier.

Interestingly, our findings also show that the equilibrium between drainage efficiency and cavitation state continues to hold even if the required cavitation state lies in the rate weakening range () (see Figure S7 in Supporting Information S1). In addition to the existence of this range being questioned for real bed topography (Helanow et al., 2021; Schoof, 2005), our results show that rate‐weakening, even if existing, would not be effective due to the hydrological feedback stabilizing the sliding speed. This is presumably a property of hard‐bedded glaciers where substantial cavitation can occur and the drainage efficiency can significantly increase with sliding speed. This would, in turn, explain the lack of surge observations of hard‐bedded glaciers (Cuffey & Paterson, 2010), such that opening of subglacial cavities would prevent sustaining high water pressure during the surge.

Our findings emphasize that glacier sliding is not only determined by a friction law but also by a “transmissivity law” that needs to be carefully determined and related to the friction state. We believe that the combination of these two laws, including an unstable friction branch, as shown here, is able to explain a wide range of glacier flow behaviors such as seasonal speed up, surge, or catastrophic detachment. Further study should examine how the parameters , , , and rely on bedrock geometry or the presence of sediment (soft bed) and investigate how much they can potentially vary from one place to another. A direct dependency of the transient friction law on effective pressure should also be introduced in order to study short‐term velocity changes due to diurnal surface melt variations, strong rain events, or lake drainages.

Supporting Information S1

Click here for additional data file.