Rate-and-state friction explains glacier surge propagation

The incomplete understanding of glacier dynamics is a major source of uncertainty in assessments of sea-level rise from land-based ice. Through increased ice discharge into the oceans, accelerating glacier flow has the potential to considerably enhance expected sea-level change, well ahead of scenarios considered by the IPCC. Central in our incomplete understanding is the motion at the glacier bed, responsible for flow transients and instabilities involving switches from slow to fast flow. We introduce a rate-and-state framework for the transient evolution of basal shear stress, which we incorporate in glacier simulations. We demonstrate that a velocity-strengthening-weakening transition combined with a characteristic length scale for the opening of subglacial cavities is sufficient to reproduce several previously unexplained features of glacier surges. The rate-and-state framework opens for new ways to analyze, understand and predict transient glacier dynamics as well as to assess the stability of glaciers and ice caps.


SI Materials and Methods
We solve the ice flow using open source finite element code Elmer/Ice [1]. Including lateral friction, conservation of momentum yields where ρ is the density, g is the gravitational acceleration and σ is the Cauchy stress tensor Here, p is the pressure, I is the identity matrix,˙ is the strain rate tensor and ν is the effective viscosity where n = 3 is Glen's exponent,˙ e = tr(˙ 2 ) 2 and A is the Arrhenius factor. We assume that the ice is incompressible so that conservation of mass is given by In addition, we solve for the evolution of the free top surface elevation h through where SM B is the prescribed surface mass balance shown in Supplementary Figure 2.
For systematic studies of surge propagation we use a simplified approach for the basal hydrology. Instead of manipulating both the water pressure and the the peak friction coefficient C, we have modified C directly. This is mathematically equivalent to assuming that the water pressure p w is proportional to the normal stress σ n . Even though the value of the maximum basal shear stress C is usually not known, we can set up the system close to the stability threshold by modifying the water pressure. This can be written as p w = f σ n , we have that σ N = (1−f )σ n so that Here, we choose C * so that the system is close to the critical threshold, and so that the frictional characteristics are qualitatively similar to the inversion of basal friction by Jay-Allemand et al. [2]; C * (x) = 0.25e −x/2000m +0.08. We further neglect seasonal variations in basal water pressure, so that a surge is triggered by changes in τ b /σ N due to topography changes alone. The simulation is stopped when the surge reaches the terminus. The parameters used in the simulations are given in table 1.
We performed systematic simulations with q ∈ [1.5, 2.0] and d c ∈ [0.1, 10]m. The systematic runs with changed d c were restarted from simulations with d c = 1m just after surge onset to reduce computational cost.To measure the front velocity we used θ(x, t) and traced the treshold value θ threshold = 0.25. Front velocity could in principle be measured from most of the variables at the glacier bed, but θ exhibits the sharpest transition at the front tip (Supplementary Figure 6 and 7).
We also perform a second set of simulations using the double-continuum hydrology model for subglacial water pressure by de Fleurian et al. [3]. We refer to their paper for a detailed description of the model, but give a brief overview of the model here. The model assumes an inefficient drainage system (IDS) and an efficient equivalent porous layer (EPL). Initially the EPL is inactive. The model solves the vertically integrated Darcy's equations for a confined aquifer for both layers where Q is the source flux per unit surface, T is the transmissivity, S is the storage coefficient of the porous media, and h w is the water head. The storage coefficient is given by where β w is compressibility of water and β s is the compressibility of the sediment. The water pressure can be found from where ρ w is the water density. The EPL is activate once the water head in the IDS reaches a the ice overburden Then, the water transfer from the IDS to the EPL is then governed by a transfer function where d IDS is the thickness of the IDS layer, and l is a leakage length scale. Depending on wether the transfer is from the EPL to the IDS or from the IDS to the EPL, S and D with subscript IDS and EPL are used respectively. The parameters used are given in table 1. The system is initialized with transmissivity and source flux yielding parabolic water pressure distribution which is similar to what is found in Jay-Allemand et al. [2]. A steady state system with a constant source 0.01m yr −1 , and the surge is triggered by an increase the source flux to 0.03m yr −1 lasting three months. For the simulations we use a domain size of 25km and a minimum ice thickness of 1m to avoid zero thickness elements at the glacier front. The grid resolution is 10m in the horizontal direction and we extrude the mesh with 20 layers for the main figures and 50m and 10 layers for the systematic restarted runs for Supplementary Figure 6. The system is integrated forward in time with an adaptive forward Euler time-integration scheme.

