Regularized Coulomb Friction Laws for Ice Sheet Sliding: Application to Pine Island Glacier, Antarctica.

The choice of the best basal friction law to use in ice-sheet models remains a source of uncertainty in projections of sea level. The parameters in commonly used friction laws can produce a broad range of behavior and are poorly constrained. Here we use a time series of elevation and speed data to examine the simulated transient response of Pine Island Glacier, Antarctica, to a loss of basal traction as its grounding line retreats. We evaluate a variety of friction laws, which produces a diversity of responses, to determine which best reproduces the observed speedup when forced with the observed thinning. Forms of the commonly used power law friction provide much larger model-data disagreement than less commonly used regularized Coulomb friction in which cavitation effects yield an upper bound on basal friction. Thus, adoption of such friction laws could substantially improve the fidelity of large-scale simulations to determine future sea level.

Introduction

Models of ice sheet evolution depend strongly on their selection of the basal friction law (Brondex et al., 2017, 2019; Joughin, Smith, & Holland, 2010; Nias et al., 2018), in which basal shear stress is expressed as ,N), where is the speed at which ice slides over its bed. Here denotes the effective pressure at the ice‐bed interface, defined as the difference between ice‐overburden and basal‐water pressure. When this relation is rearranged to instead express is a function of , the result is often referred to as a sliding law (Cuffey & Paterson, 2010).

Many models employ a general power friction law of the form, , which provides no upper limit on basal resistance for finite . A special case is Weertman (1957) friction ( = 2–3; Figure 1), which is often assumed for ice flowing over hard bedrock. The theoretical underpinnings of the Weertman power law, however, do not include cases in which the basal water pressure is sufficient to cause widespread cavitation, where pockets of water occupy a substantial fraction of the bed. When this effect is included, initially increases with speed but eventually attains a peak value, after which it effectively remains constant (Figure 1) or declines with faster flow (Gagliardini et al., 2007; Iken, 1981; Schoof, 2005), which in effect acts as a “regularized” Coulomb friction law. For regions where the ice sheet rests on a weak bed, tests on recovered till samples indicate also Coulomb‐plastic behavior () (Tulaczyk et al., 2000). Thus, with appropriate selection of parameters, a regularized Coulomb friction law can be applied to both hard bedrock and weak till.

Figure 1

Basal shear stress, as a function of speed for the friction laws used in the simulations. For the purpose of illustration, the friction coefficients have been selected to yield for a speed of 300 m/year.

Ice‐sheet models often use power law exponents that range from m = 1 (“linear viscous”) to “Coulomb‐plastic”; Nowicki et al., 2013). Narrowing this range is critical to reducing uncertainty in sea‐level projections (Stocker et al., 2013). In practice, a friction law is usually chosen (e.g., Weertman with fixed m) and then tuned by adjusting the unknown friction coefficient, to yield a modeled flow velocity that best matches the observed velocity data (MacAyeal et al., 1995). In such cases, the dependence on is often subsumed into the sliding coefficient, which is solved for directly. While such methods applied to a single snapshot of velocity can determine , they provide no real insight into which is the optimal choice of friction law to apply. The friction law, however, helps regulate an ice stream's transient response to ungrounding. Thus, comparison of simulated transient behavior with the observed evolution can help select the appropriate friction law. For this study, we simulated the transient response of Pine Island Glacier (PIG) constrained by observations to evaluate a variety of friction laws.

Friction Law Parameterizations

We used a standard shallow‐shelf model with several friction law parameterizations (see supporting information Text S1; Joughin, Smith, & Holland, 2010). For a general power friction law that relates the basal shear stress components ( and ) to the components of the sliding velocity ( and ), we determined basal friction as To include the effects of cavitation as a regularized Coulomb friction law, following Schoof (2005), we used the relation In both sets of relations, we have subsumed the potentially nonlinear dependence on effective pressure, N, into the friction coefficients, which also include the properties of the bed. In the original formulation for regularized Coulomb friction, (Schoof, 2005), where Λ is a parameter with unknown value determined by the local properties of the bed. As with the friction coefficient, we have subsumed the role of effective pressure into . Unlike the spatially varying friction coefficient, for simplifying the implementation, we assume a constant value for  over the whole domain. In all of the experiments presented here, we used a value of = 3 for regularized Coulomb friction.

The decision to subsume the role of into model parameters is dictated in part by our lack of knowledge of basal water pressure. In the case of the sliding coefficient, the effect of on or is determined through inversion as described below. In regions without seasonal input of meltwater to the bed, we assume does not change with time. While this assumption may hold well in the interior where the ice overburden changes slowly, it is likely to break down as the ice column approaches flotation. To account for this effect, we applied a reduction so that and evolved smoothly toward zero at the point of flotation (Joughin, Smith, & Holland, 2010). This linear correction was only applied once the height‐above flotation, , dropped below a threshold, , which limited its effect to the region immediately above the grounding line. The approximate extent over which this weakening applies for a given can be assessed by examining Figure S2 to determine where in Figure S2. In addition to providing ad hoc correction for weakening of the bed as flotation was approached, this correction may have compensated for errors in ungrounding imposed by errors in the surface‐ and bed‐elevation data (Figure S2; Joughin, Smith, & Holland, 2010).

