Subglacial hydrology and the formation of ice streams

Antarctic ice streams are associated with pressurized subglacial meltwater but the role this water plays in the dynamics of the streams is not known. To address this, we present a model of subglacial water flow below ice sheets, and particularly below ice streams. The base-level flow is fed by subglacial melting and is presumed to take the form of a rough-bedded film, in which the ice is supported by larger clasts, but there is a millimetric water film which submerges the smaller particles. A model for the film is given by two coupled partial differential equations, representing mass conservation of water and ice closure. We assume that there is no sediment transport and solve for water film depth and effective pressure. This is coupled to a vertically integrated, higher order model for ice-sheet dynamics. If there is a sufficiently small amount of meltwater produced (e.g. if ice flux is low), the distributed film and ice sheet are stable, whereas for larger amounts of melt the ice–water system can become unstable, and ice streams form spontaneously as a consequence. We show that this can be explained in terms of a multi-valued sliding law, which arises from a simplified, one-dimensional analysis of the coupled model.


Introduction
Variations in sliding velocity at the ice-bed interface can cause extreme spatial and temporal differentials in flow speed throughout an ice sheet. In particular, while there are many locations in Antarctica where the ice has zero velocity, velocities of up to thousands of metres per year are observed in some coastal areas [1]. In the fast-flowing regions, the rapid basal sliding is attributed to the presence of meltwater at the bed [2]. Frictional resistance at the bed is intrinsically linked to the water pressure in the meltwater drainage system [3][4][5][6]. This effect is best considered in terms of the effective pressure (the difference between ice overburden pressure and water pressure); sliding at the bed can occur when the effective pressure is small and the weight of the ice is largely supported by the underlying water. The dependence of the basal sliding velocity on both the basal shear stress and the effective pressure at the bed is typically modelled with a sliding law [7][8][9]. To improve our understanding of variations in basal velocity, we must investigate the physical basis and consequences of the sliding law.
It is the behaviour of the subglacial hydrological system that governs the effective pressure at the ice-bed interface. There are two contrasting types of subglacial drainage systems beneath ice sheets that are commonly discussed in the literature. If there is little meltwater production, the bed can be effectively drained by a distributed system of cavities [10,11] or through thin, patchy films [12]. With this style of drainage, increasing meltwater production increases the water pressure beneath the ice, which lowers the effective pressure, and hence results in more rapid sliding [10,[13][14][15]. By contrast, subglacial water transport may occur through larger scale Röthlisberger channels. These provide efficient drainage for meltwater, resulting in low water pressure at the bed, corresponding to high effective pressure, and hence slower sliding velocities [16]. Beneath Antarctica, however, there is more evidence of high water pressure, distributed drainage systems [17].
In this work, we are particularly interested in ice streams and their subglacial hydrology. A significant proportion of all ice discharge occurs through ice streams, despite the fact that they are only found over a small portion of an ice sheet. In Antarctica, for example, it is estimated that ice streams cover approximately 10% of the ice-sheet's surface, but transmit up to 90% of drainage [1,18]. There is evidence suggesting the presence of an active subglacial water system in West Antarctica [19] and it is therefore important to understand the effect of this on the ice flow. Previous work has shown that multi-valued sliding laws can result in ice-stream flow [11,13,[20][21][22][23]. More specifically, Fowler & Johnson [24,25] demonstrated that a 'hydraulic runaway' mechanism is suggestive of a triple-valued sliding law, based on thermomechanical feedbacks that arise as a result of having a distributed drainage system at the bed. Sayag & Tziperman [26] invoked a switch in drainage system as the motivation for a triple-valued sliding law, suggesting that each of two stable branches correspond to different drainage patterns at the bed.
While implementing a multi-valued flux law as a boundary condition may produce ice-stream behaviour, using knowledge of the subglacial drainage system beneath Antarctica more directly to explain ice-stream emergence would further our understanding of the system. Evidence from Antarctica suggests that the rapid basal velocities in ice streams are enabled by the presence of a layer of till at the base of the sheet [27][28][29][30]. For likely permeabilities and till thickness, Darcian porous flow through the sediment would be too slow to evacuate all the meltwater present at the bed [31,32] and the till is therefore water saturated and is thought to deform with Coulomb-plastic behaviour [33][34][35][36][37]. Knowledge of the yield strength of the till is important, but more information than this is required about the hydrological system and how its evolution affects the ice above. Recent work looking at ice-stream formation over till either disregards water transport altogether [38] or uses a simple diffusion equation to describe the evolution of a water layer over saturated till, assuming the water 'diffuses' from high to low effective thickness [39][40][41]. In this work, we wish to develop a more physically motivated model for the water layer and investigate the potential link between the hydrological system and multi-valued flux laws. This will extend the work done by Fowler & Johnson [24,25], who introduced the concept of hydraulic runaway and considered the feedbacks in a simplified one-dimensional model. In our more detailed hydrological model, we consider meltwater flowing over saturated, relatively impermeable till. A water film exists only if it is shallow enough not to submerge all surrounding bed protrusions [42]. If the water depth increases beyond some critical value (that at which all surrounding clasts are drowned by the water), the water film becomes unstable to the formation of local water streams incised in the sediment [32]. There is uncertainty surrounding some of the parametrizations made in the formulation of the model; the purpose of this work is to investigate the qualitative implications and feedbacks that arise from the coupling of the meltwater and ice. The paper is organized as follows. In §2, we present the ice-sheet model and derive governing equations for the water film depth and the effective pressure at the bed, before nondimensionalizing the model in §3. Section 4 considers behaviour of the water layer itself in more detail, specifically considering what happens in the model once the water layer becomes deep enough to submerge all surrounding clasts. In §5, we present numerical solutions to the fully coupled ice-water model, which, under certain conditions, results in ice-stream flow, and we discuss how this can be explained in relation to a multi-valued flux law and the chosen parametrizations in §6.

Governing equations (a) Ice flow
Given that ice is an incompressible, viscous fluid, an ice sheet moving under gravity is accurately described by Stokes flow. The low aspect ratio of ice sheets (approx. 10 −3 ) allows simplifications to be made to the three-dimensional nonlinear equations, with two end-member models being the shallow ice approximation (SIA) and the shallow shelf approximation (SSA). The former is a classical lubrication approximation, valid for flow frozen to the bed, and so dominated by vertical shear stresses [43][44][45]. The SSA, on the other hand, does not allow for internal deformation, and so provides an accurate flow description where horizontal ice velocities do not vary with depth [46,47], for example in ice streams or on ice shelves.
In this work, we use a model that is a vertically integrated hybrid of the SIA and SSA. It takes into account both vertical shear stresses and membrane stresses, so providing a valid flow description for all flow regimes within a shallow ice sheet. This is particularly important when modelling the emergence of ice streams. The force balance includes basal stresses, driving stresses and membrane stresses. A complete description of the model is presented in previous work [48], along with a comparison between this model and similar higher order approaches [49,50]; here just the non-dimensional equations are outlined at the end of §3.

(b) Subglacial water flow
We conceive of water flowing at the base of an ice stream as illustrated in figure 1. The ice is underlain by sediments that are saturated by water. Typically, the water is at high pressure and the till is deformable. We presume that the till is not very permeable and expect that the water generated at the base might flow in some sort of distributed system. The simplest such flow is a thin film at the interface between the ice and the till. With a clean separation between the ice and the bed, as imagined by Weertman [51], the resulting flow is an unstable configuration [52]. However, we can suppose such a film could exist stably if it is thinner than the supporting clast size [42], which is reasonable if the water film thickness is of the order of millimetres. Of course, this film is still subject to instability if it grows deep enough to separate the ice and bed, and we can expect local-scale water streams to form once the water layer depth reaches some critical value at which surrounding bed protrusions are submerged [32]. Note that these water streams form over much smaller length scales than ice streams; in this work, we consider the larger scale averaged behaviour of these small-scale features at the bed and postpone the study of a more detailed process of the formation of water streams for future work.
We consider Cartesian axes (x, y, z) where x points in the ice flow direction, y is across-flow and z is vertically upwards. As shown in figure 1, the till surface is denoted by z = b, the lower ice surface by z = s w and the resulting depth of the water film is The upper ice surface is denoted by z = s i . The hydraulic head driving water flow at the ice-till interface is where p w is water pressure, p a is atmospheric pressure, ρ w and ρ i are the densities of water and ice, respectively, and g is gravitational acceleration. ρ wi = ρ w − ρ i and N = p i − p w is the effective pressure. p i = p a + ρ i g(s i − s w ) is the basal ice pressure. Note that (2.2) reflects the well-known fact that basal water flow is driven primarily by the ice surface slope, and the basal slope only contributes approximately 1/10 to the flow direction.
Given that the water flows in a film of depth H, mass conservation of water takes the form where Γ is the melt rate of the ice (with units kg m −2 s −1 ). It includes contributions from local heat sources as well as from frictional water dissipation in the subglacial water flow. More specifically, where G is the geothermal heat flux, u b is the basal ice velocity, τ b is the basal ice stress, q T is the sensible heat flux into the overlying ice and L is the latent heat. It is important to note that the frictional heat source u b · τ b arises from an integration of the ice viscous dissipation term over the basal sliding region, and τ b and u b here refer to conditions near the base of the far-field ice-sheet flow, but far from the actual interface [53]; they are the quantities used in the ice-sheet sliding law. At this stage, we are also neglecting dissipation in the ice owing to lateral shear-an effect that could be significant in the regions around ice-stream shear margins [54,55]. Assuming a local Poiseuille flow in the water film implies that where η w is the viscosity of water.
(c) Ice closure relation While mass conservation of water tells us how the water moves in response to a hydraulic potential gradient, we also require an equation to describe the evolution of the subglacial hydraulic system. By analogy with the closure equation of [16], we consider a balance between the opening and closure rates of the system; we take  ice as being supported by load-bearing clasts, and W c represents the viscous creep of ice around these supporting clasts. Under the assumption that clasts have spacing l * (figure 1), we take the rate of closure owing to viscous creep as [42] W C = AN n l * .

(2.7)
A is the rate factor and n the exponent from Glen's flow law, which describes the rheological behaviour of ice [8]. N is the effective pressure; as the difference between ice overburden pressure and water pressure increases, the closure rate will also increase. Finally, the quantity l * is analogous to the channel width in the closure equation of Röthlisberger [16], assuming a single (wide) channel. Moreover, as the film thickness increases, the spacing between protruding clasts will increase as more of them become submerged; l * will therefore be an increasing function of H.
There is little to constrain our choice for the functional dependence of l * (H). We assume that there is some film depth H c at which all the clasts would become submerged and the ice would become locally clear of the bed. To be more specific, we will define where l 0 is a length scale that represents the typical clast spacing in the absence of water, and will be chosen to be consistent with observed effective pressures on the Whillans ice stream (B). We are then left to choose a function for Λ(H). The most simple choice, for illustrative purposes, would be a function such as providing support to the ice. More specifically, we presume that, when H approaches H c , the Walder instability is no longer suppressed, and the system becomes unstable to the formation of some kind of locally channelized drainage system at the ice-till interface [32]. We therefore consider the case where Λ(H) decreases to some small limit as H approaches H c , and specifically we take Λ(H) to be a function of the form where both Λ ∞ and δ are small. The solid black line in figure 2 illustrates this, with Λ ∞ = δ = 0.1. From (2.10), we have that N ∼ Λ 1/n . This finite limit of Λ(H) therefore now means that the effective pressure does not reach zero in our model. Locally, we expect that there are small areas where the effective pressure is zero and H > H c , but on the larger scale (over which we are interested in modelling when investigating ice streams) there are regions of finite support, now on a length scale of several metres. Our choice of Λ(H) presented here is meant to represent this globally averaged behaviour in the simplest possible way. It is also worth noting that, as we are considering water flow over sediment rather than hard bedrock, there should also be a description of bed evolution in the model. This would require a second closure relation for the bed elevation b, based on an Exner equation incorporating sediment transport. In the present work, however, we assume that there is no sediment transport so that b is fixed and H = s w − b. We therefore require only one closure equation. This is not an unreasonable assumption so long as there is no net change in bed elevation. There can still be some transport of sediment by virtue of the shear induced by the sliding of basal ice, but if the sliding is a uniform motion in the downflow direction, we would expect the deposition and erosion rates to be in equilibrium. Other work specifically considers bed evolution in an effort to explain the evolution of subglacial bedforms [56], but we postpone including details of these processes for future work.

(d) Mechanical coupling of the water and ice
It is necessary to prescribe a boundary condition at the interface between the ice and water in order to couple the subglacial hydraulic flow to the ice flow. This is done through a basal friction law, which relates the basal shear stress, τ b , to the hydrology, through the general relationship is the sliding velocity of the ice, and p and q are commonly taken as 1/3 [3,7]. We furthermore take c = 6.8 × 10 4 Pa 2/3 s 2/3 m −1/3 so as to give a sliding speed of the order of 60 m yr −1 for typical driving stresses of 10 4 Pa [57].
Coupling also occurs through feedbacks of the ice on the drainage system, in particular through the melting caused by basal friction and the creep closure of the ice over the water film. The prescription of the cooling rate, q T , in the melt rate expression (2.4) also depends on the ice velocity. To prescribe this, we first consider a simplified temperature equation. In a thermal boundary layer, the dominant balance is expected between heat advection and vertical diffusion where κ T is the thermal diffusivity of ice (and lateral heat advection is ignored because we expect u b v b in the region of ice streams). With T = 0 on z = b and T → − T far from the bed, we obtain an error function as a similarity solution of (2.13), which yields the average basal heat flux in the form     is large. For velocities u 0 ∼ 100 m yr −1 , length x 0 ∼ 500 km and depth d 0 ∼ 10 3 m, Pe = 5, and it is larger for ice-streaming speeds. The boundary layer approximation is not, therefore, a highly accurate approximation, but it probably gives a sufficiently accurate representation for the present purpose, and indicates correctly the increase of basal heat flux with basal velocity. Furthermore, the q T boundary layer formula is inappropriate at large x ∼ u 0 d 2 0 /4κ T , when the temperature profile becomes more conductive, or for small x near divides where u b → 0 and the x dependence becomes important. The Robin solution which includes vertical but not horizontal advection is often used here, but a better approach uses the von Mises transformation [58] which includes the Robin solution, and allows the similarity solution to approach divides by allowing for the x dependence of u b . However, we omit such subtleties here.
From our simplified considerations, we therefore have that the sensible heat flux increases with u 1/2 b . The frictional heat, however, increases with u b , so, for large u b , Γ (2.4) is positive. It is possible for intermediate u b that Γ could reach zero and then become negative; our model still applies until the water film thickness vanishes, at which point sliding law (2.12) would cease to apply and the model would need modification. However, in our simulations this never occurs because G and τ are sufficiently large to prevent refreezing.

Non-dimensionalization and reduction
Together with the equations for ice flow [48], the coupled system of equations includes the mass conservation of water equation (2.3) and the creep closure equation (2.6). Values of the dimensional constants are given in table 1. To couple the ice and water systems at the basal boundary, we have the sliding law (2.12).

(a) Non-dimensionalization of the subglacial water equations
To make the subglacial water model dimensionless, we guess approximate scales for some variables based on observations/physical intuition, and then use these to derive characteristic scales for other variables, as given in table 2. We assume a characteristic value for the horizontal extent of the ice x 0 ∼ 500 km based on length scales in ice-streaming regions, for example the Siple   Coast [8]. We also define an effective pressure scale approximately 4 × 10 4 Pa [59] and choose other scales by defining and The ice flow scales are furthermore defined by and where the ice flux scale Q i is either prescribed (from upstream supply) or related to the typical accumulation rate a 0 by Equations (3.8) and (3.9), together with our prescribed sliding law (2.12), determine τ 0 , u 0 and d 0 . However, as the sliding law is not well constrained, we will use observed values for these scales.     Specifically, in the Siple Coast region of Antarctica, typical ice depths are 1000 m [60]. Given a velocity scale u 0 ∼ 100 m yr −1 , we take the accumulation rate scale as Finally, as the water flows on a much shorter time scale than the ice, we define two separate time scales and (3.13) Table 2 gives values of these model variable scales. We non-dimensionalize all equations with the larger of the time scales, t i , and so introduce a non-dimensional parameter δ T = t w /t i . Using these scalings, the dimensionless form of the subglacial flow equations becomes and

(b) Simplification of the subglacial system
The parameters σ , ν and χ are all small. We thus neglect the terms of order σ and χ , representing gradients in water depth and melting due to meltwater dissipation, respectively, but retain the term in ν, as it represents a singular perturbation. From (3.15), we have (3.18) and the mass conservation of water equation (3.14) can then be written as an equation for H Furthermore, to allow consideration of the behaviour of the hydraulic system on its own, we temporarily assume that the ice has a constant, uniform surface slope 20) and zero basal slope, ∇b = 0. Equation (3.19) then reduces to Given that Λ(H) is a decreasing function of H, the essential structure of the model can be described by which is closely related to that studied by Benjamin et al. [61] as a model for long waves in shallow water; we may infer that the present equation is well posed.

