Oscillatory subglacial drainage in the absence of surface melt

The presence of strong diurnal cycling in basal wa- ter pressure records obtained during the melt season is well established for many glaciers. The behaviour of the drainage system outside the melt season is less well understood. Here we present borehole observations from a surge-type valley glacier in the St Elias Mountains, Yukon Territory, Canada. Our data indicate the onset of strongly correlated multi-day oscillations in water pressure in multiple boreholes strad- dling a main drainage axis, starting several weeks after the disappearance of a dominant diurnal mode in August 2011 and persisting until at least January 2012, when multiple data loggers suffered power failure. Jokulhlaups provide a tem- plate for understanding spontaneous water pressure oscilla- tions not driven by external supply variability. Using a sub- glacial drainage model, we show that water pressure oscilla- tions can also be driven on a much smaller scale by the inter- action between conduit growth and distributed water storage in smaller water pockets, basal crevasses and moulins, and that oscillations can be triggered when water supply drops below a critical value. We suggest this in combination with a steady background supply of water from ground water or englacial drainage as a possible explanation for the observed wintertime pressure oscillations.


Introduction
The drainage of melt water along the glacier bed can have a significant effect on ice flow velocities by modulating basal sliding. This has been documented for valley glaciers (e.g. Iken and Bindschadler, 1986;Sugiyama and Gudmundsson, 2004) as well as the marginal areas of the Greenland Ice Sheet (Zwally et al., 2002;Shepherd et al., 2009;Bartholomew et al., 2010). For systems fed by surface melt, the effect of subglacial drainage is generally strongest in summer. In Greenland, summer melt has been associated with a net speed-up as well as a net slow-down (Bartholomew et al., 2010;Palmer et al., 2011). The seasonal evolution of the drainage system is often inferred to play a dominant role: a less efficient drainage system in the early melt season can impede water flow and lead to high basal water pressures and faster sliding velocities, while a more highly developed, channelized system may lead to low water pressures and consequently to suppressed sliding velocities (see also Schoof, 2010;Hewitt, 2013).
The presence of strong diurnal cycles in water input is a key characteristic of many subglacial drainage systems during summer, leading to corresponding diurnal cycles in basal water pressure (Hubbard et al., 1995;Fudge et al., 2008;Shepherd et al., 2009) and possibly playing a role in ice flow speed-up (Schoof, 2010;Hewitt, 2013). The termination of this strong diurnal signal often occurs while there is still measurable surface melting (e.g. Fudge et al., 2005). The mechanisms behind the termination and the subsequent evolution of the drainage system are poorly understood. Different areas of the bed may simply become hydraulically disconnected from the surface and one another, or a remnant active drainage system may remain.
The wintertime behaviour of drainage and basal sliding is however crucial to long-term ice flow dynamics: the summer melt season occupies only a small part of the annual cycle, and even a significant increase in summer velocities compared with winter velocities may have a negligible effect on the annually averaged ice flux. Wintertime velocities therefore play a key role in determining the annually averaged flux. They need not be constant but can potentially be affected by the evolution of the drainage system: recent evidence from Greenland Tedstone et al., Published by Copernicus Publications on behalf of the European Geosciences Union. indicates that a hot summer and a presumably highly developed drainage system may lead to suppressed wintertime velocities, possibly as the result of continued, efficient drainage past the termination of surface melt input. In addition, glacier surges are among the most dramatic examples of changes in sliding behaviour triggered by changes in basal water, with many known to start in winter (Kamb et al., 1985;Lingle and Fatland, 2003).
Here we present direct evidence for an active subglacial drainage system persisting until at least mid-winter under a small surge-type valley glacier in the St Elias mountains, Yukon Territory, Canada. Borehole pressure records from the field site show that three out of six boreholes that straddle a probable drainage axis exhibit highly correlated pressure variations long past the termination of diurnal variations in pressure and of surface melt. Starting in September, these pressure variations take the form of oscillations that last until the end of the observation period in January, when surface temperatures are far below the melting point.
We interpret the observed variations as indicating the presence of a remnant drainage system fed by low levels of background water supply (for instance in the form of ground water), and infer that the system can undergo internally driven pressure oscillations. To explore what could trigger such oscillations, we use a one-dimensional drainage model that couples water flow to distributed englacial storage. Our results show that oscillations can be triggered spontaneously as an instability analogous to that driving jökulhlaups (Nye, 1976) when water input drops below a critical level, and that no oscillations in water supply are required to explain them.
In addition, our observations also suggest that the drainage system during summer may have a nontrivial response to diurnally varying forcing, with basal water pressures exhibiting an apparent two-day periodicity. We use our model to show that water storage may also cause this behaviour, indicating that storage capacity may play a prominent role in modulating drainage year-round.

