3.1 Effective pressure-dependent friction law

The relationship between basal velocity u and shear stress τ is strongly influenced by the presence of pressurized water at the bed, either via lowering the yield strength of subglacial till or increasing the extent of cavities between ice and bed (Cuffey and Paterson, Reference Cuffey and Paterson2010). When the effective pressure N (normal ice pressure p _i minus basal water pressure p _w) is large, velocity increases with stress, but as N decreases, high velocities can occur when basal stress is small (Lliboutry, Reference Lliboutry1968; Fowler, Reference Fowler1987b; Schoof, Reference Schoof2005). The possibility that effective pressure-dependent friction laws could lead to flow instabilities was first raised by Lliboutry (Reference Lliboutry1968), and this idea lies at the heart of most theories of surging, in which increases in basal water pressure trigger glacier acceleration (e.g. Budd, Reference Budd1975; Fowler, Reference Fowler1987a; Fowler and others, Reference Fowler, Murray and Ng2001; Kyrke-Smith and others, Reference Kyrke-Smith, Katz and Fowler2014).

3.2 Thermodynamic feedbacks

Feedbacks between ice velocity and frictional heating are a form of instability characteristic of stress-driven flows with temperature-dependent viscosity, and have been invoked as a possible cause of surges for well over half a century (e.g. Ahlmann, Reference Ahlmann1953; Robin, Reference Robin1955; Clarke, Reference Clarke1976; Clarke and others, Reference Clarke, Nitsan and Paterson1977; Murray and others, Reference Murray, Strozzi, Luckman, Jiskoot and Christakos2003). Velocity-heating feedbacks may manifest as warming and softening of cold ice (thermal runaway), water accumulation and enhanced basal motion at the beds of temperate glaciers (hydraulic runaway), or thawing and build-up of water at the beds of polythermal glaciers (a combination of the two). There is strong evidence that thermodynamic feedbacks triggered surges of small, relatively thin Arctic glaciers during the Little Ice Age (Lovell and others, Reference Lovell2015; Sevestre and others, Reference Sevestre and Benn2015).

Both thermal and hydraulic runaway mechanisms are incorporated in the surge model of Fowler and others (Reference Fowler, Murray and Ng2001), which exhibits cycling between cold-bedded (quiescent) and warm-bedded (surging) states. In addition, thermodynamic feedbacks are key elements of several models of ice flow instabilities (e.g. Yuen and Schubert, Reference Yuen and Schubert1979; MacAyeal, Reference MacAyeal1993; Fowler and Johnson, Reference Fowler and Johnson1995; Fowler and Schiavi, Reference Fowler and Schiavi1998; Tulaczyk and others, Reference Tulaczyk, Kamb and Engelhardt2000; Calov and others, Reference Calov, Ganopolski, Petoukhov, Claussen and Greve2002; Jay-Allemand and others, Reference Jay-Allemand, Gillet-Chaulet, Gagliardini and Nodet2011; van Pelt and Oerlemans, Reference Van Pelt and Oerlemans2012; Kyrke-Smith and others, Reference Kyrke-Smith, Katz and Fowler2014; Feldmann and Levermann, Reference Feldmann and Levermann2017).

3.3 Hydrologic switching

Drawing upon an unparalleled set of observations of the 1982–1983 surge of Variegated Glacier, Alaska, Kamb and others (Reference Kamb1985) argued that surges of temperate glaciers could occur due to switching of the basal drainage system from efficient conduits to an inefficient ‘linked cavity’ system. Inefficient drainage encourages high basal water pressures (low N) whereas efficient drainage reduces water storage (high N). Thus, switches in drainage system type could underpin oscillations between ‘slow’ and ‘fast’ modes of flow. This mechanism is at the core of Fowler's (Reference Fowler1987a) model of surging, which exhibits switches in drainage system type in concert with cyclic changes in ice thickness and velocity. The model assumed constant subglacial discharge, and did not consider variations in meltwater production via frictional heating or other sources.

Variations in meltwater production and discharge were incorporated in the surge model of Fowler and others (Reference Fowler, Murray and Ng2001). Thermal and hydraulic runaway processes lead to thawing and water accumulation in an inefficient Darcian-type drainage system during late quiescence, then this water is instantaneously removed from the system when N approaches zero. Rapid drainage at low N was justified by invoking efficient drainage through connected ‘blisters’, or regions of ice–bed decoupling (‘blister flux’), following Stone and Clarke (Reference Stone and Clarke1996). The model of Fowler and others (Reference Fowler, Murray and Ng2001) did not, however, include an explicit representation of such an efficient drainage system.

3.4 Balance flux vs flow mechanisms

The balance flux (the integrated surface mass-balance upstream of a specified flux gate on the glacier) determines the ice discharge required for a glacier to remain in dynamic equilibrium with climate. Budd (Reference Budd1975) argued that combinations of balance flux and velocity–friction relationships underpin the distinction between ‘ordinary’, ‘fast’ and ‘surging’ glaciers. On ‘ordinary’ glaciers, the balance flux is consistent with a slow flow mode, in which basal velocity and stress increase together, allowing glacier geometry to evolve to achieve a steady state. ‘Fast’ glaciers, in contrast, have a sufficiently high balance flux to sustain the flow in the fast mode, with high velocities and low effective pressures and basal stress. Surging glaciers oscillate between fast and slow modes, to discharge the balance flux over long periods of time.