(c) Full set of governing non-dimensional equations
We now have a full set of non-dimensional governing equations for the water system. From Kyrke-Smith et al. [48], the non-dimensional mass conservation and force balance for the ice are given by ∂h ∂t 3.23) and

25)
S is the resistive stress tensor [62,63], and a is the accumulation rate. τ is the deviatoric stress tensor, related to strain rate through Glen's flow law for ice [8]. The full derivation and non-dimensionalization are provided in the appendix of a previous publication [48].
Together with the non-dimensional sliding law as the boundary condition at the water-ice interface, (3.27) and the non-dimensional mass conservation of water equation (3.19), these equations provide a complete non-dimensionalized description of the ice-water system. Variable scales are in

Results for hydrology with constant ice (a) Hydraulic runaway
Hydraulic runaway is a physical phenomenon that models suggest may occur as a result of a positive feedback between ice velocity and melt rate [24]. As ice velocities increase, more meltwater is produced by frictional heating, which results in further ice-bed lubrication, and hence faster velocities. We consider whether the basal hydrology model presented here results in such behaviour. The melting term, Γ in (3.21), depends on the ice flow. It therefore also has indirect dependence on the water depth, H, as the sliding velocity u b depends on H through its dependence on N. The ice depth, h, evolves over a long time scale, t i , through the solution of (3.23) and (3.24). We are initially interested in changes in behaviour on a short time scale, so we assume that the ice is of constant depth, with constant surface slope (3.20). Taking the shallow-ice version of basal shear stress as an approximate sliding law, together with the water equations, we then have and which gives three equations for the three unknowns, u b , Γ and N. These can be solved to give Γ = Γ (H, H t , h). We seek to eliminate u and N from (4.1) to (4.3). Equation (4.2) gives us an expression for N(Γ , H) and substituting (4.2) into (4.1) we have an expression for u b (Γ , H) which is This can then be substituted into (4.3) so that we have Γ (H, H t , h).
In figure 3, we plot the quasi-static instance with ∂H/∂t ≈ 0 and Λ(H) given by (2.11). As H increases, the melt rate also increases owing to an increase in frictional heating. This in turn comes about from an increase in basal velocity by (4.4). From this, it is evident that some kind of hydraulic runaway is inherent in the model, with melt rate increasing rapidly as H → H c . However, complete runaway is avoided as we do not allow the model to reach zero effective pressure in the subglacial system, which would correspond to no basal stress, and so allow the ice to increase in velocity without limit.