Field site and methods
The glacier at the study site, named "South Glacier" for consistency with previous work, is a 5 km long polythermal glacier in the southern Donjek Range, located at 60 • 49 N, 139 • 8 W (see Fig. 1). The glacier is described in greater detail in De Paoli andFlowers (2009), Flowers et al. (2011) and Flowers et al. (2014). Aerial photographs taken during the 1980s indicate that the glacier last underwent a surge in 1986 (Peter Johnson, pers. comm., 2005; see also Johnson and Kasper , 1992;Flowers et al., 2014). It currently ranges in elevation from 1900 to 2900 m a.s.l., with an approximate equilibrium line altitude of 2550 m (Wheler, 2009), and is confined in a valley of L-shaped planform, with two tributary valleys. Digital elevation models (DEMs) for both surface and bedrock have been derived from extensive radar and Global Positioning System (GPS) surveys (De Paoli, 2009;Wilson, 2012).
Exposed bedrock in the valley consists mainly of highly fractured granodiorite, as well as some metasedimentary units near the glacier snout. Videos taken in boreholes through the glacier have shown the presence of granodiorite cobbles in basal ice, as well as of highly turbid water in the lowermost 3-5 m of freshly drilled, undrained boreholes. A glacier bed consisting of boulders and subglacial till has been exposed in two collapse features that have formed since 2011 in the lower ablation area of the glacier. There are extensive lateral moraines, some of them ice-cored, bordering the lower third of the glacier.
Between 2008 and 2011, 76 boreholes were drilled to the glacier bed in the upper ablation area, as shown in Fig. 1. At present, ice velocities in the study area measured in situ using Global Positioning System receivers are around 10-30 m yr −1 , with ice thicknesses ranging from about 45 to 100 m. There are numerous crevasse fields both inside and upstream of the study area, but no moulins. Two major surface streams originate within the study area, which has significant surface topography, but none enter it from above.
The hot-water drill used was previously employed in the Trapridge Glacier project (e.g. Murray and Clarke, 1995;Waddington and Clarke, 1995;Stone et al., 1997). All but three of the boreholes drilled on South Glacier were instrumented with pressure transducers. The data used in the present study were gathered using Barksdale 422H2-06A transducers, attached to data loggers via 24-gauge copper signal cable and cast in epoxy to ensure waterproofing. The transducers were calibrated in water-filled boreholes, and generally conformed to factory calibration (within a few percent) except for a change in offset that appeared to be specific to the length of signal cable used. The transducers were installed about 20 cm above the bottom of the boreholes, and connected to either Campbell Scientific CR10X or CR1000 data loggers. These were set up to log continuously, with logging intervals of 1 min for CR1000s and 2 min for CR10Xs during July and August, changing to 20 min from September to June. For the boreholes described in this paper, contact with the bed was deemed to have been established if water samples taken from the bottom of the holes showed significant turbidity. Six of the holes (A1-A3, B, C and D in Fig. 1) also drained partially near the end of the drilling process, indicating a connection to the bed. Boreholes typically froze shut within 1 to 3 days, after which they became isolated from the glacier surface. One borehole (marked A6) was additionally instrumented with Maxim DS18B20 digital thermometers on a custom-made cable supplied by Beases Stream LLC (Wilson et al., 2013).
Meteorological data were recorded by an automatic weather station on the glacier just downstream of the study area (MacDougall et al., 2011). In this paper, we utilize the temperature record from a HMP45C212 TRH air temperature probe as well as surface height soundings recorded by a and (c). The coordinate system used is UTM zone 7 north. The superimposed stars in (b) and (c) indicate the location of boreholes, with labelled stars (displayed as larger symbols) corresponding to boreholes discussed in the text. "AWS" marks the location of the automatic weather station. (c) shows the same study area, but now displays surface contours spaced at 25 m intervals, and an "upstream area" distribution based on Eq.
(1) with f = 0.5. We display this as a specific upstream area, which we have defined here as the upstream area of a single pixel in the image divided by the edge length of that pixel (20 m).
Campbell SR50 sonic ranger. The surface height record was converted to an inferred specific surface mass balance record (given in water height equivalent per day) by assuming an ice density of 900 kg m −3 , a snow density of 200 kg m −3 for new snow, and 300 kg m −3 for old snow (unpublished data, Simon Fraser University Glaciology Group). In addition to melt, observed height changes can also result from snowfall and snow redistribution. An RM Young 05103-10 wind monitor installed at the weather station provides an indication of plausible redistribution events, when significant surface height changes occur during periods of cold temperatures accompanied by strong winds, while a Kipp and Zonen CMA6 albedometer gives an indication of likely snowfall events. A stream gauging site was also established seasonally during July of each year at the main outlet stream from South Glacier, approximately 2 km downstream of the study site. The discharge of this outlet stream has multiple contributions in addition to subglacial flow through the study area, as multiple glacier surface streams and at least one major stream draining from a valley side wall are routed to the glacier bed through moulins and lateral channel portals below the study area. As a result, we have not made use of the stream gauging data for the present study.

