Glacier surges controlled by the close interplay between subglacial friction and 2 drainage 3

The fast flow of glaciers and ice sheets is largely influenced by friction at the ice8 bedrock interface, and our imperfect understanding of subglacial friction accounts for 9 one of the largest uncertainties in predictions of future sea-level rise [1]. Glacier mo10 tion ranges from slow creep to cyclic surge instabilities [2] and devastating glacier 11 collapse [3, 4] as well as continuously fast-flowing ice-streams. Glaciers dynamics also 12 exhibits seasonal velocity variations, up-glacier and down-glacier propagation of surges, 13 interrupted surges, as well as short-duration speed up events that do not develop into 14 surges [2, 5, 6]. Several aspects of this wide range of glacier dynamical behavior re15 main elusive. This knowledge gap highlights the crucial need of developing improved 16 descriptions of the physical processes that occur at the glacier bed as well as their 17 couplings. Here, we show that this wide range of sliding behavior can be understood 18 from the transient evolution of subglacial cavities and till porosity which results from 19 a feedback loop between subglacial drainage efficiency and friction. We find potential 20 for surging behavior at glaciers that exhibit low hydraulic conductivity at the base, 21 together with a weak increase in hydraulic conductivity with sliding, and where the 22 frictional response contains a transition to velocity weakening friction. This poten23 tial materializes if the local topography and surface mass balance create sufficiently 24 thick glaciers to shut down the conduit drainage system, and where the water input to 25 the glacier base is sufficiently high. Accurately accounting for feedback loops between 26 friction and drainage has the potential to improve surge understanding and future 27 assessments of hazard potential of glaciers. 28

FIG. 1. Sketch of the couplings between subglacial drainage and subglacial friction. Subglacial friction and distributed drainage evolve simultaneously through a common state parameter representing the degree of cavitation or till porosity. In addition to distributed drainage which has a strong coupling with friction, the drainage system also contains localized low pressure conduits. A number of feedback loops between friction and drainage ensure the emergence of cyclic surge instabilities. The methods section contains a detailed mathematical description of these couplings, with a visual description of the friction-drainage coupling in  of an efficient conduit system, and if the water input is sufficiently high, the only way for the glacier to evacuate excess 78 water is to increase the hydraulic gradient, which in turn leads to an increase in both hydraulic conductivity and an 79 increased sliding velocity. Since the change in state of the interface affects both drainage and friction simultaneously, 80 the way the glacier modifies its sliding speed can differ drastically depending on the precise interplay between water 81 transport at the base, the frictional response, the glacier shape, the distribution of water pressure, topography, and

87
This large variety of glacier velocity variations somewhat challenges the definition of glacier surges. In the following, 88 we define glacier surges as velocity increases that show a well-defined nucleation region as well as crack-like propagation 89 (i.e. increasing extent of the surging region) accompanied by a surface bulge. In our model, these surges are triggered 90 by a frictional instability at the glacier bed once the ratio of basal shear stress τ b to effective normal stress σ N exceeds 91 a critical value C in a sufficiently large region equivalent to a nucleation length L c in rate-and-state friction [27, 31- where the main mechanism is an increase in water pressure at the bed due to a gradual lowering of conduit drainage 97 efficiency with increasing glacier thickness (FIG. 2B, FIG. S3).

98
Surges that fall into category (1) are triggered by an increase in basal shear stress with time until the glacier surface 99 slope is sufficiently steep in a sufficiently large area of the glacier. This can in principle occur even at low water 100 pressures with an active conduit system, given that C is sufficiently low and the friction law contains a transition  Surges in category (2) require inefficient subglacial conduits, which we control by decreasing the channel constant 113 compared to scenario (1). Subglacial conduits open by melting and close by viscous creep. Since they are low pressure systems, conduits become progressively less efficient once the glacier is sufficiently thick. This gradual shut-down can manifest as a gradual increase in water pressure over several years, and the degree of subglacial cavitation changes 116 until the hydraulic conductivity in the distributed system increases to a value that allows the water to drain. In such 117 a case, frictional instability can be reached even for fairly large values of C, as along as water pressure is able to build the surge reaches the glacier front. This mechanism is similar to the hydraulic switch mechanism, although we stress 131 that we do not find a hard switch but rather a gradual conduit shutdown with increasing glacier thickness.

263
The SIA velocity is found assuming no sliding at the base and constant temperature. The stress balance where the 264 driving force due to gravity is balanced by shear deformation.
Here, y is the vertical coordinate ρ is the ice density, g is the gravitational acceleration, n is Glen's exponent, h is the 266 altitude, A is the ice rheology constant, which we consider to be constant, and b is the bedrock topography. Below 267 we will use the definition Sliding is accounted for by the SSA solution where τ b is the frictional stress explained below, and H is the thickness. In addition to the friction force at the base where d c is a length scale on the order of the characteristic cavity size, and t c is a characteristic time-scale expected 279 to be related to ice rheology and the typical cavity size, with a lower limit set by the Maxwell time of ice. The where A s is a friction prefactor and m is a creep exponent. In steady state, at the pressure melting point, the basal 284 shear stress is given by the parameterization by Gagliardini et al. [37] 285 where C determines the maximum friction threshold in steady state (Iken's bound), σ N is the effective normal stress, q 286 is an exponent that determines the degree of velocity weakening friction, χ = v C m σ m N As , and α = (q−1) q−1 q q . A transition 287 from velocity strengthening to velocity-weakening friction for hard bed glaciers has also been verified experimentally 288 [38]. The state parameter in steady state can be found by combining equation 9 and 10 and is given by For the normal stress we assume that the effective normal stress is positive, i.e. we do not solve for additional uplift 290 due to water pressure exceeding the overburden pressure. In determining the frictional stress, this approximation has 291 no practical significance as the limit of zero effective normal is τ b = 0.