(b) Solutions with a stable water film
We consider solutions of the decoupled water layer equation, to see its basic state behaviour. We take the ice to have constant surface slope and zero basal slope, as in the simplification given by equation (3.21). The coupling to the ice occurs through the melt rate term, which, from above, we know increases with both ice depth h and water layer depth H. We therefore here take an illustrative Γ function, which results in an analytic steady-state solution to (3.21) at leading order. This steady-state solution is given by is the non-dimensional position downstream where H reaches H c . This corresponds to the point at which the bed becomes submerged and the thin film flow becomes a stream. If the bed is sufficiently rough (large H c ) or the ice depth or flux is sufficiently small, then this point would be downstream of the grounding line, and the film flow can exist everywhere. We ran a suite of simulations solving (3.21) with a melt rate of form (4.5). As expected, so long as H c > (2h) 1/6 then the solutions evolve into the steady state given by (4.6) (e.g. figure 4). However, if H c < (2h) 1/6 then the water depth reaches its critical value before the grounding line. The melt rate is no longer defined and an analytic solution does not exist. Our interest therefore now turns to the case where H may reach H c , by considering the fully coupled water-ice system.

Results for coupled hydrology and ice dynamics
Having formulated a description of the subglacial water layer, we now couple it with the ice to explore the combined behaviour of the system. The water film equation takes the same form, given by (3.19), and the ice flow equations are given by (3.23) and (3.24).