Sevestre and Benn (Reference Sevestre and Benn2015) built upon this idea, arguing that both mass flux and enthalpy balances have to be satisfied for a glacier to remain in a stable steady state. For an incompressible ice mass, enthalpy relates to the temperature and liquid water content of the ice (Aschwanden and others, Reference Aschwanden, Bueler, Khroulev and Blatter2012). Enthalpy is gained at the glacier bed from the geothermal heat flux plus frictional heating (expenditure of potential energy) associated with the ice flux. Enthalpy losses occur by the dissipation of heat by conduction and loss of meltwater from the system. Ice flow rates depend sensitively on ice temperature and water storage at the bed, so for a glacier to maintain a steady flow, there must be a broad equality between enthalpy gains associated with the mass flux and enthalpy losses via conductive cooling and runoff. Sevestre and Benn (Reference Sevestre and Benn2015) postulated that under certain combinations of climate, geometry and bed geology, glaciers are unable to balance their mass and enthalpy budgets, and thus undergo mass and enthalpy cycles that manifest as quiescent and surge phases.

4.1 Overview

We now construct a general model of surging, combining elements described above. Following Fowler (Reference Fowler1987a) and Fowler and others (Reference Fowler, Murray and Ng2001), we adopt a non-spatial lumped model. We choose a lumped approach in preference to investigating specific geometries in 2-D or 3-D, because we seek to identify and explore the principles underlying the behaviour of glaciers in general, independently of the details of particular cases. In addition, a lumped approach allows a thorough exploration of model parameter space, in keeping with our aim of understanding the combined influence of climatic and geometric parameters at the glacier population scale. We represent the ice as an inclined, parallel-sided slab with slope θ and characteristic thickness H measured normal to the bed (Fig. 2). Variables can be interpreted as averages over the accumulation zone of the glacier.

Fig. 2.

Geometry of the lumped model, showing enthalpy sources (red) and sinks (blue). w is the depth of basal water and d is the thickness of the basal zone for which enthalpy is calculated. See text for other definitions.

Repeated surges represent periodic oscillations of the underlying dynamics of the glacier flow (Fowler, Reference Fowler1987a). The simplest system that can produce self-sustained oscillations is a pair of ordinary differential equations, here chosen to represent the evolution of ice thickness H and the enthalpy E in a thin layer spanning the basal ice and the ice–bed interface. This pair of ODEs represents the key principles at the heart of glacier dynamics: mass conservation and enthalpy conservation.

(1)$$\displaystyle{{{\rm d}H} \over {{\rm d}t}} = \lpar {\dot{a}-\dot{m}} \rpar -\displaystyle{{Q_{\rm i}} \over l},$$

(2)$$\displaystyle{{{\rm d}E} \over {{\rm d}t}} = \tau u + G-q_{\rm i}-\displaystyle{{\rho LQ_{\rm w}} \over l} + \rho L\beta \dot{m},$$

where $\dot{a}$ is the prescribed accumulation rate, $\dot{m}$ is the net ablation rate (i.e. meltwater that is not refrozen or otherwise stored in firn), Q _i is the ice flux, and l is the length (of the accumulation area strictly, though we assume this is proportional to the glacier length). τu is the frictional heating from sliding (τ is the basal shear stress and u the sliding speed), G is the geothermal heat flux, q _i is the conductive cooling flux, Q _w is the discharge rate of basal water, β is the proportion of surface melt that reaches the bed, ρ is the density and L is the latent heat of fusion. To simplify the calculations, we ignore the density difference between ice and water. This does not affect model physics, and given the approximate lumped nature of the model has a negligible impact on the results. The first, second and fifth terms on the right-hand side of (2) are enthalpy sources, and the third and fourth terms are enthalpy sinks. The final term on the right-hand side is water derived from the surface, which is an enthalpy source. For simplicity, we neglect viscous heating from deformation in the basal ice, which is generally very small compared to other terms (typically two orders of magnitude less than the geothermal heat flux for a unit thickness of basal ice).

We build the components of (1) and (2) from the following ingredients.

4.2 Mass and enthalpy fluxes

Net ablation is assumed to be linearly proportional to mean annual air temperature T _a:

(3)$$\dot{m} = {\rm DDF}(T_{\rm a}-T_{{\rm offset}})_ +, $$

where DDF is a degree-day factor and T _offset is the difference between mean annual temperature and the mean annual temperature below which there is no melting. We ignore seasonal cycles, so $\dot{m}$ is the mean annual melt rate.

Enthalpy E is related to basal temperature and water content according to:

(4)$$E = \rho c_{\rm p}d(T-T_{\rm m}) + \rho Lw,$$

where c _p is the specific heat capacity of ice, d is the thickness of the basal layer (taken to be d = 10 m), T is the basal ice temperature, T _m is the melting temperature (taken as constant at 0°C) and w is the average depth of water stored at the bed.

