Articles | Volume 17, issue 11
Research article
15 Nov 2023
Research article |  | 15 Nov 2023

The evolution of isolated cavities and hydraulic connection at the glacier bed – Part 1: Steady states and friction laws

Christian Schoof

Models of subglacial drainage and of cavity formation generally assume that the glacier bed is pervasively hydraulically connected. A growing body of field observations indicates that this assumption is frequently violated in practice. In this paper, I use an extension of existing models of steady-state cavitation to study the formation of hydraulically isolated, uncavitated, low-pressure regions of the bed, which would become flooded if they had access to the subglacial drainage system. I also study their natural counterpart, hydraulically isolated cavities that would drain if they had access to the subglacial drainage system. I show that connections to the drainage system are made at two different sets of critical effective pressure, a lower one at which uncavitated low-pressure regions connect to the drainage system and a higher one at which isolated cavities do the same. I also show that the extent of cavitation, determined by the history of connections made at the bed, has a dominant effect on basal drag while remaining outside the realm of previously employed basal friction laws: changes in basal effective pressure alone may have a minor effect on basal drag until a connection between a cavity and an uncavitated low-pressure region of the bed is made, at which point a drastic and irreversible drop in drag occurs. These results point to the need to expand basal friction and drainage models to include a description of basal connectivity.

1 Introduction

Subglacial drainage is often assumed to occur in part through a “distributed” drainage system: connected conduits that are not arborescent in their geometry (Fountain and Walder1998) and therefore do not localize drainage into a few large channels (Hewitt2010; Schoof et al.2012; Hewitt2013; Werder et al.2013; Rada and Schoof2018; Flowers2015). A frequently used paradigm for a distributed drainage model is that of linked cavities (Lliboutry1968; Kamb1987; Fowler1987): localized areas of ice–bed separation in the lee of bed bumps.

Large-scale models for subglacial drainage systems typically assume that the bed as a whole always remains hydraulically connected. Existing process-scale models for the evolution of subglacial cavities generally make the same assumption. In large-scale drainage models, cavities are represented by a water sheet thickness: a cavity depth averaged over a representative small area of the bed (that is, an area of the bed that is much larger than an individual cavity but much smaller than the glacier as a whole). The assumption of a connected bed here simply means that water can flow as soon as the sheet thickness exceeds zero (e.g. Werder et al.2013; Sommers et al.2018). As a result, local variations in water pressure at the scale of individual cavities are small, since they would otherwise lead to excessive water fluxes, and water pressure is a well-defined, smoothly varying variable in the large-scale model.

In process-scale models, hydraulic connectedness typically occurs through the bed itself: the bed is highly permeable. Water sourced from an ambient drainage system at some given water pressure can force its way between the ice and bed as soon as compressive normal stress at the base of the ice drops to the water pressure in the ambient drainage system, causing a cavity to form (Schoof2005; Gagliardini et al.2007; Helanow et al.2020, 2021; Stubblefield et al.2021; de Diego et al.2022, 2023).

These assumptions are at odds with a growing set of observations (Hodge1979; Murray and Clarke1995; Andrews et al.2014; Lefeuvre et al.2015; Rada and Schoof2018) indicating that hydraulic connections at the glacier bed are often patchy and evolve in time: while the bed itself may be somewhat permeable, that permeability is too low to allow significant water transport on the timescales over which the drainage system evolves. On these timescales, water must then flow predominantly along the ice–bed interface, and the topology of the conduit network present there (consisting of subglacial cavities and other forms of void space like R-channels) may not provide a connection to all parts of the bed.

Recent work in large-scale drainage modelling has attempted to address this issue (Hoffman et al.2016; Rada and Schoof2018), albeit in fairly crude form: for instance, one possibility is to assume that water can only flow when sheet thickness exceeds a critical value. The aim of the present work is to study the evolution of cavities in more detail for an effectively impermeable bed at the process scale to better understand how an ambient active drainage system can access other parts of the bed through the evolution of basal cavities. By contrast with most studies of subglacial cavity formation, my focus is mostly on the evolution of subglacial connectivity rather than on the computation of a sliding law. As a by-product, I also show that connectivity plays a major role in controlling friction at the glacier bed.

If only part of the glacier bed has access to the ambient drainage system, then isolated, uncavitated, low-pressure regions can form elsewhere, at normal stresses that would lead to ice–bed separation if water from the ambient drainage system had access. Conversely, these distant parts of the glacier bed can become flooded with water when connected cavities grow at low effective pressure. If the effective pressure in the drainage system increases again after that flooding, the intervening connections can become closed, leaving isolated cavities of fixed volume. These isolated cavities will generally be at different effective pressures than the connected drainage system.

In the present work, I have used a modified mathematical model for cavity formation to explore the physics involved. The basic physics of ice flow over an undulating bed, allowing for the possibility of ice–bed separation as water forces its way between the two, are the same as in existing models for subglacial cavity formation. However, only a pre-defined, highly permeable part of the bed, denoted by P, is assumed to be directly connected to the ambient drainage system: as in the existing models of de Diego et al. (2022, 2023), Gagliardini et al. (2007), Helanow et al. (2020, 2021), Schoof (2005), and Stubblefield et al. (2021), water is assumed to force its way between the ice and bed if compressive normal stress on P drops to the value of the water pressure in the ambient drainage system. The remainder of the bed is assumed to be completely impermeable. Water can access these other parts of the bed interface (outside of P) only if there is a hydraulic connection to P along the ice–bed interface. Moreover, if water has previously accessed some impermeable part of the bed and the hydraulic connection has subsequently been closed, then an isolated cavity is formed. The water pressure in that isolated cavity can differ from the water pressure in the ambient drainage system, but the volume of the cavity will remain fixed.

The model comes in two flavours: first, a two-dimensional, purely viscous flow model for the ice assumes that the cavity roof is in steady state and that water pressure in each separate cavity is spatially uniform. Where a cavity is in contact with the permeable part P of the bed, water pressure equals that in the ambient drainage system, while water pressure in isolated cavities is dictated by their volume. Second, a more general dynamic model assumes viscoelastic ice flow and explicitly considers how water is redistributed within the cavities by water pressure gradients, in a manner analogous to hydrofracture models for pre-existing cracks. The hydraulic conductivity that controls water flow is large within cavities (ensuring rapid equilibration) but vanishes when the ice–bed gap is zero, thereby allowing the model to capture the formation of isolated cavities and of isolated but uncavitated low-pressure regions in a dynamic framework.

The two versions of the model are susceptible to solution by different methods, making the simpler, purely viscous, steady-state version a useful test case for the more complicated dynamic version. To make the presentation more manageable, I have split these two model versions across two separate papers, focusing here on the purely viscous steady-state model. The dynamic model is presented in a companion paper (Schoof2023), which I will refer to below as Part 2. The present paper is structured as follows: first, I describe the mathematical model formulation in Sect. 2, with various technical aspects of the solution relegated to the appendices. In Sect. 3.1, I investigate how cavity extent depends on the effective pressure in the ambient drainage system, on the location of the permeable part P of the bed directly connected to the ambient drainage system, and on the past history of cavity formation across the bed. Subsequently, I use these solutions for cavity geometry in Sect. 3.2 to compute friction laws: that is, the corresponding amount of basal drag as a function of sliding velocity and effective pressure. I then investigate in Sect. 3.4 whether changes in bed geometry qualitatively affect the results. Implications for large-scale models of subglacial hydrology and glacier dynamics are discussed in Sect. 4.

2 A two-dimensional viscous steady-state model

Consider the possibility of isolated cavities in the two-dimensional, purely viscous, steady-state model of subglacial cavitation in Fowler (1986) and Schoof (2005). Based on the approximation of small bed slopes pioneered in Nye (1969) and Kamb (1970), the model can be written as follows: ice occupies the half-space y>0 in the Cartesian (x,y) plane. In that domain, ice flow satisfies Stokes' equations,

(1) η 2 u - p = 0 , u = 0 .

Here, u=(u,v) is the perturbation in ice velocity around a mean (ub,0) introduced by flow over bed topography, while p is the reduced pressure (that is, the actual pressure minus the cryostatic overburden), is the usual two-dimensional gradient operator, and η is the viscosity of ice, assumed to be constant here, and I assume that ub>0 so mean flow is to the right in Fig. 1.

Figure 1Definitions used in the model. The upstream and downstream cavity end points of the jth cavity are denoted by bj and cj, respectively. a is the width of the periodic domain. h(x) is the cavity roof height and b(x) the local bed elevation. I use beige colouring throughout the paper to indicate the permeable part P of the bed and grey for the impermeable part. The blue curve -σnn=p-2ηv/y shows that compressive normal stress against x. σnn must exceed the negative effective pressure Nj locally around any given cavity j as shown, but not globally, allowing low-pressure contact areas (with normal stress below N) to exist as shown. Cavity j=2 here is an isolated cavity, with fixed volume V2, while cavity j=1 is connected as it overlaps with P.


To be definite, I also assume the domain to be periodic in x with period a (Fig. 1). At the base of the ice y=0, let the set of points at which there is contact between ice and bed be denoted by C, and let the complement C denote cavities, or regions of ice–bed separation. Everywhere at the bed, xCC, we assume vanishing shear stress

(2) u y + v x = 0

For xC, the normal component of velocity vanishes, leading to the boundary condition

(3) v = u b b x ,

where b(x) is the elevation of the bed about a mean. Conversely, let C be composed of a set of disjoint intervals Cj=(bj,cj), each representing a separate cavity. On each Cj, normal stress is prescribed in the form

(4) p - 2 η v y = - N j ,

where Nj is the effective pressure in the jth cavity, defined as the difference between overburden and water pressure in the cavity. The cavity roof elevation hC satisfies the steady-state kinematic boundary condition

(5) v = u b h C x