(a) Model set-up and boundary conditions
We consider a 500 × 500 km domain with a uniform bedslope. The upstream end x = 0 is a coldtemperate transition line with no water at the bed (H = 0). The ice depth satisfies ∂h/∂x = 0 across the boundary, and the inflow is at a constant velocity into the domain (u b = u in at x = 0). Free slip is also allowed along this inflow boundary (τ 12 = 0). To avoid considering grounding-line behaviour, we treat x = 1 as an open outflow boundary [26,64], where the upper surface is free to evolve and ∂u/∂x = 0. Furthermore, a suitable condition for the water layer thickness, H, which also avoids the possibility of boundary layers, is ∂H/∂x = 0. The domain is taken to be periodic in the y-direction. With these boundary conditions, we solve (3.21) coupled to the hybrid ice flow model outlined in §3c. The method for the numerical solution, using the Portable Extensible Toolkit for Scientific Computation [65][66][67], is described in previous work [48]. In that work, we used a heuristic triple-valued sliding law together with an equation to describe the relaxation time for the ice to adjust to changes at the bed; these have been replaced with a sliding law that directly depends on the subglacial conditions (3.27) and the physics-based water layer equation (3.19).
The simulations are initialized with a downstream ice flow that is uniform in the cross-flow direction, i.e. (u b , v b ) = (u(x), 0). The velocity field u(x) is computed from the x-component of equation (3.24) with (3.27). There is no water at the bed and the ice has a constant uniform thickness, h = 1.5. A constant uniform accumulation rate is applied over the domain. Average accumulation rates in the Siple Coast region of Antarctica are approximately 140 kg m −2 yr −1 [68], and so we take a ∼ 140 kg m −2 yr −1 /917 kg m −3 = 0.15 m yr −1 as a suitable value. We expect the accumulation rate to be an important parameter in governing the flow behaviour, but do not present a detailed parameter study of the system in this paper; here, we discuss results that illustrate the basic flow regimes that occur as a result of a realistic accumulation rate applied over the domain.