We scale the enthalpy variable such that negative values correspond to ‘cold content’ and positive values refer to water content. Since T is fixed at T _m if w > 0 (temperate condition), and w = 0 if T < T _m (cold condition), we can invert (4) to write

(5)$$T = \displaystyle{{E_-} \over {\rho c_{\rm p}d}}\quad w = \displaystyle{{E_ +} \over {\rho L}},$$

where E_ = min(E; 0) and $\displaystyle{{E_ +}}$ = max(E; 0). Because of the non-spatial nature of the model, the bed exists in either a cold or warm state, never a mixture of the two.

The ice flux Q _i is taken to be

(6)$$Q_{\rm i} = Hu + \displaystyle{{2A{\lpar {\rho g\sin \theta} \rpar }^n} \over {n + 2}}H^{n + 2},$$

with the first term being the contribution from sliding, and the second term being the contribution from internal deformation using the shallow ice approximation. A and n are the coefficients in Glen's flow law.

The sliding velocity u is determined from a generalized Weertman friction law:

(7)$$\tau = Ru^pN^q,$$

and we take basal shear stress to be equal to the driving stress τ = ρgHsinθ. The ice surface slope is set equal to the bed slope θ, which is a free input parameter in the model. N is the effective pressure and R is a constant roughness coefficient. This can be rearranged to write

(8)$$u = \left( {\displaystyle{{\rho g\sin \theta} \over R}} \right)^{1/p}H^{1/p}N^{-q/p}.$$

Equation (8) adds effective pressure-dependent sliding to the model (section ‘Effective pressure-dependent friction law’). We choose this relationship in preference to more sophisticated approaches because it captures the essential relationships between u, τ and N while simplifying the algebraic structure of the model.

The conductive cooling term q _i is approximated as a linear function of the temperature difference between the ice surface and the bed:

(9)$$q_{\rm i} = k\displaystyle{{\lpar {{\lpar {T-T_{\rm m}} \rpar }_--{\lpar {T_{\rm a}-T_{\rm m}} \rpar }_-} \rpar } \over H},$$

where k is the thermal conductivity of ice and T _a is the mean annual air temperature (the subtraction of T _m and the minus sign are to account for the fact that neither the base nor the surface of the ice may be warmer than the melting temperature, even if the air temperature is larger than that). Equation (9) assumes steady-state thermal conditions, which is reasonable because T _a is held fixed in model runs and basal temperature T evolves very slowly.

4.3 The drainage system

Framing a satisfactory description of the basal drainage system is not entirely trivial. Previous models of surging have tended to circumvent the problem or adopt steady-state approximations of the system. However, a large body of observations shows that temporal variations in subglacial water recharge, storage and discharge play a key role in the dynamics of both surge-type and non-surge-type glaciers (e.g. Cuffey and Paterson, Reference Cuffey and Paterson2010), and this demands that any general theory of surging must address the issue of drainage system evolution. We approach the problem in two steps. First, we define a single-component drainage system in which water flux is a function of mean water depth w at the bed. This can be visualized as an inefficient, distributed system when w is small (e.g. flow of porewater in till or other subglacial aquifers, or along the ice–bed interface in a system of small, poorly-connected cavities), becoming more efficient and better connected as the amount of stored water increases (e.g. in a system of canals or well-connected cavities at the ice–bed interface). Discharge through this efficient distributed system is analogous to the ‘blister flux’ of Fowler and others (Reference Fowler, Murray and Ng2001).

Second, we define a two-component system that can include both a distributed component and subglacial R-channels or conduits. R-channels can potentially draw down water storage when w is high but also allow high water flux when w is small if supply is large. In the two-component system, both elements compete to discharge water, and typically one will dominate at the expense of the other. Our aim is not to provide an explicit picture of the drainage system, but to investigate the degree to which surging behaviour depends on how subglacial drainage adapts under different circumstances.

We assume that there is an inverse relationship between the average depth of stored water w and effective pressure N. For convenience we set:

(10)$$N = \displaystyle{{\tilde{C}} \over w},$$

where $\tilde{C}$ is a prescribed constant. The water depth w can represent water storage within a compressible layer of till, or within a film or system of cavities at the ice–bed interface; in all of these cases, one expects a (steady-state) relation of broadly this form (Hewitt, Reference Hewitt2011, Reference Hewitt2013). The effective pressure cannot exceed the ice pressure ρgH, so in terms of water depth and enthalpy, respectively, we replace (10) with

(11)$$N = \min \left( {\rho gH,\displaystyle{{\tilde{C}} \over w}} \right) = \min \left( {\rho gH,\displaystyle{C \over {E_ +}}} \right),$$

where $C = \rho L\tilde{C}.$

The single-component, distributed drainage system is defined by:

(12)$$Q_{\rm w} = \tilde{K}w^\alpha \rho g\sin \theta = KE_ + ^\alpha, $$