on C with hC=b at cavity end points, so the lower boundary of the ice is continuous. These boundary conditions are combined with far-field conditions:

(6) p , u 0 as  y .

The previous work in Schoof (2005) assumed that the water pressure in each cavity is the same, implicitly requiring a highly permeable bed and allowing a universal effective pressure to be defined as N=Nj for all j. Taking the implied permeability of the bed further, Schoof (2005) added the inequality constraints

(7)p-2ηvy-Nfor xC(8)hC>bfor xC

in order to determine the extent of cavities. Physically, these inequalities represent the idea that normal stress cannot be less than the (assumed uniform) water pressure anywhere at the bed, since water will force its way between the ice and bed in that case, forming a new cavity, and that a cavity only exists if the cavity roof is indeed above the bed.

Here I abandon the assumption of a fully permeable bed. If parts of the bed are instead impermeable, there is no universally defined water pressure, and water will not force its way between the ice and bed simply because the normal stress drops locally to the water pressure in a distant drainage system. Water pressure is still assumed to be constant in each cavity while potentially differing between cavities, so the Nj values are constants but need not be equal to one another. As a result, the constraint (7) also need no longer hold across the bed.

To be more specific, I assume that only a part P of the bed is permeable and connected to a drainage system at prescribed effective pressure N. Hence the condition (Eq. 7) holds for xP, and any cavity straddling a part of P will be “connected” at the drainage pressure N. (P is a part of the bed, but specified here only in terms of the horizontal coordinates of points in P at the ice–bed interface, since no depth-dependent physics in the bed are resolved by the model.) Any cavity not straddling P will be “isolated” and required to hold a prescribed volume of water,

(9) V j = b j c j ( h C - b ) d x ,

which, if a solution exists, determines the effective pressure Nj. The constraint (8) still holds, but the inequality (7) is instead replaced by the weaker requirement that

(10) p - 2 η v y > - N j

in some finite intervals (bj-δ,bj) and (cj,cj+δ) (that is, there is some δ>0 such that the constraint 10 holds), ensuring that the cavity remains sealed.

Outside of the intervals (bj-δ,bj) and (cj,cj+δ), the inequality (10) can, and in general will, be violated somewhere as indicated in Fig. 1. The possibility of such underpressurized regions is the primary difference between the permeable and impermeable bed models. By not bounding compressive normal stress everywhere at the bed, however, the model does not allow a vapour-filled cavity to form if the normal stress drops to the triple-point pressure of water. In order to incorporate the latter effect, I would need to add the constraint that p-2ηv/y>-pi in C, where pi is ice overburden, and set p-2ηv/y=-pi in any cavity that does not straddle P and in which the prescribed water volume Vj (potentially equal to zero) would lead to an effective pressure Nj>pi if the volume constraint (9) were imposed. I omit this complication here on the basis that I expect overburden pi to be large compared with the typical normal stress variations caused by ice flow over bed undulations; suffice it to point out that the model described in Part 2 can in principle describe vapour-filled cavities.

Also note that Stubblefield et al. (2021) employ a similar but ultimately distinct volume constraint to Eq. (9): theirs is a global constraint, in which the bed is fully permeable (equivalent to P=(0,a) here) and all cavities are at the same effective pressure, but the latter is not prescribed. Instead, the total cavity volume is prescribed through initial conditions, the constraint itself being imposed on normal velocity so as to conserve that initial volume. Equation (9) here is a local constraint instead, prescribing the volume of an individual cavity.

The specification of a permeable bed portion P may be awkward but is realistically the only way to model partial access of the drainage system to the bed in two dimensions. Strictly speaking, water here is assumed to flow through the permeable bed in P in order to access connected cavities, but P can also be thought of as locations where an ambient drainage system is able to access the flowline being modelled laterally along the ice–bed interface, with the lateral dimension remaining unresolved. Below, I will typically consider either the entire bed permeable with P=(0,a), or I will consider a small permeable patch around a single location, which I will denote by xP. I will typically choose xP to be the location of a local minimum of compressive normal stress for an uncavitated bed, since that is where cavities first form for a permeable bed. In addition, in Sect. 3.3 I consider larger permeable bed portions P that do not align with these normal stress minima.

In any case, the modified steady cavity problem can be solved by a slight modification to the complex variable method in Schoof (2005), whose numerical method I also adapt. The technical detail is relegated to the Appendix. A steady-state solution to the model is likely to be highly non-unique, since the placement of prescribed water volumes Vj in isolated cavities is history-dependent and quite arbitrary in a steady-state model (by contrast, the dynamic model described in Part 2 self-consistently determines the volume of isolated cavities, precisely because it tracks the evolution in time of cavities).

In the next subsection, I consider a system of cavities that is in quasi-equilibrium, forced by a very slowly changing effective pressure N in the ambient drainage system. I also assume that the bed starts with no cavities. The latter initially form around the permeable parts P of the bed when N is made sufficiently small. The cavities at first remain trapped between prominent protrusions but can drown these bed protrusions abruptly when N is decreased to some critical values; I describe the method by which I compute the enlarged cavity in detail in Appendix A4. If N is increased again, the extended cavity roof can then make contact again with the drowned bed protrusion, thereby (in two dimensions) sealing the lee side of that protrusion and forming an isolated cavity. The volume of that isolated cavity is dictated by cavity size at the point where the cavity roof re-contacts, making the solution unique for a sequence of slow changes in N. Again, Appendix A4 provides further detail.

3 Results

3.1 Cavity geometry

Figure 2 shows the evolution of cavity geometry for the double-bumped periodic geometry,

(11) b ( x ) = h 0 sin 2 π x a + sin 4 π x a ,

with h0 and a constant. I focus first on the reference case of a fully permeable bed, as previously considered in Fowler (1986), Schoof (2005), Gagliardini et al. (2007), Helanow et al. (2020, 2021), Stubblefield et al. (2021), and de Diego et al. (2022, 2023).

Note that, when expressed as functions of a scaled position x*=2πx/a along the bed, cavity size and shape depend only on the following dimensionless effective pressure (Fowler1981, 1986):

(12) N * = N a 2 4 π 2 h 0 η u b ,

and I adopt this here to reduce the parameter space to be explored. Similarly, a dimensionless compressive normal stress defined by

(13) - σ n n * = a 2 4 π 2 h 0 η u b p - 2 η v x y = 0

also depends only on N* when expressed in terms of the scaled position x*. I use σnn* to visualize normal stresses at the bed. In the same vein, I use b*=b/h0 and hC*=hC/h0 as scaled bed and cavity roof elevations, and I use P*={x*:x*a/(2π)P} as the scaled version of the permeable bed.

Figure 2Cavity roof shape hC*(x*) and bed elevation b*(x*) for the bed given by Eq. (11) with P*=(0,2π) and (a1) N*=7.6, (b1) N*=4.02, (c1) N*=1.19, and (d1) N*=0.91. The corresponding compressive normal stresses -σnn* are plotted in panels (a2)(d2).


With the bed geometry given by Eq. (11), the basal compressive normal stress -σnn* has two equally deep minima around x*=1.64 and x*=4.65 prior to cavity formation. Two cavities per bed period form simultaneously around these locations in the lee of the two bed protrusions when effective pressure N* drops below a critical value Ninit*=8.06 (panel a1 and a2 of Fig. 2). The cavity roof hC* remains very close to the bed b* in the cavities initially, which are therefore easier to discern in the normal stress distribution -σnn* (panel a2). The pattern of normal stress shown here is common to the steady-state cavity solutions computed elsewhere (Fowler1986; Schoof2005; Gagliardini et al.2007; Stubblefield et al.2021; de Diego et al.2022, 2023): compressive normal stress is continuous at the upstream end of the cavity, with larger values immediately outside the cavity than inside acting to contain the water in the cavity, and normal stress has a positive singularity at the downstream cavity end. I show in Appendix A5 that this stress pattern necessarily follows from the inequalities (8) and (10).

The cavities expand continuously as N* is lowered further until they merge at a second critical value Ndisconnect*=1.19 (panels c1 and c2), and the merged cavity then continues to expand further. If N* is raised again, cavity evolution is completely reversible: for instance, the merged cavity once more separates in two at N*=Ndisconnect*.

Figure 3Panel (a) shows the effective pressure N* against cavity end-point positions bj* and cj* for a fully permeable bed of the form in Eq. (11). Ninit* and Ndisconnect* are defined in the main text, and “contact” and “cavity” mark the sides of the black curve occupied by contact areas and cavities. Panel (b) shows the corresponding upper surface elevation b*(x*) of the bed against x*.


The dependence of cavity size on N* can be visualized by plotting cavity end-point positions bj*=2πbj/a and cj*=2πcj/a against N* (Fig. 3, where the two critical values are marked with broken horizontal lines). Note that there is a unique solution for every N* here, corresponding to either a single merged cavity or two separate cavities. The labels “contact” and “cavity” indicate which side of the black curves corresponds to a contact area and a cavity, respectively. A second key feature of Fig. 3 is that the contact areas disappear at N*=0 and no solution exists for negative effective pressures N*<0: naturally, when water pressure is above overburden, force balance can no longer be maintained.

Figure 4Cavity roof shape hC*(x*) and bed elevation b*(x*) for the bed given by Eq. (11) with P*={1.64} and (a1) N*=4.01, (b1) N*=0.92, (c1) N*=0.91, (d1) N*=1.19, and (e1) N*=4.02. The permeable and impermeable portions of the bed are rendered in beige and grey, respectively. The corresponding normal stresses -σnn* are plotted in panels (a2)(e2); note that the isolated cavity in (e2) is at a different constant pressure than the connected cavity around the permeable bed portion P*.