(b) Flow regimes (i) Laterally uniform water and ice flow
We first consider the case where the water film does not get thick enough to submerge all surrounding clasts (i.e. the film depth remains less than the critical depth H c ). If this is the case, the coupled system evolves into a stable steady state. This is illustrated in figure 5, where the critical water depth is H c = 2, corresponding to 6 mm.
The water film evolves rapidly into a quasi-steady-state, similar to that seen in §4b, when the water layer is considered alone with a prescribed melt rate. Meanwhile, the ice depth and basal velocities also evolve into a steady state over a longer time scale. Ice depth remains uniform, increasing by approximately 200 m over the time of the simulation run. This is due to the downstream ice velocity initially being insufficient to balance the incoming ice flux and accumulation. However, the downstream ice velocity increases and over approximately 5000 years evolves to put the system in equilibrium. The velocity profile becomes a monotonically increasing function downstream. Over this time, the downstream profile of the water layer remains very similar to the initial profile it reaches within the first year. It does increase in depth by approximately 0.25 mm downstream owing to increases in ice velocity and ice depth that cause a larger melt rate from frictional heating. The water layer rapidly adjusts for this as the melt rate increase occurs.
Overall the model is well behaved, and the coupled behaviour is straightforward for the case where there is insufficient meltwater at the bed to submerge all supporting clasts.