where ρgsinθ is the hydraulic gradient arising from the bed slope. $\tilde{K}w^\alpha$ is the hydraulic transmissivity, which varies with water film depth w. $\tilde{K}$ thus represents the invariant components of conductivity (e.g. small-scale bed roughness, till porosity) and fluid viscosity. α is set to a relatively large value (5) to account for the idea that at low water volumes discharge is low, but at high volumes the discharge can increase rapidly. The water pressure component of the hydraulic potential gradient is neglected because it is a spatially variant property not readily incorporated in the present lumped model formulation.

The two-component system includes a second term representing water discharge in a conduit:

(13)$$Q_{\rm w} = \tilde{K}\lpar {\rho g\sin \theta} \rpar w^\alpha + \phi \displaystyle{{K_{\rm c}} \over {W_{\rm c}}}\lpar {\rho g\sin \theta} \rpar ^{1/2}S^{4/3},$$

in which K _c is a constant related to Manning's roughness coefficient, S is the cross-section area of conduits, spaced at width W _c. The half power of the hydraulic gradient derives from the assumption of turbulent flow in a semicircular channel. The factor ϕ is introduced as a ‘fill-fraction’ to allow for the possibility of channels being only partially filled at times of falling discharge. According to (11), this happens when $\rho gHw \lt \tilde{C}$ (corresponding to the drainage system being at atmospheric pressure), and we therefore take