Study Area and Model Setup

In West Antarctica, PIG has sped up by nearly 75% since the 1970s (Joughin et al., 2003; Mouginot et al., 2014; Rignot et al., 2014), making it an ideal “natural laboratory” with which to evaluate friction laws. In particular, we have detailed velocity and geometry observations from 2002 to 2017, which includes its ~40% speedup from 2006 to 2009 (Joughin, Smith, & Holland, 2010; Mouginot et al., 2014). Figures 2 and S1 show the progression of thinning as the glacier evolved over this time period. As the glacier sped up moderately from 2002 to 2005 (Lee et al., 2012; Mouginot et al., 2014), thinning increased in a narrow band near the grounding line (the transition from grounded to floating ice) from ~2 to ~5 m/year (Figure 2a). As the speedup increased from 2006 to 2009 (Figure 3), the rapidly thinning region broadened, with rates of up to ~7.5 m/year. Once peak speeds leveled off at around 4,000 m/year from 2009 onward, the strong thinning near the grounding line declined rapidly as thinning diffused inland with rates of ~1 m/year observed at distances of 200 km from the grounding line (Figures 2 and S1). This fast reduction of near‐grounding‐line thinning agrees well with earlier model simulations and is important because it slows the rate of grounding‐line retreat and speedup, reducing the rate of further ice loss (Joughin, Smith, & Holland, 2010).

Figure 2

Model domain and geometry. (a) Elevation change rates, d, and (b) annual surface elevation, flotation height, bed elevation (Rignot et al., 2014), and ice‐shelf thickness. The color plotted along bed elevation profile shows . The black bar shows the area of grounding‐line retreat. The lower inset shows the profile location for both plots (white line) and the red (1996) and black (2009) lines show grounding‐line position. The lower inset speed contours show the outline of the ice stream, and the red box in the upper inset shows the location of the model domain. GL = grounding line.

Figure 3

Modeled response to ungrounding from 2002 to 2017 for various friction laws. Results with (a) m = 1 and h  = 140, (b) m = 3 and h  = 75, (c) m = 8 and h  = 43, and (d) regularized Coulomb friction with m = 3, u  = 300 m/year, and h  = 30 along the profile shown in Figure 2. Observed velocities (V2006, V2007, V2008, and V2017) as well quasi 2002 (Q2002) velocity are also shown. Each panel includes the model and data differences for 2017, for the full domain, the slow regions (<300 m/year), and the fast trunk (v > 300 m/year). The black bar indicates the extent of ungrounding over the simulation.

To better evaluate potential friction laws, we began by using a control‐method algorithm (MacAyeal et al., 1995) constrained by estimates (Text S1) of the 2002 velocity (Q2002 in Figure 2a) and surface (Figure 2b) and bed (Rignot et al., 2014) elevations to invert for the spatially varying coefficient, , for a linear‐viscous friction (m = 1) law in a depth‐averaged (shallow‐shelf) ice‐sheet model (MacAyeal, 1989). We used the resulting value of to determine the friction coefficients (e.g., ) for the other friction laws via a simple rescaling, which initializes all models to 2002 with only minor variation among initial states (Figure 3 and Text S1).

After initialization to the 2002 reference state, for each friction law, we stepped the model forward in time by updating it with the observed annual surface elevation (Figures 2 and S1). An advantage of this method over coevolving the velocity and ice‐sheet geometry is that it lessened the model‐data differences that might mask sliding‐related processes. At each time step, we ungrounded (set to zero) the points where the surface and bed elevations indicate flotation was reached (Text S1). Because the sensitivity of  to varies as a function of its corresponding exponent, a different time‐invariant scalar value of was used to best match each simulated 2017 velocity to the corresponding observation. Throughout each simulation, was held constant.

Results From Model Experiments

Figure 3 shows the simulated progression of speedup from 2002 to 2017 for power law friction with three different exponents and with regularized Coulomb friction. The largest model‐data disagreement occurs for the linear‐viscous case (m = 1), which is a result of this model's basal friction having the greatest sensitivity to speed so that most of the additional resistance occurs near the grounding line, limiting inland propagation of the speedup. For larger m values (Figures 3b and 3c), the reduced sensitivity to speed allows the speedup to migrate farther inland, producing the best agreement on the trunk for a power law with m = 8 (approaching Coulomb friction), which is consistent with other findings (Gillet‐Chaulet et al., 2016). Although the larger exponents provide better agreement on the fast‐flowing trunk (speeds >300 m/year), they produce poorer results for the slower‐moving regions. Regularized Coulomb friction (Figure 3d) produces the best results overall, including both the slow and fast regions.