The behaviour is somewhat different if I restrict the permeable portion P* of the bed to a small region around the upstream stress minimum at x*=xP*:=1.64 (Fig. 4). With only this small portion of the bed being permeable, a cavity starts to form at xP* when N*=Ninit*=8.04. As effective pressure in the drainage system is lowered, the cavity grows first on the lee side of the large bed protrusion to the left (to which the cavity is “attached”), while the lee of the smaller protrusion to the right remains uncavitated as shown in Fig. 4a1. (Note that I will use “large” or “prominent” protrusion to describe the protrusion that has the largest difference in height between the local maximum at its top and the local minimum on its upstream side.) This contrasts with the fully permeable bed case, where the lee sides of both bed protrusions become cavitated at the same N* (see also Fig. 3).

As before, the normal stress around the cavity is continuous at the detachment point at the upstream end of the cavity and singular at the reattachment point at the downstream end (Fig. 4a2). Normal stress exceeds -N* at both ends as required by the constraint (10). Note, however, that normal stresses on the lee side of the smaller protrusion are lower than -N*, and inequality (7) is violated there, away from the permeable bed portion P: an isolated underpressurized region forms here, separated from the cavity by the high-normal-stress region in the lee of the cavity.

As N* is decreased, the cavity expands, while the size of the high-stress region isolating the lee of the smaller bed protrusion shrinks. Eventually, the confinement of the cavity at its downstream end becomes marginal (Fig. 4b2) at N*=Nconnect*=0.92. A further reduction in N* causes the cavity to expand abruptly across the top of the smaller bed protrusion (Fig. 4c).

The newly expanded cavity roof now has a finite size gap above the smaller bed protrusion. It expands further, but now continuously, if effective pressure is lowered again. The expanded cavity is in fact identical in shape to the single merged cavity that forms for a fully permeable bed at the same effective pressure. More significantly, if N* is increased again from the critical value of Nconnect*, the cavity roof does not immediately re-contact the bed again. In order for the enlarged cavity roof to re-contact the smaller bed protrusion, N* has to increase by a finite amount to N*=Ndisconnect*=1.19>Nconnect* (Fig. 4d). That higher critical value is equal to the effective pressure at the merger of the two cavities that form independently in the lee of both bed protrusions when the entire bed is permeable, P*=(0,2π) (Fig. 3), and I use the same symbol Ndisconnect* deliberately.

In the present two-dimensional model, re-contact with a limited permeable bed portion immediately leads to the formation of a second isolated cavity downstream of the right-hand bed protrusion, which I treat as retaining a constant volume V2=1.062ah0/(2π) after reattachment (this being the volume at reattachment). A further increase in effective pressure N* in the permeable portion P* of the bed leads to the original cavity in the lee of the left-hand bed protrusion shrinking again, eventually disappearing at a critical value of Nshrink*=7.99, slightly less than the value of Ninit*=8.06 at which the cavity was originally formed. Meanwhile, the effective pressure N2* in the isolated cavity typically differs from N1*=N* in the connected cavity (Fig. 4e). Note that the solution is non-unique here: panels (a) and (e) of Fig. 4 correspond to (nearly) the same effective pressure N*.

Conversely, if N* is lowered again, the cavity that is attached to the larger bed protrusion on the left will reconnect to the isolated cavity that is attached to the smaller bed protrusion on the right at the same critical value Ndisconnect at which the isolated cavity originally formed: changes in cavity geometry become reversible once the lee sides of both bed protrusions have become cavitated.

Figure 5Panel (a) shows effective pressure N* against cavity end-point positions bj* and cj* for a bed of the form in Eq. (11) with P* concentrated around 1.64. Ninit*, Nshrink*, Ndisconnect*, and Nconnect* are defined in the main text, and “contact” and “cavity” mark contact areas and cavities on either side of the black curves (solutions obtained by starting with an uncavitated bed at N*=Ninit* and lowering N*). The red curves show solutions obtained when N* is lowered below Nconnect* and raised above Ndisconnect* subsequently. The newly formed isolated cavity is marked in red. Panel (b) shows effective pressure N* against cavity end-point positions bj* and cj* for a bed of the form in Eq. (11) with P* concentrated around xP*=4.65. Ndrown* is defined in the main text. Panel (c) shows the corresponding elevation b*(x*) of the upper surface of the bed against x*. The beige strips labelled Pa and Pb indicate the permeable bed portions used in panels (a) and (b), respectively.


The dependence of cavity end points on N* is again plotted systematically in Fig. 5a, which is analogous to Fig. 3. The black curves show cavity end-point positions that result if I start with an uncavitated bed and only lower N*, with the abrupt cavity enlargement at N*=Nconnect* clearly visible as a discontinuity in the downstream cavity end-point position c1. The upstream cavity end point in fact shifts discontinuously too, but by an amount that may be too small to discern. As in Fig. 3, the contact area again vanishes at N*=0 and no solution exists for negative N*: in fact, the solutions in Figs. 3 and 5a are identical for N*<Nconnect*. If, on the other hand, N* is first lowered below Nconnect* and then raised again, the cavity end-point solution follows the red curve above the disconnection value Ndisconnect*. Note that the isolated cavity that forms (indicated by red lettering) initially shifts slightly upstream as N* is increased above Ndisconnect*, but then remains relatively unaltered as the connected cavity shrinks and disappears.

Figure 6In red is the effective pressure N2*=-σnn* in the isolated cavity formed as in Fig. 4d against effective pressure N* in the connected cavity around the permeable bed portion P. The effective pressure N*=N1* in the connected cavity is plotted as a black dashed line, terminated at N*=Nshrink*, where the connected cavity disappears. The isolated cavity exists past Nshrink*, but not for N*<Ndisconnect*.


In addition, I have plotted the effective pressure N2* in the isolated cavity against the forcing effective pressure N* in Fig. 6. The effective pressure N2* mostly increases as N* does, implying a drop in water pressure in the isolated cavity as water pressure in the connected drainage system drops, albeit at a slower rate. This may be surprising given observations of anticorrelated water pressures between connected and unconnected parts of the bed (Murray and Clarke1995; Lefeuvre et al.2015; Rada and Schoof2018). There are two important differences here: first, the water pressure variations being considered are not transient, and consequently the size of both cavities has fully adjusted to steady-state conditions after a change in effective pressure N*. Second, in a flowline model, the redistribution of normal stress considered by Murray and Clarke (1995) and Lefeuvre et al. (2018) is modulated by flow over bed topography and by changes in the extent to which bed topography is drowned by cavities. For the bed geometry in Eq. (11) under consideration, an increase in N* in the connected cavity leads to more of the upstream face of the right-hand protrusion being covered by ice. The need to flow up and over that protrusion leads to a reduction in normal stress in its lee and hence to a drop in the water pressure required to maintain a cavity of fixed volume in its lee. This explains the increase in N2* with increases in N* here.

Figure 7Cavity roof shape hC*(x*) and bed elevation b*(x*) for the bed given by Eq. (11) with P*={4.65} and (a1) N*=7.60, (b1) N*=4.02, (c1) N*=0, and (d1) N*=-0.79. The permeable and impermeable portions of the bed are rendered in beige and grey, respectively. The corresponding normal stresses -σnn* are plotted in panels (a2)(d2); note that positive values of -σnn* as in panel (d2) correspond to negative effective pressure.


The ability of a cavity to expand across bed protrusions and subsequently create isolated cavities as described above depends on the position of the permeable portion of bed relative to prominent bed protrusions. Consider the same bed given by Eq. (11), but move the permeable portion of the bed to xP*=4.65. In that case, a cavity initiates here at the same initial value Ninit*=8.06. Now, however, the cavity is attached to a smaller bed protrusion and remains confined in its lee for all positive values of N, separated from the low-pressure region downstream of the more prominent protrusion by high normal stresses on either side of the cavity. This confinement in fact persists all the way to a negative effective pressure Ndrown*=-0.79 (Fig. 7). Beyond this critical effective pressure, the ice fully detaches from the bed, and vertical force balance is once more violated.

Note that the cavity is not able to expand upstream to the lee side of the bigger bed protrusion and only expands downstream past that bigger protrusion at the negative effective pressure Ndrown* when the confinement at the downstream end disappears: the normal stress upstream of the cavity remains in excess of water pressure even then. As with the previous example, I have plotted the position of cavity end points against N* in Fig. 5b. As in the fully permeable bed case in Fig. 3, there is now a unique solution, though it no longer disappears at N*=0. Note that the ability of cavities to remain contained at negative effective pressure due to uneven stresses induced by ice flow over topography may also be significant observationally, since sustained negative effective pressures are a frequent feature of borehole water pressure records (e.g. Rada and Schoof2018).

3.2 Basal drag

We can also ask how the formation of isolated cavities, and confinement of cavities, affects basal drag defined through (Fowler1986; Schoof2005)

(14) τ b = 1 a 0 a p - 2 η v x y = 0 h C x d x ,

where I treat hC=b in the contact areas C. As above, this can be cast in dimensionless form, now defining

(15) τ b * = τ b a 2 π h 0 N ,

which is then a function of N* only (Fowler1986). For consistency with Fowler (1986), Schoof (2005), and Gagliardini et al. (2007), I plot τb* against 1/N*=4π2h0ηub/(a2N) to visualize the resulting friction law, with 1/N* effectively being a proxy for ice velocity ub. Results for the double-humped bed given by Eq. (11) are shown in Fig. 8.

Figure 8Friction law for the bed given by Eq. (11): τb*=τba/(2πN) against 1/N*=4π2h0ηub/(aN) for the solutions shown in Figs. 3 and 5. The solid black curve (consisting of multiple segments) corresponds to the solution shown as a black curve in Fig. 5a (single cavity connected to a permeable bed P* around x*=1.64), while the red curve here also corresponds to the red curve in Fig. 5a (connected cavity around x*=1.64 and an isolated cavity around x*=4.65). The dashed blue curve corresponds to the solution in Fig. 5b (single connected cavity around x*=4.65); because the latter solution extends to negative N*, the friction law can likewise be extended to values of 1/N*<1/Ndrown*<0 as shown in the inset. The continuous dashed black curve (partly obscured by the solid black curve for 1/N*>1/Nconnect*) corresponds to the solution for a fully permeable bed in Fig. 3. The line labelled “Iken's bound” is at τb*=max(b*/x*).