Data
In the present study, we focus on a set of six boreholes, labelled A1-A6, spaced at 15 m intervals along a line that traverses a small plateau on the glacier surface. An "upstream area" calculation (Sharp et al., 1993) using the glacier surface and bed DEMs in Wilson (2012) suggests that this line of boreholes may intersect one of the main subglacial drainage axes through the study area (Fig. 2). The upstream area calcu-lation is based on the "D ∞ " method in Tarboton (1997) and the algorithm for handling local potential minima in Wang and Liu (2006), and assumes a hydraulic potential of the form where ρ w = 1000 kg m −3 is the density of water, ρ i = 910 kg m −3 is the density of ice and g = 9.8 m s −2 is acceleration due to gravity. f is a factor intended to mimic the effect of overburden in routing water flow. f = 0 would correspond to a spatially uniform water pressure, for instance equal to atmospheric pressure (e.g. Hooke, 1984), and f = 1 to a spatially uniform effective pressure at the bed (that is, a constant difference between overburden and water pressure over the domain; see e.g. Shreve, 1972;Sharp et al., 1993). The boreholes were drilled and instrumented between 17 and 23 July 2011. Pressure in borehole A1 was logged by a CR10X data logger until 27 February 2012, while A2-A6 were logged by a CR1000 data logger until 17 January 2012. Both loggers then suffered power failure, most likely due to their solar panels being covered by snow. We also compare the water pressure time series from these boreholes to that from a further borehole, B, that lies about 150 m west-southwest of A1 and is connected to a third data logger, and four further comparison holes, C-F. C lies about 150 m downslope from B and was instrumented in 2010, while D was drilled about 400 m to the south-west of the location of holes A1-A6 during 2009. Hole D may also lie near a main subglacial drainage axis as indicated by Fig. 1. Holes E and F lie roughly to the south-south-east of hole C, and were drilled in 2010 and 2011, respectively.
All water pressure time series are plotted in units of metres as p w /(ρ w g), where p w is the measured water pressure. In other words, we plot the contribution of water pressure to hydraulic head as a function of time or, equally, the water level above the bed that would correspond to the measured pressure. Each pressure record was resampled at 3 h intervals for plotting purposes. Alongside water pressures, we also plot surface temperature and specific surface mass balance estimated from SR50 measurements, shown only when negative and expressed in millimetres of water. Specific surface mass balance is calculated once per day from the output of the sonic ranger at the automatic weather station. Figure 3 shows data for boreholes A1-A6 and B. Each panel is separated into two columns, corresponding to different time intervals. In order to resolve diurnal variability, late July and August 2011 are shown in the left-hand column at greater temporal resolution than September 2011-January 2012 in the right-hand column. Air temperature and specific surface mass balance are shown as a solid red line and a blue outline in panel a, while the pressure records are shown in panels b-d, with a legend identifying each time series shown in the corresponding panel. The horizontal dashed lines show overburden for each borehole, computed from the depth of the hole as measured using the signal cable attached to the pressure transducer.
The water pressure time series in boreholes A1-A3 can be divided into four distinct phases (see Fig. 3), some of which also describe the behaviour in holes A4 and A5. Phase 1 lasted from the installation of the sensors until about 5 August, and corresponds to strong diurnal pressure cycles during July and early August. Boreholes A1-A3 drained when the drill reached the bed. All three boreholes immediately displayed a strong and almost identical diurnal cycle, except for a constant offset (Fig. 3). This continues until there is a hiatus in diurnal cycling between 24 and 29 July, after which water pressure once more exhibits pronounced daily peaks. One intriguing feature of the diurnal cycle in holes A1-A3, especially following the 24-29 July hiatus, is that the signal appears to consist of "doublets", with a higher water pressure peak one day followed by a lower peak the next. This corresponds to a two day-period overprinted on the diurnal signal.
Boreholes A4-A6 initially showed either no obvious diurnal signal or only a weak diurnal signal in anti-phase with A1-A3. Holes A4-A5 then switched abruptly to showing the same signal as A1-A3 from 2 August, again with a constant offset, while A6 showed a still weak signal, which was however now more clearly anti-correlated with those in A1-A5 (note for instance the pressure peaks in holes A1-A5 on 3 and 5 August corresponding to distinct pressure troughs in hole A6). Hole B also exhibited noticeable diurnal cycling but with smaller amplitude. The variations in diurnal cycle amplitude are also different from holes A1-A6, and obvious diurnal cycles remain visible in the time series for hole B until 10 August.
Phase 2 (approximately August 5-10) corresponds to the final disappearance of the dominant diurnal signal, and a gradual drop in water pressure over several days in early August. In holes A1-A5, the diurnal cycle amplitude is reduced abruptly after 4 August, with only very slight diurnal pressure variations, accompanied by a gradual decrease in mean water pressure. The pressure signals remain highly correlated at the beginning of this drop-off phase, until about 8 August. Hole B by contrast ceased to experience strong diurnal cycles only on 10 August, with water pressure remaining relatively steady thereafter. Phase 3 (approximately 10 August-4 September) is a period of gradual increase in water pressure over a period of several days, followed by stabilization. Not all holes in the line A1-A5 exhibited this behaviour at the same time, and some of the time series are qualitatively different. The water pressures in holes A4 and A5 experienced an earlier increase, starting on 7 August, which resulted in water pressure rising to values near or just above overburden. In the process, the pressure records from A4 to A6 diverged noticeably not only from those measured in A1-A3, but also from each other. By contrast, water pressure in holes A1-A3 did not begin to rise until 15 August. Around 22 August, holes A1-A3 reached a water pressure that then remained nearly constant until 5 September. This level is slightly lower than mean water pressure recorded during phase 1, and remains significantly below overburden. Throughout phase 3, A1-A3 produced almost identical records, again with constant offsets.
Phase 4 began with an abrupt drop in water pressure on 5 September in holes A1-A3, equivalent to about 15 m in hydraulic head over one day. Initially, all three holes continued to show the same signal during the subsequent, more gradual drop in pressure, but then began to diverge. Over several days following 16 September, the record generated by hole A1 undergoes several abrupt jumps between consecutive measurements, equivalent to several metres of head, not reproduced by any of the other transducers. These jumps are shown as gaps in the time series in Fig. 3. One possible instrumental cause of such abrupt jumps is damage to the transducer diaphragm due to transient pressure spikes as discussed in Kavanaugh and Clarke (2001), but we cannot reconstruct with certainty whether this occurred in hole A1.
Phase 4 continued with the onset of oscillations in water pressure in hole A3, which occurred on 16 September. Similar oscillations were subsequently recorded in A1 and A2, starting on 27 September. Several prominent, in-phase oscillations occurred in boreholes A2 and A3 in early October (dotted vertical lines in Fig. 3). The amplitude of these oscillations differed between these two boreholes. Close inspection of the record from hole A1 reveals that this hole also experienced pressure maxima at the same time as A2 and A3 during this period, though these maxima are much less pronounced. The oscillations were not strictly periodic, but every cycle lasted multiple days. Each pressure oscillation consisted of a slowly rising limb, followed by a more rapid drop-off lasting 12 h or less. During this part of phase 4, the time-averaged water pressure rose slowly. Around 18 October, water pressure once more reached levels similar to those seen at the end of phase 3. There was a hiatus in the oscillations between 18 October and 8 November in all three boreholes. Pressure oscillations subsequently resumed in holes A2 and A3, with both holes producing very similar pressure records, but these oscillations have longer periods of around 12 days, as well as a higher mean pressure and a smaller amplitude than those recorded in October. December, with water pressure approaching overburden as in holes A4-A5. The pressure signal from A1 rose sharply to tens of metres head equivalent above overburden, starting on 30 November, which may again indicate instrument damage. It is instructive to compare these water pressure records with the surface temperature record and melt estimates ( Fig. 3a). Note that the temperature record was measured at the lower end of the study area, and is therefore likely to overestimate temperatures in the higher reaches of the glacier. The initial hiatus in diurnal cycling during phase 1 (24 to 29 July) occurs during a period of reduced surface melt and less-pronounced daily temperature variations, but temperatures remain above the melting point until the night of 25 July, after which significant diurnal temperature variations recommence. By contrast, the termination of strong diurnal cycling for holes A1-A5 occurred following the night from 4 to 5 August. Five August had positive maximum temperature, but this was followed by pronounced low temperatures around −4 • C on the night from 5 to 6 August. The albedometer record also indicates that snow fell during that night. That said, daytime highs subsequently continued to reach positive values similar to the period immediately preceding 6 August. The termination of strong diurnal cycling in hole B on 11 August corresponds to a longer period of sustained negative surface temperatures at night. There are several multi-day episodes of noticeable surface lowering as recorded by the SR50 sonic ranger at the weather station after 11 August, though these are of lesser magnitude than surface melt prior to 11 August. Surface temperatures reach positive values only on a few days in September, and no positive surface temperatures are recorded after 4 October. Measured surface melt is also minimal after the end of September. Additional data from the weather station show that the spike in surface lowering observed on 2 November coincides with high winds and low temperatures, and is therefore likely to correspond snow redistribution recorded by the SR50, rather than to actual melting.
The behaviour exhibited by the holes in the line A1-A5 during phases 1, 2 and 3 is qualitatively similar to that seen in some other records collected on South Glacier, except that we have found no clear evidence for two-day doublets elsewhere. Figure 4 shows a pressure time series recorded in another hole, marked C in Fig. 1, in 2010. A similar sequence of diurnal cycles (phase 1) followed by a drop-off (phase 2) and subsequent increase to pressures near overburden (phase 3) can be seen here, though phases 3 appears to include a long initial interval of very low and slowly rising water pressures. Note that the termination of diurnal cycling in this case corresponds straightforwardly to the onset of night-time freezing. Figure 5 shows data recorded for hole D in 2009, and plotted the same way as in Fig. 3. The summertime behaviour is similar to that seen in holes A1-A5 and hole C, but diurnal oscillations persisted in this borehole until mid-September, during a period of time in which night-time temperatures regularly fell below zero. Note that holes C and D are drilled into thinner ice (approximately 67 and 45 m, respectively) than holes A1-A6, and lie in areas of much steeper surface slope. Hole D also appears to be connected to an efficient drainage system in which water pressure often drops to nearatmospheric values. Consequently there is no analogue for phase 2 here, and, instead, only an increase in water pressure analogous to phase 3 is seen. Note that the records for holes C and D do not correspond to the same year as those for holes A1-A6, and the early summer of 2009 in particular was unusually warm at the field site, which may explain some of the qualitatively different behaviour of hole D. Figure 6a shows a record of diurnal cycles throughout an entire melt season at the hole marked E in Fig. 1. As in hole B, the amplitude of water pressure variations in the melt season was more limited than in holes A1-A5, C and D, and there was no obvious phase 2 or 3 following the termination of diurnal cycling. Instead, the water pressure simply remained close to its summertime mean. The diurnal pressure variations in hole E in late July and early August may be suggestive of a doublet structure, but this is not as pronounced as in holes A1-A3.
It is unclear whether the qualitative behaviour of phase 4 is replicated in any of the other time series recorded at South Glacier. The strongest indication of similar behaviour in another borehole comes from hole D during the 2009-2010 winter. During the equivalent of phase 3, the record from hole D includes two prolonged water pressure spikes in mid-to late October 2009 that have no equivalent in the A1-A3 record. These pressure spikes occurred immediately after a prolonged period of above-melting surface temperatures. The spikes were followed by period of relatively steady water pressure from late October 2009 to the end of December. Subsequently, hole D exhibited an abrupt drop in water pressure similar to holes A1-A3 on 8 January 2010, followed by an increase in water pressure and the subsequent onset of oscillations lasting about a week each. While these features are similar to phase 4 as seen in holes A1-A3 in 2011-2012, the record from hole D comes from a single instrument and is not replicated elsewhere, so we cannot rule out the possibility of the transducer having been corrupted.
An oscillatory behaviour similar to phase 4 is also evident in some boreholes during summer. Figure 6b shows an example from hole F. Water pressure initially underwent weak diurnal cycling between 13 and 21 July, followed by an interval of aperiodic pressure oscillations, each having a more gently rising limb followed by an abrupt pressure drop. These oscillations lasted between a few hours to two days, so considerably shorter than the wintertime oscillations of phase 4. It is unclear whether all these oscillations are due to the same process. Note that summertime oscillations similar to those in hole F were previously reported on Bench Glacier (see Fig. 2 in Fudge et al., 2005).