Note that the depth‐averaged model does not include vertical shear, so in the slow‐moving interior, simulated sliding may emulate some portion of the motion that is actually due to internal deformation. In such cases, regularized Coulomb or m = 3 power law friction yields results that better match the data because their sensitivity to driving stress is similar to that of a model that includes simple shear in the ice column. While a higher‐order model might improve the overall results, it likely would not drastically alter the relative performance of the various friction laws.

Using our elevation time series, we can examine the individual contributions to the speedup from loss of basal traction and evolving geometry (Figure 4). We began with a simulation that maintains the 2002 surface and its influence on driving stress but otherwise allows to evolve as it would have due to the observed thinning and the friction law's dependence on speed. For this case, the speedup on the lower trunk and shelf more than doubled (cf. Figure 4a with Figure 3d) but fell of rapidly with distance from the grounding line region where the loss of traction occurred. By contrast, when we held the friction parameters constant (i.e., no reduction in due to thinning or flotation) but allowed the surface to evolve (Figure 4b), the speedup migrated inland due to steepening of interior slopes (Figure 4b; Joughin et al., 2003). Interestingly, in this experiment, although the ice‐shelf thinning over this time period would be expected to reduce buttressing and increase speed (Alley et al., 2005), there is instead a slowdown above and below the grounding line (Figure 4b). While thinning does reduce buttressing, the effect is more than offset by a reduction of the slope just upstream of the grounding line, causing the net slowdown. These results indicate that, rather than considering the effect of ice‐shelf melt and loss of buttressing in isolation (e.g., Christianson et al., 2016), it is important to consider such changes in the context of the full evolution of the coupled ice‐sheet‐ice‐shelf system.

Figure 4

Modeled response to ungrounding from 2002 to 2017 for regularized Coulomb friction parameters used in Figure 3d. Cases with (a) surface elevation maintained at the 2002 level (holdThick) but with ungrounding as determined by the observed elevation change and (b) no change in friction parameters (holdtaub) as ice near the grounding line approaches and reaches flotation. The black bar indicates the extent of ungrounding over the simulation. Model data differences are computed as described as in Figure 3.

Role of Basal Sliding in Controlling Ice Stream Flow

Although our study is confined to PIG, it is useful to consider the results in the context of other glaciers. For instance, a recent study suggests that exerts no control on glacier speed, finding instead that ~ with (Stearns & Van der Veen, 2018). In actuality, what their results show is that is insensitive to for fast flowing glaciers in Greenland, which is fully consistent with regularized Coulomb friction. This lack of sensitivity, however, does not indicate that exerts no influence on flow. On the contrary, while can be relatively insensitive to locally (Figure 1), this property means that a localized loss of traction must be balanced by redistributing stresses over a wider area. For example, the reduction in concentrated near the grounding line due to thinning and ungrounding represents just ~0.5% of the total traction integrated over the PIG basin but produces a large speedup that extends over a broad area, indicating a strong sensitivity of the overall system to . Because cavitation reduces the sensitivity to speed in fast‐flowing regions, the speedup tends to extend farther inland to slower flow areas where  increases more readily with speed (e.g., slower speeds in Figure 1).

The notion that speed depends solely on is based on several flawed arguments (Stearns & Van der Veen, 2018). First, in a regression of the speeds of 140 Greenlandic glaciers against , they found no apparent relationship between the two, which was taken to indicate that exerts no control. But since slope, thickness, bed depth, trough width, and presence or absence of sediments are contributing factors that vary from glacier to glacier, there is no reason to expect a clean relation between speed and . Second, is assumed to be a direct proxy for . Although limited in number, borehole observations generally contradict this assumption, indicating low (water pressure near overburden) and little correspondence to (Engelhardt & Kamb, 1997; Luthi et al., 2002; Meierbachtol et al., 2013; Ryser et al., 2014). In addition, for PIG (Figure 2b) and three major glaciers in Greenland (Shapero et al., 2016), inversions indicate extensive patches of weak bed that are indicative of low (<250 kPa; Tulaczyk et al., 2000), in locations where the implies much higher (>1 MPa). Finally, the purely ‐dependent model is inferred from an apparent relationship between and and, by extension, . For reasons related to continuity, however, this relation should exist regardless of the friction law or other process governing sliding (Text S1 and Figure S3). Thus, while such a relation exists, it does not form the basis for eliminating the friction law in ice‐flow models and assuming instead that speed is only a function of or , which also would imply sliding that speed was independent of driving stress (Minchew et al., 2019).

Influence of Effective Pressure