The standard assumption of a fully permeable bed P*=(0,2π) gives rise to the single-valued, continuous black dashed curve (partly obscured by the solid black curve as indicated by the arrow marked “fully permeable bed”). It corresponds to relatively small values of τb* that satisfy Iken's bound τb*max(b*/x*) (Schoof2005): the maximum possible basal drag that can be attained is bounded by bed slope, where, with the bed shape given by Eq. (11), max(b*/x*)=3. The shape of the dashed black curve mirrors some of those in Schoof (2005).

With a small P* centred around xP*=1.64, the friction law changes significantly: the relationship of τb* and 1/N* come in multiple branches, depending on the presence of isolated cavities. When there is only a cavity in the lee of the prominent bed protrusion on the left, basal drag is quite high and can exceed Iken's bound (whose derivation in Schoof2005, is based on a permeable bed). Drag τb* drops abruptly when Nconnect* is reached and the cavity expands to drown out not only the second smaller bed protrusion, but also a significant part of the steeper slope behind it (solid black curve). Once the cavity has expanded and N* is increased again, an isolated cavity forms, leading to values of τb* that are generally comparable to those for a fully permeable bed (solid red curve). Below 1/Nshrink*0.125, τb* then simply becomes linear in 1/N*: this implies that τbub independent of N, as is familiar from theories of basal sliding in the absence of cavitation (Nye1969; Kamb1970). In the absence of expanding or shrinking connected cavities, an isolated cavity simply adopts a constant shape and changes its internal water pressure to keep that shape. The effect of such a constant-shape isolated cavity on steady-state basal drag is the same as for a rigid bed, since the shape of the base of the ice remains constant (although different from the uncavitated bed).

For the alternative case of P* centred around xP*=4.65 (dashed blue curve), the formation of a single confined cavity means there is only a single branch of the relationship between τb* and 1/N*. By contrast with the other cases considered above, τb* now increases without bound as 1/N* increases and, in fact, does so linearly in 1/N* for large 1/N*. The reason is simply that a finite cavity size is approached as N*0, and simultaneously a finite τb is approached so that τb/N must increase linearly in ub/N.

We, however, can also view 1/N* as the limit of a large velocity ub rather than the limit of a vanishing effective pressure. Once more, I find linear behaviour analogous to that in Nye (1969) and Kamb (1970) precisely because the confined cavity adopts a constant steady-state shape in the limit of large ub and therefore has the same effect as a rigid bed in the sense that the base of the ice retains its shape when ub changes, provided ub remains large. That shape differs from the case of an isolated cavity discussed above, which explains why the slope of the dashed blue curve for large 1/N* differs from the solid red curve at small 1/N*: even though the cavity shape becomes independent of 1/N* in both cases, those cavity shapes and locations differ from one another.

An oddity of the solution with xP*=4.65 is that it also exists with negative values 1/N*<1/Ndrown*-1.27 (see inset in Fig. 8); this is not to be interpreted as a valid solution for negative ub and positive N (which would give negative N*) but arises because although ub>0 is assumed throughout here, N* can be negative for xP*=4.65 (Fig. 5b).

As a further caveat, note that for a fixed N, unbounded τb as shown in Fig. 8 results from the ability to generate arbitrarily large compressive normal stresses on the upstream side of the smaller bed protrusion, balanced by correspondingly negative compressive normal stresses on the downstream side in the hydraulically isolated low-pressure region on the downstream side of the larger protrusion (Fig. 7c2, where -σnn* is scaled with 1/ub, so the actual stress is the pattern shown multiplied by a coefficient proportional to ub). As described in Sect. 2 immediately after Eq. (10), arbitrarily negative normal stresses cannot be generated since a vapour-filled cavity will eventually form, and this should ultimately lead to a bounded basal drag satisfying an amended version of Iken's bound, τbmax(b/x)pi, where pi is once more overburden. The model here ignores that possibility, effectively treating pi as infinite for the purposes of bounded basal drag.

3.3 More complicated permeable bed portions

The results above were computed either for completely permeable beds or for beds that had permeable sections located at normal stress minima prior to cavity formation. As pointed out, I view these permeable bed portions P as potential proxies for lateral access from a three-dimensional ambient drainage system along an unmodelled part of the ice–bed interface, to one side of the flowline that the model describes. In that case it may make sense for that lateral access to reach the modelled flowline in places where compressive normal stress has local minima. Locating the permeable bed where cavities form at the highest possible values of N is also convenient as it reduces the number of additional parameters that describe the bed in the absence of a more sophisticated three-dimensional model.

In order to investigate the effect of choosing different permeable bed portions P, I plot in Fig. 9 the dependence of cavity end-point positions on effective pressure N* for the same bed geometry (Eq. 11) as before, but for two alternative choices of P*: in panel (a), P* is the union of the intervals Pa=(5.890,6.283) and Pb=(2.316,2.749), while in panel (b), P*=Pb. In both cases, if I start with an uncavitated bed, a cavity first forms around the permeable patch Pb at a critical value N*=Ninit2*2.00; this the normal stress -σnn* at the upstream end x*=2.316 of the interval Pb when the bed is uncavitated.

Once formed, the cavity immediately has a finite size that extends beyond Pb and is identical to the single-cavity solution shown by the black curve in Fig. 5a. The black curve in Fig. 9a and b traces the growth of the cavity as N* is lowered below the initiation value Ninit2* and remains identical to the corresponding portion of the black curve in Fig. 5a, with the cavity expanding past the second smaller bed protrusion at the same critical value Nconnect* as in Fig. 5a.

Conversely, if N* is increased after initial cavity formation Ninit2*, the single cavity in the lee of the larger bed protrusion remains connected up to a much higher critical pressure Ndisconnect2*6.70: because the cavity that is established at Ninit2* extends far beyond the permeably patch, a significant increase in N* is needed to shrink it to the point at which the connection is lost to the permeable part of the bed patch. This is shown by the red curves in Fig. 9a and b, which remain identical to the black curve in Fig. 5a up to Ndisconnect2*, after which an isolated cavity forms as the downstream contact point moves past the upstream end of the interval Pb.

The only difference between the solutions in Fig. 9a and b arises if effective pressure is lowered below the critical value Nconnect*, at which the cavity extends past the lower bed protrusion, and subsequently raised to the critical value Ndisconnect* at which the downstream portion of the enlarged cavity becomes separated again by a contact area in Figs. 3a and 5a. Solutions for that situation are shown as blue curves in Fig. 9a and b. In Fig. 9a, the permeable portion Pa in the lee of that smaller bed protrusion keeps the downstream cavity connected up to a critical effective pressure Ndisconnect3*4.00. The blue curve in that case remains identical to the solution for a fully permeable bed shown in Fig. 3a up to that point; at Ndisconnect3*, the downstream end of the cavity attached to the smaller bed protrusion moves past the upstream end of Pa. By contrast, the absence of a downstream permeable interval immediately isolates the downstream cavity when Ndisconnect* is reached in Fig. 9b. The blue curve in Fig. 9b is therefore identical to the solution for an isolated downstream cavity shown in red in Fig. 5a. This is true until the upstream cavity also disconnects again at a critical pressure Ndisconnect4* very close to Ndisconnect2*.

These results serve as an illustration of how the placement of drainage access can serve to complicate the computation of cavity extent. Note, however, that in many cases the solutions shown in Fig. 9 correspond to appropriately spliced-together segments of the simpler solutions of Figs. 3 and 5, with each segment limited by the value of N* at which a cavity loses connection to P*. The key takeaway is probably that the location of the permeable patch P* may make irreversibility under changes in N* more pronounced: if P* is not centred around the location of the lowest normal stress for an uncavitated bed, then a fairly low effective pressure N* may be required to initiate cavity formation (given by Ninit2* in Fig. 9), but the cavity that is then formed can remain connected to the drainage system up to much larger effective pressures (Ndisconnect2* in Fig. 9).

Figure 9The effect of permeable bed patch location, using the same plotting scheme as Fig. 3. Panel (a) shows the cavity end-point location against N* for P*=PaPb, Pa=(5.890,6.283), and Pb=(2.316,2.749). Black shows cavity end points obtained when decreasing N* from the initiation value Ninit2*, and red shows cavity end points obtained when increasing N* after first establishing a cavity at Ninit2*. Blue (partly obscured by red) shows cavity end points obtained when increasing N* after first lowering effective pressure below Nconnect*. The black curve in Fig. 5a is shown as a dashed grey curve, and the black curve in Fig. 3 is shown as a dot-dashed grey curve. Panel (b) shows cavity end points for P*=Pb. The black curve and red curves are identical to panel (a), and the blue curve is analogous to panel (a), showing the cavity end-point position against N* if N* has previously been lowered below Nconnect*. The black and red curves in Fig. 5a are shown as dashed and dot-dashed grey curves, the latter almost completely obscured by the blue curve. Panel (c) shows the bed and location of Pa and Pb.


3.4 A more complicated bed shape

The results I have found above for the double-humped bed given by Eq. (11) translate qualitatively to other more complicated bed geometries. Below, I use the following triple-humped periodic bed profile as an illustration.

(16) b ( x ) = h 0 sin 2 π x a + 1 2 cos 4 π x a - sin 4 π x a + sin 8 π x a