(ii) Laterally unstable stream flow
Our interest now turns to the case where the amount of meltwater is sufficient to locally drown all supporting clasts; we are interested in what effect the decrease in support provided by the bed will have on the ice. Analysis in §4 suggested that melt rate will increase rapidly as the water  depth reaches its critical value. In that analysis, we assumed that the ice was stationary as we were only interested in behaviour over a short time scale. As a result of that assumption, there is nothing to prevent water depth increasing indefinitely, in turn also allowing the ice velocity to increase (because there is little support at the bed). However, there is a long time-scale feedback that was not taken into account; over this longer time scale the ice depth decreases as a result of increased outflow from the domain. This causes the basal stress, and thus Γ , to decrease. The water depth will start to decrease as a result, therefore decreasing the ice velocity and the total ice flux. This cycle could then repeat once the water depth builds up again to reach the critical depth. In a narrow, confined flow, we might expect this to lead to cyclic surging behaviour (similar to [69]). However, in a wide flow, a lateral instability that leads to ice streaming occurs; we illustrate this with numerical solutions to the governing equations. Figure 6 shows maps of basal velocity, water depth and surface elevation at three different times for a simulation run with H c = 1.5, corresponding to a critical water depth of 4.5 mm. It is evident that this small change in the critical water depth results in dramatic changes in the behaviour of both the water layer and the ice, when compared with the results of the previous section. Initially, the flow is stable and uniform cross-stream, but, as more water is generated at the bed, H approaches H c and an instability emerges, as shown in the plots at t = 82 years. There is a sudden increase in water depth near the boundary, which results in an increase in ice velocity, because the ice can slip more easily over the bed. A consequential surface lowering of the ice is also evident in the ice depth field. A non-uniformity in the water depth and the ice velocity develops across the domain at this time, as a result of lateral instability (there has been no imposed perturbation to the system). Patches of faster ice flow form, corresponding to the positioning of the discrete areas of deeper water. There is a coarsening of this instability as the system evolves to create distinct stream-like features with a longer wavelength than the initial instability (t = 105 years). The streamlines of the water flow show that water from upstream is drawn towards these water-rich areas, resulting in a positive feedback. These unstable features propagate up the domain, and evolve to create pronounced ice streams, flowing at hundreds of metres per year and of width approximately 60 km (t = 360 years). The values are of the same order of magnitude as observed velocity and width values of ice streams in Antarctica [8], and in other numerical studies [38,70]. A discussion of the stress balances that are necessary to maintain ice-stream flow for different basal conditions is in previous work [48] and demonstrates the  importance of membrane mechanics in determining ice-stream width. Again, the channelizing of water towards these deep-water patches is evident from the streamlines of water flow. There are depressions in surface elevation at the onset region of the ice streams, as a result of the rapid increase in discharge from the region. This is consistent with some observations [71]. At this point, the simulation has reached a quasi-steady state. The distinct streams do not grow further. There are, however, periodic wave-like features at the edge of the ice streams that propagate along the boundaries (particularly clear in the plots of water depth at 360 years and in the movie of the velocity field provided in the electronic supplementary material). These could be the result of travelling wave solutions that exist for the water layer equation (3.19), although this is purely speculative at present.
Plots of effective pressure and melt rate at t = 360 years are shown in figure 7. The finescale structures at the edge of the ice streams are evident in these fields as well. As would be expected from the relationship between N and H (equation (4.2)), we see that the areas of low effective pressure correspond directly to where there is a lot of meltwater. Low effective pressure furthermore allows the ice to move more rapidly because there is little resistance provided by the bed in these areas (from the sliding law (3.27)).
Moreover, the melt rate field illustrates the positive feedback that occurs with more melt being produced underneath the ice streams owing to frictional heating from the rapid velocities. The melt rate is, in fact, the greatest at the onset zone of the ice streams. This is due to the large driving stress here from the gradient in surface elevation, combined with a higher effective pressure, and therefore higher basal stress, which results in a larger melt rate from frictional dissipation. There is therefore enhanced melting at the onset zone-a feedback that will help to maintain ice-stream flow.