As noted above, we reduce the friction coefficient where the bed approaches flotation (linearly with for ). In some instances, this decline is similar to that incorporated in other models (see Text S1). For the simulation with regularized Coulomb friction, this scaling accounts for about two thirds of the reduction in , with the rest of the loss of traction due to full ungrounding. It is important to note that the regularized Coulomb model requires the smallest correction (see in Figure 3). While this simulated bed weakening compensates for uncertainty in the ungrounded extent, some portion likely accounts for actual bed weakening. This loss of traction may be related to changes in as the ice column nears flotation, which is physically plausible but not rigorously modeled in our simulation. Alternatively, bed weakening at high speeds occurs for some regularized Coulomb models (Gagliardini et al., 2007; Schoof, 2005), although not for the model we used. More data from speedups on other glaciers would better constrain the relative roles of each of these effects. Nonetheless, our present ad hoc weakening scheme works well for PIG, and it likely should generalize to other glaciers, perhaps with some degree of tuning.

Away from the grounding‐line region, the time‐invariant friction coefficient () in our model implicitly holds fixed with time, and at least at the decadal scale, it reproduces ice‐sheet behavior over the course of a 40% speedup. This finding along with similar results (Gillet‐Chaulet et al., 2016) indicates that the assumption of fixed and may be sufficient for century‐scale simulations. The evolution of is difficult to model as it is governed by melt volume and timing as well as poorly known bed properties (Schoof, 2010). Thus, the ability to represent by values fixed in time rather than a poorly constrained time‐varying model greatly simplifies the ability to make sea‐level projections for societally relevant scales.

The regularized Coulomb friction law ((2a)) is parameterized by a velocity threshold, , which governs the transition between sliding with and without cavitation. As described above, , but since we have no reliable knowledge of either or Λ, we have instead assumed an empirically selected fixed  m/year over the entire model domain. The transition to Coulomb friction demarcated by generally is well removed from the grounding line (see 300 m/year contour in Figure S1), making it relatively insensitive to changes in as the ice column approaches flotation.

While all the results presented here use 300 (see contour in Figure S1), we tested a range of values (e.g., 200–500 m/year; not shown) and found similar performance. The lack of sensitivity to the value of likely occurs because this range of speeds often tends to be restricted to narrow shear margins (Joughin, Smith, Howat, et al., 2010; Rignot et al., 2011), so that a change in would not drastically relocate the transition to Coulomb friction. This relative insensitivity suggests the parameters will transfer well to other glaciers. Further work with other glaciers is needed to determine the robustness of this conclusion.

Comparison With Other Regularized Friction Results

As originally formulated (Gagliardini et al., 2007; Schoof, 2005), regularized sliding laws depend on . Lacking reliable knowledge of , however, we subsume its effect into our choice of and by solving for . This approach relies on underlying assumptions that (i) Coulomb friction occurs in the model domain where the motion is fast () and (ii) that is relatively low in these regions to facilitate cavitation or till dilation. Our numerical experiments support our decision to use the threshold to separate the model domain in to Coulomb and Weertman regimes, and they are consistent with earlier work indicating more Coulomb‐like behavior for the fast regions of PIG (Gillet‐Chaulet et al., 2016).

Some other works that employ a regularized Coulomb friction law use an explicit representation of based on an equivalence with (Brondex et al., 2017, 2019; Nias et al., 2018), which may be problematic as described above (see also Text S1). When parameterized this way, the regions where Coulomb friction occurs are limited by geometrical constraints to the regions within ~1–15 km of the grounding line (e.g., < ~100 m in Figure S2) as described in Text S1. Thus, while the same basic friction law is used in both cases, these two implementations likely would yield substantively different results when used for simulations of future behavior.

Conclusions

While friction laws that include the effect of cavitation on sliding are not new (Gagliardini et al., 2007; Schoof, 2005), ice‐sheet models used for sea‐level projection have been slow to adopt them (Nowicki et al., 2013). The PIG results indicate that cavitation effects are important and a relatively simple friction law with cavitation outperforms more widely used models. To the extent that regularized Coulomb laws have been applied, as parameterized, they limit Coulomb friction effects to regions in the immediate vicinity of the grounding line (Brondex et al., 2019; Nias et al., 2018). While Coulomb behavior has been demonstrated empirically for till (Tulaczyk et al., 2000), the form of the regularized Coulomb friction law we have used was developed for sliding with cavitation over a hard bed. In the case of PIG, there are extensive regions of both strong and weak bed (Joughin et al., 2009). Thus, an advantage of regularized Coulomb models is that a single friction law can reproduce behavior on a glacier over the entire domain irrespective of whether sliding occurs with weak till or hard bedrock. While further research is needed to improve regularized Coulomb friction laws, including better understanding of near‐grounding line behavior, further validation and adoption of such models should lead to substantially improved estimates of future sea level.

Supporting Information S1

Click here for additional data file.