Figure 10 shows N* against the location of cavity end points as in Figs. 3 and 5. I see similar behaviour as for the double-bumped bed: with spatially limited drainage access P, cavities can expand to drown bed protrusions in their lee, but not on their upstream side (panel b). In order to drown a lee-side bed protrusion at a positive effective pressure, the cavity in question needs to be attached to a larger bed protrusion than that being drowned (panel b). That drowning is also irreversible, leaving isolated cavities in place if N is increased again by a sufficient amount (red and blue solution curves in panel b). Where a cavity is attached to a small bed protrusion upstream of a larger one, it typically remains confined even at small negative effective pressures up to a critical value beyond which force balance is violated and no solution exists (panels c and d).

Figure 10Panel (a) shows effective pressure N* against cavity end-point positions for a fully permeable bed with shape given by Eq. (16) as solid black curves. Note that the solution is unique. Panel (b) shows the cavity end-point positions for the same bed with a small P*=Pb centred around xP*=3.23 (in the lee of the large bed protrusion; see panel e). Black shows the solution for a single cavity initiated around xP*. Red shows the solution with a single isolated cavity and blue with two isolated cavities. The dashed black line shows values of N* at which the single cavity expands abruptly, and the dashed red and blue curves show the formation of isolated cavities and the closing of the connected cavity in the presence of one or two isolated cavities. See the inset for details on cavity expansion and formation of an isolated cavity. Panel (c) shows the cavity end-point positions for the same bed with a small P*=Pc centred around xP*=5.25 (in the lee of the smallest bed protrusion as shown in panel e). The dashed line shows the negative value of N* at which the cavity no longer remains confined and the ice detaches from the bed. Panel (d) shows the cavity end-point positions for the same bed with a small P*=Pd centred around xP*=1.03 (the medium bed protrusion; see panel e). Panel (e) shows the corresponding bed surface elevation b*(x*) defined by Eq. (16) against x*. The beige strips show the permeable areas Pb, Pc, and Pd used in panels (b)(d), respectively.


The critical effective pressure at which a cavity extends abruptly across a smaller protrusion in its lee is marked by dotted black lines in Fig. 10b (this is equivalent to Nconnect in Fig. 5a, although there are two such critical values in Fig. 10b as there are two smaller bed protrusions in the lee of the largest protrusion). Once the critical effective pressure has been reached and the cavity has extended, contact with the cavity roof can only be re-established by increasing effective pressure to a somewhat different higher effective pressure shown as blue and red dotted lines in Fig. 10b (equivalent to Ndisconnect in Fig. 5a; see also the inset in panel b of Fig. 10). At the point of re-contact, an isolated cavity is created behind the lee-side bed protrusion. That isolated cavity will reconnect again if the effective pressure is lowered back to the same critical value at which it first formed. In other words, once a cavity has expanded past a lee-side bed protrusion, it will do so again more readily, facilitated by the presence of an isolated cavity behind that protrusion, instead of the original isolated region of ice–bed contact at low normal stress. The disconnection of isolated cavities is reversible, unlike the flooding of low-pressure contact areas.

Figure 11Friction law – the equivalent of Fig. 8 for the bed given by Eq. (16): τb* against 1/N* for the solutions shown in Fig. 10. Dashed black (partially obscured by solid black as indicated by arrows) corresponds to the permeable bed solution in Fig. 10a. Solid black (multiple segments), red, and blue correspond to the solutions shown in black, red, and blue, respectively, in Fig. 10b. The dashed red curve corresponds to the solution in Fig. 10c and dashed blue to the solution in Fig. 10d. The latter two do extend to some negative values of 1/N* (not shown).


The friction law for the triple-humped bed (Fig. 11) is more complicated than for the double-humped bed on account of the fact that different numbers of isolated cavities can form, but it again retains similar features: high levels of basal drag τb* are favoured when smaller lee-side bed protrusions have not been drowned yet or when cavities remain confined in the lee of small bed protrusions. For the latter case, basal drag is again unbounded as 1/N*. The abrupt expansion of a cavity corresponds to an abrupt drop in drag, as it does in Fig. 8. The lowest levels of basal drag are typically generated for permeable beds and for fully cavitated beds.

Figure 12Effective pressure in the isolated cavities shown as the blue solution in Fig. 10b against the corresponding effective pressure N*=N2* in the connected cavity (the equivalent of Fig. 6 for the triple-bumped bed given by Eq. 16). Blue shows the effective pressure N1* in the isolated cavity around x*=1.03, and red shows effective pressure N3* in the isolated cavity around x*=5.25, while the black dashed line shows N2*=N* for the range of values of N* for which the j=2 cavity is open. Note that N1* decreases slightly with forcing effective pressure N* once the second smaller isolated cavity around x*=5.25 has formed.


One behaviour that differs subtly between the two bed geometries considered here is the dependence of effective pressure in isolated cavities on the effective pressure in connected cavities. For the triple-humped bed (Fig. 12), I see that effective pressure in an isolated cavity directly downstream of the connected cavity increases with forcing effective pressure N* as in Fig. 6 (with the increase again being rapid when the cavity first forms and then much less than linear in N*). This corresponds to the upward slope of both the blue and red curves near their left-hand starting points, which mark the effective pressures at which the corresponding cavities first become isolated. However, once the larger isolated cavity becomes separated from the connected cavity upstream by an additional isolated cavity in the lee of the smallest bed protrusion, then effective pressure in that larger isolated cavity actually decreases with N*: the blue curve in Fig. 12, representing effective pressure in the isolated cavity in the lee of the second-tallest bed protrusion, actually slopes downwards very slightly underneath the red curve (that is, once there is an isolated cavity in the lee of the smallest bed protrusion, which separates the second-tallest from the tallest protrusion as shown in Fig. 10b).

4 Discussion

4.1 Steady-state subglacial hydrology

The steady-state solutions in Sect. 3 point to three primary insights: first, if the bed is forced by slow changes in drainage system effective pressure N and is therefore always in steady state except during brief transients, then connections to previously uncavitated parts of the bed are made at critical values of N/ub. These critical values depend on the geometry of the bed and on the locations of the parts of the bed that are permeable and therefore intrinsically connected to the ambient drainage system. The model denotes these parts by P, and they are indicated by beige colouring throughout the paper.

Second, when such connections occur, they invariably extend the existing cavity in the downstream direction and never upstream. This has major implications for the evolution of connectedness of the bed and for the effective pressures that can be sustained. For cavities that are caused by drainage system access P immediately in the lee of prominent bed bumps, downstream connections occur at positive effective pressures, and smaller bed bumps are submerged by expanding cavities first, as might be expected. If drainage system access P is located in the lee of less prominent bed bumps, then (perhaps counterintuitively) connections are made only once sufficiently negative effective pressures are reached and result in complete ice–bed detachment. Importantly, this implies that sustained negative effective pressures at the glacier bed are possible, as has been inferred from observations (Rada and Schoof2018).

Third, once a connection has been made and the lee of a smaller bed protrusion has become submerged, the cavity space on that lee side can subsequently become isolated due to an increase in effective pressure (or decrease in sliding velocity), which causes the cavity roof to be lowered. The critical value for the disconnection between the upstream cavity and newly isolated, cavity however, occurs at a higher critical value N/ub than the original connection (Fig. 5). Importantly, connection and disconnection become reversible at this point: once the downstream side of a smaller bed bump becomes cavitated, connection and disconnection happen at the same critical value of N/ub. A corollary of this third point is that it is easier to create connections once there are isolated cavities in place, in the sense of that connection happening at a higher value of N/ub than in the absence of those isolated cavities.

The reader may wonder at this point why one would bother with considering isolated, low-pressure contact areas at the bed at all: since their flooding is irreversible, are they irrelevant, since they will connect sooner or later and henceforth remain flooded, even if they become hydraulically isolated again? The point here is that treating the bed as fully impermeable outside of the region P is likely to be an idealization: in reality, there is almost certainly slow leakage through the “impermeable” bed portions as also envisaged in Hoffman et al. (2016) and Rada and Schoof (2018). If there are lengthy periods outside of the active drainage season (with the latter often occupying a relatively short part of the annual cycle) during which that leakage can drain isolated cavities, then it is possible that the bed starts each season in an uncavitated state. In that case, the expansion of cavities initially confined to locations with access to the drainage system occurs seasonally.

A second point that needs to be addressed here is the limitation imposed by using a two-dimensional domain. True hydraulic connections over distances longer than a single bed wavelength a are clearly only possible in two dimensions if the ice becomes fully detached from the bed, which is clearly not the object of the present study. In reality, hydraulic connections have to be made by connected cavity space that goes around rather than over prominent bed bumps in three dimensions. I anticipate that the results obtained here are still relevant to individual connections between cavities in three dimensions, with those cavities being extended laterally and connecting further downstream or upstream at a lateral offset. Studying these more complicated geometries requires a three-dimensional model (see also Helanow et al.2020, 2021) that can capture the dynamics of hydraulically isolated cavities and of isolated, uncavitated, low-pressure regions. The model presented in Part 2 is in principle capable of doing that, although in practice I have not been able to run it in a three-dimensional configuration due to computational constraints: three-dimensional cavity dynamics with hydraulic isolation remain an obvious area of future research.

4.2 Steady-state friction law

For a fully permeable bed, the ratio τb/N of basal drag to effective pressure is a single-valued function of the ratio of sliding velocity to effective pressure ub/N or, more generally, of ub/Nn for a power-law Glen's law rheology with exponent n (Fowler1986; Schoof2005; Gagliardini et al.2007; Helanow et al.2021). That function behaves roughly as a regularized Coulomb friction law, at least for highly irregular beds (Schoof2005; Helanow et al.2021). By contrast, partial permeability of the bed has a major effect on basal friction: basal drag τb/N now depends not only on ub/N, but also critically on where along the bed the region of drainage access P is located and on whether isolated cavities have previously been formed (Figs. 8 and 11).

The first qualitative difference between a permeable and impermeable bed is that Iken's bound τbNmax(b/x) need not hold for the latter: the derivation of that bound (Schoof2005) specifically relies on there being no compressive normal stresses at the bed below the water pressure in the ambient drainage system.