Discussion (a) Relationship to multi-valued flux laws
We have seen in the previous section that, under certain parameter combinations, ice flow can become laterally unstable; this leads to two distinct modes of flow-fast ice streams and slowmoving ice. This is similar to the behaviour seen in numerical models that use a triple-valued sliding law as a basal boundary condition [26,48]. Can the streaming results from the coupled ice-hydrology model be understood in terms of an emergent triple-valued sliding law, as with the more simplistic model of Fowler & Johnson [25]?
The water film evolves on a much faster time scale than the ice. As it is the ice flux we are interested in here, we assume now that the water film is in equilibrium, i.e. we take  where ∂s i /∂x = −Ω is the ice surface slope, a constant in the downflow direction (we also assume there is no cross-flow variation in H or Ω). This is in contrast to the analysis of the water film in §4a when we were interested in the behaviour of the water film over a short time scale and so assumed the ice was in equilibrium.
Assuming that the melt rate is constant with x in this steady-state water case, and integrating (6.1) over x gives under the assumption that H = 0 at x = 0, i.e. the inflow boundary is taken to be a divide or a cold-temperature transition line. We can write the non-dimensionalized sliding law (3.27) as where This gives a relationship between ice depth, h, and basal ice velocity, u b , for ice overlying a water film which is in equilibrium.
Using Λ(H) as described by (2.11), we have from (6.4) with Γ once again given by (6.5). We solve this algebraic relationship between h and u b [72] and plot the solution in figure 8a with the same parameter values that produced the ice-streaming behaviour in §5b(ii). In this case, the relationship between ice depth and ice velocity is multi-valued, suggesting that distinct velocity states can occur at the same ice thickness, h. In figure 8b, we plot this multi-valued relationship against a rescaled version of the triple-valued sliding law derived in Sayag & Tziperman [26]. In their work, the sliding law is obtained analytically from a simple onedimensional force balance, based on the cross-flow structure of an ice stream; it has been used for subsequent analysis of ice-stream behaviour by [48,73]. We see that the sliding law that arises from a basic coupling of our hydrology model and the ice flow exhibits very similar properties to the triple-valued function from Sayag & Tziperman [26]. As a result, when solving our coupled model, similar behaviour arises to that seen when using a triple-valued sliding law of the form presented in [26] as the basal boundary condition.