Phase 1
The strong diurnal pressure cycle during phase 1 indicates that surface melt reaches the bed with a delay of less than a day. Peak water pressures in holes A1-A6 are typically attained around midnight, lagging around 10 h behind peak air temperatures.
Holes A1-A3 show essentially the same signal during phase 1, indicating that they are connected to the same drainage system. They also exhibit very significant variations in water pressure, far more pronounced than hole B. There is a nearly constant offset between the measured water pressures in holes A1-A3 during phase 1. This difference in water pressure may simply account for different bed elevations at the hole bottoms, and correspond to negligible hydraulic gradients. We have no inclinometry data for the boreholes, and although we were able to measure borehole lengths and surface elevations, we are unable to locate the bottom of the boreholes exactly; the observed offsets between the pressure records are certainly consistent with the possibility that there are insignificant hydraulic gradients between the holes. In addition, there may be calibration errors due to damage to the transducers caused by pressure spikes (Kavanaugh and Clarke, 2001). This generally affects only the offset in the calibration but not the multiplier, so the measured pressure is the actual borehole pressure plus an offset that generally cannot be determined after the fact.
The occurrence of diurnal doublets is a striking feature of the record from holes A1-A3, especially from 1 to 10 August. This may be driven purely by the external forcing of the system, but there is no obvious doublet structure in the surface temperature record during the same period. It is therefore possible that the doublets are internal modes with a twoday period excited by external water supply. We expand on this further below and in Sect. 5, where we model their occurrence.
The behaviour of holes A4-A6 is consistent with the behaviour observed on Trapridge Glacier (Murray and Clarke, 1995;Stone and Clarke, 1996) as well as Haut Glacier d'Arolla (Hubbard et al., 1995) and Bench Glacier (Fudge et al., 2008). The sudden switch to showing the same pressure time series as A1-A3 on 4 August indicates that holes A4-A5 become connected to an expanding subglacial drainage system at that time. Meanwhile, the anti-phase behaviour observed in hole A6 can be explained through perturbations in stress in the ice caused by proximity to pressurized basal conduits as discussed in Murray and Clarke (1995). There is no obvious evidence of a diffusive connection between boreholes, unlike what was inferred for Haut Glacier d'Arolla in Hubbard et al. (1995).
Hole B is clearly not connected directly to the same part of the drainage system. It shows a qualitatively similar set of diurnal cycles, but with much smaller amplitudes that do not change over time in the same way as the connected holes in the A1-A6 line. This is consistent with the results of the upstream area calculation in Fig. 2, which indicate that hole B is at a significant distance from the main drainage axis straddled by A1-A6. Similar behaviour, with nearby boreholes displaying dissimilar but diurnally varying pressures, has been recorded at a number of other glaciers (e.g. Gordon et al., 1998;Fudge et al., 2005Fudge et al., , 2008. In addition to differences in their catchments, the apparently spontaneous generation of doublets in the A1-A3 record hints at an additional reason why holes A1-A6 and hole B show such marked differences in the amplitude of diurnal oscillations during phase 1, even though the time series are recorded simultaneously on the same glacier. Melt input to the drainage system could be similar for the catchments feeding into holes A1-A6 and hole B. Even then, however, the observed water pressure record is the nontrivial response of the drainage system to that input. The parts of the drainage system that run through A1-A6 and B are likely to have different modes that can be excited by surface melt forcing. The termination of strong diurnal pressure cycles in the A1-A6 line does not coincide with surface temperatures dropping below freezing in a sustained manner; instead it appears to be triggered by a single cold night. We hypothesize that this is sufficient to cause some englacial conduits to become blocked, preventing further efficient water supply to the bed. Due to the absence of moulins in or above the study area, it is likely that direct routing of melt from the glacier surface into the drainage system during phase 1 occurs primarily through crevasses, possibly with narrow cracks facilitating a more rapid closure than would be possible in the case of surface streams reaching the bed through moulins. By contrast with holes A1-A6, diurnal pressure cycles in hole B terminate several days later, during a period of several days of low air temperatures that reach barely above the melting point. The later termination suggests that water pressure in hole B may be controlled by water supply from a lower altitude, possibly from the crevasse field directly above the hole (Fig. 1), while A1-A6 lie along a possible major drainage axis that may be fed from higher elevations. Note also that a similar nonsimultaneous pattern of termination of diurnal cycling has previously been observed at Bench Glacier (Fudge et al., 2005), with surface meteorological data furnishing no obvious explanation for this behaviour.
The drop in water pressure during phase 2 can be understood as resulting from a drainage system whose capacity to discharge water suddenly exceeds water input into the system. A temporary lowering of pressure gradients is required to account for the reduced discharge through the system, leading to reduced water pressure (see for instance the seasonal cycle simulations in Schoof, 2010). The drop in water pressure then allows basal conduits to shrink in size, reducing the capacity of the drainage system.

Phase 3
The subsequent increase in water pressure in holes A1-A3 during phase 3 is also consistent with the conduits reequilibrating to a smaller size. As they do so, water pressure can rise again. Boreholes A4 and A5 show evidence of becoming disconnected from each other during this phase, as their pressures begin to evolve apparently independently, and approach overburden. We suggest that this is the result of these boreholes turning into isolated water pockets, in which changes in water pressure result from the response of the presumably slowly deforming pocket to normal stresses in the ice.
Boreholes A1-A3 remain highly correlated with each other during phase 3, which we take to be evidence that they remain physically well-connected. Water pressure reaches an almost constant level that remains significantly below overburden. This suggests that the boreholes may remain connected to a drainage system during phase 3, and that wa-ter throughput in this system is relatively steady. As we expect these boreholes to straddle a major drainage axis, continued drainage is a likely consequence of ongoing groundwater flow feeding a subglacial conduit, or of slow englacial drainage. Note that the 2009 water pressure record for hole D, which may also be near a major drainage axis, shows qualitatively similar behaviour, with water pressure reaching an almost constant level below overburden in November; the pressure spike seen in late October 2009 can possibly be attributed to surface melt during a warm spell immediately preceding the pressure spikes.

Phase 4
Phase 4 is unique to the boreholes A1-A3, at least among the time series recorded during the 2011-2012 winter. We stress that the pressure signal during this period appears to be robust. Aside from the sharp jumps in the record shown by A1 following the initial pressure drop, there is little reason to believe that the observed signal is the result of faulty instrumentation. Three different pressure transducers connected to two independent data loggers all show the same initial pressure drop, and exhibit several in-phase oscillations later on.
The temperature record clearly indicates that the oscillatory behaviour shown subsequent to the initial pressure drop is unlikely to be the result of surface melt, and neither the initial pressure drop nor the subsequent oscillations can be explained in a simple way by direct surface forcing. The lengthy period of steady water pressures at the end of period 3, during which positive surface temperatures are still reached on most days (22 August-5 September), also suggests that the pressure variations during phase 4 are unlikely to be the result of direct surface forcing. Nonetheless, as the initial pressure drop and some of the subsequent pressure oscillations are common to these three boreholes, which are located in what is likely to be the main drainage axis of the glacier, we surmise that they are the result of ongoing drainage activity. With no obvious external driver for these oscillations, we interpret them as being the result of an unforced, internal mechanism. Note that we also expect hole D to lie near a different major drainage axis, and it exhibits similar behaviour. Although the observations were made during a different year, this suggests a common physical process at work.
The most plausible explanation we are able to give is that the pressure oscillations are the result of a quasi-periodic build-up and release of subglacial or englacial water storage. This may be fed by ongoing, low-level but relatively steady water supply to the glacier bed from englacial sources, geothermal or frictional heating or, possibly more likely, from ground-water flow after the end of the melt season. Such water supply is known to occur for instance under Storglaciären (Stenborg, 1965), and continued englacial water supply has also been observed in Glacier de Tête Rousse (O. Gagliardini, personal communication, 2012 (Lingle and Fatland, 2003).
In the presence of a low-level water supply, a possible mechanism to explain spontaneous oscillatory behaviour could be a variation of Nye's (1976) jökulhlaup instability. Nye's instability, strictly speaking, requires a glacierdammed lake or other large water reservoir drained by a subglacial Röthlisberger-type channel. While mostly familiar as the mechanism by which a flood from a lake occurs through unstable enlargement of the channel, an entire jökulhlaup cycle of flooding and refilling can be seen as a pressure oscillation in the lake, where the oscillations do not require any temporal variations in water supply (see also Fowler, 1999;Kingslake and Ng, 2013).
Radar surveys have not shown any obvious evidence for large water bodies within South Glacier. Although such water bodies could exist in the areas between recorded radar tracks, we consider it more likely that a substantial water reservoir in the present case would have to be distributed, consisting of water stored in basal crevasses (Harper et al., 2010) or other smaller water pockets. We present a set of idealized numerical model calculations in the next section to show that such distributed storage can indeed drive oscillatory drainage. We also show that oscillatory behaviour in that case does not necessarily require a Röthlisberger-type channel but is also possible in a linked-cavity drainage system, and contrast this with the case of localized storage in a larger water body, for which oscillations require a channel-like conduit. Given that internally driven oscillations can occur, an obvious question is then why these oscillations should only commence some time after the melt season. To address this, we also show that the onset of oscillatory drainage corresponds to water supply to the conduit dropping below a critical value.
In addition, we also address the observed doublet structure in diurnal pressure peaks using the same model. We demonstrate that such doublets, with a two-day period, can result from the nonlinear coupling between water storage at or near the bed and a surface water input that varies purely diurnally. Similar ideas have recently been developed by Kingslake (2013, chapter 3). We show that a doublet response can be generated with water input to the system remaining above the upper threshold for unforced oscillatory drainage. The modelled doublets are therefore not directly the result of the instability mechanism that we propose for wintertime oscillations in water pressure, but are linked to it through the role played by basal water storage.

The model
To illustrate the possibility of initiating pressure oscillations due to mechanisms purely internal to the drainage system and not driven directly by variability in surface melt, we use  (2010), building on similar models used elsewhere (e.g. Ng, 1998;Hewitt and Fowler, 2008;Schuler and Fischer, 2009;Hewitt et al., 2012). Denoting conduit cross section by S(x, t), effective pressure by N (x, t) and discharge by Q(x, t), where x is downstream distance and t is time, we put Here c 1 , c 2 , c 3 , n, α, v o , S 0 and v p are positive constants with α > 1 (see Appendix A). 0 is a geometrically determined background gradient, given by Note that effective pressure N is linked to water pressures p w , which is the primary observable in the field, through where p i is overburden (or, more precisely, spatially averaged normal stress in the ice at the bed). Physically, the model is based on a conduit whose size evolves due to a combination of dissipation-driven wall melt at a rate c 1 Q , opening due to ice sliding over bed roughness at rate v o (1 − S/S 0 ) (where S 0 is a cut-off due to bed roughness being drowned out; see Schoof et al., 2012) and creep closure at rate c 2 SN n , with discharge in the conduit given by a Darcy-Weisbach or Manning friction law through Eq. (2b). Water in the drainage system is stored in void space connected to the conduit; the filling level of this void space depends on water pressure, and hence on N. v p is a measure of the storage capacity of these voids. For instance, if there are parallel-sided cracks in basal ice connected to the conduit, and these cracks have base area a 0 per unit length of the flow path, then stored water volume per unit length of the flow path is v = a 0 h w = a 0 p w /(ρ w g), where h w is the filling level in the crevasses. In terms of effective pressure, water storage per unit length of the flow path becomes v = a 0 (p i − N )/ρ w g if overburden p i can be treated as constant. Then dv/dt = −a 0 /(ρ w g) dN/dt, which becomes the first term in Eq. (2d) if v p = a 0 /(ρ w g). m is a water supply rate per unit length of the conduit.
As boundary conditions, we impose zero flux Q = 0 at the upstream end of the domain x = 0, and zero effective pressure N = 0 at the glacier snout x = L 0 . The latter corresponds to having both zero water pressure and zero overburden at the snout. For computational tractability and in keeping with many other similar drainage models (e.g. Werder et al., 2013), we do not impose the upper and lower bounds on water pressure considered in Schoof et al. (2012) and Hewitt et al. (2012).
Note that similar results are obtained if we set distributed water supply m to zero and instead input an equivalent flux Q at the upstream end of the domain, but we focus on the distributed supply case here for simplicity. The computations were carried out with the parameter values in Table 1. Numerically, solutions are computed by first discretizing the equations above to yield the one-dimensional version of the discrete model in Schoof (2010) (see also Appendix B), and time-stepping is done by means of a backward Euler step.

Instability and oscillatory drainage
To gain an understanding of possible self-sustained oscillations in the system, we first impose a fixed, spatially uniform water supply rate m. We compute the corresponding steady state. This can be done by recognizing that (see also Hewitt et al., 2012) Q(x) = x 0 m(x ) dx = mx in steady state, and rewriting the steady-state version of Eq. (2) as a first-order nonlinear differential equation for N. This is done by solving for S in Eq. (2a), Together with Eq. (2b), Eq. (4a) defines implicitly as a function of the known discharge Q(x) as well as the unknown effective pressure N(x): Solving this for (N, Q), we then have from Eq. (2c) which can be integrated backwards from the glacier snout, where N = 0. To test whether the steady state is stable to . Both panels show mean effective pressure over the domain in steady state plotted against net water input mL into the system. Rhombus-shaped markers correspond to cavity-like conduits (see main text), while circles correspond to channel-like conduits. An empty marker corresponds to a steady state that is unstable in the presence of distributed storage; a solid marker corresponds to a stable solution. The two panels differ in the amount of englacial storage: (a) has v p = 5.1 × 10 −7 m 2 Pa −1 ; (b) has v p = 5.1 × 10 −6 m 2 Pa −1 . Note that there are unstable cavity-like solutions. The large markers in (a) and (b) correspond to the initial conditions in Figs. 8 and 9, respectively.
perturbations or whether it spontaneously evolves into an oscillatory solution when perturbed, we perform a numerical linear stability analysis of the discretized model about the steady-state solution (see Appendix B). In addition, we also perturb the steady-state solution by adding a small amount of noise to S and N, and subsequently compute the evolution of the solution to the full nonlinear model to confirm the results of the stability analysis and understand the evolution of S and N once perturbations have grown to finite size. Our results are shown in Fig. 7. Here we plot the mean effective pressure N in the steady state as a function of total water supply rate mL. A solid symbol indicates a stable steady state, while an empty symbol denotes an unstable steady state. Invariably, the onset instability corresponds to at least one pair of complex conjugate eigenvalues with positive real part in the linear stability analysis, corresponding to an oscillatory instability, and this is confirmed by direct numerical solution of the full nonlinear model. We see that instability occurs for an intermediate range of water supply rates mL, with mL lying between a lower and an upper critical value. This unstable range is shifted to lower water supply rates for smaller storage capacities v p . Two sample time series of mean water pressure over the domain for an unstable solution evolving away from the steady state are given in Figs. 8 and 9, along with several snapshots of N as a function of position during the ensuing limit-cycle pressure oscillations. shows snapshots of effective pressure against position x along the flow path during one limit cycle. Note that the pressure oscillations occur over the entire length of the flow path.

Distributed versus localized water storage
For the steady states shown in Fig. 7, we also tested whether they are overall channel-or cavity-like in the sense of Schoof (2010) by placing two of these conduits side by side and permitting them to exchange water freely, so that both conduits are at the same effective pressure. Each conduit, labelled by i = 1 and 2, then satisfies Eqs. (2a) and (2b) separately, with S and Q replaced by S i and Q i . Equation (2c) still holds for the single effective pressure N and hydraulic gradient in the system. In order to suppress storage-driven instabilities, we do not allow water storage, so that Eq. (2d) is replaced by In order to classify the conduit overall as channel-like we tested whether the system of two parallel steady-state conduits is stable to perturbations. This was again done by performing a numerical linear stability analysis on the discretized model, confirmed by direct solution of the full nonlinear model starting from a slightly perturbed steady state (through adding a small random amount of noise to S 1 and S 2 ). In each instance of instability, the eigenvalues are now purely real and positive, and the result of solving the corresponding full nonlinear model is channelization with one conduit growing at the expense of the other. We classify the steady state as cavity-like if it remains stable.
As in Schoof (2010), the channelizing instability is purely the result of competition between neighbouring conduits and is not oscillatory but results in a new steady state. By contrast, the oscillations shown in Figs. 8 and 9 and indicated by empty symbols in Fig. 7 result from the interaction between a single conduit and distributed water storage at the bed. Note however that the definition of cavity-and channellike behaviour used above differs slightly from that proposed in Schoof (2010); the analytical criterion developed there is best applied in idealized situations where discharge Q and hydraulic gradient are constant along the conduit, which is not the case here.
Based on the above definition, we find that the change from cavity-to channel-like in our conduit occurs between mL = 10 3 m 3 day −1 and mL = 2 × 10 3 m 3 day −1 for the parameter values in Table 1, while oscillatory behaviour persists down to mL = 250 m 3 day −1 (note that the transition from channel-like to cavity-like corresponds to the change from circular to rhombus-shaped markers, while oscillations occur for empty markers). This contrasts with the behaviour of an analogous drainage system in which distributed storage is replaced by a single reservoir at the upstream end of the conduit: the results in Schoof (2010, Supplement) indicate that a jökulhlaup instability is then only expected to occur for a channel-like conduit.
To illustrate this further, we also solve Eq. (2) with v p and m set to zero, but with the upstream boundary condition Q = 0 at x = 0 replaced by where Q 0 is a fixed water input into the reservoir at the head of the channel, which we assume to have a base area A and vertical sides. In that case, V p can be identified analogously to v p as A/(ρ w g). These modifications turn our model into a more standard jökulhlaup-type model (e.g. Fowler, 1999;Kingslake and Ng, 2013), although without a "seal" in the sense of Fowler (1999) but with the ability to describe the transition of the conduit from cavity-like to channel-like. To create a roughly equivalent set-up to that used for distributed storage, we set Q 0 = mL and V p = v p L for the values of m and v p used to simulate flow with distributed drainage. Figure 10 shows steady-state results for mean effective pressure over the domain against water supply Q 0 . The solutions shown are slightly different from those in Fig. 7 because water input is localized at the upstream end of the conduit rather than being distributed along its length. Rhombusshaped and circular markers once more depict cavity-like and Q 0 (10 3 m 3 / day) Figure 10. Mean effective pressure over the domain in steady state plotted against water input Q 0 into the system with localized upstream water input and storage. Rhombus-shaped markers correspond to cavity-like conduits (see main text), while circles correspond to channel-like conduits. An empty marker corresponds to a steady state that is unstable in the presence of storage, with V p = 2.55×10 −3 m 3 Pa −1 . The large marker in panels corresponds to the initial condition in Fig. 11. channel-like steady states, respectively, with empty markers corresponding to states that are unstable to perturbations with V p = 2.55 × 10 −3 m 3 Pa −1 . This storage capacity is equivalent to the distributed storage capacity v p = 5.1 × 10 −7 m 2 Pa −1 used in Fig. 7a. The main difference between the case of distributed and localized (lake-like) storage is clearly that the localized storage case only leads to instability when the conduit is channel-like (there are no empty rhombus-shaped markers in Fig. 10). Figure 11 shows a sample time series for the nonlinear evolution of the localized storage system into a limit cycle solution; note that the steady-state solution for the chosen parameter values is only marginally unstable. The growth rate in amplitude is slow, which explains why the system only approaches the limit cycle very slowly.

Water pressure doublets
Finally, we show by way of example that the drainage model described above in Eq.
(2) with distributed storage and water input can generate a doublet structure analogous to that seen in Fig. 3 during summer. We force the model with diurnally varying water input, of the form and choose parameter values that ensure that m always remains above the threshold for the onset of the self-sustaining pressure oscillations that we have modelled above. Figure 12 shows results for a set of parameter choices that generates a response with a one-day period with v p = 0 but has a two-day period with sufficient water storage capacity.

Discussion
There are four important qualitative results from our model: first, distributed storage is capable of generating spontaneous oscillations in a subglacial drainage system through a positive feedback mechanism, as is also true for similar drainage  Figure 11. Evolution of the instability for localized water input and storage, with Q 0 = 4 × 10 3 m 3 day −1 , V p = 2.55 × 10 −3 m 3 Pa −1 . (a) shows mean effective pressure over the domain plotted against time. The solution evolves very slowly into a limit cycle, with a period around 20 days. (b) shows snapshots of effective pressure against position x along the flow path during one limit cycle. Note that these effective pressure profiles are qualitatively similar to those for distributed storage with a channel-like conduit as shown in Fig. 8. systems fed by an upstream reservoir (Nye, 1976;Fowler, 1999). Second, the onset of these oscillations occurs only below an upper critical water supply rate: at sufficiently high discharge rates through the conduit, the effect of water storage is insufficient to generate self-sustaining pressure oscillations. This critical water supply rate depends on the amount of distributed storage available, and increases with v p . Third, oscillations occur only above a lower critical value for water supply rate. For distributed water storage, that lower threshold for oscillatory behaviour is small enough that the conduit can be cavity-like in the sense of Schoof (2010). By contrast, Nye's instability for drainage fed by a single reservoir only occurs in channel-like conduits. Lastly, even if water input to the system is large enough for it not to exhibit spontaneous oscillations, we find that forcing the system through water supply at a single diurnal frequency can give rise to a doublet structure with a two-day period in the pressure response. This results from storage capacity of the drainage system retaining excess water injected during the first day of the cycle and releasing it on the second. The first three qualitative observations made above can be explained in greater theoretical detail, which we will do elsewhere (paper in preparation).
The most relevant of these results to our field data are that pressure oscillations can be driven by distributed, englacial or subglacial water storage such as cracks in basal ice, and that the onset of oscillations corresponds to water supply mL dropping below a critical threshold. This is consistent with the observed onset of oscillations on South Glacier after the end of the melt season, when water supply rates are likely to decrease as residual water stored englacially or as ground water drains away.
There are some noticeable discrepancies between the model results above and the observed time series that we do not wish to trivialize. While our simple model is able to predict the onset of oscillations at reduced water supply rates, it does not account for the abrupt drop-off in water pressure  6), m 0 = 6.4 × 10 4 m 3 day −1 , m = 8 × 10 3 m 3 day −1 , ω = 2π day −1 . (a) shows mean effective pressure over the domain plotted against time for v p = 0; (b) shows the same but for v p = 5.1 × 10 −7 m 2 Pa −1 . These parameter values ensure that water input does not drop below the threshold for instability depicted in Fig. 7 at any point in the cycle. Note that the amplitude of pressure oscillations is unrealistically large in both panels, though smaller when there is water storage. A damping effect due to water pressure is to be expected as variations in discharge along the channel can be reduced by allowing water to be stored temporarily. In turn, this requires smaller variations in pressure gradient and hence in pressure at the bed during the cycle.
at the start of phase 4, nor does it account for the apparently equally abrupt onset of oscillations shown in Fig. 3. Instead, Figs. 8 and 9 show an initial period of exponential growth in oscillation amplitude predicted by the model. Another feature of the observed oscillations is a lengthening of the period of oscillation, accompanied by a drop in oscillation amplitude and increase in mean as well as minimum and maximum water pressure during each oscillation cycle. The lengthening of the oscillation cycle may conceivably result from a decreasing water supply, which is known to increase the interval between jökulhlaups in drainage systems with a single reservoir (e.g. Fowler, 1999). The concurrent rise in mean water pressure and reduced oscillation amplitude are harder to explain, and we have thus far not been able to replicate this behaviour conclusively with our simple model (although the onset of oscillations can lead to a drop in mean effective pressure relative to the unstable steady state; see Fig. 8).
Our simulations have been based on a particular set of parameter choices that are motivated by the geometry of our field site. Ultimately, many of the parameters in the model are not well constrained, and our flowline model is highly idealized. It would be unwise to attach too much importance to the precise range of unstable water supplies in Fig. 7 as our model results are parameter-dependent and Fig. 7 only explores changes in a a single parameter, namely the storage capacity v p . If we do take the unstable water supply rates at face value, we find values in Fig. 7a in the range 200-30 000 m 3 day −1 . The upper end of that range is very unlikely to be realized under wintertime conditions, while the lower end is plausible for wintertime discharge if there is ground-water discharge under the glacier. The borehole line A1-A6 has a total upstream area of about 0.5 km 3 for f = 0 and f = 0.5 in figure 2 if we include only glacier-covered areas, and about twice that if ice-free valley slopes are included. Assuming a geothermal heat flux of 0.04 Wm −2 and with measured ice velocities around 10 m a −1 (De Paoli and Flowers, 2009;Flowers et al., 2014) and driving stresses around 10 5 Pa, we can estimate basal melt rates of around 2 × 10 −5 m day −1 , giving an integrated water supply of 10 3 m 3 day −1 , below the unstable range. To explain higher supply rates in excess of 200 m 3 day −1 requires ongoing ground-water drainage, with area-averaged discharge rates (and therefore long-term recharge rates) around 2×10 −4 m day −1 , or 0.07 m a −1 . This amounts to about 10 % of the glacier-wide annual accumulation rate estimated in Wheler et al. (2014), and is at least not grossly unrealistic.
The model assumes that the necessary water storage consists primarily of cracks or other storage bodies that can fill and drain easily with water, and in which filling level, and therefore storage, is an increasing function of water pressure at the bed. It is possible but not very likely that the boreholes themselves constitute significant storage volume, in which case the act of observation could interfere with the observed system. Two reasons make this unlikely. The first is that the boreholes freeze shut relatively within a few days of being drilled. Limited borehole thermometry evidence (Wilson, 2012) indicates that the top 50 m of one borehole (A6) reached temperatures below the melting point within 25 days. Second, even if they remain open, the boreholes have extremely small volume storage. With a cross section of about 0.01 m 2 , each borehole would give a storage capacity of V p = 1.6 × 10 −3 m 3 Pa −1 in the model for lakelike drainage above, three orders of magnitude smaller than the storage volumes we have used in our simulations. With this limited storage capacity, a variation in hydraulic head of 10 m typical of the later pressure oscillations in holes A2-A3 would correspond to a water volume of 0.1 m 3 stored in or discharged from each borehole over an approximately 12-day period. This contrasts with the approximately 200 m 3 we expect to be produced over the same time period by friction and geothermal heating over the upstream area for the borehole line A1-A6. The storage capacity of boreholes is therefore unlikely to be the cause of the observed oscillations if their drainage patterns indicated in Fig. 2 have any bearing on wintertime drainage. Further, boreholes by themselves constitute localized storage, and the proposed instability mechanism would require them to connect to a channelized drainage system in order to account for spontaneously generated oscillations. At the low discharge rates we have just estimated, it is unlikely that this would be the case. season in a group of boreholes that are likely to straddle a major subglacial drainage axis. We interpret these as the result of a remnant drainage system that stays active due to a background input of water not derived directly from the glacier surface. This input could take the form of ground water, slow englacial drainage through the temperate ice that exists in the upper glacier (Wilson et al., 2013), or melt generated by geothermal heat or frictional dissipation. We have also shown that distributed water storage connected to a subglacial drainage system can cause an instability leading spontaneously to pressure oscillations at sufficiently low but steady water input. The instability is similar to that described in Nye (1976) for concentrated water storage in the form of a lake, but importantly can occur even when water supply is low enough for drainage to occur purely through cavity-like conduits. We suggest this mechanism as a possible explanation for the pressure oscillations observed during winter at South Glacier. Water storage connected to the subglacial drainage system may also account for a nontrivial drainage response to diurnal forcing, with a two-day periodicity resulting from a forcing that has a purely diurnal period. This has been explored previously by Kingslake (2013).

Conclusions
It is tempting to link the observed drainage phenomenon with the dynamics of the glacier, which has a history of surging. In particular, water pressure spikes during the pressure oscillations might be expected to facilitate faster sliding. However, these oscillations only occur in a spatially confined area, and we have at present no evidence that they have an effect on the flow of the glacier. Further work, both observational and modelling, is required to understand both the physics behind the observed drainage oscillations and their effect on glacier dynamics. The discretized version of the single conduit model solves for conduit size S i between nodes labelled i and i + 1 at which effective pressures N i and N i+1 are solved for. With a uniform node spacing x and the i = 1 node marking the inflow boundary and i = n 0 being the glacier snout, these are determined through where Denoting steady-state solutions by a superscript 0, S 0 i and N 0 i can be computed from Eq. (B1) analogously to Eq. (4a). We have Q 0 i = i j =1 m i x from Eq. (B2d). In addition, Eq. (B2a) becomes which together with Eqs. (B2b) and (B2c) gives an implicit relationship between N 0 i and N 0 i+1 : This can be solved recursively for N 0 i , starting with i = n 0 − 1, in which case N 0 i+1 = N 0 n 0 = 0. The numerical linear stability analysis linearizes the dynamical system Eq. (B1) about the steady state, solving for exponentially growing or shrinking solutions N i = N 0 i + N i exp(λt), and S i = S 0 i + S i exp(λt). λ is then the solution of the generalized eigenvalue problem