If drainage system access P is located in the lee of one of the smaller bed bumps, then the resulting cavities remain confined and do not lead to widespread ice–bed separation until effective pressure becomes negative as discussed above (see also Figs. 5 and 10). In that case, τb/N increases without bound in ub/N, with the relationship becoming linear at large ub/N so that τbub approximately (see the blue and red dashed curves in Figs. 8 and 11). This result is familiar from Nye–Kamb sliding theory (Nye1969; Kamb1970) for ice of constant viscosity (as is assumed here) sliding over a rigid bed in the absence of cavitation: the confined cavity modifies the shape of the lower boundary of sliding ice, but because cavities do not expand to cover the entire bed as ub/N (as would be the case for a fully permeable bed; see Fig. 3 or 10a), that modification approaches a finite limit for large ub/N, explaining why behaviour analogous to Nye–Kamb sliding is obtained. Importantly, the modification of the lower boundary of the ice depends on the precise location of the confined cavity, and the approximate constant of proportionality relating τb to ub depends on the location of P: this explains, for instance, why there are distinct dashed red and blue lines in Fig. 11.

The most dramatic changes in basal friction occur when P is immediately in the lee of the largest bed bump. In that case, τb/N will increase approximately linearly in ub/N until the cavity connects with the remainder of the bed (see the solid black curves in Figs. 8 and 11, with the discontinuity that corresponds to the connection point marked as 1/Nconnect* in Fig. 8). Iken's bound may be significantly exceeded during that initial increase in ub/N. Once the connection with the remainder of the bed occurs, basal drag τb/N drops dramatically by factors of approximately 3 and 10 in Figs. 8 and 11, respectively. This is not particularly surprising, as the extension of the cavity drowns out much of the previously uncavitated bed topography, forcing the ice to flow over fewer bed obstacles and thereby reducing form drag (that is, drag caused by flow over basal topography).

Once connection has occurred, the friction law mimics the friction law for a fully permeable bed. This remains the case even if ub/N decreases again to the point where isolated cavities form in the lee of some of the smaller cavities (compare the dashed black curve for a fully permeable bed with the solid red curve for a single isolated cavity in Fig. 8 and with the solid blue curve for isolated cavities in Fig. 11): the smaller obstacles remain drowned once these isolated cavities form, and form drag remains low.

Computation of steady-state friction τb (the dynamic case being even more complicated; see e.g. de Diego et al.2022 and also Gilbert et al.2022) therefore requires not only knowledge of ub and N, but also of the prior history of the bed and of hydraulic connections that have been made. This suggests that at least one additional state variable may need to be included in the formulation of steady-state basal friction laws, possibly the cavitation ratio θ of Thøgersen et al. (2019). The latter is defined as the fraction of the bed that has become cavitated. In fact, the results here suggest that changes in cavitation ratio may have a dominant effect on basal friction: a significant and abrupt increase in cavitation ratio occurs when a cavity extends or “connects” downstream (Figs. 5 and 10), and that increase in cavitation ratio corresponds to an equally abrupt, large drop in basal drag as discussed above.

In fact, a prototype parameterization for the friction laws shown in Figs. 8 and 11 is

(17) τ b = C ( θ ) u b 1 / n + C 0 N u b 1 / n Λ 0 N n + u b 1 / n ,

where the second term on the right-hand side is the regularized Coulomb friction law of Schoof (2005) and Gagliardini et al. (2007). Here n is the exponent in Glen's law (Paterson1994), while C0 and Λ0 are constants determined by the geometry of bed roughness. The first term on the right reflects the fact that the friction laws in Figs. 8 and 11 behave as effectively linear relationships when there is an isolated low-pressure region that has not become cavitated yet. In Eq. (17), I propose capturing this behaviour through a linear term in ub1/n (recall that n=1 in Figs. 8 and 11), with the slope of that linear term depending on the cavitation ratio θ. In order to emulate the behaviour of a fully permeable bed, C(θ) needs to approach zero for cavitation ratios sufficiently close to unity.

A friction law of this form in turn implies that subglacial drainage models may need to incorporate a description of the evolution of cavitation ratio. As I will show in Part 2, cavitation ratio and mean cavity depth (the variable commonly used to define cavity geometry in large-scale drainage models) are not simple proxies for each other, implying that the introduction of cavity ratio into friction laws and drainage parameterizations would indeed imply an increase in model complexity.

There is a second complication in the definition of a friction law that deserves to be stressed for an impermeable bed: the quantity that is commonly understood as “effective pressure”, overburden minus water pressure at the bed, is not uniquely defined but potentially varies from cavity to cavity. That is, effective pressure varies over length scales that are treated as microscopic in typical subglacial drainage models because water pressure differs between cavities. In the idealized model I use here, I define a unique “ambient drainage system effective pressure” N in the permeable bed portions P and am able to express a friction law in terms of N and ub (albeit in the form of a multi-valued friction law) as is done in Figs. 8 and 11.

The effective pressure in the connected portion of the drainage system is likely to be the only useful effective pressure that can be defined, as it will in general vary smoothly in space and can therefore be modelled at the large scale, at least in principle. That observation does underline, however, the need to include additional degrees of freedom that capture the degree of cavitation in friction laws, since effective pressure is then meaningless in a part of the bed that is fully hydraulically isolated, with no drainage system access at all: there may still be isolated cavities in that case, and their presence will affect basal friction as discussed above. To compound matters, this situation also significantly complicates any attempts to constrain such a friction law observationally: while effective pressure in a connected drainage system can in principle be measured by borehole access to the bed, the presence and extent of isolated cavities at the bed are much harder to determine.

5 Conclusions

Using a simple extension of an existing purely viscous model for steady-state basal cavities in two dimensions, I have shown that uncavitated regions of the bed can persist indefinitely at low normal stress provided there is no drainage pathway along which water can reach them. Such drainage pathways are created under slow changes in forcing effective pressure N when that effective pressure reaches a critical value. The creation of such connections is not reversible by simply raising N back above its critical value, but requires a greater increase in N and leaves behind an isolated cavity. The formation of connections also leads to a significant drop in basal friction that is likewise irreversible, since the isolated cavity that is left behind by a subsequent increase in N significantly reduces contact between the ice and bed even when the hydraulic connection is closed again. To the best of my knowledge, few if any of these phenomena are included in current large-scale subglacial drainage models or basal friction laws.

The main limitations to the work presented here derive from its assumption of quasi-steady conditions and its restriction to two dimensions. Dynamic cavity connections have significantly richer behaviour than the quasi-steady solution in the present paper suggests and are investigated in detail in a companion paper. Three-dimensional bed topography by contrast remains an open problem and holds the key to a more complete understanding of hydraulic connectivity. Connections at the bed are presumably more likely to occur when bed topography is three-dimensional: in a two-dimensional setting, connectivity along the entire model domain is only possible when ice–bed contact is lost completely, whereas this is not the case in three dimensions. Similarly, contact of the ice roof between two cavities in three dimensions does not necessarily make them disconnected, whereas it does in two dimensions.

Appendix A: Complex variable solution of the viscous steady-state problem

A1 Complex variable formulation

The construction in Fowler (1986) and Schoof (2002, pp. 51–54) allows the problem consisting of Eqs. (1), (3), (4), and (6) to be written in the following form: let z=x+iy, and find an analytic function Ω(z) in the complex plane cut along the real axis, satisfying

(A1a)Ω(z)=Ω(z),(A1b)-2iΩ+(x)-Ω-(x)=-Njfor xCj,(A1c)Ω+(x)+Ω-(x)=ηubb′′(x)for xC,(A1d)Ω(z)0as I(z)±,

where a prime indicates differentiation (in this case, with respect to x), an overbar signifies complex conjugation, and superscripts + and denote limits taken from above and below the real axis. The constraints (8) and (10) become

(A1e)hC(x)>b(x)for xC, whereΩ+(x)+Ω-(x)=ηubhC′′(x)(A1f)and hC(bj)=b(bj),(A1g)hC(cj)=b(cj);(A1h)-2iΩ+(x)-Ω-(x)>Njfor bj-δ<x<bj(A1i) and cj<x<cj+δ

and some finite δ.

Let ζ=expi2πz/a and ξ=expi2πx/a. The assumed periodicity of the solution ensures that Ω(z) can be mapped one to one to G(ζ)=Ω(z), and similarly b2(ξ)=b′′(x) and h2C(ξ)=hC′′(x) are one-to-one mappings. The functions G, b2, and h2C satisfy

(A2a)G(ζ)=G(1/ζ),(A2b)G()=0,(A2c)-2iG+(ξ)-G-(ξ)=-Njfor ξΓj,(A2d)G+(ξ)+G-(ξ)=ηubb2(ξ)for ξΓ,(A2e)G+(ξ)+G-(ξ)=ηubh2C(ξ)for ξΓ,

where Γj, Γ, and Γ are Cj, C, and C mapped into the complex ζ plane (where they become subsets of the unit circle), and + and now indicate limits taken from within and without the unit circle in the ζ plane.

The solution method followed here is that of Schoof (2002, 2005), slightly modified to account for cavities at different effective pressures. I outline the procedure in full below, adding detail omitted in the original account by Schoof (2002, 2005).

A2 Cavity roof re-contact constraints

As in Schoof (2005), it is possible to conclude that the cavity roof must disconnect and reattach tangentially and that it suffices to impose this on n−1 of n cavities since any valid solution to Eq. (A2) ensures that re-contact is then also tangential for the nth cavity. Consider the integral

(A3) I = 0 a Ω + ( x ) + Ω - ( x ) d x = η u b j = 1 n h C ( c j ) - b ( c j ) - j = 1 n h C ( b j ) - b ( b j ) ,

