# Glacier Surges Controlled by the Close Interplay Between Subglacial Friction and Drainage

Kjetil Thøgersen, Adrien Gilbert, Coline Bouchayer, Thomas Vikhamar

Schuler

► To cite this version:

Kjetil Thøgersen, Adrien Gilbert, Coline Bouchayer, Thomas Vikhamar Schuler. Glacier Surges Controlled by the Close Interplay Between Subglacial Friction and Drainage. Journal of Geophysical Research: Earth Surface, 2024, 129 (10), 10.1029/2023JF007441. [https://doi.org/10.1029/2023JF007441](https://doi.org/10.1029/2023JF007441)

HAL Id: hal-04945336

https://hal.science/hal-04945336v1

Submitted on 21 Mar 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# JGR Earth Surface

## RESEARCH ARTICLE

10.1029/2023JF007441

**Key Points:**

* We present a model combining rate and state friction with subglacial hydrology processes that spontaneously generates glacier surges

* A wide range of glacier behavior emerges from the model, including periodic and quasi-periodic glacier surges

* Surges arise if the drainage system is unable to sufficiently drain water from the glacier base, triggering a frictional instability

**Supporting Information:**

Supporting Information may be found in the online version of this article.

**Correspondence to:**

T. V. Schuler,
t.v.schuler@geo.uio.no

**Citation:**

Thøgersen, K., Gilbert, A., Bouchayer, C. & Schuler, T. V. (2024). Glacier surges controlled by the close interplay between subglacial friction and drainage. *Journal of Geophysical Research: Earth Surface*, 129, e2023JF007441. https://doi.org/10.1029/2023JF007441

Received 13 SEP 2023
Accepted 9 SEP 2024

**Author Contributions:**

**Conceptualization:** Kjetil Thøgersen, Adrien Gilbert, Coline Bouchayer, Thomas Vikhamar Schuler

**Funding acquisition:** Thomas Vikhamar Schuler

**Investigation:** Kjetil Thøgersen, Adrien Gilbert, Coline Bouchayer, Thomas Vikhamar Schuler

**Methodology:** Kjetil Thøgersen, Adrien Gilbert, Coline Bouchayer

**Project administration:** Thomas Vikhamar Schuler

**Software:** Kjetil Thøgersen, Adrien Gilbert, Coline Bouchayer

**Supervision:** Thomas Vikhamar Schuler

**Validation:** Kjetil Thøgersen, Adrien Gilbert

© 2024. The Author(s).
This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited.

<figure><img src="image_2.png" /></figure>

## Glacier Surges Controlled by the Close Interplay Between Subglacial Friction and Drainage

Kjetil Thøgersen¹, Adrien Gilbert² ID, Coline Bouchayer¹ ID, and Thomas Vikhamar Schuler¹ ID

¹Departments of Physics and Geosciences, The Njord Centre, University of Oslo, Oslo, Norway, ²Université Grenoble Alpes, CNRS, IRD, Grenoble INP, IGE, Grenoble, France

**Abstract** The flow of glaciers and ice sheets is largely influenced by friction at the ice-bed interface that can trigger rapid changes in glacier motion ranging from seasonal velocity variations to cyclic surge instabilities or even devastating glacier collapse. This wide range of transient glacier dynamics is currently not captured by models, and its implications for long-term glacier evolution are uncertain. This highlights the need of developing improved descriptions for processes that occur at the glacier bed. Here, we present a model that describes the evolution of basal friction inspired by a "rate and state" approach, coupled to models of subglacial drainage and glacier flow, and investigate how these couplings affect the dynamics of glaciers. We show that a wide range of sliding behavior results from a feedback loop between subglacial drainage efficiency and friction which depends on the evolution of a frictional state that can be interpreted as the degree of cavitation or till porosity for hard and soft beds, respectively. In our simulations, we find that glaciers are susceptible to surging if they exhibit a transition to velocity weakening friction associated with a poor sensitivity of the drainage capacity to the frictional state. This potential materializes if the local topography and mass balance create the conditions for high water pressure to build up in an area sufficiently large to exceed a critical length. We advocate accounting for feedback loops between friction and drainage as a promising avenue for better understanding dynamical instabilities of glaciers and ice sheets.

**Plain Language Summary** Glaciers move through a combination of basal sliding and viscous deformation of the ice. Observed velocities range widely from glaciers creeping at a few meters per year, to continuously fast-flowing ice-streams at velocities of kilometers per year, all the way to devastating ice avalanches. Some glaciers can alter their velocity in quasi-periodic intervals, known as surge-type glaciers. These glaciers can move steadily for decades, then speed up for a short period of time by up to two orders of magnitude or more, before returning to quiescence. In this study, we show that a wide range of glacier behavior can be understood from the interaction between friction and the subglacial drainage system. Glacier surges are possible if the efficiency of the drainage system is not high enough to drain the water that enters and is produced at the glacier base. This allows water pressure to increase over time, which can potentially trigger a frictional instability. Depending on the configuration of the subglacial drainage and how it evolves in response to sliding, glaciers can be stable or surge-type, in addition to exhibiting several different classes of velocity variations.

## 1. Introduction

Glacier surges are characterized by cyclic velocity variations over orders of magnitude, with fast flow typically restricted to a few years followed by a longer period of quiescence (e.g., Truffer et al., 2021). These events feature propagation of surge fronts down-glacier, up-glacier or in both directions with an accompanying change in glacier thickness and advance of the glacier terminus (Altena et al., 2019; Cuffey & Paterson, 2010; Dunse et al., 2015; Truffer et al., 2021). The intervals between surges are often fairly regular (e.g., Eisen et al., 2001; Meier & Post, 1969), but highly irregular intervals have also been reported (e.g., Bhambri et al., 2017; Björnsson et al., 2003). Surging glaciers show strong geographical clustering and are most frequently found within an optimal climatic envelope in terms of temperature and precipitation, ranging from cold-dry to warm-humid conditions (Sevestre & Benn, 2015). Glacier surges are typically assumed to occur in glaciers based on deformable till, often referred to as soft-bedded glaciers (e.g., Björnsson et al., 2003; Hamilton & Dowdeswell, 1996; Jiskoot et al., 2000). Several studies (e.g., Harrison et al., 1986; Iken et al., 1983) also highlighted the diversity of transient glacier dynamics that are typically not characterized as glacier surges. These range from multi-year pulsing with smaller velocity increases than ordinary surges, to events that share characteristics with surge nucleation but are interrupted before they affect a large portion of the glacier (e.g., Cuffey & Paterson, 2010;
**Visualization:** Kjetil Thøgersen, Adrien Gilbert, Coline Bouchayer, Thomas Vikhamar Schuler

**Writing – original draft:** Kjetil Thøgersen, Adrien Gilbert, Coline Bouchayer, Thomas Vikhamar Schuler

**Writing – review & editing:** Adrien Gilbert, Coline Bouchayer, Thomas Vikhamar Schuler

Herreid & Truffer, 2016; Turrin et al., 2014). In light of this wide spectrum of velocity variations, even obtaining a robust definition of glacier surges can be difficult (Herreid & Truffer, 2016).

Historically, theories of glacier surges and other flow instabilities have been based on switches in the subglacial environment which can alter glacier velocity, and several surge theories exist (Benn et al., 2019; Fowler, 1987; Fowler et al., 2001; Minchew & Meyer, 2020; Thøgersen et al., 2019). The two basic ideas are that temperate glaciers can surge due to changes in the drainage system from efficient conduits to less efficient distributed drainage, increasing water pressure and causing a surge (e.g., Benn et al., 2019; Flowers, 2015; Fowler, 1987; Kamb, 1987; Kamb et al., 1985) and polythermal glaciers can surge due to changes in the thermal regime at the glacier base from cold to thermal allowing elevated sliding velocities (e.g., Benn et al., 2019; Fowler et al., 2001; Jiskoot et al., 2000). A recent non-spatial theoretical approach combining both previous theories focuses on the balance between enthalpy sources and sinks (Benn et al., 2019), and produces cyclic glacier surges. However, by design, non-spatial models do not predict spatial characteristics of glacier surges.

A usual assumption made in previous studies is that friction couples to drainage through the effective normal stress (Flowers, 2015; Iken, 1981; Kamb, 1987; Walder, 1986). However, there are strong indications that considering subglacial friction and drainage as two separate systems, that only interact through the effective normal stress, is at least in some cases not sufficient (Gilbert et al., 2022). Hard-bedded glaciers develop subglacial cavities with increasing sliding velocity, which is expected to significantly increase drainage efficiency (Gilbert et al., 2022; Hoffman & Price, 2014). Soft-bedded glaciers experience dilation of sediments with increasing basal motion that can alter the drainage efficiency through changes in porosity (Clarke, 1987). Decoupling at the ice-till interface has been suggested as a mechanism for fast glacier motion (Kjær et al., 2006) which equally would have a huge impact on drainage. Recent theoretical and numerical developments suggest a strong influence of till dilation at the onset of glacier surges (Minchew & Meyer, 2020), and that velocity-weakening friction is a key ingredient in capturing the propagation of glacier surges (Thøgersen et al., 2019).

Here, we present a glacier flow model integrating a consistent, two-way coupling between subglacial drainage and friction (Gilbert et al., 2022) for a physically based description of the spatio-temporal evolution of the frictional state. Our approach enables simulations of both normal behavior (i.e., slowly varying flow) as well as instabilities such as surge nucleation, propagation and rapid termination, associated glacier advance and geometrical changes, as well as cyclic surges, but also phenomena such as short-duration speed up events, pulsing events and interrupted surges.

## 2. Model Description

To investigate the origin of the large variety of observed glacier flow variations, we build on the rate-and-state formulation of Thøgersen et al. (2019), which was further developed by Gilbert et al. (2022), and introduce a fully coupled rate-and-state friction-drainage description for temperate glaciers, where both non-monotonic friction and multi-component subglacial drainage are interlinked through a shared state parameter. The mathematical framework presented here is based on experimental (Zoet & Iverson, 2015, 2016) and numerical evidence (Gagliardini et al., 2007) of temperate ice sliding over rigid bedrock. The model components and their interactions are illustrated in Figure 1. The coupling introduces multiple feedback loops between friction and drainage, diverging from the conventional approach that relies on unique and strictly increasing friction laws based on sliding velocity and effective normal stress.

Our goal is to set up a model with a minimum number of ingredients that can reproduce the large variety of observed glacier velocity variations. Exploiting simulations of synthetic glacier geometries over thousands of years, we then investigate glacier surges and speed-up events. Decreasing model complexity reduces the computational cost and allows us to perform a large number of simulations over long periods of model time. Here, we consider situations where the dynamical behavior is largely dominated by sliding, rendering accurate descriptions of ice deformation subordinate. Therefore, we accept a reduced accuracy in the description of ice deformation for the sake of computational efficiency.

### 2.1. Glacier Flow

To solve for ice deformation we use a flowline model and follow Bueler and Brown (2009) in combining the integrated shallow shelf approximation (SSA) and the shallow ice approximation (SIA). The SIA velocity is
<figure><img src="image_3.png" /><figcaption>**Figure 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 the till porosity (hard-bed or soft-bed scenarios, respectively). 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 model description section contains a detailed mathematical description of these couplings.</figcaption></figure>

found assuming no sliding at the base and constant temperature and the driving stress due to gravity being balanced by shear deformation. The vertically averaged SIA velocity reads:

$$u_{\text{SIA}} = -\frac{2A}{n+2}(\rho_{\text{i}}g)^n \left|\frac{\partial h}{\partial x}\right|^{n-1} H^{n+1} \frac{\partial h}{\partial x}. \quad (1)$$

Here, $\rho_i$ is the ice density, $g$ is the gravitational acceleration, $n$ is the exponent in Glen's flow law, $A$ is the ice rheology factor, which is assumed constant for temperate conditions, $h$ represents surface elevation, $H$ is the ice thickness and $x$ is the horizontal coordinate along the glacier. The mass continuity equation can then take the form of a diffusion equation using a non-linear diffusivity:

$$D_{\text{SIA}} = \frac{2A}{n+2} (\rho_{\text{i}} g)^n \left| \frac{\partial h}{\partial x} \right|^{n-1} H^{n+2}. \quad (2)$$

Sliding is accounted for by the SSA momentum equation

$$\frac{\partial}{\partial x} \left( 2A^{-1/n} H \left| \frac{\partial}{\partial x} u_{\text{SSA}} \right|^{\frac{1}{n}-1} \frac{\partial}{\partial x} u_{\text{SSA}} \right) - \tau_b - \tau_m - \rho_{\text{i}} g H \frac{\partial h}{\partial x} = 0 \quad (3)$$

where $u_{\text{SSA}}$ is the velocity of the shelf, $\tau_b$ is the basal shear stress (Equation 9) and $\tau_m$ is an additional frictional stress to approximate the drag from the margins (Equation 7).

The horizontal glacier velocity $U$ is obtained from a weighted combination of the SIA and SSA solutions:

$$U = f(|u_{\text{SSA}}|) u_{\text{SIA}} + (1 - f(|u_{\text{SSA}}|)) u_{\text{SSA}} \quad (4)$$

where we adopt the same weighting function as Bueler and Brown (2009)

$$f(|u|) = \frac{2}{\pi} \arctan \left( \frac{u^2}{100 \left[ \text{m/yr} \right]^2} \right). \quad (5)$$

For mass continuity we solve the advection-diffusion equation
$$ \frac{\partial h}{\partial t} = \frac{\partial}{\partial x} \left( f(|u_{SSA}|) D_{SIA} \frac{\partial h}{\partial x} \right) + \frac{\partial}{\partial x} \left( (1 - f(|u_{SSA}|)) u_{SSA} H \right) + \dot{b} \quad (6) $$

where $\dot{b}$ is the surface mass balance rate.

The lateral drag $\tau_m$ scales with ice thickness $H$ and is computed as a function of the glacier width $W$ and the glacier speed (taken as the weighted SSA velocity) using

$$ \tau_m = \frac{2H}{W} \left( \frac{5(1 - f(|u_{SSA}|))u_{SSA}}{AW} \right)^{\frac{1}{n}}, \quad (7) $$

a relation that has been derived by Nick et al. (2010) based on Van der Veen (1999, section 5.5).

## 2.2. Basal Friction

We base the friction law on the newly developed rate-and-state friction formulation for temperate beds (Thøgersen et al., 2019), coupled to the ice flow model described above (Section 2.1) by identifying the sliding velocity $v$ with $(1 - f(|u_{SSA}|))u_{SSA}$. We then couple basal friction with the subglacial drainage through the dimensionless state parameter $\theta$. The state of the interface, $\theta$, determines both the frictional stress and the hydraulic conductivity through transient opening of subglacial cavities and can be understood as the fraction of the interface area that is occupied by cavities. Consistent with the notation by Gilbert et al. (2022), it takes values between 0 and 1, corresponding to full contact and full cavitation, respectively. The state evolution with respect to a steady state $\theta^\dagger$ can be written as

$$ \frac{\partial \theta}{\partial t} = \frac{v}{d_c}(\theta^\dagger - \theta) + \frac{1}{t_c}(\theta^\dagger - \theta) \quad (8) $$

where $d_c$ is a length scale on the order of the characteristic cavity size, and $t_c$ is a characteristic time-scale expected to be related to ice rheology and cavity shape, with a lower limit set by the Maxwell time of ice. The length scale was used by Thøgersen et al. (2019), but in order to account for cavity opening and closure due to water pressure variations that could occur at low sliding speeds, a time-scale is needed. First, the basal shear stress is found through the state parameter:

$$ \tau_b = (1 - \theta) \left( \frac{v}{A_s} \right)^{1/m}. \quad (9) $$

where $A_s$ is a friction prefactor and the creep exponent is taken $m = 3$ (Gilbert et al., 2023). In steady state, at the pressure melting point, the basal shear stress is given by the parameterization by Gagliardini et al. (2007)

$$ \tau_{b,ss} = \sigma_N C \left( \frac{\chi}{1 + \alpha \chi^q} \right)^{1/m}, \quad (10) $$

where $C$ determines the maximum friction threshold in steady state, Iken's bound (Iken, 1981), $\sigma_N$ is the effective normal stress, defined as ice overburden minus water pressure, $q$ is an exponent that determines the degree of velocity weakening friction, $\chi = \frac{v}{C^{m} \sigma_{N}^{n} A_{s}}$, and $\alpha = \frac{(q-1)^{q-1}}{q^{q}}$. A transition from velocity strengthening to velocity-weakening friction for hard bed glaciers has been proposed by Fowler and Nye (1986) for a Newtonian rheology and has been experimentally confirmed by Zoet and Iverson (2015). The state parameter in steady state can be found by combining Equations 9 and 10 and is given by:
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Response of the rate-and-state friction law to step changes in effective normal stress. Markers are consistent throughout the panels. (a) Effective normal stress $\sigma_N$ as a function of time. (b) Basal shear stress $\tau_b$ as function of time. (c) Normalized hydraulic conductivity $K/K_0$ as a function of time. (d) Steady state basal shear stress $\tau_{b,ss}$ as a function of effective normal stress $\sigma_N$. (e) Steady state, normalized hydraulic conductivity $K_{ss}/K_0$ as a function of effective normal stress $\sigma_N$. Parameters used: $v = 50$ m yr$^{-1}$, $q = 2.5$, $m = 3$, $C = 0.4$, $d_c = 1$ m, $t_c = 10^6$ s, $A_s = 10^{-23}$ s$^{-1}$ Pa$^{-3}$, $K_0 = 10^{-8}$ m s$^{-1}$, $K_s = 5 \times 10^{-8}$ m s$^{-1}$.</figcaption></figure>

$$ \theta^* = 1 - \left( \frac{1}{1 + \alpha \chi^q} \right)^{1/m} \qquad (11) $$

For the normal stress we assume that the effective normal stress is positive, that is, we do not solve for additional uplift due to water pressure exceeding the overburden pressure. In determining the frictional stress, this approximation has no practical significance as the limit of zero effective normal stress is $\tau_b = 0$.

In Equation 9, the basal shear stress depends only on sliding velocity and the frictional state of the interface. Any instantaneous dependency on normal stress (Iken, 1981; Schoof, 2005) is therefore implicitly neglected and the basal shear stress reacts to a change in normal stress through its effect on the state variable only (Figure 2). This assumption has been shown to be valid at time scales longer than a few days (Gilbert et al., 2022) which is within the range of the time scales considered in this study. The instantaneous response to change in sliding velocity is strain hardening even though the steady-state may be in the rate-weakening regime (Figures 3b and 3d). This has been directly observed experimentally (Zoet & Iverson, 2016) and supports the form of Equation 9. The state evolution (Equation 8) is somewhat arbitrarily chosen but physical arguments can give us approximate values for the characteristic length and time-scales involved. In their experiments, Zoet et al. (2022) found $d_c$ corresponding

<figure><img src="image_5.png" /><figcaption>**Figure 3.** Behavior of the rate-and-state friction law in response to step changes in sliding speed. Markers are consistently colored throughout the panels. (a) Sliding speed $v$ as a function of time. (b) Basal shear stress, $\tau_b$ as function of time. (c) Normalized hydraulic conductivity, $K/K_0$ as a function of time. (d) Steady state basal shear stress, $\tau_{b,ss}$ as a function of sliding speed at steady state, $v_{ss}$. (e) Steady-state relation between normalized hydraulic conductivity $K_{ss}/K_0$ and sliding speed $v_{ss}$. Parameters used: $\sigma_N = 0.5$ MPa, $q = 2.5$, $m = 3$, $C = 0.4$, $d_c = 1$ m, $t_c = 10^6$ s, $A_s = 10^{-23}$ s$^{-1}$ Pa$^{-3}$, $K_0 = 10^{-8}$ m s$^{-1}$, $K_s = 5 \times 10^{-8}$ m s$^{-1}$.</figcaption></figure>
to the wavelength of bed obstacles, and Woodard et al. (2023) showed that there is ample power at wavelengths around 1–10 m when mapping previous glacier beds; our choice of values is well supported by these findings. Even though we consider it unlikely that $d_c$ and $t_c$ are truly constant across the entire ranges of $\theta$, $\sigma_N$ and $\nu$, treating them as constant is likely acceptable as a first order estimate.

## 2.3. Basal Drainage System

We picture drainage of water along the glacier bed in two communicating systems, a distributed system of interconnected cavities and another of arterial, low-pressure conduits. In our model, only the first one is coupled to the friction, whereas the subglacial conduit system is a highly localized feature that does not directly affect friction. Assuming laminar flow in the sheet, the water flux is given by

$$q_s(\phi, \theta) = -\frac{K(\theta)}{\rho_w g} \frac{\partial \phi}{\partial x} \quad (12)$$

where $\phi = \rho_w g z_b + p_w$ is the hydraulic potential with $z_b$ denoting bed elevation, $\rho_w$ the density of water and $p_w$ the water pressure (e.g., Flowers, 2015; Hewitt, 2011). The hydraulic sheet factor $K$ includes conductivity, the conversion from hydraulic potential to hydraulic head and an assumed sheet thickness (corresponding to $h_0$ used in Equation 17). $K$ is determined by the amount of cavity opening through the relation

$$K(\theta) = K_0 + K_s \theta^a f_{\text{perc}}(\theta). \quad (13)$$

Here, we use $a = 3$, $K_0$ is the background value in the limit of zero cavity opening, $K_s$ is the maximum when the cavities are fully open. $f_{\text{perc}}(\theta)$ is a function accounting for drainage percolation that takes values between 0 and 1, that is, the hydraulic conductivity increases rapidly when the cavities start to connect. The exponent $a = 3$ arises if we assume $\theta$ is proportional to the often used sheet thickness (e.g., Flowers, 2015). This function is subject to uncertainty, but we expect it can be better constrained in the future through experiments and simulations simultaneously measuring friction and distribution of cavities.

In order to perform simulations over thousands of years we need to use an approximate model of the subglacial conduit system. Solving for the conduit system water pressure requires rather small time-steps, so to keep the time-step sufficiently large we want to avoid this. In principle, this can be achieved through any assumption of how water pressure in the conduits relates to the water pressure in the distributed system, where the end members are equal water pressure in the two systems on one hand, and zero water pressure in the conduit system on the other hand. Using these assumptions, we simplify the problem by avoiding the need to solve for separate pressure in the conduit and the sheet. Instead, we focus on solving for the water exchange between the conduit and the distributed drainage system. This allows for a larger time-step in the simulations. While the reason for making this assumption is numerical efficiency, the conduits still exhibit the usual characteristics: (a) their geometry is subject to variations, (b) they grow with increasing water discharge, and (c) close by viscous creep. These two assumptions are referred to as sheet water pressure and zero water pressure (Table 1), and differ in the following way: using the sheet pressure assumption, the hydraulic potential in the conduit is $\phi_c = \rho_w g z_b + p_w$ and the effective normal stress is $\sigma_N = \rho_i g H - p_w$ while for the zero pressure assumption $\phi_c = \rho_w g z_b$ and $\sigma_N = \rho_i g H$. Discharge through conduits is given by Werder et al. (2013):

$$Q_c = -k_c S^{\alpha_c} \left| \frac{\partial \phi}{\partial s} \right|^{\beta_c - 2} \frac{\partial \phi}{\partial s} \quad (14)$$

where $k_c$ is a conduit constant, $S$ is the cross-sectional area, $\frac{\partial \phi}{\partial s}$ is the hydraulic gradient along the conduit axis $s$; taken as $\frac{\partial \phi}{\partial x}$ for our one-dimensional model and equal to the bed slope when assuming zero pressure. $\alpha_c$ and $\beta_c$ are exponents accounting for turbulent flow (Werder et al., 2013). The conduit cross-sectional area $S$ increases by melting, and closes by viscous creep (Flowers, 2015). Its evolution is given by

$$\frac{\partial S}{\partial t} = \frac{l_c q_s + Q_c}{\rho_i L} \left| \frac{\partial \phi_c}{\partial x} \right| - 2S \left( \frac{\sigma_N}{nB} \right)^n \quad (15)$$
<table><caption>**Table 1**
List of Parameters That Are Varying for Different Simulations, Organized in Sets Which Relate to Presented Simulations</caption><thead><tr><th>Parameter</th><th>Description</th><th>Units</th><th>1</th><th>2</th><th>3</th><th>4</th><th>5</th></tr></thead><tbody><tr><td>z<sub>b</sub>(x)</td><td>Bed topography</td><td>[m]</td><td>Equation 18</td><td>Equation 19</td><td>Equation 18</td><td>Equation 18</td><td>Equation 18</td></tr><tr><td>k<sub>c</sub></td><td>Channel conductivity</td><td>[m<sup>3/2</sup> kg<sup>-1/2</sup>]</td><td>0.01</td><td>0.15</td><td>0.01</td><td>0.001</td><td>0.003</td></tr><tr><td>C</td><td>Maximum friction threshold</td><td>–</td><td>0.3</td><td>0.14</td><td>0.4</td><td>0.4</td><td>0.4</td></tr><tr><td>Ψ<sub>source</sub></td><td>Water source term</td><td>[m s<sup>-1</sup>]</td><td>10<sup>-8</sup> sin<sup>2</sup>(πt/yr)</td><td>10<sup>-8</sup> sin<sup>2</sup>(πt/yr)</td><td>10<sup>-8</sup> sin<sup>2</sup>(πt/yr)</td><td>10<sup>-8</sup> sin<sup>2</sup>(πt/yr)</td><td>Ψ<sub>0</sub> sin<sup>2</sup>(πt/yr)</td></tr><tr><td>˙b(h)</td><td>Surface mass balance rate</td><td>[m yr<sup>-1</sup>]</td><td>4 × 10<sup>-3</sup>(h − 1600)</td><td>4 × 10<sup>-3</sup>(h − 1600)</td><td>4 × 10<sup>-3</sup>(h − 1600)</td><td>4 × 10<sup>-3</sup>(h − 1600)</td><td>∂b/∂z (h − 1600)</td></tr><tr><td>K<sub>s</sub></td><td>Sheet conductivity constant</td><td>[m<sup>4</sup> s kg<sup>-1</sup>]</td><td>5 × 10<sup>-8</sup></td><td>2 × 10<sup>-7</sup></td><td>1 × 10<sup>-6</sup></td><td>1 × 10<sup>-6</sup></td><td>5 × 10<sup>-8</sup></td></tr><tr><td>h<sub>0</sub></td><td>Characteristic drainage thickness</td><td>[m]</td><td>0.1</td><td>0.01</td><td>1.0</td><td>1.0</td><td>0.1</td></tr><tr><td>W</td><td>Glacier width</td><td>[m]</td><td>1,200</td><td>1,500</td><td>1,500</td><td>1,500</td><td>1,200</td></tr><tr><td>l<sub>c</sub></td><td>Sheet width contributing to channel growth</td><td>[m]</td><td>2</td><td>20</td><td>2</td><td>2</td><td>2</td></tr><tr><td>B</td><td>Conduit closure parameter</td><td>[s<sup>1/3</sup> Pa]</td><td>5.3 × 10<sup>7</sup></td><td>1.0 × 10<sup>8</sup></td><td>5.3 × 10<sup>7</sup></td><td>5.3 × 10<sup>7</sup></td><td>1.0 × 10<sup>8</sup></td></tr><tr><td>–</td><td>Conduit pressure assumption</td><td>–</td><td>Sheet pressure</td><td>Zero pressure</td><td>Sheet pressure</td><td>Sheet pressure</td><td>Zero pressure</td></tr><tr><td>–</td><td>Illustrated in Figure</td><td>–</td><td>4, 5B, 6, 7</td><td>5A, 6</td><td>8B</td><td>8A</td><td>9</td></tr></tbody></table>

where $l_c$ is the width of the sheet contributing dissipation energy to the conduit (Hewitt, 2011), $n$ is Glen's exponent, $L$ is the latent heat of fusion and $B$ is an ice rheology constant which determines how fast conduits close. We note that for the zero pressure assumption, the conduit geometry may become unstable, the consequences of which are discussed later (Section 4).

Assuming a completely water-filled drainage system, we determine the water exchange $\Psi^{szc}$ between the two subsystems from continuity

$$l_s \Psi^{szc} = \frac{\partial S}{\partial t} + \frac{\partial Q_c}{\partial x} - \frac{l_c q_s + Q_c}{\rho_w L} \left| \frac{\partial \phi_c}{\partial x} \right| \quad (16)$$

The direction of $\Psi^{szc}$ is from the sheet toward the channel and $l_s$ denotes the lateral conduit spacing to provide an appropriate length scale. Following Werder et al. (2013), we solve for the hydraulic potential in a distributed drainage system

$$\frac{e_v}{\rho_w g} \frac{\partial \phi}{\partial t} + \frac{\partial q_s}{\partial x} (\phi, \theta) + h_0 \frac{\partial \theta}{\partial t} - \Psi_{\text{source}} - \Psi_{\text{melt}} + \Psi^{szc} = 0 \quad (17)$$

where $e_v$ is an aquifer void ratio, $h_0$ is a thickness constant relating the state parameter to a physical sheet thickness via $h = h_0 \theta$, $\Psi_{\text{source}}$ is a combined water source term from rainfall and surface melt reaching the bed, as well as melting due to geothermal heat. $\Psi_{\text{melt}}$ is a frictional melt rate $\Psi_{\text{melt}} = \frac{v \tau_b}{\rho_w L}$.

## 2.4. Model Implementation

We employ finite differences to represent spatial derivatives along a constant grid of 20 km length using a spatial increment $dx = 20$ m. We integrate the combined set of equations for a prescribed synthetic glacier geometry using an implicit time-stepping algorithm with adaptive time-stepping (maximum time-step of $10^3$ s). We first solve the system with an SIA solver assuming no sliding for 500 years, and then run it for a further 1,000 years with the coupled SIA/SSA solver to reach steady conditions. The glacier is defined where $H > 0$ m. Parameter sets used in the different simulations are given in Tables 1 and 2.

Most simulations use the same bed topography given by

$$z_b(x) = 2000 \text{ m } e^{-x/15000\text{m}}, \quad (18)$$
<table><caption>**Table 2**
List of Parameters That Are Kept Constant for All Simulations Presented in the Manuscript</caption><thead><tr><th>Parameter</th><th>Description</th><th>Units</th><th>Value</th></tr></thead><tbody><tr><td>ρ<sub>i</sub></td><td>Ice density</td><td>[kg m<sup>-3</sup>]</td><td>900</td></tr><tr><td>ρ<sub>w</sub></td><td>Water density</td><td>[kg m<sup>-3</sup>]</td><td>1,000</td></tr><tr><td>g</td><td>Gravitational acceleration</td><td>[m s<sup>-2</sup>]</td><td>9.8</td></tr><tr><td>L</td><td>Latent heat of fusion</td><td>[J kg<sup>-1</sup>]</td><td>3.0 × 10<sup>5</sup></td></tr><tr><td>m</td><td>Friction law rheology exponent</td><td>[]</td><td>3</td></tr><tr><td>n</td><td>Glen's exponent</td><td>[]</td><td>3</td></tr><tr><td>β<sub>c</sub></td><td>Conduit drainage exponent</td><td>[]</td><td>3/2</td></tr><tr><td>α<sub>c</sub></td><td>Conduit drainage exponent</td><td>[]</td><td>5/4</td></tr><tr><td>A</td><td>Ice rheology constant</td><td>[s<sup>-1</sup> Pa<sup>-3</sup>]</td><td>2.4 × 10<sup>-24</sup></td></tr><tr><td>l<sub>s</sub></td><td>Lateral conduit spacing</td><td>[m]</td><td>1,000</td></tr><tr><td>d<sub>c</sub></td><td>Characteristic cavity opening length scale</td><td>[m]</td><td>1.0</td></tr><tr><td>t<sub>c</sub></td><td>Characteristic cavity closure time scale</td><td>[s]</td><td>1.0 × 10<sup>6</sup></td></tr><tr><td>A<sub>s</sub></td><td>Basal shear stress prefactor</td><td>[m s<sup>-1</sup> Pa<sup>-3</sup>]</td><td>1.0 × 10<sup>-23</sup></td></tr><tr><td>f<sub>perc</sub></td><td>Percolation function</td><td>[]</td><td><sup>1</sup>&frasl;<sub>2</sub>(tanh(50(θ - <sup>1</sup>&frasl;<sub>2</sub>)) + 1)</td></tr><tr><td>e<sub>v</sub></td><td>Aquifer void ratio</td><td>[]</td><td>1.0 × 10<sup>-2</sup></td></tr><tr><td>q</td><td>Friction decay exponent</td><td>[]</td><td>2.5</td></tr><tr><td>K<sub>0</sub></td><td>Sheet conductivity background</td><td>[m<sup>4</sup> s kg<sup>-1</sup>]</td><td>1.0 × 10<sup>-8</sup></td></tr></tbody></table>

but the frictional triggered surge in Figure 5A has a slight modification in topography compared to the rest of the simulations

$$z_b(x) = 2000 \text{ m } e^{-x/15000\text{m}} - 10 \text{ atan}\left(\frac{x - 7000\text{m}}{500\text{m}}\right). \qquad (19)$$

The reason for this is to ensure that the largest slope of the glacier bed approximately coincides with the maximum surface slope, and that this occurs in the central part of the glacier. This design allows us to set the nucleation region of the surge to the central part of the glacier.

In our model, the climatic forcing on the glacier happens both through surface mass balance as well as through melt water input to the bed. The surface mass balance acts directly on the glacier geometry through adding or withdrawing mass at the surface (see Equation 6), and is computed as a linear function of surface elevation using prescribed values for the equilibrium line altitude (1,600 m a.s.l.) and mass balance gradient (4 × 10⁻³ yr⁻¹, see Table 2). In the present study, we consider these parameters as constant in time, but we systematically investigate the sensitivity of our results for varying mass balance gradient later, the results of which are presented in Figure 9. We note that in nature, most of the surface meltwater would run off in discrete channels either supra-, en- or subglacially, and only a fraction of this would directly affect a larger portion of the glacier base. Instead of resolving water exchange between these drainage structures and the entire ice-bed interface, we describe a distributed water supply that has a seasonal variability but a lower magnitude than surface meltwaters, similar to the approach by Benn et al. (2019). We take Ψ<sub>source</sub> as uniform in space but varying periodically over 1 year, following a sin² function with a peak value of 10⁻⁸ m s⁻¹ (corresponding to ≈0.9 mm day⁻¹). Furthermore, frictional melt water, generated as a consequence of sliding, is added to the melt water input.

## 3. Results

Forced by simple and regular mass balance and water-supply forcing, our model spontaneously generates a wide range of different and complex surges and other velocity variations. As an example, Figure 4 shows snapshots of glacier geometry, glacier velocity as well as effective normal stress, basal shear stress and stress from the glacier margins before, during and after a surge. The surge is triggered when τ<sub>b</sub> > Cσ<sub>N</sub> in a sufficiently large part of the
<figure><img src="image_6.png" /><figcaption>**Figure 4.** Glacier geometry (a, d, g), velocity *U* (b, e, h), shear stress bound *C*σ<sub>N</sub>, basal shear stress τ<sub>b</sub> and margin stress τ<sub>m</sub> (c, f, i) at nucleation (a, b, c), during (d, e, f) and after (g, h, i) a surge in the model (corresponding to the simulation presented in Figure 5B).</figcaption></figure>

base (a–c). In the unstable region, glacier velocity increases accompanied by a sharp reduction in τ<sub>b</sub> that is compensated for by an increase in τ<sub>m</sub> (d–f) but also by longitudinal stress transfer. The surge ends when τ<sub>b</sub> < Cσ<sub>N</sub> is reestablished (g–i). Animated versions illustrating the simulations presented in Figure 5 are available as Supporting Information S1.

Both sub-units of the drainage system, that is, the distributed sheet and the conduit system, have efficient and non-efficient modes. However, each is controlled by different mechanisms governing the transition between effective and ineffective drainage. (a) The low-pressure conduits open by melting and close by viscous creep, hence they adjust to variations in conduit discharge. They can efficiently evacuate water from the glacier bed as long as the effective normal stress is low (Equation 15). This is typically the case for thin ice, and conduit drainage efficiency decreases with increasing glacier thickness (Meierbachtol et al., 2013). (b) For the distributed sheet, we adopt the approach by Gilbert et al. (2022) and let the hydraulic conductivity vary with sliding speed, thereby accounting for the transient opening of interconnected subglacial cavities (Equation 13). This results in a strong feedback between friction and drainage which can reduce water pressure at higher sliding speeds and thereby stabilize basal motion. In the absence of an efficient conduit system, increases in water supply will result in increased water pressure coinciding with increased sliding. The concomitant increase in cavitation will also increase the conductivity of the sheet, which has a reducing effect on the water pressure. Since the change in state of the interface affects both drainage and friction simultaneously, the way the glacier modifies its sliding speed can differ drastically depending on the precise interplay between water transport at the base, the frictional response, the glacier geometry, the distribution of water pressure, topography, and the interplay with subglacial conduits.

The range of dynamical behaviors that we have observed in our simulations includes cyclic surge instabilities with close to constant recurrence intervals (Figures 5A and 6), irregular surges where both the magnitude and the recurrence interval of the surges vary (Figures 5B and 6), short duration speed-up events that do not turn into glacier-wide surges (Figure 7), stable slow velocities (Figure 8B), as well as pulse-like speed-up events (Figure 8A) reminiscent of self-healing slip pulses in frictional rupture (Shlomai & Fineberg, 2016).

This large variety of glacier velocity variations challenges the definition of glacier surges. In the following, we define glacier surges as velocity increases that show a well-defined nucleation region as well as an increasing extent of the surging region accompanied by a surface bulge. In our model, these surges are triggered by a
<figure><img src="image_7.png" /><figcaption>**Figure 5.** A: Surges cycle governed by increasing the basal shear stress $\tau_b$ through mass balance-induced thickening and steepening of the glacier surface (parameter set 2 in Table 1). B: Surge cycles governed by increasing the basal water pressure due to gradual shut-down of the conduit system (parameter set 1 in Table 1). In each part (A and B) the panels show the following: (a) Average effective normal stress ($\sigma_N$). (b) Average velocity $U$. (c) Spatio-temporal evolution of glacier thickness $H$ against the glacier geometry (from the top of the glacier to the front). (d) Spatio-temporal evolution of the shear stress bound $C\sigma_N$. (e) Spatio-temporal evolution of the basal shear stress $\tau_b$, using the same color coding as in (d). In A, (f) illustrates water fluxes from the sheet to the conduit system $\Psi_{\text{source}}$ (blue), and water supply from the glacier surface $\Psi_{\text{con}}$ (red). The values represent averages over the center region $x \in [2.5, 7.5]\,\text{km}$. The faded lines show the raw data, while the deep lines show a moving Gaussian average with standard deviation 2.6 years. In B, (g) illustrates the spatio-temporal evolution of conduit cross-sectional area $S$. The purple shading in (d)–(h) highlights the surge phase, and the hatched area just prior to model year 285, indicates $\tau_b > C\sigma_N$. (h) Spatio-temporal evolution of velocity $U$ during a surge. (i) Spatio-temporal evolution of the state parameter $\theta$ during a surge.</figcaption></figure>

frictional instability at the glacier bed, once the ratio of basal shear stress to effective normal stress $\tau_b/\sigma_N$ exceeds a critical value $C$ (Equation 10) in a sufficiently large region, equivalent to a nucleation length in rate-and-state friction (Handwerger et al., 2016; Rice & Ruina, 1983; Scholz, 1998; Thøgersen et al., 2019). The instability criterion can be reached either through a combination of increased $\tau_b$ through an increase of the glacier surface slope and thickness with time, through an increase in water pressure at the glacier bed, or a combination thereof. In Figure 5, we present examples of the two end-member scenarios; (a) surges caused by increased basal shear stress $\tau_b$ due to changes in glacier geometry (Figure 5A), and (b) surges driven by an increase in water pressure at the bed due to a gradual lowering of conduit drainage efficiency with increasing glacier thickness (Figures 5B and 4). We refer to these categories as shear-stress driven and hydraulically driven surges, respectively.

### 3.1. Shear-Stress Driven Surges

Shear-stress driven surges are caused by an increase in basal shear stress over time until a sufficiently high value is reached in a sufficiently large area of the glacier. This can in principle occur even at low basal water pressures, given that the critical shear stress $C$ is sufficiently low and that the friction law contains a transition from velocity strengthening to velocity weakening friction. This end-member is shown in Figure 5A, where we have selected a large value for the conduit conductivity to ensure efficient conduit drainage and assume no water pressure in the
<figure><img src="image_8.png" /><figcaption>**Figure 6.** Surge intervals corresponding to the surges in Figure 5. The surge intervals are determined using a surge criterion of $\langle U \rangle > 400$ m/yr. The data are derived from 9,000 years of simulation time and correspond to the cases illustrated in Figure 5A (red, parameter set 2) and B (blue, parameter set 1).</figcaption></figure>

conduit (Table 1). For this type of surge, we can distinguish two phases. In the quiescent phase, the glacier retreats and surface slope increases gradually. This corresponds to increasing basal shear stress $\tau_b$, while the effective stress $\sigma_N$ also increases due to thickness change (Figures 5Ad and 5Ae). The transport of water from the distributed system to the conduit system $\Psi^{sc}$ closely matches the water source term at the base $\Psi_{source}$ (Figure 5Af). A surge is triggered once the frictional stability threshold is surpassed. During the surge phase (highlighted in Figure 5A), the average velocity $\langle U \rangle$ increases by around two orders of magnitude (Figure 5Ab), and the surge propagates both up- and down-glacier, starting about 2 km upstream from the terminus at $x = 5$ km (Figures 5Ag and 5Ah). In our model, this type of surge is relatively insensitive to the hydraulic conductivity at the base, and entirely relies on whether or not the surface mass balance can build a sufficiently high basal shear stress to trigger a frictional instability. For this surge-type, the onset occurs in a large area close to the glacier front and could also be related to tidewater glacier surges with frontal destabilization, triggered by an increase in driving stress upstream of a frozen margin (Nuth et al., 2019; Sevestre et al., 2018).

## 3.2. Hydraulically Driven Surges

Hydraulically driven surges require high basal water pressures, for instance caused by low subglacial drainage efficiency. The model produces such behavior for low values of the channel conductivity $k_c$, relative to the scenario described above (Section 3.1, Table 1). As mentioned, conduit drainage becomes progressively less efficient with increasing glacier thickness requiring a larger hydraulic gradient. Steep glacier beds can provide a hydraulic gradient sufficient to drive the water through the drainage system, but at shallower slopes, water pressure has to build up to provide a sufficient hydraulic gradient. Hydraulically driven surges are therefore more likely at lower bed slopes.

A gradual decrease of the drainage efficiency during the quiescence will lead to a gradual increase in water pressure, corresponding to decreasing effective stress $\sigma_N$ (Figure 5Bd); the concomitant increase in sliding will also increase the degree of subglacial cavitation until eventually cavities suddenly become interconnected and the enhanced hydraulic conductivity in the distributed system allows the water to drain. In such a case, a frictional

instability cannot develop as the water pressure would drop as soon as velocity starts to increase, this is probably the case of most hard-bedded glaciers (Gilbert et al., 2022). This case is illustrated in Figure 8. In the absence of strong feedback between sliding speed and hydraulic conductivity, a frictional instability can be reached even for fairly high values of $C$ (hatched area in Figures 5Bd–5Bh), as long as water pressure is able to build up high enough over several years. In the quiescent phase, subglacial conduits dominate the drainage and $\sigma_N$ is high (Figure 5Bd) but gradually decreases, due to a gradual decrease of the drainage efficiency of the subglacial conduit system caused by ice thickening. Once effective pressure is sufficiently low, a stability criterion is eventually surpassed, triggering a surge (shading in Figures 5Bd–5Bh). The velocity increases by orders of magnitude (Figure 5Bb), and the instability quickly propagates up-glacier from the initiation region close to the front (Figures 5Bg and 5Bh). The rapid redistribution of ice mass during the surge causes reductions of glacier slope and ice thickness, paralleled by increasing drainage efficiency. The latter is accomplished by rapid growth of conduits (Figure 5Bf) and subsequent arrest of the surge due to restoration of low water pressure (high values of $\sigma_N$ in Figure 5Bd). This mechanism is similar to the hydraulic switch mechanism, although we stress that we do not find a hard switch but rather a gradual conduit shutdown with increasing glacier thickness. Figures 5Bd–5Be also display a similar evolution of the drainage system ~5 years after the surge. However, this time, the reorganization of subglacial

<figure><img src="image_9.png" /><figcaption>**Figure 7.** The model predicts short-duration speed up events that do not develop into full surges. The hydraulic conductivity of the distributed drainage system increases when the sliding velocity increases through the concomitant increase in $\theta$. The temporal evolution of the drainage system then alters the distribution of effective normal stress, terminating the glacier acceleration. The figure shows spatio-temporal plots of the glacier velocity $U$ (a), the effective normal stress $\sigma_N$ (b), the glacier thickness $H$ (c), and the state parameter $\theta$ (d). The figure illustrates a subset of the simulation using parameter set 2, entirely shown in Figure 5B.</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>**Figure 8.** Two classes of non-surging glaciers that are non-surge type because of an efficient distributed drainage system (A, parameter set 4 Table 1), or efficient conduit system (B, parameter set 3 Table 1). (a) Average effective normal stress $\langle\sigma_N\rangle$ as a function of time. (b) Average velocity $\langle U \rangle$ as a function of time. (c) Spatio-temporal evolution of glacier thickness $H$. (d) Spatio-temporal evolution of velocity $U$. (e) Spatio-temporal evolution of the state parameter $\theta$.</figcaption></figure>

hydraulic conditions is not accompanied by a surge because low slope and ice thickness of the post-surge geometry do not exert high enough basal shear stress to exceed the friction threshold. These two events illustrate how different behavior occurs depending on the competing evolution of drainage efficiency and glacier geometry.

### 3.3. Inter-Surge Intervals

Our model simulations exhibit different inter-surge intervals, and both regular as well as irregular length is observed (Figures 5 and 6). Due to the interaction between friction and the drainage system, irregular intervals can occur even under constant mass balance conditions where water input is only subject to seasonal variations. We have not carried out a detailed analysis of the origin of irregular surge cycles. We also stress that inter-surge intervals are not a characteristic of the different surge types but regular as well as irregular intervals can occur for both surge types.

### 3.4. Non-Surging Velocity Variations

Traversing the parameter space, our model simulations also revealed different classes of non-surging glaciers. Glaciers do not surge

1. If the background hydraulic conductivity is high or the water supply term $\Psi_{source}$ is small, such that water pressure does not build up.

2. If the drainage efficiency of subglacial conduits is high enough to drain all incoming water. This criterion is met if the glacier stays sufficiently thin, and links to surface mass balance conditions.

3. If the hydraulic conductivity of the distributed drainage system allows for drainage of all incoming water through increased sliding velocity, but without surpassing a frictional stability threshold.

Even though these non-surging cases fall beyond the surge definition we have adopted here, they can show strong velocity fluctuations, seasonal variations, and pulse-like events (Figure 7). Similar speed-up events also precede surges, and are noticeable in Figure 5Bb as minor spikes. These speed-ups occur before most but not all of the surges, and they are not systematic either.

Sometimes the frequency of occurrence of such events appears to increase toward the critical point (for instance just prior to model year 1,000) but most of the time they remain regular until the surge happens. These speed-ups are purely frictional as it is a low basal shear stress zone that propagates. Presumably, such events cannot be regarded as reliable precursors to predict a surge, due to their irregular occurrence and different appearance. Figure 8 illustrates details of two examples of simulations that did not develop surges. Figure 8A shows strong velocity variations (Figures 8Ab and 8Ad) that do not affect the entire glacier and the glacier geometry remains constant in terms of thickness and length (Figure 8Ac). This behavior is caused by the feedback between sliding velocity and distributed drainage that stabilizes variations in $\theta$ (Figure 8Ae). Figure 8B exhibits constant sliding velocity $U$ (Figures 8Bb and 8Bd) because the water pressure is low at all times, corresponding to almost constant, high values of $\langle\sigma_N\rangle$ (Figure 8Ba). This behavior does not trigger variations in the state parameter $\theta$ (Figure 8Be) and accordingly, the glacier geometry remains constant in terms of thickness and length (Figure 8Bc).

These criteria highlight the role of basal conditions in glacier surges. In particular, in our model, the hydraulic conductivity at the bed and its link to friction directly controls whether a glacier can be surge-type or not. Glaciers where the hydraulic conductivity is either sufficiently high at low sliding speeds, or where the hydraulic conductivity can increase drastically with relatively small increases in sliding speed without triggering frictional instability will be non-surge type (given that $C$ is sufficiently large compared to the surface slope). Surges can occur for glaciers where a frictional stability threshold is surpassed in the process of increasing the hydraulic conductivity of the distributed drainage system. This means that understanding the precise interplay between friction and drainage, and how it is affected by glacier geometry and climatic variables, is key to determine if a glacier can be surge type or not. When it comes to basal conditions, we identify low hydraulic conductivity,
<figure><img src="image_11.png" /><figcaption>**Figure 9.** Sensitivity of simulated surge behavior to variations in climate conditions, both in terms of the water supply $\Psi_{source}$ and mass balance gradient $\partial b/\partial z$. All other parameters are kept constant (Parameter set 5 in Table 1). Values of $\max((U))/\text{median}((U)) > 10$ are indicative of surge behavior. Each point is found from a 3,000 years simulation where the violins show the probability distribution of the maxima in glacier-wide velocities $\langle U \rangle_m$ (prominence 1% of maximum velocity and width 3.65 days) where we have excluded all minor variations defined by $\langle U \rangle_m < (\max((U)) - \min((U)))/2$. The solid lines connect the mean values of the selected velocity maxima.</figcaption></figure>

combined with a limited ability to increase the hydraulic conductivity through increasing sliding speed as key factors. If combined with a transition from velocity-strengthening friction at low sliding velocities to velocity-weakening friction at higher velocities this can produce surges that propagate both up- and down-glacier with well-defined surge fronts.

Climate conditions are also key factors in determining whether a glacier is of surge-type or not. Figure 9 shows how glacier velocity variations are affected by the mass balance gradient (which controls the glacier thickness) and the melt water supply to the bed, keeping all other parameters constant (parameter set 5 in Table 1). We find that:

1. A low water input to the bed inhibits large velocity fluctuations.

2. The limiting water input for the onset of large velocity fluctuations decreases with increasing surface mass balance gradient.

3. The limiting surface mass balance gradient for the onset of large velocity fluctuations increases with decreasing water input.

These findings imply that non-surge type glaciers, given favorable basal conditions, can switch to surge-type through an increasing surface mass balance gradient, through an increased water input to the bed, or a combination thereof. Changes in surging behavior of temperate glaciers are thus expected as a result of changes in the duration and intensity of the melt-season and precipitation patterns.

**4. Discussion**

The mathematical framework presented here is based on experimental and numerical evidence on temperate ice sliding over rigid bedrock. Despite the ability of the model to capture a large variety of glacier velocity variations, several parameters of the friction-drainage coupling are poorly constrained. First, the values of hydraulic conductivity needed to produce surges are on the low side of the expected range and therefore may be representative for soft-bedded glaciers. This raises a number of questions about the precise nature of the basal conditions of surge type glaciers, in particular how till influences the parameterization presented here. One way forward is through studies following up on recent work by Minchew and Meyer (2020), where porosity and the state of the interface for a soft-bedded glacier were described through two separate state evolution laws. Another possibility is the existence of mixed beds, where cavities partially filled with till could cause lower hydraulic conductivity than
what is expected for a hard-bedded glacier. Further, the direct coupling between friction and drainage has so far received little experimental attention. Experiments simultaneously measuring hydraulic conductivity and sliding speed, including transient changes, can provide very useful information on constraining the coupling between subglacial friction and drainage, for both soft-bedded and hard-bedded glaciers. If combined with experimental and observational constraints on empirical parameters and functional dependencies, we believe the physical framework adopted here has the potential to significantly improve our ability to predict glacier velocity variations. Unlocking this potential opens new possibilities to improve assessments of hazard potential of glaciers, but it also offers a largely unexplored path to reduce uncertainties related to basal friction in ice sheet models that form the basis for predictions of sea-level rise.

In our model, oscillatory solutions of glacier velocity can occur for purely velocity strengthening friction, but the propagation of surges with a sharp rupture front requires a transition from velocity strengthening to velocity-weakening friction. If, and under which conditions, velocity weakening friction can occur is still debated. While simulations and experiments on hard-bedded periodic topographies have reported on velocity weakening friction (Gagliardini et al., 2007; Zoet & Iverson, 2016), recent 3D simulations using real bedrock topographies in front of retreating glaciers report weak or no velocity weakening at high sliding speeds (Helanow et al., 2021; Woodard et al., 2023). However, recent studies by Gimbert et al. (2021) and Gilbert et al. (2022) fit Gagliardini et al. (2007)'s friction law to observations from Argentière glacier. Their findings suggest velocity weakening at high sliding speeds can exist for hard-bedded glaciers (Gimbert et al., 2021), but they suggest that a large increase in hydraulic conductivity due to opening of cavities can prevent the glacier from entering the unstable regime (Gilbert et al., 2022).

Soft-bedded glaciers also have a maximum friction threshold often called the Coulomb criterion. The existence of velocity-weakening friction for soft bedded glaciers is also debatable. So far, experiments of clean ice sliding over subglacial till at the pressure melting point have not found substantial velocity weakening friction (e.g., Zoet & Iverson, 2020). However, plowing clasts (Iverson, 2010; Thomason & Iverson, 2008), ice sliding over frozen till (Saltiel et al., 2021), and entrained debris (Lipovsky et al., 2019) have all been suggested as possible mechanisms causing velocity weakening friction, similar to the unstable behavior proposed by Tulaczyk et al. (2000a, 2000b). Furthermore, velocity weakening may arise from deformation in the till itself, though available data is limited. Iverson et al. (1998) showed slight velocity weakening with a till in sediment ring shear experiments, Tika et al. (1996) showed velocity weakening in fine-grained sediments analogous to till. Rathbun et al. (2008), however, showed mostly velocity strengthening for rate and state tests on two tills.

In this study, we introduced a number of simplifications in order to make long time-series simulations tractable. We used a coupled shallow ice/shallow shelf approximation combined with a parameterization of glacier width which was assumed constant. The shallow shelf approximation is activated at high sliding speeds so that longitudinal stress gradients are accounted for during high velocity events in the simulations. These assumptions have an impact on model accuracy which is likely not crucial for the purpose of studying synthetic glaciers as synthetic geometries and dynamical behavior are largely dominated by evolution of friction, rendering accurate descriptions of ice deformation subordinate. Nevertheless, we stress that for real-world applications with more complex geometries, the flow-band formulation of our ice-flow model should be extended to ensure its suitability, especially with respect to lateral shear (Sergienko, 2012).

We also made the assumption of laminar water flow in the distributed drainage system. This is most accurate in the limit of low hydraulic conductivity associated with low degree of cavitation, but we do expect a transition to turbulent flow at higher degrees of subglacial cavitation, and if comparing directly with observations this refinement will be necessary. In addition, the conduit system water pressure is not solved for. Instead, we solve for the water exchange term between the two sub-systems required to keep the conduit water pressure at the prescribed value (*sheet pressure* or *zero pressure*, Table 1).

For the *zero pressure* assumption, the conduit geometry may become unstable (Equation 15); the opening rate (first term on the right-hand side of Equation 15) varies with water flow whereas the closure rate (second term) slowly varies with ice thickness. This gives rise to unstable behavior with either collapse or unbounded growth of the conduit, once a delicate balance between conduit discharge and ice thickness is perturbed. In consequence, the drainage system displays binary behavior with complete dominance of either sheet or conduit drainage. The latter results in negative water pressure in the sheet, when all available water is attracted by the conduit system, and provokes a violation of mass conservation. To set a limit to such unwanted behavior, we artificially increase the
closure rate by a factor 10³, once water pressure in the sheet turns negative; this ensures fast termination of the runaway. We admit that this ad hoc formulation may be unsatisfactory and should be addressed in future work. As such, it removes the possibility of efficient conduit drainage underneath thick glaciers; in nature this situation may exist for conduits with high water pressure. This is clearly a weakness of our study as the effect of glacier thickening on reducing drainage efficiency may be overestimated. However, the existence of conduits at high water pressure would not necessarily prevent a frictional instability. For a given discharge, in steady state, the effective pressure of the conduits tends to remain the same through time with water pressure increase compensating the ice thickness increase (Gimbert et al., 2021). Meanwhile, the basal shear stress generally increases with glacier thickening, so $\tau_b$ may eventually exceed $C\sigma_N$, triggering instability. Thus, the mechanism of hydraulically driven surges presented in this study would remain similar, although the water pressure regulation from conduits would act differently.

## 5. Conclusions

Despite its numerous simplifications, our model reproduces a wide range of complex glacier behavior, including steadily sliding glaciers, surges propagating both down- and up-glacier, pulse-like velocity increases, speed-up events that do not turn into full-scale surges and surges that terminate before affecting the entire glacier. We show that our model spontaneously produces unstable glacier behavior such as surges. Such complex behavior arises from the consistent, two-way coupling between rate-and-state friction and subglacial hydrology, allowing for mutual interaction and feedbacks. Glacier surges are favored by low hydraulic conductivity that is largely insensitive to sliding speed variations and by a frictional response that contains a velocity-weakening transition. This study can be seen as a first step toward a physically based identification of the surge conditions. In contrast to the prevalent perception of glacier surging as an intrinsic phenomenon, we have also shown that changing climatic conditions have the potential to cause surges of previously stably flowing glaciers, as previously proposed by Sevestre and Benn (2015).

We hope that future studies will help refining some of the poorly constrained parameters and equations. In particular, open questions remain concerning hydraulic conductivity as a function of cavitation, as well as the time and length scales of the transient response of friction to normal stress. The role of velocity weakening friction in the surge behavior as well as the effects of sub-freezing ice temperatures should also be investigated within the proposed framework. Finally, implementing our approach in more complex models taking into account the three dimensional geometry and a more realistic subglacial hydrology would be technically possible although the computational cost remains challenging. In general, we believe that accurately accounting for transient feedback loops between friction and drainage has the potential not only to improve our understanding of glacier stability, but also to improve predictions of future glacier changes.

## Data Availability Statement

Software—The presented simulations and figures have been generated using version 7 of Thøgersen et al. (2024), publicly accessible at Zenodo https://zenodo.org/records/11109612.

## Acknowledgments

This study has been supported by Earth Flows, a strategic research initiative of the Faculty of Mathematics and Natural Sciences, University of Oslo, and the project MAMMAMIA, funded by the Research Council of Norway (project no. 301837). The authors appreciate the feedback from the editor Olga Sergienko, an associate editor, Martin Truffer, Luke Zoet and three further anonymous reviewers, that helped improving previous versions of the manuscript. TVS acknowledges support from the project FricFrac funded by the Center for Advanced Study (CAS) at the Norwegian Academy of Science and Letters during the academic year 2023–2024.

## References

Altena, B., Scambos, T., Fahnestock, M., & Kääb, A. (2019). Extracting recent short-term glacier velocity evolution over southern Alaska and the Yukon from a large collection of Landsat data. *The Cryosphere*, 13(3), 795–814. https://doi.org/10.5194/tc-13-795-2019

Benn, D., Fowler, A. C., Hewitt, I., & Sevestre, H. (2019). A general theory of glacier surges. *Journal of Glaciology*, 65(253), 701–716. https://doi.org/10.1017/jog.2019.62

Bhambri, R., Hewitt, K., Kawishwar, P., & Pratap, B. (2017). Surge-type and surge-modified glaciers in the Karakoram. *Scientific Reports*, 7(1), 1–14. https://doi.org/10.1038/s41598-017-15473-8

Björnsson, H., Pálsson, F., Sigurðsson, O., & Flowers, G. E. (2003). Surges of glaciers in Iceland. *Annals of Glaciology*, 36, 82–90. https://doi.org/10.3189/172756403781816365

Bueler, E., & Brown, J. (2009). Shallow shelf approximation as a “sliding law” in a thermomechanically coupled ice sheet model. *Journal of Geophysical Research*, 114(F3), F03008. https://doi.org/10.1029/2008jf001179

Clarke, G. K. (1987). Subglacial till: A physical framework for its properties and processes. *Journal of Geophysical Research*, 92(B9), 9023–9036. https://doi.org/10.1029/jb092ib09p09023

Cuffey, K. M., & Paterson, W. S. B. (2010). *The physics of glaciers*. Academic Press.

Dunse, T., Schellenberger, T., Hagen, J., Kääb, A., Schuler, T., & Reijmer, C. (2015). Glacier-surge mechanisms promoted by a hydro-thermodynamic feedback to summer melt. *The Cryosphere*, 9(1), 197–215. https://doi.org/10.5194/tc-9-197-2015

Eisen, O., Harrison, W., & Raymond, C. (2001). The surges of Variegated Glacier, Alaska, USA, and their connection to climate and mass balance. *Journal of Glaciology*, 47(158), 351–358. https://doi.org/10.3189/172756501781832179
Flowers, G. E. (2015). Modelling water flow under glaciers and ice sheets. *Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences*, **471**(2176), 20140907. https://doi.org/10.1098/rspa.2014.0907

Fowler, A. (1987). A theory of glacier surges. *Journal of Geophysical Research*, **92**(B9), 9111–9120. https://doi.org/10.1029/jb092ib09p09111

Fowler, A., Murray, T., & Ng, F. (2001). Thermally controlled glacier surging. *Journal of Glaciology*, **47**(159), 527–538. https://doi.org/10.3189/172756501781831792

Fowler, A., & Nye, J. F. (1986). A sliding law for glaciers of constant viscosity in the presence of subglacial cavitation. *Proceedings of the Royal Society of London. A. Mathematical and Physical Sciences*, **407**(1832), 147–170. https://doi.org/10.1098/rspa.1986.0090

Gagliardini, O., Cohen, D., Råback, P., & Zwinger, T. (2007). Finite-element modeling of subglacial cavities and related friction law. *Journal of Geophysical Research*, **112**(F2), F02027. https://doi.org/10.1029/2006jf000576

Gilbert, A., Gimbert, F., Gagliardini, O., & Vincent, C. (2023). Inferring the basal friction law from long term changes of glacier length, thickness and velocity on an alpine glacier. *Geophysical Research Letters*, **50**(16), e2023GL104503. https://doi.org/10.1029/2023GL104503

Gilbert, A., Gimbert, F., Thøgersen, K., Schuler, T. V., & Kääb, A. (2022). A consistent framework for coupling basal friction with subglacial hydrology on hard-bedded glaciers. *Geophysical Research Letters*, **49**(13), e2021GL097507. https://doi.org/10.1029/2021GL097507

Gimbert, F., Gilbert, A., Gagliardini, O., Vincent, C., & Moreau, L. (2021). Do existing theories explain seasonal to multi-decadal changes in glacier basal sliding speed? *Geophysical Research Letters*, **48**(15), e2021GL092858. https://doi.org/10.1029/2021GL092858

Hamilton, G. S., & Dowdeswell, J. A. (1996). Controls on glacier surging in Svalbard. *Journal of Glaciology*, **42**(140), 157–168. https://doi.org/10.1017/s0022143000030616

Handwerger, A. L., Rempel, A. W., Skarbek, R. M., Roering, J. J., & Hilley, G. E. (2016). Rate-weakening friction characterizes both slow sliding and catastrophic failure of landslides. *Proceedings of the National Academy of Sciences of the United States of America*, **113**(37), 10281–10286. https://doi.org/10.1073/pnas.1607009113

Harrison, W., Raymond, C., & MacKeith, P. (1986). Short period motion events on Variegated Glacier as observed by automatic photography and seismic methods. *Annals of Glaciology*, **8**, 82–89. https://doi.org/10.3189/S0260305500001191

Helanow, C., Iverson, N. R., Woodard, J. B., & Zoet, L. K. (2021). A slip law for hard-bedded glaciers derived from observed bed topography. *Science Advances*, **7**(20), eabe7798. https://doi.org/10.1126/sciadv.abe7798

Herreid, S., & Truffer, M. (2016). Automated detection of unstable glacier flow and a spectrum of speedup behavior in the Alaska Range. *Journal of Geophysical Research: Earth Surface*, **121**(1), 64–81. https://doi.org/10.1002/2015jf003502

Hewitt, I. J. (2011). Modelling distributed and channelized subglacial drainage: The spacing of channels. *Journal of Glaciology*, **57**(202), 302–314. https://doi.org/10.3189/002214311796405951

Hoffman, M., & Price, S. (2014). Feedbacks between coupled subglacial hydrology and glacier dynamics. *Journal of Geophysical Research: Earth Surface*, **119**(3), 414–436. https://doi.org/10.1002/2013jf002943

Iken, A. (1981). The effect of the subglacial water pressure on the sliding velocity of a glacier in an idealized numerical model. *Journal of Glaciology*, **27**(97), 407–421. https://doi.org/10.1017/s0022143000011448

Iken, A., Röthlisberger, H., Flotron, A., & Haeberli, W. (1983). The uplift of Unteraargletscher at the beginning of the melt season—A consequence of water storage at the bed? *Journal of Glaciology*, **29**(101), 28–47. https://doi.org/10.3189/S0022143000005128

Iverson, N. R. (2010). Shear resistance and continuity of subglacial till: Hydrology rules. *Journal of Glaciology*, **56**(200), 1104–1114. https://doi.org/10.3189/002214311796406220

Iverson, N. R., Hooyer, T. S., & Baker, R. W. (1998). Ring-shear studies of till deformation: Coulomb-plastic behavior and distributed strain in glacier beds. *Journal of Glaciology*, **44**(148), 634–642. https://doi.org/10.3189/s0022143000002136

Jiskoot, H., Murray, T., & Boyle, P. (2000). Controls on the distribution of surge-type glaciers in Svalbard. *Journal of Glaciology*, **46**(154), 412–422. https://doi.org/10.3189/172756500781833115

Kamb, B. (1987). Glacier surge mechanism based on linked cavity configuration of the basal water conduit system. *Journal of Geophysical Research*, **92**(B9), 9083–9100. https://doi.org/10.1029/jb092ib09p09083

Kamb, B., Raymond, C., Harrison, W., Engelhardt, H., Echelmeyer, K., Humphrey, N., et al. (1985). Glacier surge mechanism: 1982–1983 surge of Variegated Glacier, Alaska. *Science*, **227**(4686), 469–479. https://doi.org/10.1126/science.227.4686.469

Kjær, K. H., Larsen, E., van der Meer, J., Ingólfsson, Ó., Krüger, J., Benediktsson, Í., Ó., et al. (2006). Subglacial decoupling at the sediment/bedrock interface: A new mechanism for rapid flowing ice. *Quaternary Science Reviews*, **25**(21–22), 2704–2712. https://doi.org/10.1016/j.quascirev.2006.06.010

Lipovsky, B. P., Meyer, C. R., Zoet, L. K., McCarthy, C., Hansen, D. D., Rempel, A. W., & Gimbert, F. (2019). Glacier sliding, seismicity and sediment entrainment. *Annals of Glaciology*, **60**(79), 182–192. https://doi.org/10.1017/aog.2019.24

Meier, M. F., & Post, A. (1969). What are glacier surges? *Canadian Journal of Earth Sciences*, **6**(4), 807–817. https://doi.org/10.1139/e69-081

Meierbachtol, T., Harper, J., & Humphrey, N. (2013). Basal drainage system response to increasing surface melt on the Greenland ice sheet. *Science*, **341**(6147), 777–779. https://doi.org/10.1126/science.1235905

Minchew, B., & Meyer, C. R. (2020). Dilation of subglacial sediment governs incipient surge motion in glaciers with deformable beds. *Proceedings of the Royal Society A*, **476**(2238), 20200033. https://doi.org/10.1098/rspa.2020.0033

Nick, F. M., Van der Veen, C. J., Vieli, A., & Benn, D. I. (2010). A physically based calving model applied to marine outlet glaciers and implications for the glacier dynamics. *Journal of Glaciology*, **56**(199), 781–794. https://doi.org/10.3189/002214310794457344

Nuth, C., Gilbert, A., Köhler, A., McNabb, R., Schellenberger, T., Sevestre, H., et al. (2019). Dynamic vulnerability revealed in the collapse of an Arctic tidewater glacier. *Scientific Reports*, **9**(1), 5541. https://doi.org/10.1038/s41598-019-41117-0

Rathbun, A. P., Marone, C., Alley, R. B., & Anandakrishnan, S. (2008). Laboratory study of the frictional rheology of sheared till. *Journal of Geophysical Research*, **113**, F02020. https://doi.org/10.1029/2007jf000815

Rice, J., & Ruina, A. L. (1983). Stability of steady frictional slipping. *Journal of Applied Mechanics*, **50**(2), 343–349. https://doi.org/10.1115/1.3167042

Saltiel, S., McCarthy, C., Creyts, T. T., & Savage, H. M. (2021). Experimental evidence of velocity-weakening friction during ice slip over frozen till: Implications for basal seismicity in fast moving, soft-bed glaciers and ice streams. *Seismological Research Letters*, **92**(5), 2793–2810. https://doi.org/10.1785/0220200480

Scholz, C. H. (1998). Earthquakes and friction laws. *Nature*, **391**(6662), 37–42. https://doi.org/10.1038/34097

Schoof, C. (2005). The effect of cavitation on glacier sliding. *Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences*, **461**(2055), 609–627. https://doi.org/10.1098/rspa.2004.1350

Sergienko, O. V. (2012). The effects of transverse bed topography variations in ice-flow models. *Journal of Geophysical Research*, **117**(F3), F03011. https://doi.org/10.1029/2011JF002203

Sevestre, H., & Benn, D. I. (2015). Climatic and geometric controls on the global distribution of surge-type glaciers: Implications for a unifying model of surging. *Journal of Glaciology*, **61**(228), 646–662. https://doi.org/10.3189/2015jog14j136
Sevestre, H., Benn, D. I., Luckman, A., Nuth, C., Kohler, J., Lindbäck, K., & Pettersson, R. (2018). Tidewater glacier surges initiated at the terminus. *Journal of Geophysical Research: Earth Surface, 123*(5), 1035–1051. https://doi.org/10.1029/2017JF004358

Shlomai, H., & Fineberg, J. (2016). The structure of slip-pulses and supershear ruptures driving slip in bimaterial friction. *Nature Communications*, 7(1), 1–7. https://doi.org/10.1038/ncomms11787

Thøgersen, K., Gilbert, A., Bouchayer, C., & Schuler Vikhamar, T. (2024). Simulations and scripts for “Glacier surges controlled by the close interplay between subglacial friction and drainage” [Software]. Zenodo. https://doi.org/10.5281/zenodo.11109612

Thøgersen, K., Gilbert, A., Schuler, T. V., & Malthe-Sørenssen, A. (2019). Rate-and-state friction explains glacier surge propagation. *Nature Communications*, 10(1), 1–8. https://doi.org/10.1038/s41467-019-10506-4

Thomason, J. F., & Iverson, N. R. (2008). A laboratory study of particle ploughing and pore-pressure feedback: A velocity-weakening mechanism for soft glacier beds. *Journal of Glaciology*, 54(184), 169–181. https://doi.org/10.3189/002214308784409008

Tika, T. E., Vaughan, P. R., & Lemos, L. J. L. J. (1996). Fast shearing of pre-existing shear zones in soil. *Géotechnique*, 46(2), 197–233. https://doi.org/10.1680/geot.1996.46.2.197

Truffer, M., Kääb, A., Harrison, W. D., Osipova, G. B., Nosenko, G. A., Espizua, L., et al. (2021). Glacier surges. In *Snow and ice-related hazards, risks, and disasters* (pp. 417–466). Elsevier.

Tulaczyk, S., Kamb, W. B., & Engelhardt, H. F. (2000a). Basal mechanics of ice stream B, west Antarctica: 1. Till mechanics. *Journal of Geophysical Research*, 105(B1), 463–481. https://doi.org/10.1029/1999JB900329

Tulaczyk, S., Kamb, W. B., & Engelhardt, H. F. (2000b). Basal mechanics of ice stream B, west Antarctica: 2. Undrained plastic bed model. *Journal of Geophysical Research*, 105(B1), 483–494. https://doi.org/10.1029/1999JB900328

Turrin, J. B., Forster, R. R., Sauber, J. M., Hall, D. K., & Bruhn, R. L. (2014). Effects of bedrock lithology and subglacial till on the motion of Ruth Glacier, Alaska, deduced from five pulses from 1973 to 2012. *Journal of Glaciology*, 60(222), 771–781. https://doi.org/10.3189/2014jog13j182

Van der Veen, C. (1999). *Fundamentals of glacier dynamics*. A.A. Balkema.

Walder, J. S. (1986). Hydraulics of subglacial cavities. *Journal of Glaciology*, 32(112), 439–445. https://doi.org/10.1017/s0022143000012156

Werder, M. A., Hewitt, I. J., Schoof, C. G., & Flowers, G. E. (2013). Modeling channelized and distributed subglacial drainage in two dimensions. *Journal of Geophysical Research: Earth Surface*, 118(4), 2140–2158. https://doi.org/10.1002/jgrf.20146

Woodard, J. B., Zoet, L. K., Iverson, N. R., & Helanow, C. (2023). Inferring forms of glacier slip laws from estimates of ice-bed separation during glacier slip. *Journal of Glaciology*, 69(274), 324–332. https://doi.org/10.1017/jog.2022.63

Zoet, L. K., & Iverson, N. R. (2015). Experimental determination of a double-valued drag relationship for glacier sliding. *Journal of Glaciology*, 61(225), 1–7. https://doi.org/10.3189/2015jog14j174

Zoet, L. K., & Iverson, N. R. (2016). Rate-weakening drag during glacier sliding. *Journal of Geophysical Research: Earth Surface*, 121(7), 1206–1217. https://doi.org/10.1002/2016jf003909

Zoet, L. K., & Iverson, N. R. (2020). A slip law for glaciers on deformable beds. *Science*, 368(6486), 76–78. https://doi.org/10.1126/science.aaz1183

Zoet, L. K., Iverson, N. R., Andrews, L., & Helanow, C. (2022). Transient evolution of basal drag during glacier slip. *Journal of Glaciology*, 68(270), 741–750. https://doi.org/10.1017/jog.2021.131