(14)$$\phi (w,H) = \left\{ {\matrix{ {w\left( {\displaystyle{{\rho gH} \over {\tilde{C}}}} \right)} & {w \lt \left( {\displaystyle{{\tilde{C}} \over {\rho gH}}} \right)} \cr 1 & {w \ge \left( {\displaystyle{{\tilde{C}} \over {\rho gH}}} \right)} \cr}} \right.,$$

which makes flow in partially-filled channels proportional to the global water thickness w. Without this modification, a previously enlarged channel would unphysically continue to drain a large flux of water even when the supply drops off.

Channel cross-section area S evolves through time according to the kinematic condition:

(15)$$\displaystyle{{{\rm d}S} \over {{\rm d}t}} = \displaystyle{{K_{\rm c}{\lpar {\rho g\sin \theta} \rpar }^{3/2}S^{4/3}} \over {\rho L}}-\tilde{A}SN^n + \dot{S}_0.$$

$\tilde{A}$ is a rate factor analogous to A in Glen's flow law. The first term on the right-hand side represents channel enlargement by wall melting due to turbulent heat dissipation, the second term represents creep closure, and the third is a small opening term associated with wall melting driven by other heat sources (geothermal and frictional heat) or ice sliding over bed roughness elements (Schoof, Reference Schoof2010; Hewitt, Reference Hewitt2011; Werder and others, Reference Werder, Hewitt, Schoof and Flowers2013). We choose this approach over setting a minimum channel cross-section area since it allows the same Eqn (15) to hold all the time. Where two-component drainage is included in the model, (15) adds a third ordinary differential equation to the coupled system (1) and (2) (see Appendix).

The proportion β of surface melt that reaches the bed (Eqn (2)) can be set to zero or varied to represent surface-to-bed drainage through moulins nucleated on crevasses (cf. van der Veen, Reference Van Der Veen1998; Benn and others, Reference Benn, Kristensen and Gulley2009). The amount of surface-to-bed drainage depends on the extent of surface crevasses, which in turn is related to spatial phenomena (e.g. velocity gradients) that are not explicitly simulated. We assume that β is proportional to the ice velocity:

(16)$$\beta = \left\{ {\matrix{ 0 \hfill & {u = u_1} \hfill \cr {(u-u_1)/(u_2-u_1)} \hfill & {u_1 \lt u \lt u_2} \hfill \cr 1 \hfill & {u \ge u_2} \hfill \cr}} \right.,$$

i.e. ice with a velocity less than u ₁ is assumed to be uncrevassed and impermeable, whereas if velocity equals u ₂ or above, the ice is pervasively crevassed and all net surface meltwater reaches the bed. This is the simplest possible formulation that allows the model to explicitly exhibit feedbacks between surface-to-bed drainage and ice dynamics.

Variables were scaled and non-dimensionalized prior to model implementation. Details are given in the Appendix and parameter values are listed in Table 1. The model is structurally similar to that of Fowler and others (Reference Fowler, Murray and Ng2001), which solved coupled ordinary differential equations for ice thickness H and effective pressure N, rather than H and E. Indeed, Fowler and others (Reference Fowler, Murray and Ng2001) defined N in a way that closely matches our definition of E. The main innovations in our model are: (1) partitioning of surface mass balance into accumulation and melt components to allow investigation of their influence on surging behaviour; (2) more detailed treatment of basal hydrology, particularly the introduction of a time-evolving two-component drainage system; and (3) incorporation of surface-to-bed drainage. In our model, subglacial discharge and drainage efficiency evolve in response to changing water sources and sinks, unlike the original hydrological switch model of Fowler (Reference Fowler1987a) in which discharge was assumed to be constant.

Table 1.

Default parameter values, scales and dimensionless parameters

Due to the approximate, lumped nature of the model, the values of some constants have been rounded. The flow-law parameter A is taken to be a representative depth-averaged value (across the full ice thickness, since the place it enters is in the shallow-ice flux) and therefore smaller than the textbook value at 0°C.

6.1 Predictions of enthalpy balance theory

The model results can be summarized as a set of hypotheses that can be tested against existing or future observations. For the simplest implementation of the model (single-component drainage system and no surface-to-bed drainage), enthalpy balance theory predicts:

1. (1) Non-surge-type or ‘normal’ glaciers have stable steady states with respect to their mass and enthalpy budgets. In terms of E–H phase portraits, the nullcline intersections are point attractors and lie on the upper or lower branches of the E-nullcline (Fig. 3).

2. (2) Surge-type glaciers have unstable steady states with respect to their mass and enthalpy budgets. These steady states are unstable to small perturbations, which set in train feedbacks between enthalpy balance and dynamics, which in turn lead to periodic mass, enthalpy and velocity cycles. In phase portraits, solutions follow limit cycles, orbiting around nullcline intersections that lie on the middle branch of the E-nullcline.

3. (3) Climate, glacier geometry and bed properties determine whether glaciers have stable or unstable steady states.

4. (4) In cold, dry climates, glaciers have low balance fluxes and are low-enthalpy producers; enthalpy production can be balanced by conductive heat losses due to low mean annual air temperatures, particularly where ice is thin.

5. (5) In warm, humid climates, glaciers have high balance fluxes and are high-enthalpy producers. High enthalpy at the bed allows enthalpy production to be offset by high subglacial water discharge. In the single-component model, high discharge occurs when water storage w is high (12).

6. (6) In intermediate climates, neither conduction nor water discharge can evacuate enthalpy at steady rates that match the production from geothermal heat flux plus frictional heating (expenditure of potential energy). The resulting steady-state conditions are unstable, and ice thickness, enthalpy and velocity undergo periodic oscillations.

7. (7) The model predicts the following sequence of processes over the course of surge cycles. During quiescence, velocities are less than balance velocities, and mass accumulates. Enthalpy at the bed is low at the start of the quiescent period, and gradually increases because geothermal heating exceeds enthalpy losses. Once enthalpy rises high enough to initiate sliding, feedback between frictional heating and sliding speed causes rapid flow acceleration and transition to the surge phase. The surge terminates when basal water is lost faster than it is produced. Depending on the environmental conditions, the glacier bed may remain entirely temperate during the surge cycle, or enter a cold-based condition during quiescence. In both cases, the physical processes underlying the surge cycle are essentially the same.

8. (8) Glacier length and slope influence the distribution of unstable states with respect to accumulation and mean annual air temperature. In our model, length affects dynamics through its impact on the balance flux and subglacial discharge, and slope affects dynamics through its impact on the hydraulic gradient and basal shear stress. In cold, dry climates, long, gently sloping glaciers are more likely to be surge-type than short, steep glaciers.

9. (9) Basal hydraulic conductivity influences the distribution of surging glaciers by altering drainage system efficiency (strictly, this statement applies to our hydraulic conductivity parameter $\tilde{K}$, which has different units to hydraulic conductivity as normally defined; here we use ‘hydraulic conductivity’ for brevity, with this caveat). Low hydraulic conductivity (inefficient drainage) encourages surging behaviour in warm, humid climates, whereas high conductivity (efficient drainage) encourages stable steady states. Variations in length, slope and hydraulic conductivity in a population of glaciers mean that the distribution of surge-type glaciers and non-surge-type glaciers will overlap.

10. (10) Climate change can alter the geographical distribution of surge-type glaciers.

When surface water is allowed to access the bed via crevasses as flow speeds increase, additional behaviour emerges:

1. (11) Surface water input can influence surge phase dynamics by increasing basal enthalpy and ice velocities more rapidly than is possible by velocity-heating feedbacks alone. The influence on quiescent phase duration is complex, however, and the surge cycle may be shortened or lengthened depending on input parameters.

2. (12) Inputs of surface water eliminate unstable steady states in part of the warm–humid end of the climate envelope. Where surface-to-bed drainage occurs, stable flow speeds can be maintained by surface meltwater inputs rather than through frictional heating alone. This statement applies for the single-component drainage model, where high water throughput corresponds to high basal water storage.

Finally, the addition of a two-component drainage system with both distributed and channelized elements allows hydraulic switching to influence the dynamics:

1. (13) Switching between distributed and channelized modes of the drainage system can occur over the course of surge cycles. A low-enthalpy distributed system exists during quiescence, which gradually transitions to a high-enthalpy distributed system during late quiescence and the surge phase. When the effective pressure N becomes very low, reduced creep closure allows rapid development of an efficient channel system, terminating the surge. Thus, hydrologic switching (from a distributed to a channelized system) triggers surge termination, and the opposing switch (channels to a distributed system) occurs when the bed runs out of water. The switch to a distributed drainage system does not trigger a surge, but sets up the initial conditions necessary for the build-up of mass and enthalpy during quiescence.

2. (14) It is possible that stable steady states exist with low water storage and high throughput, with efficient R-channels maintained by surface-to-bed drainage. In the present model, these high-throughput states have almost constant ice thickness and low-velocity throughout, but the subglacial channel undergoes spontaneous oscillations in size.

6.2 Comparison with observations

Some of the hypotheses listed above can be tested with existing data, while others will require new, targeted observations before rigorous testing is possible.

Predictions of the theory are consistent with global-scale data on the distribution of surge-type glaciers with respect to mean annual precipitation and air temperature (Sevestre and Benn, Reference Sevestre and Benn2015; Figs 1, 6, 8). Importantly, we have made no attempt to tune model parameters to optimize the match between model results and observations; such an attempt would have little value given the idealized, lumped model geometry, the approximate nature of key equations, and the indicative character of the ERA-I climatic data. Given this lack of tuning, the similarity between model predictions and observations is striking, particularly the concentration of surge-type glaciers in climates intermediate between cold–dry and warm–humid environments, and the overlapping distributions of surge-type and non-surge-type glaciers with different geometries.

The model results are consistent with the observed relationships between glacier geometry and surging behaviour, and indicate that length and slope are distinct influences and not simply two perspectives on a single dependency (cf. Clarke, Reference Clarke1991). Notably, the model correctly predicts that in cold, dry environments (i.e. Arctic Canada: Fig. 1), only long, low-gradient glaciers will be of surge type. In the model, short glaciers in dry climates have low balance fluxes and are low-enthalpy producers; steeper glaciers are thin and lose enthalpy efficiently by conduction to the cold atmosphere. New observations are needed to further test this model prediction.

In the model, low hydraulic conductivity encourages surging behaviour because it reduces the ability of the subglacial drainage system to evacuate enthalpy from the bed. This result offers an explanation for the observed association between surge-type glaciers and weak sedimentary rocks (e.g. Jiskoot and others, Reference Jiskoot, Murray and Boyle2000). Glaciers eroding weak lithologies (e.g. shales) tend to produce thick, fine-grained tills with lower hydraulic conductivity, likely to impede water flow at the bed (Clark and Walder, Reference Clark and Walder1994). However, substrate type likely influences glacier dynamics in complex ways, and detailed observations are required to further illuminate the relationships between enthalpy balance and substrate type.

Surge cycle lengths are very long (hundreds of years) in the examples illustrated in Figures 4, 10, 12, reflecting the low net accumulation prescribed in these cases. Examples with higher accumulation (not shown) exhibit shorter cycle lengths (decades), consistent with the idea that accumulation is a key control on quiescent phase duration (e.g. Dowdeswell and others, Reference Dowdeswell, Hodgkins, Nuttall, Hagen and Hamilton1995; Eisen and others, Reference Eisen, Harrison and Raymond2001). However, cycle length must also reflect the amount of ice drawn down during the surge phase, which determines the total amount of build-up that must occur before the next surge is initiated. Surge-phase draw-down is a function of ice velocity and surge duration, which in turn reflect a complex web of dynamic and hydrological processes. Closer correspondence between modelled and observed surge periods could be achieved by tuning several possible combinations of the parameters governing relationships between stress, effective pressure, ice velocity, surface-to-bed drainage, water storage and subglacial discharge (e.g. R, p, q, $\tilde{C}$, α, β, u ₁ and u ₂ in Eqns (7), (10), (12), (13) and (16)). We have not attempted this, because any parameter combination choice would be arbitrary and the exercise would not yield additional insights into how these processes interact to control cycle length in real-world cases. Careful work with fully spatial models will be needed to disentangle the roles of dynamic and hydrological processes in controlling surge cycle length.

Our model runs focus on glaciers under constant climate conditions, but the predicted climate dependence of surging implies that the spatial distribution of surge-type glaciers could shift in response to climate change. This is consistent with evidence from Svalbard, where many small glaciers are known or inferred to have surged during the Little Ice Age, but now have strongly negative mass balance and can no longer build up the necessary mass to maintain the surge cycle (James and others, Reference James2012; Sevestre and others, Reference Sevestre and Benn2015). Indeed, it may be said that these glaciers are no longer in quiescence, but a climatically induced senescence. In terms of the climatic parameter space shown in Figure 6, this is equivalent to moving a glacier from the surge regime into the region of glacier unviability in the top left of the diagram. A similar regime change may have affected glaciers in the Pyrenees that appear to have been surge-type during the Little Ice Age (Cañadas and Moreno, Reference Cañadas and Moreno2018). In the Austrian Alps, Vernagtferner repeatedly surged during the Little Ice Age, but ceased to do so in the warmer climate of the 20th century (Hoinkes, Reference Hoinkes1969; Kruss and Smith, Reference Kruss and Smith1982). Evidence from lake sediments indicates that a change in the opposite direction affected Eyjabakkajökull in Iceland, which switched from non-surge type to surge type ~2200 years ago (Striberger and others, Reference Striberger2011).

Several studies have attributed changes in surge magnitude and frequency to climate change. For example, Frappé and Clarke (Reference Frappé and Clarke2007) showed that the ‘slow’ 1980–2000 surge of Trapridge Glacier lacked the fast flow phase that occurred during the previous surge in the 1940s, probably due to diminished mass build-up. In contrast, Copland and others (Reference Copland2011) argued that recent temporal clustering of surges in the Karakoram occurred because surge cycle length has been shortened by an increase in precipitation. Dunse and others (Reference Dunse2015) proposed that a major surge of the Basin 3 outlet glacier of Austfonna, Svalbard, was triggered by increased amounts of surface meltwater penetrating crevasses and thawing formerly cold ice. Less convincingly, Nuth and others (Reference Nuth2019) suggested that some recent surges in Svalbard were triggered by near-margin thinning and freezing. In all these cases, climatic factors appear to have influenced the timing and/or magnitude of surges, rather than initiated surging behaviour in hitherto non-surge-type glaciers. As such, these cases are qualitatively different from the climatically driven regime changes implied by the current suite of model runs. Implications of enthalpy balance theory for understanding glacier dynamic response to transient climatic states is the focus of ongoing work.

Enthalpy balance theory makes the important prediction that polythermal and temperate glacier surges occur by essentially the same mechanism. In both cases, an excess of enthalpy gains relative to losses results in the accumulation of water in a distributed drainage system. In polythermal glaciers, the bed may undergo an intermediate cycle of freezing and thawing during quiescence, although it should be noted that larger polythermal glaciers may remain largely warm-based throughout the surge cycle (Sevestre and others, Reference Sevestre and Benn2015). The prediction that polythermal and temperate glacier surges have a common dynamical basis is consistent with population-scale data on the distribution of surge-type glaciers, and the observation by Frappé and Clarke (Reference Frappé and Clarke2007) that thermal regime is ‘collateral to the surging phenomena rather than essential’. At the individual glacier scale, detailed mass and enthalpy budget analysis of the polythermal glacier Morsnevbreen, Svalbard, has shown that enthalpy sources and sinks underwent the predicted sequence of changes over the course of its most recent surge cycle (Benn and others, Reference Benn2019). Comparative analyses will be conducted for temperate surge-type glaciers to further test the theory.

The theory also predicts that temperate glacier surges can occur in the absence of switches in drainage system configuration. Introduction of a two-component drainage system affects some of the details of surge evolution, but not the underlying relationships between mass and enthalpy budgets for glaciers with unstable steady states. Importantly, and in contrast with previous models of hydrologic switching (e.g. Jay-Allemand and others, Reference Jay-Allemand, Gillet-Chaulet, Gagliardini and Nodet2011; Mayer and others, Reference Mayer, Fowler, Lambrecht and Scharrer2011), in our model, a transition from R-channels to a distributed system does not act as a surge trigger. Rather, a channel-to-distributed transition occurs after surge termination when basal water storage is drawn down, and simply resets the system in a low-storage state at the start of quiescence. This offers a solution to a long-standing difficulty with hydrological switch theory, namely, exactly what mechanisms could trigger a switch from channelized drainage to a high-storage distributed system, and why this should occur on some glaciers but not others (Harrison and Post, Reference Harrison and Post2003). Moreover, it is important not to view distributed and channelized drainage systems in terms of a simple binary switch. Both drainage types can exist simultaneously on different parts of the same glacier, and in principle water could accumulate in a distributed system under some parts of a glacier system while other parts continue to be drained by surface-fed conduits (Fischer and Clarke, Reference Fischer and Clarke2001; Benn and others, Reference Benn, Kristensen and Gulley2009; Schoof and others, Reference Schoof, Rada, Wilson, Flowers and Haseloff2014). If a surge is initiated, conduits could then be destroyed following surge propagation, in which case a switch from channelized to distributed drainage would be a consequence rather than the cause of a surge. The difficulty of determining drainage system configurations means that testing hydrological aspects of surge theory presents major challenges.

In all cases, our model predicts abrupt termination of surges due to the discharge of water from the bed, either via a high-storage distributed system (‘blister flux’) or a switch to an efficient channelized system. This is consistent with the previous work on the hydrology of surge-type glaciers that has shown that abrupt surge terminations are triggered by step changes in drainage efficiency (e.g. Kamb and others, Reference Kamb1985; Kamb, Reference Kamb1987). Not all surges terminate abruptly, however, and some instead undergo gradual slowdowns spread over multiple seasons (e.g. Luckman and others, Reference Luckman, Murray and Strozzi2002; Sund and others, Reference Sund, Lauknes and Eiken2014; Sevestre and others, Reference Sevestre2018; Benn and others, Reference Benn2019). Gradual slowdowns are typical of polythermal glaciers in Svalbard (though not all, see Kristensen and Benn, Reference Kristensen and Benn2012), and may reflect the removal of hydraulic barriers during surge-front advance rather than a drainage system switch. Simulation of these effects will require a spatial treatment of basal hydrology.

For combinations of input parameter values that allow a large volume of surface meltwater to access the bed, the model predicts almost stable steady states with respect to mass and velocity but with oscillations of channel size and basal water storage. Stable steady states with respect to both mass and water storage were not found with the present model. This reflects the simplified representation of hydrology in the lumped model, particularly the assumptions that the effective pressure in the channel and distributed system are the same, and that the hydraulic gradient is constant and determined by the bed slope. Indeed, a slightly modified version of the present model (not shown) can allow for the possibility of stable channels at high discharge. On real glaciers, surface meltwater production (and hence surface-to-bed drainage) fluctuates on diurnal, seasonal and intermediate timescales, so examples of such steady states are unlikely to occur in nature. It is also notable that spontaneous subglacial water pressure oscillations unrelated to water supply variations have been observed (Schoof and others, Reference Schoof, Rada, Wilson, Flowers and Haseloff2014), and modelling indicates that these are due to conduit growth and decay cycles analogous to the jökulhlaup instabilities described by Nye (Reference Nye1976). This raises the possibility that unforced drainage system oscillations can occur across a wide range of timescales, with a corresponding range of impacts on ice dynamics. This may account for some short-term dynamic instabilities such as ‘mini-surges’ and pulses (e.g. Humphrey and others, Reference Humphrey, Raymond and Harrison1986; Kamb and Engelhardt, Reference Kamb and Engelhardt1987; Turrin and others, Reference Turrin, Forster, Sauber, Hall and Bruhn2014; Herreid and Truffer, Reference Herreid and Truffer2016).

Our model experiments do not include seasonal cycles, but the predicted impact of surface water input is consistent with the observed timing and location of rapid speed-up events during surges. In Svalbard, several studies have shown that expansion of crevasse fields during the early stages of surges has created templates for areas of flow acceleration at ablation season onset (Dunse and others, Reference Dunse2015; Sevestre and others, Reference Sevestre2018; Gong and others, Reference Gong2018; Benn and others, Reference Benn2019). The role of surface-to-bed drainage in stabilizing the flow is illustrated by Kronebreen, a perennially fast-flowing tidewater glacier in the midst of the Svalbard surge cluster. Convergent flow on the upper glacier focuses accumulation from a large catchment area, giving Kronebreen a higher balance flux than its neighbours. The effect of convergent flow is equivalent to shifting the glacier into a higher-precipitation environment than the adjacent glaciers. High velocities are maintained by water reaching the bed via the intensely crevassed glacier surface (How and others, Reference How2017; Vallot and others, Reference Vallot2017).

6.3 Model formulation: limitations and future prospects

In keeping with our aim of exploring fundamental principles underlying glacier dynamics, we have adopted a lumped model approach incorporating simple equations to represent system components. Inevitably, this approach glosses over much fine detail, and representations of processes such as conduction of heat, basal sliding and hydrology are necessarily approximate. The model results, however, are robust for plausible values of most input parameters. Model sensitivity is greatest for the exponent α in the drainage law (12). To achieve the stable upper branch of the E-nullcline, α must be larger than 1/p (the exponent in the sliding-law: Eqn (8)). In physical terms, this means that for stable steady states to occur under warm, humid conditions, drainage must be efficient enough to cope with the energetic consequences of faster ice flow (i.e. frictional heating). Since non-surge-type glaciers are common in relatively warm, humid environments, this feature of the model likely reflects an important principle, although of course only in an approximate, generalized way.

The present implementation of the model does not allow free evolution of the ice-surface slope and omits the pressure term in the hydraulic potential gradient. These approximations likely limit the range of behaviour exhibited by the model, but were unavoidable because of its non-spatial character. More sophisticated treatment of these system components will require a fully spatial model, either in 2-D (flowline) or 3-D. This will allow particular processes and geometries to be examined in more detail, but at the expense of generality.

The present lumped model is based on the assumption that surge initiation occurs in a reservoir zone where accumulation exceeds net ablation, whereas some surges initiate in glacier ablation zones (e.g. Björnsson and others, Reference Björnsson, Pálsson, Sigurđsson and Flowers2003; Dowdeswell and Benham, Reference Dowdeswell and Benham2003; Murray and others, Reference Murray2012; Sevestre and others, Reference Sevestre and Benn2015; Abe and others, Reference Abe, Furuya and Sakakibara2016). Potentially, this type of surge behaviour could be accommodated in the model by allowing ice inflow (i.e. modifying the flux divergence term in Eqn (2)), thus permitting thickening to occur even if the surface mass balance is negative. This possibility will be explored in future work.

Many important phenomena associated with glacier surges are inherently spatial. It is widely accepted that spatial variations in basal thermal regime strongly influence the trajectories of polythermal glacier surges (e.g. Clarke and others, Reference Clarke, Collins and Thompson1984; Fowler and others, Reference Fowler, Murray and Ng2001; Murray and Porter, Reference Murray and Porter2001). For example, cold-based ice near glacier termini may impede water discharge and exacerbate instabilities associated with basal water accumulation in warm-bedded regions up-glacier. Similarly, spatial variations in basal drainage systems (including the location of water sources and sinks) likely play a key role in surge propagation (e.g. Jay-Allemand and others, Reference Jay-Allemand, Gillet-Chaulet, Gagliardini and Nodet2011; Sevestre and others, Reference Sevestre2018). Spatial phenomena can be considered in terms of enthalpy gradients, and the advection of enthalpy between sources and sinks. A fully spatial development of enthalpy balance theory has the potential to shed light on individual cases of surge initiation and propagation, and to explore the ways that local conditions manifest general principles. For further discussion of ways in which spatial phenomena can modulate enthalpy and mass fluxes, see Benn and others (Reference Benn2019).