where I have used Eqs. (A1c) and (A1g). Enforcing the contact condition (A1e) combined with the constraint that hC(bj)=b(bj), hC(cj)=b(cj) implies that hC(cj)b(cj), hC(bj)b(bj), and hence I≤0. On the other hand, transforming to the ζ plane,

(A4) I = a 2 π i Γ Γ G + ( ξ ) + G - ( ξ ) ξ d ξ = 0

on account of Cauchy's theorem, since ΓΓ is the unit circle and therefore a closed contour, and G(0)=G()=0. I=0 in turn implies that the cavity roof detaches and re-contacts tangentially, so

(A5) h C ( c j ) = b ( c j ) , h C ( b j ) = b ( b j )

for j=1,,n.

In fact, tangential cavity roof detachment and re-contact are required not only by Eq. (A4), but also by the original construction of the model in Eq. (A1), which requires differentiation of the original normal velocity condition v=ubb or v=ubhC (Schoof2002, p. 44); recovery of the original boundary condition in terms of antiderivatives of Ω confirms that no discontinuity between hC and b can appear if Ω is sectionally holomorphic in the sense of Muskhelishvili (1992) (meaning it gives rise to an integrable stress field).

The point here is really to account for the independent number of constraints on the solution that arise from the tangential re-contact. In integrating Eq. (A1g) (or Eq. A2e), the relevant continuity constraints can always be imposed on one cavity end point (say, the upstream end), and integration forward to the other cavity end point then creates a constraint on the solution. Thus, integrating once, I obtain n equations of the form

(A6) b ( c j ) = b ( b j ) + b j c j Ω + ( x ) + Ω - ( x ) d x ,

where Ω±(x)=G±(exp(i2πx/a)). Integrating twice, I obtain another n constraints

(A7) b ( c j ) = b ( b j ) + b ( b j ) ( c j - b j ) + b j c j ( c j - x ) Ω + ( x ) + Ω - ( x ) d x .

Note, however, that one of the n constraints (A6) is redundant for a valid solution G satisfying G(0)=G()=0, since this ensures that I=0 and the remaining equation of the form in Eq. (A6) is automatically satisfied.

A3 Solution

Armed with this result, I can again follow the same solution procedure as in Schoof (2002). G can be written in the form (Muskhelishvili1992)


where P is a polynomial and χ is a Plemelj function, holomorphic in the complex plane cut along Γ, on which it satisfies χ+(ξ)+χ-(ξ)=0. There are multiple choices of χ that give rise to a sectionally holomorphic solution G, differing in the number and location of singularities at the cavity end points x=bj and x=cj. As in Fowler (1986) and Schoof (2005), I default to the choice

(A9) χ ( ζ ) = j = 1 n ζ - ξ b j ζ - ξ c j 1 / 2 ,

behaving as χ→1 as ζ→∞, with ξbj=exp(i2πbj/a) and ξcj=exp(i2πcj/a). This choice of χ generally places a stress singularity at cavity re-contact points x=cj but ensures that stress is continuous at detachment points x=bj. That choice is not arbitrary: in Sect. A5 I confirm that stress at x=bj must be continuous in order to simultaneously satisfy Eqs. (A1e) and (A1i) and that, in general, the stress field at x=cj will be singular when the same constraints are satisfied locally near the re-contact point.

In order for G to satisfy G(∞)=0 with χ given by Eq. (A9), P≡0 is necessary and sufficient. The remaining constraint on G is that G(ζ)=G(1/ζ); when the latter is satisfied, G(0)=G()=0 follows automatically. Again as in Schoof (2005, p. 618), it is possible to show that

(A10) χ + ( ξ ) = - χ + ( ξ ) / χ ( 0 ) on  Γ , χ + ( ξ ) / χ ( 0 ) on  Γ , , ξ = 1 ξ , d ξ = - 1 ξ 2 d ξ .

Using these, it follows that

(A11) G ( 1 / ζ ) = G ( ζ ) - 1 2 π i Γ η u b b 2 ( ξ ) χ + ( ξ ) ξ d ξ + j = 1 n Γ j - i N j / 2 χ + ( ξ ) ξ d ξ χ ( ζ ) ,

and the required constraint is to set the term in square brackets to zero,

(A12) J := Γ η u b b 2 ( ξ ) χ + ( ξ ) ξ d ξ + j = 1 n Γ j - i N j / 2 χ + ( ξ ) ξ d ξ = 0 .

Suppose that the Nj is prescribed. With P≡0, the solution G in Eq. (A8) contains 2n unknown parameters in the form of the cavity end-point locations ξbj and ξcj. Assuming that G(ζ)=G(1/ζ) so G++G- is real, I have 2n−1 real constraints through Eqs. (A6) and (A7). This leaves a single real constraint to close the system, and it therefore remains to show that Eq. (A12) constitutes that single real equation. Taking the complex conjugate of the left-hand side of Eq. (A12) and using Eq. (A10), it is possible to show that J=χ(0)J. Since χ(0)=exp[iπj=1n(bj-cj-1)/a] (Schoof2002, p. 98) and 0<j=1n(bj-cj-1)/a<1, it follows that the real and imaginary parts of χ(0) are non-zero, and hence ℜ(J)=0 implies ℑ(J)=0 and vice versa. Equation (A12) therefore constitutes a single real constraint, and together with Eqs. (A6) and (A7) I have 2n real constraints to determine the 2n cavity end points. Prescribing cavity volume Vj rather than effective pressure Nj does not lead to further complications since putting Vj=bjCjhc(x)dx simply adds the required additional constraint to determine the corresponding Nj. The implementation of Eqs. (A12), (A6), and (A7) (combined with additional constraints on Nj when cavity volume is prescribed) follows the same numerical method as in Schoof (2002).

A4 Arc length continuation

In practice, I introduce the smallest new cavity possible when the inequality (7) is violated somewhere on P (note that this is generally simple to do when P is a small region around the location xP where normal stress has a local minimum in the absence of cavitation). I then use an arc length continuation to solve the system of equations in Eqs. (A6), (A7), and (A12) while decreasing the effective pressure N, forcing cavity end points to change continuously where they can.

Neighbouring cavities j and j+1 can merge when cj=bj+1 for some critical value of N, in which case I simply deleted cj and bj+1 and created a single enlarged cavity with end points bj and cj+1. Abrupt enlargement of cavities into a previously uncavitated low-pressure region occurs when the solution computed by arc length continuation violates the local constraints in Eqs. (8) and (10) near a cavity end point. This generally corresponds to a fold bifurcation along the solution curve (plotting cavity end-point locations against N). N begins to increase again along the curve at such a fold, signalling that the actual solution under a further decrease in N is not continuous. I use arc length continuation to extend the solution further until I reach another solution with the value of N at the fold bifurcation, but for which the inequalities (8) and (10) are satisfied. I treat that solution as representing the enlarged cavity that results from decreasing N past the fold bifurcation and discard solutions computed by arc length continuation that do violate the inequalities in Eqs. (8) and (10).

In order to capture the effect of cavity isolation, I compute solutions by arc length continuation under increases in N, checking whether inequalities in Eqs. (8) and (10) are satisfied. An isolated cavity forms when the cavity roof contact constraint (8) is violated in the interior of a cavity. In that case, I introduce new contact points where re-contact occurs and check whether either of the two new cavities created in the process no longer straddles P. Such a cavity is then isolated. I compute its volume and restart a computation by arc length continuation while imposing Eq. (9) for that cavity.

A5 Cavity end-point singularities revisited

Here I show that continuous stress at cavity detachment points and a stress singularity at reattachment points are natural consequences of the inequality constraints (A1e) and (A1i). I use the complex variable formulation deployed above, but note that the same result could be obtained by looking for a stream function solution of the Stokes flow problem (Eq. 1) in terms of local polar coordinates centred at x=bj or x=cj (see e.g. Fontelos and Muñoz2007).

Consider the original Hilbert problem (Eq. A1) locally, in a neighbourhood of a cavity end point z=bj or z=cj. Consider first the detachment point bj, and let

(A13) Ω ( z ) = - i N j / 4 + η u b b ′′ ( b j ) + F ( z ) for  I ( z ) > 0 i N j / 4 + η u b b ′′ ( b j ) + F ( z ) for  I ( z ) < 0 .

Then, in some sufficiently small open disc D around z=bj, F(z)=F(z) is holomorphic with a branch cut L along the intersection of D with the half-line L0 given by y=0, x<bj. On that branch cut

(A14) F + ( x ) + F - ( x ) = u b [ b ′′ ( x ) - b ′′ ( b j ) ] .

Assuming again that F is sectionally holomorphic in the disc to ensure integrable stresses, the solutions in the disc take the form (Muskhelishvili1992)

(A15) F ( z ) = Φ ( z ) + 1 2 π i L u b b ′′ ( x ) - b ′′ ( b j ) χ + ( x ) ( x - z ) d x χ ( z ) .

Here, χ(z)=(z-bj)-1/2 is analytic in the plane cut along L0, behaving as χ(x)=1/x-bj for x>bj along the real axis, and Φ is holomorphic in D. Assume b′′ is continuously differentiable. Then the limiting values as y→0 of the integral in the curly brackets behave as a constant plus a term of O(|z-bj|3/2log|z-bj|) (by a straightforward adaptation of the derivation in Muskhelishvili1992, pp. 45–49), while the analytic function Φ can be expanded as a Taylor series around z=bj as Φ(z)=a0+a1(z-bj)+O(|z-bj|2). To ensure that F(z)=F(z), a0 and a1 must be real. To a quadratic error in (zbj), I simply have F(z)a0+a1(z-bj)(z-bj)1/2, and I can evaluate

ηubhC′′=Ω+(x)+Ω-(x)ηubb′′(bj)+2[a0+a1(x-bj)]/x-bj(A16)for x>bj,p-2ηvx=-2iΩ+(x)-Ω-(x)-Nj-4a0+a1(x-bj)/bj-x-Nj(A17)for x<bj.