(b) Parameter dependence
The analysis above allows us to investigate what parameters are critical in governing multivalued flux behaviour. Throughout the paper, we have used Λ ∞ = δ = 0.01 and have shown in figure 8 that this results in a multi-valued flux law when H c = 1.5 (corresponding to 4.5 mm). We were motivated to plot this as we observed multi-valued behaviour from numerical simulations run with this parameter combination. The dependence of this behaviour on H c is of interest; in figure 9, we therefore plot the relationship between ice depth and ice velocity for three different values of H c with Λ ∞ = δ = 0.01. We see that while for very small H c the relationship is not multi-valued, as H c increases and the function becomes multi-valued the effect of further increasing H c is that the turning points of the function occur at larger values of h and u b . This is because more meltwater is required for the water film depth to come close to the critical depth; a larger ice depth results in a larger driving stress that in turn causes more frictional heating.
For larger values of H c (e.g. figure 5), it is necessary to have thicker ice if there is to be sufficient meltwater to initiate the instability. This is in contrast to the behaviour we would expect if H c were reduced too low; the flux law is no longer multi-valued, so we would not expect multiple velocity states to coexist in steady state, however large the driving stress. As well as considering variations in H c , in figure 10 we plot a phase diagram showing the distinct regions of Λ ∞ -δ parameter space where the relationship is triple-valued and singlevalued for the case with H c = 1.5. It is clear that if both Λ ∞ and δ are sufficiently small then the relationship (6.6) is multi-valued and there is therefore the potential for different velocity states to coexist.

Summary and conclusion
The theory developed in §2 provides a simplified model of subglacial water flow at the base of an ice sheet. We assume the meltwater flows in a Weertman-Creyts-Schoof style rough-bedded film [42,51] and show that this results in hydraulic runaway of the water film when the depth of the water layer reaches some critical depth H c at which all supporting clasts are submerged. In reality, a water sheet with depth greater than H c is an unstable configuration [52] and we would expect the water film to form discrete streams incised into the sediment; we present a simplified representation of this idea motivated by the fact that while locally there might be areas of zero effective pressure, at a larger scale (that is still sub-grid scale) there will be regions of finite support between water streams. Numerical solutions of the fully coupled ice-water system result in distinct flow regimes dependent on whether the water depth reaches the critical depth within the domain. Under the model set-up considered in this paper, the system evolves into a uniform steady state if the bed is sufficiently rough (i.e. H c sufficiently large) that all the melt can flow in the thin film configuration. However, reducing the critical water depth (so that H reaches H c ) results in the ice flow becoming laterally unstable and non-uniform in the cross-flow direction; ice-stream-like features develop.
A basic analysis of the coupled system suggests that the hydrological system can result in a triple-valued flux relationship for the ice. It has previously been shown that a triple-valued sliding law is associated with ice-stream behaviour and such multi-valued flux laws have therefore been used to investigate ice-stream dynamics [26,48,73]. By showing that this coupling between subglacial hydrology and the ice can result in similar multi-valued flux behaviour, this work helps to contextualize analyses of ice-stream flow that result from using a triple-valued sliding law at the bed. Given that the sliding law is multi-valued, however, does not enforce instability in the system for all simulations; it still requires the ice to get sufficiently thick and the water layer to be sufficiently deep. Furthermore, some parameter combinations (we specifically consider variations of Λ ∞ , δ and H c in this paper) do not result in triple-valued flux behaviour at all. In these cases, under no circumstances would we expect a lateral instability in the system to occur.
While this paper presents theory and intriguing results associated with a new coupling between physically reasonable models of ice and hydrology, a more detailed investigation of the coupled model behaviour is needed. Future work will provide a more comprehensive parameter study, as well as analysing in detail the stress balances that occur across emergent ice streams, in a bid to explain their spacing. Furthermore, an important limitation of the current model is that it neglects energy conservation. The model instead assumes that the bed is always at the melting point. While this might not be an unfair assumption in the region of ice streams, other work with thermodynamically coupled ice dynamics has shown that fast flow may begin owing to purely thermomechanical feedbacks [70,74,75]. It is interesting to observe from these studies that similar ice-stream widths are obtained both in thermal theories and in the hydraulic theory presented here. We expect that this is due to the role of membrane mechanics in resisting ice-stream flow, as discussed in our previous work investigating the necessary stress balances to maintain ice-stream flow [48]. It will be important in the future to consider a detailed model with energy conservation included, to explore the result of combined thermal and hydraulic feedbacks in a bid to further our understanding of ice-stream dynamics.