292
A few comments on the rationale behind the rate-and-state formulation are in order. Equation 9 for θ = 1 is 293 the limit of the steady state solution when v → 0. In this limit, there is no direct normal stress dependency on 294 basal shear stress because there is no cavitation. Using equation 9 requires the assumption that the basal shear 295 stress depends on the configuration of subglacial cavities, which in steady state depends on normal stress and sliding 296 velocity. With the additional assumption that the low velocity limit is independent of normal stress because there is 297 no additional cavities forming when you increase the sliding speed, the natural conclusion is that for a step change 298 in velocity where cavities have no time to develop there will be a strain hardening following the form of the low 299 velocity limit. This will be modified by a factor in [0,1] because there is not necessarily full contact between ice and 300 bedrock in this case. The form of equation 9 has experimental support, and follows closely the observations by Zoet where φ is the hydraulic potential, and the hydraulic conductivity K is determined by the amount of cavity opening 313 through the empirical relation Here, K 0 is the background conductivity in the limit of zero cavity opening, K sheet is the maximum conductivity 315 when the cavities are fully open. f perc (θ) is a function accounting for drainage percolation that takes values between 316 0 and 1, i.e. the hydraulic conductivity increases rapidly then the cavities start to connect.

317
In order to perform simulations over thousands of years we use an approximate solution of conduits by assuming 318 that the water pressure in the conduits p w,c = 0. This assumption allows us to avoid solving for pressure variations 319 within the conduits, and instead solve for the water exchange between the conduits and the distributed drainage where k c is a channel constant, S is the cross-sectional area, φ c is the hydraulic potential in the sheet, and α c and where ∆ c is the distance between conduits, φ c = ρ w gb is the hydraulic potential in the conduits at zero water pressure, 325 n is Glen's exponent, σ n is the normal stress at the base, L is latent heat, and B is a ice rheology constant which 326 determines how fast conduits close. A correction term is to the water pressure needed because conduits contain a 327 solution where they grow indefinitely, draining water from the sheet at negative water pressures. To remove this 328 solution from the equations, we increase the closure term B of the conduits by a factor 10 3 when the water pressure 329 in the sheet is negative. The exchange source term between the two systems required to keep p w,c = 0 is then given Following [40] we solve for the hydraulic potential in a distributed drainage system where e v is an aquifer void ratio, h 0 is a thickness constant relating the state parameter to a physical sheet thickness,

333
Ψ source is a combined water source term from rainfall and surface melt reaching the bed, as well as melting due to 334 geothermal heat. Ψ melt is a frictional melt rate Ψ melt = U SSA τm Lρice .

335
It should be noted that the mathematical framework presented here is based on experimental and numerical evidence

347
Parameter sets are given in tables S2 and S1. Most simulations use the same topography given by The reason for this is to make sure the largest slope of the glacier bed more or less coincides with the maximum 351 surface slope, and that this occurs in the central part of the glacier. This design allows us to set the nucleation region 352 of the surge to the central part of the glacier.

353
The climatic forcing is set through the surface mass balance and the water input to the bed. The surface mass