Dilation of subglacial sediment governs incipient surge motion in glaciers with deformable beds

Glacier surges are quasi-periodic episodes of rapid ice flow that arise from increases in slip rate at the ice–bed interface. The mechanisms that trigger and sustain surges are not well understood. Here, we develop a new model of incipient surge motion for glaciers underlain by sediments to explore how surges may arise from slip instabilities within a thin layer of saturated, deforming subglacial till. Our model represents the evolution of internal friction, porosity and pore water pressure within the till as functions of the rate and history of shear deformation, and couples the till mechanics to a simple ice-flow model. Changes in pore water pressure govern incipient surge motion, with less permeable till facilitating surging because dilation-driven reductions in pore water pressure slow the rate at which till tends towards a new steady state, thereby allowing time for the glacier to thin dynamically. The reduction of overburden (and thus effective) pressure at the bed caused by dynamic thinning of the glacier sustains surge acceleration in our model. The need for changes in both the hydromechanical properties of the till and the thickness of the glacier creates restrictive conditions for surge motion that are consistent with the rarity of surge-type glaciers and their geographical clustering.

Glacier surges are quasi-periodic episodes of rapid ice flow that arise from increases in slip rate at the ice-bed interface. The mechanisms that trigger and sustain surges are not well understood. Here, we develop a new model of incipient surge motion for glaciers underlain by sediments to explore how surges may arise from slip instabilities within a thin layer of saturated, deforming subglacial till. Our model represents the evolution of internal friction, porosity and pore water pressure within the till as functions of the rate and history of shear deformation, and couples the till mechanics to a simple ice-flow model. Changes in pore water pressure govern incipient surge motion, with less permeable till facilitating surging because dilation-driven reductions in pore water pressure slow the rate at which till tends towards a new steady state, thereby allowing time for the glacier to thin dynamically. The reduction of overburden (and thus effective) pressure at the bed caused by dynamic thinning of the glacier sustains surge acceleration in our model. The need for changes in both the hydromechanical properties of the till and the thickness of the glacier creates restrictive conditions for surge motion that are consistent with the rarity of surge-type glaciers and their geographical clustering. are assumed to be the same in both layers. Our idealized glacier has a subglacial hydrological system that, like any glacier, evolves as a result of changes in meltwater flux and basal slip rate [68][69][70]. Here, we assume that both the state of the hydrological system and the basal water flux are accounted for in p w r , the water pressure within the hydrological system, depicted as a reservoir in the system diagram (figure 1). We assume that basal slip is due entirely to deformation of the upper till layer and there is no sliding between the ice and uppermost till layer. As a result, p w r only influences ice flow through its influence on p w . We make this simplifying assumption in spite of the fact that p w r may cause sliding of the ice relative to the bed [69,[71][72][73][74] because our focus is on how the mechanical properties of till might induce surging in the absence of meltwater flux to the bed. This assumption of nearly constant p w r is merely conceptual and is not a necessary condition in the subsequent derivation because time-varying p w r is accounted for in the model. Indeed, in future work, subglacial hydrological models could be readily bolted onto the model presented here. For simplicity, we ignore potential changes in pore water pressure caused by ploughing particles [41,45], and begin our study at the glacier bed with an exploration of till mechanics.

(a) Mechanical properties of till
We adopt a phenomenological model for the mechanical strength of till that depends on basal slip rate u b and the state of the subglacial till θ . This rate-and-state friction model accounts for instantaneous basal slip rate and, importantly, basal slip history, and was derived to explain numerous laboratory measurements of sliding on bare rock and granular interfaces. Rate-and-state friction is widely used in studies of earthquake nucleation and slow-slip events on tectonic faults, and gives the instantaneous shear strength of subglacial till as [55,56] where μ n is the coefficient of nominal internal friction (i.e. friction coefficient at steady state and u b = u b n ), d c is a characteristic slip displacement, u b n is a constant reference velocity, and the constants a and b are material parameters that define the magnitude of the direct (velocity, u b ) and evolution (state, θ ) effects, respectively. As we will discuss, b is important for this study because it encodes the effect of dilation on the bulk friction coefficient μ. In our idealized glacier geometry, the bed is horizontal and effective normal stress is equal to effective pressure N, defined as with p w the pore water pressure in the till and the ice overburden pressure p i defined as where ρ i is the mass density of ice and g is gravitational acceleration. Rate-and-state friction has received attention in studies of the ice-bed interface [38,41,44,75,76] and is widely studied for slip on tectonic faults containing gouge [52,[77][78][79], a material mechanistically similar to till [80]. Though distinct in many respects, earthquakes and glacier surges are analogous in the sense that both involve long quiescent periods and relatively short activation time scales. Slow slip on tectonic faults is particularly relevant to studying glacier surges because of their comparable slip durations and slow slip rates compared with major earthquakes [52,53]. Incipient motion in both earthquakes and glacier surges is brought on by excess applied stress relative to frictional resistance. While stresses and displacement rates are orders of magnitude higher in earthquakes than in glaciers, the experimentally verified rateand-state friction model is applicable to glacier surges as there is no known lower bound on velocity for the model to be valid [81]. Furthermore, because the rate-and-state friction model given in equation (2.1) represents friction along grain boundaries within the deforming till layer, this model is consistent with the skin-friction regime in regularized Coulomb friction rules [82,83].
When till is deformed, individual grains are mobilized by cataclastic flow (which includes grain rolling and grain boundary sliding), dilation and comminution. Under small displacements, the granular structure of the till is related to the pre-deformed structure, meaning that the till essentially remembers its prior state. Memory is represented by the state variable θ (units of time). State has been taken to represent the product of the contact area and intrinsic strength (quality) of the contact [84], but also has been interpreted as the average age of contacts between loadbearing asperities [85]. Under either interpretation, state is expected to evolve as a function of time, slip and effective normal stress [55,[85][86][87]. To represent the evolution of θ, we adopt the state evolution equation, sometimes referred to as the slip law [56]  Increasing u b beyond d c /θ -through enhanced surface meltwater flux, calving or other external forcing-will reduce θ over time. Similarly, when u b < d c /θ , θ will increase towards steady state.
In the next section, we show that changes in θ are brought about through till compaction and dilation. As such, θ accounts for the basal slip history and plays a key role in determining bed strength and the response of bed strength to shear and external forcing. Steady-state till shear strength occurs when state evolution ceases (θ = 0) and is defined aŝ The other important factor to consider in this study is d c , the slip distance over which state (and porosity) evolve. Computational and microphysical studies have concluded that d c is proportional to the thickness of the deforming layer [79,88,89], which can be expected to be of order 0.1-1 m in subglacial till and varies with permeability [59,90]. Other factors influencing d c include grain size and porosity [79].
(b) Pore water pressure then the rate of change of water mass per unit volume within the till is given aṡ where φ is the (dimensionless) effective till porosity, defined as the ratio of pore volume to total volume, and ρ w is the density of water. In this section, we seek to understand the rate of change in pore water pressure as a function of basal slip rate under the basic assumptions that water is incompressible over the range of reasonable subglacial pressures and that frictional heating at the ice-bed interface and plastic dissipation within the till are negligible.

(i) Evolution of porosity
Assuming that individual grains in the till are rigid, strain within the till will be accommodated by changes in porosity. Adopting an elastic-plastic model for the deformation of granular till, wherein the total strain is equal to the sum of the elastic and plastic strains, we separate porosity changes into an elastic componentṗ w β and a plastic componentφ p such that [78,91] φ =ṗ w β +φ p , is the till compressibility and e is the elastic compressibility coefficient, taken to be in the range e ∼ 10 −3 -10 −1 [92]. Following work by Segall & Rice [78] and Segall et al. [53] on slow-slip events on tectonic faults, we take the plastic component of porosity to have the same form as the evolution component of the rate-and-state model for till shear strength (equation (2.1)), namely where φ c is a (constant) characteristic porosity and p is a dilatancy coefficient, a dimensionless parameter hereafter assumed constant and in the range 10 −4 ≤ p ≤ 10 −2 [53]. We note that the only sensitivity in our model to the absolute value of p is to the evolution of porosity; surge behaviour, the main focus of this study, is influenced only by the ratio p /β, which represents the relative importance of each term in equation (2.7). By adopting equation (2.9), we are assuming that plastic deformation of the till is completely determined by changes in state, θ, the only variable in equation (2.9). This assumption is physically justifiable: irreversible changes in porosity necessitate a change in the average age of granular contacts and, equivalently, a change in the product of the contact area and quality, both of which are the physical interpretations of state (θ ) discussed above. Differentiating equation (2.9) in time yieldṡ an expression that indicates that shearing of the till layer causes it to compact (φ p < 0) when θ is below steady state (θ < d c /u b ) and to dilate when θ is above steady state. Such behaviour is consistent with observations of the response of over-and under-consolidated soils to shear [50]. This relationship between plastic till deformation and state gives rise to rich mechanical relationships between compaction, dilation and shearing, as is expected from sediments.
(ii) Evolution of pore water pressure Let us now consider water flux in the till in response to changes in porosity and sources outside the till shear layer.
Conservation of water mass gives ∂q w ∂z +ṁ w = 0, (2.12) where q w is the vertical water mass flux. Here, we have assumed that horizontal gradients in water pressure are negligible compared with vertical gradients and the bed slope is sufficiently shallow to allow us to consider only vertical water flux. Taking the basal ice to be impermeable then requires water flux to be entirely into and out of the deforming till layer. Under these conditions, Darcy's law is given as where γ h is the till permeability and η w is the dynamic viscosity of water. Combining equations (2.11)-(2.13) under the assumption that till permeability is spatially constant and independent of porosity givesṗ (2.14) where is the hydraulic diffusivity of the deforming till layer. Measurements of hydraulic diffusivity in till give a range for κ h of approximately 10 −9 -10 −4 m 2 s −1 , with a strong sensitivity to clay content [93,94]. We take constant effective permeability to be a reasonable first approximation given the small change in permeability under glaciologically relevant pressures and strains found in discrete-element modelling studies [90]. A more general treatment of pore water pressure evolution would include a porosity-dependent permeability in place of a constant effective permeability-for example, the Kozeny-Carman model used by Clarke [39]. We reserve this additional complexity for future work as our simple model retains the salient physical processes. Shearing in till concentrates in a thin, multi-layer zone that is typically several centimetres thick [60,[95][96][97]. We, therefore, approximate where h s is the thickness of the shear zone in the till, p w ∞ is the water pressure in the underlying permeable half-space and p w r is the water pressure in the basal hydrological system (figure 1). With this approximation, equation (2.14) becomeṡ where the first term represents Darcian flow into and out of the deforming till layer and the second term represents dynamical (dilation-driven) changes in pore water pressure. The Darcy flow component of pore water pressure evolution is inversely proportional to the characteristic diffusive time scale for pore water in the deforming till layer we can suppose, to a reasonable approximation, t h ∼ N, which should retain the same order of magnitude during incipient surge motion. Similarly for permeability, where compaction-driven reductions in permeability will induce relatively small (factor of 2) decreases in thickness h s [90]. Such small changes are unlikely to dramatically alter the dynamics of surge motion captured here, and we leave for future work a more detailed analysis involving variable t h . From the second term in equation (2.17), we can see that the sign of the dynamical (or dilationdriven) component ofṗ w is determined by the state of the till. When state θ is below (above) steady state and t h > 0, pore water pressure will increase (decrease) until steady state is achieved. These changes in pore water pressure are entirely due to changes in till porosity: compaction (φ p < 0) results in non-zero rates of change in the dynamical component of water pressure because the second term in equation (2.17) is pθ N/[ e θ (1 − φ) 2 ] = −φ p /β. Whether p w decreases or increases following step changes in basal slip rate depends on whether the ratio θu b /d c is greater than or less than unity. Equation (2.17) also shows that steady-state pore water pressure isp w = p w ∞ = p w r whenθ = 0.

(c) Basal slip acceleration
Glacier ice is an incompressible viscous fluid in laminar flow, and the momentum equation, incompressibility condition and continuity equation, respectively, take the forms where u i is the ice velocity vector,ū i is the depth-averaged ice velocity vector, τ ij is the deviatoric stress tensor, δ ij is the Kronecker delta,p is the mean isotropic ice stress (pressure),Ṁ is the total surface mass balance (which includes surface and basal mass balance and is positive for mass accumulation), and we employ the summation convention for repeated indices. To simplify our analysis, we neglect vertical shearing in the ice column and adopt a depth-integrated momentum equation (often referred to as the shallow-shelf approximation) [98] 2 where τ xx is the extensional deviatoric stress, τ xy is the lateral shear stress and we have neglected the transverse normal (deviatoric) stress τ yy . In some surge-type glaciers, vertical shearing may be the dominant flow regime during the quiescent phase, while basal slip is the dominant flow regime during the surge phase. Equation (2.22) is valid only when basal slip is dominant, and thus a model of basal slip acceleration derived from equation (2.22) may not fully detail glacier flow during incipient surge acceleration in some glaciers. Nevertheless, this simplification is reasonable because the focus of this work is on till mechanics and the flow model based on equation (2.22) will represent the salient processes of nascent surge acceleration. We reserve for future work a more detailed analysis that retains more components of the stress divergence and is able to capture the transition from vertical-shear-dominated flow to basal-slip-dominated flow. Force balance dictates that basal shear traction cannot exceed the lesser of applied stress and yield stress of the till, giving rise to the relation [19,46] where τ t = μN is the till shear strength (equation 2.1) and the gravitational driving stress is defined as where α is the ice surface slope, assumed small such that sin (α) ≈ α. Recall that we are focusing on the case in which rapid flow during the surge is accommodated primarily by deformation of the bed, giving rise to the relations τ b = τ t and u s ≈ u b , where u s is the ice-surface velocity. We note that equation (2.23) is consistent with the so-called regularized Coulomb sliding law, which has recently emerged as a candidate for a universal form of the sliding law, because, in this work, we are focused on the skin-friction regime defined by equation (2.1) [74,82,83,99].
Let us now focus only on the region where the surge is initiated and assume the areal extent of incipient surge motion is large enough to make the gradient of longitudinal stress (first term in equation (2.22)) negligible during the nascent surge phase. Taking ice to be a shear-thinning (i.e. pseudoplastic) viscous fluid, the constitutive relation, commonly known as Glen's law [100], iṡ whereε e = ε ijεij /2 is the effective strain rate, τ e = τ ij τ ij /2 is the effective deviatoric stress, the rate factor A is a scalar and the stress exponent is n = 3. Hereafter, A and n are assumed constant. Under our prior assumptions, 2ε e ≈ ∂u b /∂y and τ e ≈ τ xy . Integrating the reduced form of equation (2.22) twice along y subject to the symmetry condition τ xy = 0 at the centreline and no-slip condition at the margins gives the centreline basal slip rate [94] Taking w to be constant and differentiating equation (2.26) with respect to time yields an expression for acceleration of basal sliṗ where the rates of change in glacier geometry (ḣ andα), pore water pressure (ṗ w ) and state (θ) all contribute to the basal slip acceleration, along with instantaneous geometry (h and α), pore water pressure (p w ), state (θ ) and basal slip rate (u b ). Note that the conditions discussed and imposed in the model development-τ d > τ b (glacier is slipping at the bed), τ b = τ t (the till is deforming) and p w < p i (till has non-zero shear strength)-ensure that the denominator in equation (2.27) is always greater than zero. Equation (2.27) is the central result of this study. This formula describes the dependence of surge acceleration on glacier geometry, pore water pressure and the properties of the till. The terms in the numerator can be related to the processes of interest during the surge. Namely, the first term in the numerator (α) essentially represents the rate of change in the gravitational driving stress. The second term in the numerator captures the evolution of effective pressure (N), which governs the shear strength of the bed. The third and final term in the numerator accounts for the influence of dilation on the internal friction coefficient of the till. We spend the remainder of this study investigating the influence of the various physical processes represented in equation (2.27).

Results
Since shear strength of the till is the governing factor in surge motion and is defined by three variables (overburden pressure p i , pore water pressure p w and the internal friction coefficient μ), we present the results in three sections. In the first section, we discuss the evolution of pore water pressure following an increase in basal slip rate. Second, we consider the acceleration of basal slip for a glacier with a fixed geometry (i.e. fixed overburden pressure). Last, we explore the full model, which allows for variations in pore water pressure, glacier geometry and internal friction coefficient for till.  (a) Evolution of pore water pressure Pore water pressure in the deforming till layer evolves through dilation and compaction of the till as well as through the exchange of water between the deforming till layer, the subglacial hydrological system and the stagnant till layer that underlies the deforming layer (equation (2.17) and figure 2). In our model, the pressures in the stagnant till layer (p w ∞ ) and the subglacial hydrological system (p w r ) are assumed constant in time, and the flow of water into or out of the deforming till layer is described by Darcy's law (equation (2.13)). Using the parameter values given in the caption of figure 2, we integrate equations (2.4), (2.7) and (2.17) forwards in time from the (steady-state) initial conditions , an open-source Python toolkit [101]. The results shown in figure 2 illustrate how the evolution of pore water pressure p w following a step increase in basal slip rate (u b = 10u b 0 ) is influenced by the hydraulic diffusion time scale of the deforming till layer (t h ) and the relative values of the elastic ( e ) and plastic ( p ) compressibility coefficients. Note that, because we hold t h fixed in time, only the relative compressibility ratio e / p influences pore water pressure, not the absolute values of e and p . All cases shown in figure 2 start at steady state and indicate initial decreases in pore water pressure p w in response to till dilation followed by a return to steady state (p w =p w 0 = p w r = p w ∞ ) via Darcian flow over a time scale proportional to the diffusion time scale (cf. equation (2.17)). The minimum pore water pressure is determined by the diffusion time scale t h and the relative compressibility e / p . For a given relative compressibility, longer diffusion time scales, corresponding to lower till permeabilities, lead to a greater drop in pore water pressure (figure 2a). For a given diffusion

(b) Acceleration with fixed ice thickness
We now consider glacier acceleration. As a first step, we simplify our analysis by assuming that the time scale of interest is longer than the time scale for pore water diffusion (t > t h ) but short enough to allow us to reasonably neglect changes in glacier geometry. While it can be argued that this condition may be physically contrived in some cases, it is useful for exploring surge dynamics and the behaviour of the till in the absence of some complicating factors (in the next section, we will allow glacier geometry to evolve). After fixing glacier geometry by imposingḣ = 0 andα = 0 at all times, we solve the system of equations defined by equations (2.4), (2.7), (2.17) and (2.27). For all results discussed here, we prescribe as the initial velocity u b = 1.1u b 0 at t = 0, whereû b 0 = 10 m/yr, and set the initial values for all other variables to their respective steadystate values. The system of equations is stiff; therefore, we integrate forwards in time using an implicit Runge-Kutta method-specifically, the Radau IIA fifth-order method-implemented in In the cases shown in figure 3, we focus on the influences of a range of viable evolution effects (b values; indicated by line intensity and thickness) and different hydraulic diffusion time scales (t h ; indicated by colours). Aside from b and t h , all parameters are the same for all cases and are listed in the caption to figure 3. Note that a = 0.013, so in terms of the till friction coefficient μ, the cases shown in figure 3 are both rate weakening (a < b; solid lines) and rate strengthening (a > b; dashed lines).
The most notable feature in all cases shown in figure 3 is the lack of unstable acceleration. Steady-state speed is governed by the steady-state shear strength of till (equation (2.5)) and is therefore sensitive to the rate-and-state parameters (a − b) and μ n . Since the direct effect (a) is constant in all cases in figure 3, increasing the evolution effect (b) leads to a greater steady-state stress drop and faster steady-state basal slip rate because of the increasingly negative value (a − b). The steady-state values for all state variables are independent of the diffusion time scale t h and characteristic slip length d c . The primary influences of t h and d c are on the time the system takes to reach steady state and the peak change in pore water pressure. These results show that the system tends to steady state over a characteristic time scale that scales with the (dimensionless) hydraulic transmittance where λ is the (dimensionless) perturbation in u b (λ = 1.1 in figures 3-7). Equation (3.1) is defined as the ratio of the hydraulic diffusion time scale t h to the time scale for dilation-driven changes in pore water pressure d c e /( pûb 0 ), which follow from the coefficients in equation (2.17). The dependence on ψ 0 of the time to steady state is indicated in figure 3 by noting that the only term in ψ 0 that changes between the different cases is the t h . The time axes in figure 3 are normalized by d c e /( pûb 0 ), the time scale for dilation-driven changes in pore water pressure, to help show that model realizations in which the diffusion time scale t h is an order of magnitude longer take an order of magnitude longer time to evolve to steady state. As we show in the next section, the time required to reach steady state is a crucial factor governing whether or not a glacier surges.
The behaviour of the model in the absence of changes in glacier geometry (figure 3) provides further insight that helps to explain some of the results of the full model presented in the next section. For instance, the till dilates in all cases because of initial step and subsequent changes in basal slip rate (figure 3). The amplitude of the change in till porosity scales with the evolution parameter b, with larger values of b resulting in greater dilatancy. As seen in the previous section, higher dilatancy results in a larger drop in pore water pressure as the glacier accelerates. Dilatancy also drives a reduction in the internal friction coefficient of till, as a dilated till provides less resistance to shearing owing to reduced contact areas between grains. This drop in the internal friction coefficient commensurately reduces the shear strength of the till.  where ζ =ū/u s ,ū is the depth-averaged glacier speed and u s is the glacier surface speed. Since we have taken ice to be a non-Newtonian viscous fluid, we have (n + 1)/(n + 2) ≤ ζ ≤ 1, where n is the stress exponent in the constitutive relation for ice (equation (2.25)) [103]. In this study, we adopt the most common value for the stress exponent, n = 3, and we prescribe ζ = 1 for consistency with the reduced momentum equation in equation (2.22) (when ζ = 1, u s = u b ). We further simplify the expression for dynamical thinning by neglecting extensional strain rates (consistent with the assumptions in §2c), yieldingḣ where u * =Ṁ/(αζ ) is the balance velocity. Finally, the rate of change in surface slope becomeṡ where equation (3.4b) follows from the assumption of a parabolic surface profile for the glacier [92].  , long t h ), the latter of which we define as a period of rapid flow speeds (factor of 2 or more faster than quiescent speeds) that do not meet the definition of a surge, followed by a slowdown and evolution to steady state. To clarify the distinction: initial acceleration is unstable in surges and stable in abandoned surges.
To explore the processes that govern whether a surge develops, is abandoned or is essentially absent, let us focus on some illustrative cases shown in figure 4. We start with two prominent cases: those with the highest b values (and therefore the heaviest lines in figure 4) and different hydraulic diffusivities (i.e. t h values). The case with b = 0.05 and higher diffusivity (and, consequently, higher hydraulic permeability and shorter t h ), shown with the thick red lines in figure 4, undergoes an abandoned surge, defined by a brief acceleration phase, resulting in a maximum velocity of approximately twice the steady-state slip rate (u b /û b 0 ≈ 2), followed by deceleration and evolution to steady state. In this case, the glacier thins somewhat, but the till tends to steady state before there is any marked change in the effective pressure at the bed (N). The case with b = 0.05 and lower hydraulic diffusivity (thick blue line) surges, with muted acceleration (relative to the case with higher hydraulic diffusivity) preceding a continual reduction in state, pore water pressure, ice thickness and till internal friction coefficient. The rates of change in each of these values when the integration was terminated (at u b /û b 0 = 10) show that the glacier would continue to accelerate in the absence of contravening processes, such as increases in extensional stresses, that are not considered in our model but could manifest in a natural glacier. It is important to note that the effective pressure N continually decreases despite reductions in pore water pressure p w because of the dynamic thinning of the glacier. In other words, reductions in overburden pressure p i = ρ i gh outpace reductions in pore water pressure p w , leading to a net decrease in N = p i − p w that complements reductions in the friction coefficient μ, ensuring that basal drag (τ b = τ t = Nμ) diminishes in time. Sustained acceleration of the glacier unequivocally indicates that the decline in basal drag outpaces thinning-induced reductions in gravitational driving stress.
Other cases shown in figure 4 indicate the same basic behaviour: till with higher values of hydraulic permeability (shorter t h ) allows for faster acceleration, which causes the till to evolve to steady state before significant thinning of the glacier can occur. Rates of acceleration and evolution to steady state are slower in less permeable till, allowing rapid ice flow to persist for longer periods of time, facilitating dynamic thinning of the glacier. Longer time scales with  relatively muted acceleration allow for thinning because dynamic glacier thinning scales as the time integral of ice velocity (equation (3.3)), meaning that longer periods of moderately rapid flow can produce more thinning than much shorter periods of somewhat faster flow. These results suggest that it is the reduction in overburden pressure p i , and therefore effective pressure N, through dynamic thinning that is ultimately responsible for sustaining surge motion. The lack of unstable acceleration when glacier geometry is fixed in time (discussed in the previous section) and the manifestation of surging in cases of rate-strengthening friction coefficients (dashed lines in figure 4) both serve to highlight the importance of dynamic thinning for sustaining surge motion.
The evolution of till porosity when the glacier geometry is allowed to vary (figure 4) is markedly different from the case with fixed glacier geometry (previous section). With a fixed glacier geometry (constant overburden pressure), till consistently dilated because the effective pressure initially decreased and then returned to steady state along with water pressure. But when we allow the glacier to thin (thereby decreasing overburden), the dependence of the rate of change in porosity on the effective pressure (via β; equations (2.7) and (2.8)) results in compression of the till. As the effective pressure decreases as a result of thinning of the glacier, the rate of change in porosity becomes increasingly sensitive to changes in pore water pressure (cf. equation (2.8)). Since pore water pressure decreases in response to the evolution of till state (equation (2.17)), till compaction lags reductions in pore water pressure.
The results discussed in this section indicate that the principal factors governing the surge behaviour of a glacier are the hydraulic diffusion time scale of the deforming till layer t h , the relative compressibility e / p , and the evolution parameter b, the last of which dictates the response of the internal friction coefficient to till dilation. We explore this parameter space in figure 5; except where indicated, model parameters are the same as for figure 4, and we use the same numerical solver. The results in figure 5 show that, for any relative compressibility e / p , surge-type behaviour is favoured in glaciers with high b values and long diffusion time scales (i.e. relatively impermeable beds). Higher b values imply a greater reduction in the internal friction coefficient of till (μ) in response to changes in porosity (and, therefore, state), with rate-weakening values (b > a) resulting in a reduced steady-state friction coefficient. Positive glacier acceleration is generally expected as the friction coefficient decreases in response to state evolution, causing surges to be favoured at higher b values. As previously discussed, longer diffusion time scales (i.e. lower hydraulic permeability) diminish the rate of porosity (state) evolution, and, therefore, slow the increase in effective pressure caused by dilation-driven reductions in pore water pressure (i.e. dilatant hardening). Thus, slow diffusion of pore water enables a longer acceleration period that allows time for dynamic glacier thinning to drive a net reduction in the effective pressure. Surgetype glaciers are more likely to manifest in tills that have a high relative compressibility, e / p > 10, as these higher values imply less dilatant hardening (the reduction in pore water pressure due to shearing; cf. figure 2).
The rich dynamical behaviour illuminated in figure 5 is enhanced by the manifestation of regions (in the parameter space) of abandoned surges adjacent to the regions of surging behaviour. Abandoned surge regions are indicated in figure 5 by maximum basal slip rates greater than the initial value (u b max /û b 0 > 2, as shown in purple-to-red hues) and final basal slip rates less than the initial value (u b final /û b 0 < 1, as shown in grey tones). Abandoned surges manifest only where b values are relatively large but not large enough to produce a surge and diffusion time scales are slightly too short to allow for a full surge. According to our results, it is possible for a glacier to exhibit abandoned surges for any value of e / p , but the region in the parameter space that produces abandoned surges shrinks with increasing e / p (i.e. as dilatant hardening decreases).
Two other remarkable and persistent features of the parameter space are worth highlighting. First, abandoned surge regions are accompanied by a region in the parameter space that takes the shape of an aerofoil and contains points suitable for surge-type glaciers. In all cases, these aerofoil features are isolated from the main region of surging, are oriented at roughly the same angles in the parameter space, have long-axis lengths that scale nonlinearly with e / p , and have positions  Second, the boundary separating the surging region from the non-surging and abandoned surge regions is sharp, rather than diffuse, suggesting the existence of a supercritical Hopf bifurcation at the (approximately) linear boundary between surging and non-surging in the t h -b parameter space. As expounded on in the Discussion, this sharp boundary and possible bifurcation illuminates some potential mechanisms that cause surging to switch on and off over longer (multi-centennial) time scales in a given glacier system, and for surging glaciers to be relatively rare and geographically clustered. We reserve for future work detailed exploration of bifurcations in the system. To better understand the features in figure 5, we explore the dynamics in figure 6, which shows that small variations in b for fixed values of t h and e / p lead to a range of responses. The parameter values represented in figure 6 are shown with corresponding colours in figure 5. In order of decreasing b, we observe surging following the perturbation (blue line; b = 0.03), abandoned surging (orange line; b = 0.028), an abandoned surge followed by a surge at longer time scales (red line; b = 0.026) and slight dynamical variations (green and olive lines; b ≤ 0.024). These transitions in dynamical behaviour as a function of decreasing b can be understood in the context of changes in μ, the internal friction coefficient of the till. The sensitivity of μ to changes in state increases with b, allowing for greater and more rapid reductions in the friction coefficientand, by extension, the shear strength of the till, τ t (lowest panel of figure 6)-at higher b values. Thus, higher b values lead to unstable acceleration immediately following the perturbation by allowing dynamic glacier thinning to drive a net reduction in the effective pressure, further decreasing the shear strength of the till. Slightly smaller b values in the abandoned surge region result in slightly smaller changes in μ, which creates a situation that is unfavourable to surging because the acceleration in basal slip rate is sufficiently fast to drive till evolution but not significant dynamic thinning of the glacier. As a result, the initial acceleration is facilitated by reductions in both the effective pressure and internal friction coefficient, but decreases in pore  case, small initial decreases in μ driven by state evolution allow for basal slip acceleration, which drives the till towards steady state and ultimately increases state θ beyond the initial steady-state value as the glacier slows. Since basal slip does not stagnate as it did in the previously discussed case, the till continues to evolve, eventually leading to compaction and commensurate increase in pore water pressure. This increase in pore water pressure drives a reduction in effective pressure that leads to glacier acceleration, which eventually becomes self-sustaining as the glacier thins and effective pressure drops.
We find good agreement between our model behaviour and observations of surge motion in two natural glaciers (figure 7). Our model reproduces both the timing and order of magnitude of the speed-up with b = 0.03 and t h = 3000 days and other parameters corresponding to values used in figures 3 and 4. Note that our focus in this study has been on the incipient acceleration phase of the surges. Simplifications in the model, namely the lack of an evolving subglacial hydrological system and consideration of extensional stresses in the momentum balance, prevent the model from decelerating [10]. The agreement between our model and these data, however, is encouraging as it suggests that the dilation and glacier-thinning time scales we consider in our model may work in concert to trigger glacier surges.

Discussion
At this point, we have derived and explored the behaviour of a fundamentally new dynamical model of incipient surge motion that considers the mechanics of subglacial till and ice flow. Few comparable models exist in the literature, thus we endeavour to develop the simplest model capable of capturing the salient physical processes of ice slipping due to deformation of beds composed of water-saturated till. As detailed later in this section, natural glacier systems will, of course, be more complex than our model. Nevertheless, our model evinces rich dynamical behaviours consistent with observations, suggesting that our model strikes an appropriate balance between capturing the salient physical processes while remaining simple enough to allow for physical insight.

(a) Mechanics of incipient surge motion
Rich dynamical behaviour in our model is driven by the interactions of the three factors that define the shear strength of the till τ t = (p i − p w )μ: the overburden pressure p i = ρ i gh, pore water pressure p w , and the rate-and-state-dependent internal friction coefficient μ = μ(u b , θ). To understand surge behaviour in glaciers with till-covered beds, it is important to recognize that pore water pressure tends to decrease as a result of dilation, which strengthens till and resists surge motion, while the internal friction coefficient can increase or decrease, often by small amounts. Rate-weakening internal friction (a − b < 0) can help to facilitate surges but is not a necessary condition as surges are possible with rate-strengthening friction coefficients (a − b > 0) under conditions that allow for reductions in effective pressure ( figure 5).
The key process governing incipient surge motion is suction caused by till dilation in relatively impermeable till. In this case, pore water pressure decreases in response to shear-driven dilation, and the drop in pore water pressure diminishes the ability of till to evolve to a new steady state. If hydraulic permeability is sufficiently low (i.e. if the diffusion time of the deforming till layer t h is sufficiently long), slowing of state evolution allows the glacier to accelerate for longer periods of time. This longer acceleration phase gives the glacier time to thin dynamically, which reduces the overburden pressure (p i ). In the region of the parameter space shown in figure 5, the reduction in overburden pressure outpaces drops in pore water pressure (p w ), leading to a net reduction in the effective pressure (N = p i − p w ) and thereby the shear strength of till (τ t = μN). From equations (2.24) and (3.4), we can see that the rate of change in driving stress iṡ τ d ≈ 2ṗ i α, indicating that driving stress evolves at least an order of magnitude more slowly than changes in overburden owing to the shallow slopes of glaciers (α 1). As a result, reductions in overburden pressure facilitate sustained excess driving stress (τ d > τ b ), the key ingredient for sustained incipient surge motion. It is necessary, then, that the initial acceleration be large enough and last for long enough to generate sufficient dynamical thinning of the glacier.

(b) Implications of surge mechanics
The need for dynamic thinning to sustain surge motion gives two necessary conditions for glacier surging: till must have sufficiently low hydraulic permeability to allow for incipient surge motion to be maintained over a long enough period of time, and the velocity during the nascent surge must exceed the balance velocity to allow for dynamical thinning. The latter condition implies a third necessary condition: shear strength of the till must be less than the balance driving stress, defined as the driving stress at which the balance velocity is achieved through internal deformation of the ice column. Consequently, yielding of the till must occur at glacier velocities slower than the balance velocity to allow for continual shear loading of the till.
In the accumulation zones of surging glaciers, flow speeds must be slower than the balance velocity to build an ever-thickening reservoir of ice [15]. This condition must persist throughout the quiescent phase because, once the flow speed reaches the balance velocity, there would be no way to further increase driving stress and load the bed as ice mass would be evacuated by flow accommodated through vertical shearing of the ice column. In other words, mass balance along with the geometric and rheological properties of surge-type glaciers allow them to build a reservoir that exerts a driving stress equal to bed failure strength before flow rates reach the balance velocity. To illustrate this point, consider that the maximum shear stress a glacier can apply to its bed is given by the gravitational driving stress when the surface velocity of the ice equals the balance velocity and basal slip rate is negligible (τ b ≈ τ d ). Surface velocity due solely to vertical shearing within the ice column u v is given by assuming that stress increases linearly with depth, that ice rheology is constant with depth and that ice flow is parallel to the ice surface, yielding where A is the prefactor and n is the stress exponent in the constitutive relation for ice (equation ( where equation (4.2b) comes from recognizing that α = τ d * /(ρ i gh). The variablesṀ, A and, to a lesser extent, ρ i and n are governed by local climate [94]. Although mass density cannot vary more than 25% and n should be approximately 3,Ṁ and A can vary independently by orders of magnitude. Thus, the balance driving stress, τ d * , for an idealized glacier is determined primarily byṀ/A, the ratio of mass balance,Ṁ, to the rate factor, A, the latter of which depends on ice temperature and interstitial meltwater content, along with crystallographic fabric [104]. Equation (4.2) underpins a necessary condition for surging: at a minimum, surging glaciers must have a climate, and geometry, that allows for sufficiently high balance driving stresses (τ d * )-a combination of high, positive mass balance and stiff ice (small A)-to overcome the strength of their beds. As a result, the geographic distribution of surge-type glaciers will reflect areas that combine sufficiently high rates of snowfall, relatively low summertime melt at the surface and cold, stiff ice with beds that have yield stresses below the respective τ d * but are strong enough to allow the glacier to develop driving stresses that allow for order-of-magnitude increases in ice flow during the surge. To get a rough estimate for the pre-surge driving stress needed to produce a given speed-up, let us assume that the pre-surge surface velocity, u s pre , in the region where a surge begins is primarily due to viscous deformation in the ice column (given by equation (4.1)) and that surface velocity at peak surge speeds, u s surge , is due primarily to basal slip (given by equation (2.26)). Taking the ratio u s pre /u s surge and rearranging the terms (recalling that h w) gives where τ t surge is the shear strength of the till when the glacier is flowing at peak surge speed.
Combining equation (4.3) with the balance velocity explicitly gives the necessary condition which to a good approximation is simplyτ t < τ d * , whereτ t is the long-term average shear strength of the till in the region where surges nucleate. The range of reasonable values on ρ i g is small, so, to a good approximation, whether a glacier meets the condition in equation (4.4) is determined primarily by mass balance, ice rheology, bed strength and cross-sectional aspect ratio (h/w). The conditions defined by equations (4.2)-(4.4) yield surge conditions discussed in previous observational studies. The dependence on mass balance is consistent with observations that have shown cumulative quiescent-phase mass balance to be a reliable predictor of surging on Variegated Glacier, Alaska [105,106]. The temperature-dependent ice rheology broadly agrees with the climatic trends reported in [12] for surge-type glaciers (figure 8). In this framework, warmer climate (and ice temperatures) require higher values of surface mass balance to satisfy the condition that the bed yields before the driving stress becomes high enough to cause the glacier to flow at the balance velocity through internal deformation within the ice.
Further insight into the spatial distribution and longer term evolution of surge-type glaciers can be gleaned from the boundaries between surge-type and non-surge-type glaciers illuminated in the permeability versus evolution effect parameter space (figure 5). The sharp, diagonal boundary between surging on non-surge behaviour suggests the existence of a bifurcation in the system and lies at values that are likely to be relatively rare in nature and closely linked to local lithology and degree of weathering. In particular, our model suggests that values of hydraulic diffusivity for till in surge-type glaciers fall in the lower range of observed values (∼ 10 −9 m 2 s −1 ) for the range of b values explored in this study. Such low hydraulic diffusivities are consistent with canonical values of permeability expected for fine-grain sediments and loams [50,94]. The need for such low values of hydraulic permeability and fine-grained sediments suggests a potential role for comminution and sediment transport in activating and deactivating surging over millennial time scales, though future work is needed to elucidate these connections.
The governing role of till dilation and evolving pore water pressure in our model points to further methods for testing the model in nature. In addition to the comparisons with data similar to those given in this study (namely figure 7 and the preceding discussion of geographical  , Q c is the activation energy that increases from 60 kJ mol −1 for T ≤ T * to 115 kJ/mol for T * < T ≤ 0 • C, and R = 8.314 J/(K · mol) is the ideal gas constant [94]. The colourmap is defined to capture the range of driving stresses typically found on Earth [42]. (Online version in colour.) distribution of surge-type glaciers), we propose that passive seismic data collected during the incipient surge phase would provide valuable insight into the salient processes and could be used to test our model. Passive seismic data are routinely used to estimate the seismic moment from which estimates of the bulk shear modulus can be gleaned. The shear modulus is sensitive to both the porosity and pore water pressure, and so can be used as a means to observe till dilation and variations in pore water pressure.

(c) Model limitations and future development
Our goal with this work is to better understand basal mechanics by developing a model for incipient surge motion in glaciers with till-covered beds. We do not attempt to capture all of the processes that may be important in initiating and sustaining glacier surges. As a result, our model has some limitations that provide avenues for future work.
A notable limitation is the lack of explicit treatment for evolution of the subglacial hydrological system during any stage of the surge or the quiescent phase. The influence of basal hydrological characteristics is manifested in the model through the system water pressure p w r , but we implicitly treat this water pressure as passive in the model development. A fully passive basal hydrological system is unlikely given the rapid, extreme changes in glacier dynamics that define a surge. During surges, significant volumes of till are displaced, filling most existing cavities, basal crevasses or channels that constitute the contemporaneous hydrological system [17]. This lack of explicit treatment for changes in p w r owing to till displacement leaves open the possibility that increases in basal water pressure caused by changes in the basal hydrological system can cause surges. What we have provided in this study are proposed mechanisms of incipient surge motion in glaciers with deformable beds that are not dependent on changes in the basal hydrological system. The existence of such a mechanism, which works equally well for temperate and polythermal glaciers, and observations of surges beginning in times of the year when there is little or no additional surface meltwater available to pressurize a basal hydrological system (e.g. during winter) support the hypothesis that it is the incipient surge motion that diminishes the efficiency of any extant hydrological system rather than changes in the hydrological system that lead to surges [10].
We do not explicitly consider enhanced melting of basal ice through frictional heating or viscous dissipation, the former being a key mechanism in some existing surge models [10]. The reasons for this exclusion are twofold. First is the model set-up. We focus on instabilities that may arise from till mechanics and assume that rapid flow is due to basal slip, and basal slip is due entirely to deformation of the till. Thus the ice is in stationary contact with the top layer of the till and vertical shearing within the ice column is negligible. The latter case minimizes melting of basal ice through viscous dissipation. For frictional heating, the rate of melt scales with the product of rate of sliding along the ice-bed interface (i.e. the velocity of the ice relative to the top of the till) and the drag at the interface. As a result, there is no frictional heating in our model because we do not allow sliding at the ice-till interface. The addition of sliding at the ice-till interface is an appealing avenue for future research as the effective rheology of this interface and the resultant heating are non-trivial given that deformation of the bed is through cataclastic flow, facilitated by boundary sliding and rolling of sediment grains, and should facilitate comminution. The second reason we exclude slip-induced melting is that melting only influences ice dynamics through changes in basal and pore water pressure. Without a reliable model for subglacial hydrology on deforming sediment, there is no way to effectively link basal melt rate and water pressure. However, we reiterate that our model produces surge-like behaviour without representing frictional heating or evolving subglacial hydrology, suggesting that these mechanisms are not necessary for surge-like behaviour. The existence of unstable sliding phenomena, such as earthquakes and landslides, that do not rely on frictional heating or evolving hydrology at the slip interface supports our supposition that these mechanisms are not general requirements for surges, though they undoubtedly influence glacier surge dynamics [10]. The only general requirement for unstable acceleration is the existence of mechanisms that allow for the sustainment of excess driving forces relative to resistive forces, which our model achieves through reductions in overburden (and, thereby, effective) pressure. Finally, our model does not capture the down-glacier propagation of mechanical, kinematic or basal water pressure waves [21,107,108]. This limitation arises from the fact that our model is essentially one-dimensional (in the vertical), meaning that we neglect the gradient of extensional (along-flow normal) stresses and strain rates (equations (2.26) and (3.3)) along with horizontal gradients in water pressure. During the quiescent phase, neglecting extensional stresses is reasonable in the upper accumulation zone where surges are prone to begin. Here, surface velocities tend to be slow and relatively consistent over large spatial scales, meaning that alongflow strain rates are small relative to the effective strain rate; since ice is a viscous fluid, low strain rates mean low stresses. During the surge, the surface velocities are high, with the exception of the period when surge waves are present, and velocities can be expected to have small spatial gradients [6,29]. A more complete model of glacier surges would include more terms of the stress divergence such that it could account for the propagation of surge motion through the glacier. This more complete model would be useful for further investigating the influence of glacier length on surge behaviour [10]. However, we consider our box-model analysis to be a prerequisite to more complicated flowline and three-dimensional studies, which we reserve for future work.

Summary and conclusion
In this paper, we develop a new model of incipient surge motion in glaciers with till-covered beds. Incipient surge motion in our model occurs in the absence of enhanced water flux to the bed, changes to the basal hydrological system, frictional heating due to slip at the ice-bed interface and freeze-thaw cycles in till. Our model is based on granular mechanics of the till and focuses on processes that can lead to unstable acceleration in glaciers with deformable beds. Our model is unique among existing surge models in that it accounts for till porosity and pore water pressure, and represents the evolution of internal friction, porosity and pore water pressure within the deforming till layer as functions of the rate and history of shearing within the deforming till layer. This combination of mechanisms allows for exploration of the rich dynamics that arise from changes in the three factors that govern the shear strength of till: ice overburden pressure, pore water pressure and the internal friction coefficient. To represent these factors, we adopt the phenomenological rate-and-state model commonly used in studies of slip on tectonic faults. We link the state variable, which encodes the history of basal slip, to till porosity and derive a model in which pore water pressure evolves as a result of changes in porosity and transport of pore water (i.e. Darcy flow) into and out of the deforming till layer.
We find that till dilation, and more specifically suction caused by the reduction of pore water pressure in response to dilation, is a fundamental control on incipient surge motion. This control arises from the need for dynamic thinning of the glacier to sustain surge motion by reducing the effective pressure at the bed. Glacier thinning is necessary because, following a perturbation, till tends toward a new steady state while flow of water into and out of the deforming layer acts to equalize pore water pressure between the underlying static till layer, the deforming till layer and the subglacial hydrological system. As a result, the shear strength of the bed tends to a new steady state, leading to stable acceleration, unless the glacier thins. If the permeability of the till is sufficiently low, the evolution of the till to a new steady state is slow enough to allow accelerated surge motion to thin the glacier, as long as flow speeds during the nascent surge exceed the balance velocity of the glacier. Thinning of the glacier allows for unstable acceleration of the glacier owing to reductions in effective pressure, and, consequently, shear strength of the till, leading to order-of-magnitude increases in flow velocity that characterize surges and are consistent with observations of glacier acceleration during surges.
The hydromechanical properties of till, namely the need for low till permeability, required to induce rapid glacier thinning and surge motion give rise to restrictive conditions for glacier surges and rich dynamics. The necessary conditions for surging illuminated by our model are low hydraulic permeability in the deforming till layer, surge velocities that exceed the balance velocity and maximum shear strength of till that is less than the driving stress needed to achieve the balance velocity through vertical shearing in the ice column. These conditions are consistent with the rarity of surge-type glaciers; the geographical and climatic distribution and clustering of surge-type glaciers; and millennial time-scale evolution of surge behaviour. Furthermore, the rich dynamics produced by our model allow for abandoned surges along with a spectrum of surgelike behaviours that are consistent with kinematic observations of natural glaciers but are lacking in existing surge models.
Our model is necessarily simplified but contains important new physical processes-namely, till mechanics-that have been neglected in virtually all previous studies of glacier surges. To focus on the complex processes of water-saturated till, we deliberately ignore other processes that may be essential for a complete understanding of surge dynamics. Most notably, we neglect extensional stresses and vertical shearing in the ice column, and we treat the subglacial hydrological system as static. As a result, our model only captures the incipient surge phase and not slowdowns that terminate surges. We derive our model such that the inclusion of a dynamic subglacial hydrological system should be a relatively straightforward addition, and extension and vertical shear stresses can be included with the application of a more sophisticated flow model that accounts for more terms of the stress divergence in the momentum equations. These avenues provide numerous opportunities for future exploration of surge dynamics.