Since hC=b from Sect. A2, it follows that I must have hC′′(x)b′′(bj) for x>bj in order to ensure that hC>b inside the cavity. It follows that a0≥0, with a1≥0 if a0=0. In order for the normal stress constraint (A17) to be satisfied, ensuring the cavity remains sealed (by considering a local solution, I can dispense with the machinery of requiring a constraint only over a finite region of size δ as in Eq. A1i), I see that I must have a0≤0 and a1≥0 if a0=0. The only way that both of these constraints can be satisfied is that a0=0 and a1≥0. This immediately ensures that normal stress p-2ηv/x-Nj+a1bj-x is non-singular at the detachment point.

The same approach can be used near a re-contact point, but with different conclusions. Replacing bj by cj, I can still define F inside a small open disc D centred on z=cj through Eq. (A13). F still satisfies Eq. (A14), but now on the intersection L0 of D with the half-line L={(x,0):x>cj}. Similarly, χ(z)=(z-cj)-1/2 is holomorphic in the plane cut along L0, behaving as 1/x-cj when the branch cut is approached from the upper half-plane. F(z)=F(z) now requires that I write Φ(z)=ia0+a1(z-cj)+O(|z-cj|2) with a0 and a1 real. The equivalent of Eq. (A16) and inequality (A17) becomes

ηubhC′′=Ω+(x)+Ω-(x)ηubb′′(cj)+2[a0+a1(x-cj)]/cj-x(A18)for x<cj,p-2ηvx=-2iΩ+(x)-Ω-(x)-Nj+4a0+a1(x-cj)/x-cj-Nj(A19)for x>cj.

With hC=b at x=cj, I must still have hC′′b′′ for x<cj to ensure that hC>b, and hence a0≥0 with a1≤0 if a0=0 from Eq. (A18). Satisfying the normal stress condition (A19) requires that a0≥0, with a1≥0 if a0=0. In general I therefore expect a solution with a0>0 and a singular normal stress of the form p-2ηv/x4a0/x-cj.

Code availability

The code used is available from the author on request.

Data availability

No data sets were used in this article.

Competing interests

The author has declared that there are no competing interests.


Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.


This work was supported by NSERC Discovery Grant RGPIN-2018-04665. I would like to acknowledge the helpful comments of two anonymous referees, which improved the paper significantly.

Financial support

This research has been supported by the Natural Sciences and Engineering Research Council of Canada (grant no. RGPIN-2018-04665).

Review statement

This paper was edited by Nanna Bjørnholt Karlsson and reviewed by two anonymous referees.


Andrews, L., Catania, G., Hoffman, M., Gulley, J., Lüthi, M., Ryser, C., Hawley, R., and Neumann, T.: Direct observations of evolving subglacial drainage beneath the Greenland Ice Sheet, Nature, 514, 80–83, 2014. a

de Diego, G. G., Farrell, P. E., and Hewitt, I. J.: Numerical approximation of viscous contact problems applied to glacial sliding, J. Fluid Mech., 938, 21,, 2022. a, b, c, d, e

de Diego, G. G., Farrell, P. E., and Hewitt, I. J.: On the Finite Element Approximation of a Semicoercive Stokes Variational Inequality Arising in Glaciology, SIAM J. Numer. Anal., 61, 1–25,, 2023. a, b, c, d

Flowers, G. E.: Modelling water flow under glaciers and ice sheets, P. Roy. Soc. A, 471, 20140907,, 2015. a

Fontelos, M. and Muñoz, A.: A free boundary problem in glaciology: the motion of grounding lines, Interface. Free Bound., 9, 67–93, 2007. a

Fountain, A. and Walder, J.: Water flow through temperate glaciers, Rev. Geophys., 36, 299–328, 1998. a

Fowler, A.: A theoretical treatment of the sliding of glaciers in the absence of cavitation, Philos. T. Roy. Soc. Lond., 298, 637–685, 1981. a

Fowler, A.: A sliding law for glaciers of constant viscosity in the presence of subglacial cavitation, P. Roy. Soc. Lond. A, 407, 147–170, 1986. a, b, c, d, e, f, g, h, i, j

Fowler, A.: Sliding with cavity formation, J. Glaciol., 33, 255–267, 1987. a

Gagliardini, O., Cohen, D., Raback, P., and Zwinger, T.: Finite-element modeling of subglacial cavities and related friction law, J. Geophys. Res., 112, F02027,, 2007. a, b, c, d, e, f, g

Gilbert, A., Gimbert, F., Thøgersen, K., Schuler, T. V., and Kääb, A.: A Consistent Framework for Coupling Basal Friction With Subglacial Hydrology on Hard-Bedded Glaciers, Geophys. Res. Lett., 49, e2021GL097507,, 2022. a

Helanow, C., Iverson, N. R., Zoet, L. K., and Gagliardini, O.: Sliding Relations for Glacier Slip With Cavities Over Three-Dimensional Beds, Geophys. Res. Lett., 47, e2019GL084924,, 2020. a, b, c, d

Helanow, C., Iverson, N. R., Woodard, J. B., and Zoet, L. K.: A slip law for hard-bedded glaciers derived from observed bed topography, Sci. Adv., 7, eabe7798,, 2021. a, b, c, d, e, f

Hewitt, I.: Modelling distributed and channelized subglacial drainage: the spacing of channels and eskers, J. Glaciol., 57, 302–314,, 2010. a

Hewitt, I.: Seasonal changes in ice sheet motion due to melt water lubrication, Earth Planet. Sc. Lett., 371, 16–25, 2013. a

Hodge, S. M.: Direct Measurement of Basal Water Pressures: Progress and Problemss, J. Glaciol., 23, 309–319,, 1979. a

Hoffman, M., Andrews, L., Price, S., Catania, G., Neumann, T., Lüthi, M., Gulley, J., Ryser, C., Hawley, R., and Morris, B.: Greenland subglacial drainage evolutoin regulated by weakly connected regions of the bed, Nat. Comms., 7, 13903,, 2016. a, b

Kamb, B.: Sliding motion of glaciers: Theory and observation, Rev. Geophys., 8, 673–728, 1970. a, b, c, d

Kamb, B.: Glacier Surge Mechanism Based on Linked Cavity Configuration of the Basal Water Conduit System, J. Geophys. Res., 92, 9083–9100, 1987. a

Lefeuvre, P.-M., Jackson, M., Lappegard, G., and Hagen, J. O.: Interannual variability of glacier basal pressure from a 20 year record, Ann. Glaciol., 56, 33–44,, 2015. a, b

Lefeuvre, P.-M., Zwinger, T., Jackson, M., Gagliardini, O., Lappegard, G., and Hagen, J. O.: Stress Redistribution Explains Anti-correlated Subglacial Pressure Variations, Front. Earth Sci., 5, ISSN 2296-6463,, 2018. a

Lliboutry, L.: General Theory of Subglacial Cavitation and Sliding of Temperate Glaciers, J. Glaciol., 7, 21–58, 1968. a

Murray, T. and Clarke, G.: Black-box modeling of the subglacial water system, J. Geophys. Res., 100, 10231–10245, 1995. a, b, c

Muskhelishvili, N.: Singular Integral Equations, Dover Publications, Inc., New York, unabridged republication of 2nd Edn., P. Noordhoff, Groningen, 1953, ISBN 10 0486668932, 1992. a, b, c, d

Nye, J.: A calculation of the sliding of ice over a wavy surface using a Newtonian viscous approximation, P. Roy. Soc. Lond. A, 311, 445–467, 1969. a, b, c, d

Paterson, W.: The Physics of Glaciers, Pergamon, Oxford, 3rd Edn., ISBN 0-08037945 1, 1994. a

Rada, C. and Schoof, C.: Channelized, distributed, and disconnected: subglacial drainage under a valley glacier in the Yukon, The Cryosphere, 12, 2609–2636,, 2018. a, b, c, d, e, f, g

Schoof, C.: Mathematical Models of Glacier Sliding and Drumlin Formation, PhD thesis, Oxford University, 2002. a, b, c, d, e, f, g

Schoof, C.: The effect of cavitation on glacier sliding, P. Roy. Soc. Lond. A, 461, 609–627,, 2005. a, b, c, d, e, f, g, h, i, j, k, l, m, n, o, p, q, r, s, t, u, v

Schoof, C.: The evolution of isolated cavities and hydraulic connection at the glacier bed – Part 2: A dynamic viscoelastic model, The Cryosphere, 17, 4817–4836,, 2023. a

Schoof, C., Hewitt, I., and Werder, M.: Flotation and free surface flow in a model for subglacial drainage. Part 1. Distributed drainage, J. Fluid Mech., 702, 126–156, 2012. a

Sommers, A., Rajaram, H., and Morlighem, M.: SHAKTI: Subglacial Hydrology and Kinetic, Transient Interactions v1.0, Geosci. Model Dev., 11, 2955–2974,, 2018.  a

Stubblefield, A. G., Spiegelman, M., and Creyts, T. T.: Variational formulation of marine ice-sheet and subglacial-lake grounding-line dynamics, J. Fluid Mech., 919, A13,, 2021. a, b, c, d, e

Thøgersen, K., Gilbert, A., and Schuler, T.: Rate-and-state friction explains glacier surge propagation, Nat. Comms., 10, 2023,, 2019. a

Werder, M., Hewitt, I., Schoof, C., and Flowers, G.: Modeling channelized and distributed subglacial drainage in two dimensions, J. Geophys. Res., F118, 2140–215,, 2013. a, b

Short summary
Computational models that seek to predict the future behaviour of ice sheets and glaciers typically rely on being able to compute the rate at which a glacier slides over its bed. In this paper, I show that the degree to which the glacier bed is hydraulically connected (how easily water can flow along the glacier bed) plays a central role in determining how fast ice can slide.