Effective friction force with lateral friction
We add the friction force from lateral drag as an additional body force. Elmer/Ice has this feature available and it is documented in [4]. One has to supply a constant K resulting in the body force due to lateral friction. K can be (approximately) related to the glacier width 2W through under the assumption of no slip on the margins. Here, n is Glen's exponent, W is the glacier half width and A is the Arrhenius factor. For a glacier moving at high velocities, the velocity does not vary that much in the vertical direction, and we can treat the problem as onedimensional. Then, we can approximate the lateral drag as a friction coefficient: where φ is the glacier slope. The steady state friction coefficient is given by equation 2 With the one-dimensional approximation of the lateral friction, the effective steady state friction force is then given by where we have approximated cos(φ) 1. From twodimensional simulations, we can approximate the stress from the margins as a function of x as Linear stability The stability of a glacier within the rate-and-state framework introduced here is not simply given by a single point reaching the velocity weakening regime. For nonzero d c , a critical nucleation length appears from linear stability analysis in line with classical rate-and-state friction. However, in our case the nucleation length depends on the nonlinear viscous stresses rather than the shear modulus which is typically used in classical rateand-state theory. Ruina et al. [5] performed linear stability analysis of a generalized rate-and-state friction law on the form For such system, if we assume that a slider is pushed by a spring moving at a constant velocity with spring constant k, we have that the instability criterion is The stability criterion is most useful in terms of the spring constant where Note that for our sliding law, the instability can only occur for ∂τ b,ss ∂v < 0, i.e. in the velocity weakening regime. In addition, the spring constant has to be small enough. From k critical the usual approach is to approximate a critical nucleation length from the system shear modulus. Dietrich uses the relation k G L [6]. However, we are dealing with ice, which is best approximated as nonlinear viscous at long time scales. Then it is more appropriate to approximate k as where˙ ν is the viscous stress in a patch of size L in the velocity weakening regime. This gives a critical nucleation length The critical nucleation length depends on temperature through ice rheology, characteristic length scale d c , the velocity of the sliding patch, as well as the steady state basal shear stress parameters q, A s , σ N m and C through the velocity derivative of τ b,ss (equation 2). Although we do not explore the concept of a nucleation length in detail in this paper, it is important to bear in mind that the stability criterion for the onset of a glacier surge is not given by the sign of ∂τ b,ss ∂v alone; a patch of size L critical has to reach the velocity-weakening regime.

Scaling relation: Sliding velocity -front propagation speed
The propagation speed in a one-dimensional model is given by where the derivatives are taken at the front tip. ∂v/∂t can be found from the equation relating basal shear stress, velocity, and θ (equation 1), combined the state evolution law (equation 3).
which can be rewritten as At the tip, we can approximate τ b θ( v slip As ) 1/m . This gives us the relation for the velocity derivative We then insert for the state evolution law and find If the timescale of the velocity increase is short compared to the evolution of frictional strength, we can approximate θ 1. If we further assume that θ † is small behind the surge front, we end up with the proportionality Next, assume that the velocity gradient at the front tip can be approximated by the sliding velocity at the tip multiplied by a decaying function that f = f (x) that has dimension length −1 The scaling relation for the front propagation speed is then Even though the function f (x) is unknown, we can scale the sliding velocity with d c in order to collapse the data.

Scaling relation: Surge front ice thickness -Front propagation speed
The thickness change can be understood from a simple conservation of mass argument. During a surge, the velocity in the vertical direction is fairly constant, so we use a one-dimensional argument. Conservation of mass for a slice at position x yields where we consider the propagation to occur along the x-coordinate. Assume that dv x (x)/dx is nonzero with average value to v slip /∆x in spatiotemporal interval ∆x, ∆t where v front = ∆x/∆t, we have that which gives This means that the elevation change in the full system is governed by v slip /v front . Supplementary Figure 1: Initial glacier geometry used in the simulations in this paper. To model a realistic glacier geometry, we parametrize the Variegated glacier geometry data from 1973 [2].The bottom geometry follows y = 2000me −x/6700m while the thickness follows y = 61.5m + 0.066x − 3.6 × 10 −6 x 2 with a smoothed step function to avoid negative thickness at the front.