Antarctic subglacial lakes drain through sediment-floored canals : Theory and model testing on real and idealized domains

Over the past decade, satellite observations of ice surface height have revealed that active subglacial lake systems are widespread under the Antarctic ice sheet, including the ice streams. For some of these systems, additional observations of ice stream motion have shown that lake activity can affect ice-stream dynamics. Despite all this new information, we still have insufficient understanding of the lake-drainage process to incorporate it into ice sheet models. Process models for drainage of ice-dammed lakes based on conventional “R-channels” incised into the base of the ice through melting are unable to reproduce 5 the timing and magnitude of drainage from Antarctic subglacial lakes estimated from satellite altimetry given the low hydraulic gradients along which such lakes drain. We have developed an alternative process model, in which channels are mechanically eroded into the underlying deformable subglacial sediment. When applied to the known active lakes of the Whillans/Mercer ice stream system, the model successfully reproduced both the inferred magnitudes and recurrence intervals of lake volume changes, derived from Ice, Cloud and land Elevation Satellite (ICESat) laser altimeter data for the period 2003–2009. Water 10 pressures in our model changed as the flood evolved: during drainage, water pressures initially increased as water flowed out of the lake primarily via a distributed system, then decreased as the channelized system grew, establishing a pressure gradient that drew water away from the distributed system. This evolution of the drainage system can result in the observed internal variability of ice flow over time. If we are correct that active subglacial lakes drain through canals in the sediment, this mechanism also implies that active lakes are typically located in regions underlain by thick subglacial sediment, which may 15 explain why they are not readily observed using radio-echo sounding techniques.

Abstract.Over the past decade, satellite observations of ice surface height have revealed that active subglacial lake systems are widespread under the Antarctic Ice Sheet, including the ice streams.For some of these systems, additional observations of ice-stream motion have shown that lake activity can affect ice-stream dynamics.Despite all this new information, we still have insufficient understanding of the lake-drainage process to incorporate it into ice-sheet models.Process models for drainage of ice-dammed lakes based on conventional "R-channels" incised into the base of the ice through melting are unable to reproduce the timing and magnitude of drainage from Antarctic subglacial lakes estimated from satellite altimetry given the low hydraulic gradients along which such lakes drain.We have developed an alternative process model, in which channels are mechanically eroded into the underlying deformable subglacial sediment.When applied to the known active lakes of the Whillans-Mercer ice-stream system, the model successfully reproduced both the inferred magnitudes and recurrence intervals of lake-volume changes, derived from Ice, Cloud and land Elevation Satellite (ICESat) laser altimeter data for the period 2003-2009.Water pressures in our model changed as the flood evolved: during drainage, water pressures initially increased as water flowed out of the lake primarily via a distributed system, then decreased as the channelized system grew, establishing a pressure gradient that drew water away from the distributed system.This evolution of the drainage system can result in the observed internal variability of ice flow over time.If we are correct that active subglacial lakes drain through canals in the sediment, this mechanism also implies that active lakes are typically located in regions underlain by thick subglacial sediment, which may explain why they are not readily observed using radio-echo-sounding techniques.

Introduction
Since the initial observation of "large, flat, circular basins" in the ice surface of Antarctica by Russian pilots during the International Geophysical Year (Robinson, 1964), there has been significant interest about the role of lakes within the larger ice-sheet system.Beginning in 1972, radio-echo sounding (RES) began confirming that these surface features reflect storage of free water at the ice-sheet base (Oswald and Robin, 1973) and subsequently continued to be a primary technique to identify subglacial lakes, including Lake Vostok, one of the largest freshwater lakes in the world (Ridley et al., 1993;Kapitsa et al., 1996).With increased availability of RES data, and as our ability to observe the ice surface precisely has improved, the number of known subglacial lakes in Antarctica has increased.At the time of writing, over 300 subglacial lakes have been discovered throughout the continent using a variety of geophysical methods (Wright and Siegert, 2012).Until the mid-2000s RES was the primary technique for identifying subglacial lakes (e.g., Siegert et al., 2005;Carter et al., 2007).Most of the lakes found in RES surveys tended to be located beneath the slowmoving ice near the divides (Fig. 1a), driving initial research questions on whether lakes were open or closed systems (e.g., Bell et al., 2002;Tikku et al., 2005), with considerable speculation about their impact on local ice dynamics (e.g., Dowdeswell and Siegert, 1999;Bell et al., 2007;Thoma et al., 2012).(a) Antarctic subglacial lakes inventory showing detected "active" lakes, satellite-detected "big" lakes, and lakes detected by radio-echo-sounding (RES) (Wright and Siegert, 2012).Background grayscale corresponds to ice velocities from Rignot et al. (2011).(b) Zoomed-in map for Whillans-Mercer ice streams showing the lakes (SLM: subglacial Lake Mercer; SLC: subglacial Lake Conway; SLW: subglacial Lake Whillans; SLE: subglacial Lake Engelhardt; L78: Lake 7-8) and flow paths used in this study (see black box in panel (a) for location) Since 2005, a variety of repeat observations of the ice surface from satellite missions have revealed patterns of surface uplift and subsidence consistent with the filling and draining of subglacial water bodies (e.g., Gray et al., 2005;Wingham et al., 2006;Fricker et al., 2007).In contrast to "RES lakes", these "active lakes" have been found beneath fastflowing ice streams and outlet glaciers (Fig. 1a, b).Many of the active lakes that have been surveyed with RES (e.g., Christianson et al., 2012;Siegert et al., 2014;Wright et al., 2014) lacked the characteristic basal reflections traditionally used to identify the presence of subglacial lakes (hydraulic flatness, specularity, and brightness relative to surroundings; see Carter et al., 2007).Multiple hypotheses have been proposed to explain this discrepancy, mostly relating to unconstrained hydrodynamics of the system or data processing artifacts (Siegert et al., 2015).
Active lakes are of particular glaciological interest due to their potential impact on the dynamics of fast-flowing ice streams.Surface altimetry observations suggest that these lakes can hold back and then episodically release large volumes of water into the subglacial environment (e.g., Fricker et al., 2007;Smith et al., 2009;Carter et al., 2011), which could alter ice-stream dynamics for hundreds of kilometers downstream.Repeat measurements of ice velocity coincident with the existing record of lake dynamics are sparse, with only three published instances: Byrd Glacier, East Antarctica (Stearns et al., 2008); Crane Glacier, Antarctic Peninsula (Scambos et al., 2011); and Whillans Ice Stream, West Antarctica (Siegfried et al., 2016).In all three cases, subglacial lake activity inferred from surface height anomalies correlated to episodes of temporary ice acceleration, but understanding of the longer-term impact of lake drainage on ice dynamics is significantly limited by a short (maximum 12-year) observational window (Siegfried et al., 2016).
Critical to resolving the link between ice dynamics and lake activity is determination of the mechanism by which lake drainage occurs.While some ice-sheet models have started to incorporate primitive elements of subglacial lake dynamics (e.g., Goeller et al., 2013;Livingstone et al., 2013), they still do not have a realistic treatment of observed lake-drainage processes.Most process-based treatments of Antarctic subglacial lakes (e.g., Wingham et al., 2006;Evatt et al., 2006;Carter et al., 2009;Peters et al., 2009) have hypothesized that active lakes drain via a mechanism similar to that of ice-dammed lakes in temperate glacial environments, in which narrow, semi-circular conduits are melted into the basal ice ("R-channels").More recently, however, several complementary lines of evidence have called into question the ability of an R-channel to form (Hooke and Fastook, 2007) and close (Fowler, 2009) in subglacial conditions typical of Antarctica.Fowler (2009) suggested that channels incised into the underlying sediment may be the preferred mechanism.
In this paper we describe the development of a new model for the filling and drainage of Antarctic subglacial lakes on ideal and real domains based upon several well-studied lakes in the hydrological system of Mercer and Whillans ice streams, West Antarctica (Fricker et al., 2007;Fricker and Scambos, 2009;Siegfried et al., 2014;Tulaczyk et al., 2014;Siegfried et al., 2016).Lake drainage occurs in our model through channels in the underlying sediment (termed "canals" by Ng, 2000) and we compare the output lake volumes with those from an existing R-channel model (following Kingslake and Ng, 2013).We aim to (i) develop a lakedrainage model that reproduces the recurrence interval and inferred volume ranges for lake-drainage events on the Siple Coast, (ii) provide better context for ongoing observations of lake-volume change and lake distribution with respect to the subglacial hydrology, and (iii) move towards a consistent parameterization of subglacial lake activity in ice-sheet models.
This paper begins with background on existing subglacial water models, both for Antarctica and for ice-dammed lakes in temperate environments (Sect.2).Next, we describe the theory of the models used in this paper, including background hydropotential theory, our implementations of R-channel and canal theory, and our coupling between a background distributed water system and a channelized system.We also detail the domains over which we apply the models (Sect.3).We then highlight the results from our experiments using the R-channel and canal models, compare these results to published observations of Antarctic subglacial drainage events, and explore the sensitivity of our model to various physical parameters required as model inputs (Sect.4).We discuss the implications of these model sensitivities, suggest reinterpretations for existing observations of subglacial hydrology in Antarctica based on the viability of a canal mechanism for lake drainage, and consider issues related to the inclusion of realistic models of subglacial lake drainage into large-scale ice-sheet models (Sect.5).Finally, we summarize our findings (Sect.6).
2 Basal water models and subglacial lake drainage

Antarctic basal water models
Models for subglacial water transport and distribution include at least one of three processes: distributed sheet flow, groundwater, and channelized flow.The simplest and most common models for ice-sheet basal water flow invoke some form of distributed system that spreads water laterally (e.g., Le Brocq et al., 2009).In such systems, water pressure increases with water flux, while basal traction decreases.More sophisticated models, however, prefer to accommodate sliding by deformation of the subglacial sediment, making basal traction decrease through increasing sediment porosity (e.g., Tulaczyk et al., 2000).Given that subglacial sediment is widely understood to lack the transmissivity necessary to accommodate the water fluxes at the base of the Antarctic Ice Sheet (Alley, 1989), changes in sediment water storage have been regulated by exchange with a distributed system (e.g., Christoffersen et al., 2014;Bougamont et al., 2015).These more sophisticated distributed-groundwater exchange models show the most consistency with borehole (Engelhardt and Kamb, 1997;Christner et al., 2014;Tulaczyk et al., 2014) and seismic (Blankenship et al., 1987) observations of the basal environment.
Channelized systems are those in which water flux is concentrated in one or more discrete conduits.The type of conduit that is most commonly modeled is referred to as an "Rchannel", which is eroded thermally into the ice by turbulent heat generated by water moving down a hydraulic gradient (Röthlisberger, 1972;Nye, 1976).As the relative area of the ice-bed interface occupied by these systems is small, they can support lower water pressures.R-channels have been well studied in Greenland, where they have been directly observed at the ice-sheet margin, and are understood to be supplied primarily by surface water, much of which enters at discrete recharge points known as moulins (e.g., Zwally et al., 2002).Basal meltwater systems in Antarctica differ in two significant ways from those in Greenland: they are isolated from the atmosphere and so do not receive surface meltwater; and they are located in regions of low hydraulic slopes, such that the heat generated by water moving down gradient is likely not sufficient to erode an R-channel (Alley, 1989).Therefore, channelization has typically not been considered in large-scale basal water models for Antarctica (e.g., Le Brocq et al., 2009;Bougamont et al., 2011).In the last decade, however, increased consideration has been given to the role of channelization in the drainage of subglacial lakes (e.g., Evatt et al., 2006;Carter et al., 2009) and to the possibility of subglacial channels incised in the sediments instead of the overlying ice (e.g., van der Wel et al., 2013;Kyrke-Smith and Fowler, 2014).

Lake-drainage theory
Most of our understanding of the drainage of subglacial lakes derives from the Nye (1976) model for the drainage of icedammed lakes on temperate glaciers in alpine environments (e.g., Clarke, 2003;Werder et al., 2013), where floods descend hundreds of meters over 1's to tens of kilometers on timescales of days, and channelization is well documented (e.g., Roberts, 2005).Fowler (1999) explained the repeated drainage in these systems with the following model (see also Fig. 2): (i) when the lake is at low stand, water is trapped behind a local maximum in hydropotential, known as "the seal"; (ii) as lake water levels rise, a hydraulic connection forms over the seal initiating thermal erosion by outflow from the lake; (iii) during the early stages of lake drainage, the potential gradient is relatively steep and effective pressure is low causing melt to exceed creep closure; (iv) with ongoing drainage, the level of the lake lowers causing a decrease in the hydraulic gradient, while the effective pressure at the seal increases, allowing the channel to continue to siphon outflow despite the lower lake level.The effective pressure change also increases the rate of creep closure of the channel by the overlying ice, while the energy available for thermal erosion decreases simultaneously.(v) As the channel closure rate increases, a reduction in effective pressure ultimately reforms a hydraulic seal between the lake and points downstream.
Antarctic subglacial floods occur on larger spatial and longer temporal scales than alpine subglacial floods: water typically descends tens of meters over hundreds of kilometers and drainage sometimes persists for multiple years (Wingham et al., 2006;Fricker et al., 2007)  (a) initially the lake is filling and no outflow takes place; (b) as lake level increases, effective pressure differences between lake and seal allow for some water to escape via a distributed system (a.k.a."water sheet"); (c) as the distributed system grows, a low pressure channel is eroded that begins to draw water from the surroundings; (d) while the low pressure channel allows for the lake to drain below the level at which outflow was initiated, as pressure lowers, sediment creeps inward; and (e) the channel contracts until its flow is negligible and the lake refills.
dence from borehole observations of the basal environment along major flow paths connecting lakes indicates that water travels via distributed systems (e.g., Engelhardt and Kamb, 1997).Although the patterns of lake volume and outflow over time are qualitatively similar to those observed during the drainage of alpine ice-dammed lakes, both historical and more recent work has shown that R-channels are unlikely to play a major role in the Antarctic subglacial water system due to several issues: 1.The inability of water to melt a channel on an adverse bed slope: a substantial number of flow paths that drain known subglacial lakes appear to exist on an adverse bed slope from thicker ice into thinner ice, where the pressure melting point is higher.Alley et al. (1998) showed that once the basal slope exceeds 1.2-1.7 times the surface slope, the heat generated by turbulent dissipation is insufficient to maintain the water at the pres-sure melting point and therefore cannot melt surrounding ice to form an R-channel.
2. Thermal issues related to polar ice: most models for Rchannel formation imply that the surrounding ice is temperate and isothermal.In polar ice, where a significant temperature gradient is likely to be present in the ice immediately above the bed, more turbulent heat will be required to melt a given amount of ice than would be needed for temperate ice (Hooke and Fastook, 2007).
3. The slow rate of closure predicted for R-channels at pressures observed under Antarctica: creep closure of R-channels requires a drop in water pressure on the order of several tens of meters water equivalent, 5-10 times higher than the surface drawdown observed during the drainage of most Antarctic subglacial lakes (Fowler, 2009).Field reports suggest that subglacial Lake Whillans (SLW), for example, contained water even at low stand (Christianson et al., 2012;Horgan et al., 2012;Christner et al., 2014), suggesting that the channel may have closed before the lake was completely drained.
In our work we directly compared output from a model for lake drainage via an R-channel (Kingslake and Ng, 2013) with a model in which channels are formed via mechanical erosion of underlying sediment (canals) on a domain with geometry similar to that found for flow paths draining Antarctic subglacial lakes.By replacing a channel incised into the overlying ice with one incised into the sediment we address issues (1) and (2), as the erosion of sediment is less sensitive to temperature than the erosion of ice.This change also addresses (3) because the deformability of sediment is more sensitive to small changes in water pressure than deformability of ice (Fowler, 2009).We tested the output from these models against observations to determine which model was better able to reproduce estimates of the inferred magnitude and recurrence interval for known subglacial lake-drainage events.

Model description
Our model is a prototype for lake drainage via a single channel incised into a linearly viscous subglacial sediment.The overall principle of lake drainage via a channel where water pressures are significantly below overburden pressure (in our case 2-8 m water equivalent, or w.e.) is well established in the literature (e.g., Nye, 1976;Fowler, 1999;Evatt et al., 2006).Only recently has it been suggested that such drainage may occur via channels forming in sediment rather than the overlying ice (e.g., Fowler, 2009).In order to accommodate a channel incised into sediment rather than in ice, we adopted principles of sediment mechanics from theoretical work presented in Walder and Fowler (1994) and Ng (2000).
Given that borehole observations on the Siple Coast region have indicated a distributed system at the ice-bed interface (e.g., Engelhardt, 2004), we needed a simple model for the coupling of the channel to a distributed system.Therefore, we adapted a formulation presented in Kingslake and Ng (2013) which, although designed for a rocky-bottomed mountain glacier lake drainage, contained the four most essential elements we required for our ice-stream experiment: a lake, a distributed-flow system, a channelized system, and a means for communication between all parts of the system.We developed two models based on the Kingslake and Ng (2013) foundation: one model with channelization through conduits melted into the basal ice based on R-channel theory (henceforth called the R-channel model) and one model with channelization through canals eroded into the basal sediment (henceforth called the canal model).Although there may be more appropriate models for each of the individual components in the Antarctic subglacial environment (see Flowers (2015) for a review of existing water models), we found the ease of coupling between the lake, channelized system, and the distributed system provided by the Kingslake and Ng (2013) formulation to be the most effective for proof of concept.Input data for the model come from several sources: previous water-budget studies for lake-inflow estimates, RES for ice thickness, and satellite and airborne radar altimetry for surface height changes.By comparing the lake-volume change inferred from satellite observations with lake-volume change as simulated with R-channel and canal models, we can (i) perform a diagnostic test for the hypothesis by Fowler (2009) that drainage could occur for sediments that behave like erodible-deformable ice and (ii) develop a potential prototype for simulating Antarctic subglacial lake drainage.
Our model consists of a system of equations on a 1-D finite difference grid, with scalar values, such as ice base elevation, ice thickness, water pressure, and channel width calculated at the center of the grid cells, while fluxes of sediment and water are calculated at intermediate points.The systems of equations for the model are described in Sect.3.1-3.3.Methods of solution and time integration are described in Sect.3.4.We describe our method of obtaining input data and the size of the domain in Sect.3.5.

Theory for subglacial water flow
Subglacial water flows from areas of high hydropotential to areas of lower hydropotential.Hydropotential at the base of the ice, θ 0 (calculated in m w.e.), is the sum of the watersystem elevation, z b , and water pressure, p w , normalized by water density, ρ w , and gravity, g (Shreve, 1972): Many models (e.g., Le Brocq et al., 2009;Carter et al., 2011) calculate θ 0 assuming water pressure at the ice base is equal to the overburden pressure, p o , because the quanti-ties used to calculate p o (ice surface elevation z s , z b , and ice density ρ i ) are easily measured: In reality, p w is the difference between the overburden pressure and effective pressure, N, which can only be calculated with seismic data (e.g., Blankenship et al., 1987) or borehole data (e.g., Engelhardt, 2004).Although N does not typically affect regional water routing, modeling work has shown that temporal change in N is a critical part of the lake-drainage process (e.g., Fowler, 1999;Evatt et al., 2006;Fowler, 2009).As N decreases, it enables the formation of a temporary hydraulic divide between the lake and points downstream.As N increases, water can then overcome the divide leading to the onset (and on-going) drainage of that lake.Most of the theoretical work surrounding subglacial floods expresses hydropotential, θ , in terms of pressure units (Pa), which in this work is calculated as follows: (1c) While we also use θ (in Pa), we mostly express hydropotential as θ 0 (in m w.e.) to allow more direct comparison between it and the measurements used to calculate it.
Given that θ and N, as well as channel cross-sectional area (S) and water flux (Q), are each defined differently for the R-channel, canal, and distributed-flow systems, we will use the subscripts RC , CC , and S , respectively, to avoid confusion in the following sections.Explanations for all symbols not defined explicitly in the text can be found in Tables 1 and 2.

Subglacial channel formation
The principal governing equations for the two channelization systems in our model all come from previous work: (i) the evolution of the R-channel from Röthlisberger (1972) and Nye (1976) and (ii) the evolution of the canal model from Walder and Fowler (1994) and Ng (2000).The method for solving both models has been largely adopted from recent work by Kingslake and Ng (2013), which specifically described the drainage of an ice-dammed lake via an R-channel coupled with a distributed system and thus was already well posed for subglacial lake drainage.Exchange between channelized and distributed water systems in both models also follows Kingslake and Ng (2013) and is mediated by a term describing water transfer between the systems (T RC for the R-channel model; T CC for the canal model), which is a function of the pressure difference between the channelized and distributed systems.To adapt the Kingslake and Ng (2013) model fully to a domain containing Antarctic subglacial lakes, however, we made additional modifications, including defining the hydropotential gradient with Eq. ( 1), so that it is a function of the slope of the ice base and the effective pressure as well as the slope of the ice surface as defined in Kingslake and Ng (2013).We also allow for the change in Table 1.List of all symbols used in this paper with definitions, dimensions, and value/method used to estimate them.

Symbol
Meaning Dimensions Value/method Constant of sediment deformation -1.33 (Walder and Fowler, 1994) b Constant of sediment deformation -1.8 (Walder and Fowler, 1994) Characteristic grain size L Input (see Table 2) Glen's flow law parameter for ice L 3 T 5 M −3 10 −24 Pa −3 s −1 (Kingslake and Ng, 2013) Coefficient for partitioning of turbulent energy between heating water and melting surrounding ice -0.309 (Hooke, 2005) Coefficient for transmission efficiency between R-channels and distributed-flow systems -Input (0.05) Critical hydraulic shear stress necessary to initiate erosion Output by model pressure melting point of water with change in pressure (e.g., Alley et al., 1998;Hooke, 2005), as this effect has been reported to be significant to hydrological evolution over other flow paths following adverse bed slopes in Antarctica (e.g., Carter et al., 2009;Fricker et al., 2014).

Channels incised into the ice (R-channels)
In the classic R-channel model, transmissivity through a channel is controlled by the channel's cross-sectional area (S RC ), which is a balance of channel-wall melt rate (m RC ) and viscous-ice deformation rate (C VRC ) such that This formulation, originally presented in Röthlisberger (1972), has since been adapted for a variety of situations.The melt rate is related to the flux of water through the channel (Q RC ) by www.the-cryosphere.net/11/381/2017/The Cryosphere, 11, 381-405, 2017 while transfer between the R-channel and the distributed system (see Sect. 3.2.3) is governed by the transfer term T RC (Kingslake and Ng, 2013) using the equation where R kRC is the transfer efficiency, k is a constant for connectivity from Kingslake and Ng (2013), and N RC and N S are the effective pressures of the R-channel and distributed system, respectively.It should be noted that our value for R kRC is near the lowest end of values explored by Kingslake and Ng (2013), primarily due to model stability issues.Pressure along the channel co-evolves with water flux through an adaptation of the Manning friction formula, following Kingslake and Ng (2013): where f r is the hydraulic roughness.

Channels incised into the sediment (canals)
Our model for channelization by mechanical erosion into the sediments is adapted extensively from concepts in existing models of lake drainage via R-channels eroded into the ice (e.g., Röthlisberger, 1972;Nye, 1976;Fowler, 1999;Evatt et al., 2006;Kingslake and Ng, 2013).It includes the basic principle of semi-circular channels capable of sustaining water pressures several meters w.e.below flotation (Röthlisberger, 1972;Fowler, 1999), as well as more recent concepts regarding exchange of water between the channelized and distributed systems and evolution of pressure (Kingslake and Ng, 2013).In contrast to R-channel models, our model replaces melting of ice with sediment erosion and ice closure with deformation of sediment with linear viscous rheology, borrowing concepts and formulations from Walder and Fowler (1994) and Ng (2000).In describing channels eroded into the sediment, our work follows other efforts to model canals incised into the sediment behaving as conduits (e.g., van der Wel et al., 2013;Kyrke-Smith and Fowler, 2014) but, in contrast to these models, focuses specifically on subglacial lake drainage.Our descriptions of sediment erosion and deformation and of channel geometry are all extensively simplified; we regard the model developed below as a proof of concept that can be further refined if it is able to reproduce lake-volume change as inferred from satellite and GPS observations (e.g., Fricker and Scambos, 2009;Siegfried et al., 2014Siegfried et al., , 2016) ) with realistic parameter choices.
As with the R-channel model, transmissivity for a channelized system with canals is governed by the canal crosssectional area (S CC ).In our canal model, we ignore melting (m RC ) and deformation (C VRC ) of the ice above and assume all change to S CC is due to erosion and deformation of the sediment.Change in aperture in this model therefore is a balance of net erosion (erosion, E CC , minus deposition, D CC ) and sediment deformation (C VCC , for which we assume a linear viscous rheology): Adopted from Walder and Fowler (1994), E CC is a function of the mean sediment settling velocity (v s ), the channel geometry (α CC ), the stress exerted on the bed by the flowing The Cryosphere, 11, 381-405, 2017 www.the-cryosphere.net/11/381/2017/water (τ CC ), and the particle size (d 15 ): where K T 1 is a sediment erosion constant and "max" refers to the maximum of τ CC −τ k and zero.α CC is a dimensionless correction factor between 5500 and 233 000, which accounts for the geometric differences between the semi-circular channel geometry implied by the formulation and the actual geometry which is likely more elliptical in nature (e.g., Ng, 2000;Winberry et al., 2009;van der Wel et al., 2013).D CC , which also follows from Walder and Fowler (1994), is defined similarly: where K T 2 is a sediment deposition constant and c is the concentration of sediment in the water column.Equations ( 8a) and (8b) imply that erosion and deposition are both a function of the hydraulic shear stress (τ CC ), which is a function of water velocity.Following Walder and Fowler (1994), we use the mean water velocity (u CC ) such that where f CC is a constant describing roughness.In our model, erosion begins once τ CC exceeds a critical value, τ k , which is a function of median grain size as described by In this description, erosion of the sediments occurs only when water velocity exceeds a certain value since τ CC is only dependent on u CC .Sediment settling velocity (v s ), important for both Eq. ( 8a) and (8b), is treated as a function of sediment particle size (d 15 ), the density contrast between sediment and water (ρ CC − ρ w ), and water viscosity (µ w ), by the classic derivation from Stokes' law (Lamb, 1932): Given that the net deposition rate as defined in Eq. ( 8b) is sensitive to the concentration of sediment already present, a mechanism is required for keeping track of sediment concentration.To achieve this, we assume that the transport rate of sediment downstream comes into equilibrium instantaneously.Sediment concentration ( c) at the most-upstream cell where erosion takes place is given by For each cell downstream, c is calculated as the sum of sediment eroded within the cell locally and the concentration of sediment in the water within the next cell upstream ( cup ): Closure occurs via viscous sediment creep, C VCC , though this may in reality be a convenient continuum representation of discrete sediment-collapse events on the sides of the channel: Here N ∞ is the regional sediment effective pressure (which is closely correlated with sediment strength), N CC is the effective pressure of the water within the canal, and A CC , a, and b are flow-law constants from Walder and Fowler (1994).
A CC is similar to K RC in Eq. ( 4).Although our formulation for Eq.(8a-h) is largely adapted from work by Walder and Fowler (1994), we follow Ng (2000) in decoupling N CC from N ∞ .In this work, however, we depart from Ng (2000) and treat N ∞ as a constant and explore this simplifying assumption in more detail in Sect.4.3.
In this system of equations, channel growth via erosion is a function of water velocity (through the relationship between erosion and shear stress) and is sensitive to the sediment size and the channel geometry.At lower water velocity, channelization in the sediment will not initiate and sheet flow will dominate.Channel closure via deformation is a function of the effective pressure (N CC ) and sensitive to the chosen value for sediment effective pressure (N ∞ ).Although our formulation of Eq. (8h) theoretically allows for viscous channel growth if water pressure exceeds overburden pressure, such pressures are not expected (see also Sect. 4.1.3).
Conservation of water mass is accomplished by T CC is determined by the effective pressure difference between the canal and the distributed sheet via with R KCC behaving analogously to R KRC (Eq.6).Finally we define the propagation of effective pressure (N CC ) along flow using a version of the Manning friction formula adapted from Kingslake and Ng (2013): where we have introduced the roughness parameter f CC = f r = 0.07 m −2/3 s 2 based on the assumption of a semicircular channel geometry following Kingslake and Ng (2013).The distributed sheet-flow system, which is a component of both the R-channel and canal model, is governed by three primary equations all modified from Kingslake and Ng (2013): two concern the conservation of mass, and a third governs the evolution of N S .Water flux through the sheet (Q S ) is function of the hydraulic gradient (∂θ S /∂x) and cross-sectional area (S S ): where Equations ( 12) and (12a) assume flow-path width, y S , is constant with depth, such that water flux increases linearly with water layer thickness, h S .Evolution of S S is governed by the conservation of mass, where M S is a source term for water flowing into the system from the sides and T is the flux between the sheet and channelized system, equal either to T CC or T RC depending on whether the distributed-flow system is coupled to the canal or R-channel model.Equation ( 12) implies that transmissivity increases linearly with water storage, while Eq. ( 13) implies that local water storage will increase if outflow from a cell is less than inflow.If water pressure is equal to overburden pressure then the only way for a water layer to exit an enclosed basin in the hydropotential is to thicken until it overtops the lowest saddle.With the inclusion of the effective pressure term (N S ) in the equation for hydropotential (Eq.1b), distributed water layers have another way over small obstacles.In a distributed system, N S tends to increase as water thickness decreases, by the following relation: Our parameterization for N S is adapted from the Kingslake and Ng (2013) description of a linked-cavity system, but substitutes till cavities as observed by Engelhardt (2004) for hard bed cavities.This formulation assumes N S reaches steady state instantaneously.Thinner water layers (and therefore higher values of N S ) are maintained over hydropotential maxima, while thicker water layers (and therefore lower values of N S ) are maintained over hydropotential minima.This coupling allows for a monotonically decreasing pathway in θ S in spite of undulations in overburden pressure, which would inhibit the flow of water downstream without taking variable N S into account.

Modeling the lake
For a lake, we assume that the hydropotential is the sum of the ice base elevation (z b ) and overburden pressure (derived from h i ) at the center of the lake (Fig. 3).While hydropotential over the lake can then rise and fall in response to filling and drainage, we assume that the variation of hydropotential with time is spatially uniform across the lake.With the subscript L referring to conditions (hydropotential, elevation change, etc.) across the lake, the change in lake level with time is defined as Change in lake volume is calculated as and where Q in is inflow, Q out is outflow, A L is lake area (which is held constant in both the R-channel and canal models), and F L is a dimensionless number that varies between 1 and 2 depending on the degree to which the overlying ice is supported by flexural bridging stresses (effectively a simplification of the parameterization of surface deformation in Evatt and Fowler, 2007).For large lakes that have a large section of freely floating ice, where the zone of flexural deformation is less than 20 % of the total lake area, F L = 1; subglacial Lake Engelhardt (SLE) and subglacial Lake Conway (SLC) fall into this category (for lake locations, see Fig. 1b).For smaller lakes under thicker ice where the flexure zone comprises the most or all of the lake area and surface deformation can be approximated by a parabola, F L = 2; subglacial Lake Mercer (SLM) and SLW fall into this category.We assume that there is no change in area as the lake state changes, and the implications of this assumption are explored in Sect.4.2.2.We also assume that all surface displacement over the lake results directly from changes in lake volume and neglect any temporal changes in ice thickness.We also assume that there is only a single channel.Based on earlier runs of the R-channel and canal models as well as assertions in Shoemaker (2003), we assume that the destination lake is the nearest major low in the hydraulic potential, empirically defined to be an enclosed basin 3 m w.e.deeper than its surrounding.This assumption generally agrees with observations (e.g., Fricker et al., 2007).We neglect changes in the destination lake's level and assume that N = 0 at its center.For the canal model, we also assume (i) no change in sediment water content or rheology over time; (ii) negligible bridging stresses; (iii) near-instantaneous advection of sediment downstream.These systems of Eqs.(2-17) are solved on a 1-D finite difference domain, consisting of a source lake, intermediate points, and a destination lake.The point at which θ 0 reaches a local maximum downstream of the lake is termed "the seal".This point is controlled by ice-sheet geometry and is treated as immobile, in contrast to its definition in some higher-order models which define the seal as the location where transmissivity is lowest at any given time (e.g., Fowler, 1999;Clarke, 2003).

Model spinup and initialization
Given our initial uncertainties about the onset of channelization, it is simpler to run our model when the source lake is filling, such that outflow from the lake is minimized, and the lake level is well below discharge level.We initialize the model assuming that water pressure equals overburden pressure and that there is a constant supply of meltwater, M C , along the flow path downstream of the lake, such that initial water flux in the distributed system (Q b ) increases with distance downstream of the lake.With these initial values, we calculate an initial water layer thickness, assuming N S = 0, reworking Eq. ( 12) such that We then use Eq. ( 14) to calculate a new value for N S , based on S S and Eq. ( 13).For each subsequent time step we calculate Q S from Eq. ( 12), ∂S S /∂t from Eq. ( 13), as well as change in lake properties from Eqs. (15)(16)(17).With the resulting new water layer thickness, we finally recalculate N S and S S .

Criteria for onset and shutdown of channelization
Initially water will flow from the seal towards the lake; as the lake level and lake hydropotential increase, however, water will begin flow out of the lake over the seal.Due to the inverse relationship between N S and Q S , Q out will initially be much lower than Q in .We assume that the initial outflow does not erode a channel at the seal but remains as sheet flow.Walder (1982) suggested that water velocity in a distributed system would tend to be higher in places where irregularities in the overlying ice base and underlying bed lead to a thicker water layer.When the bed consists of sediments, τ CC would exceed a critical threshold (τ k ) in these locations www.the-cryosphere.net/11/381/2017/The Cryosphere, 11, 381-405, 2017 before doing so elsewhere.As the water layer thickness increases, erosion of the bed and/or melting of the ice above would be concentrated in these locations, leading to a positive feedback between local water layer thickness and water flow until a channel develops either in the bed below (a canal) or ice above (an R-channel).
Our model tries to simplify these complex dynamics by using a threshold criterion for channel initiation in which a channel carrying 0.5 m 3 s −1 appears spontaneously once Q S exceeds a critical threshold value (Q onset ).Similar threshold methods have been employed in previous works, such as Flowers and Clarke (2002), Peters et al. (2009), de Fleurian et al. ( 2014), and Hoffman and Price (2014).The value for flux through a newly formed channel is based empirically on work by Peters et al. (2009), which indicated that higher initial values would lead to an unrealistically rapid growth of outflow.As this value is also well below Q onset , it allows for a gradual transition between a distributed and channelized system.
Effective pressure is initially set equal to the effective pressure in the distributed system so that N RC = N S for the Rchannel model (or N CC = N S for the canal model).These initial values for water pressure then allow us to calculate an initial S RC (or S CC ) through Eq. ( 7) (or Eq. 11).Cessation of channelization occurs once Q RC (or Q CC ) falls below a threshold value, Q shutdown , at which point it becomes computationally simpler to eliminate the channelization and revert back to sheet flow.Although it is possible that small incipient channels are always present, such as is parameterized by Flowers et al. (2004), ignoring them improves computational speed.In Sect.4.2, we explore in more detail how variations in Q onset and Q shutdown affect model performance.

Evolution of channelized flow
Once the channel is initiated we calculate the geometry of a proto-channel (Eqs.7, 11), assuming Q RC = Q onset (or Q CC = Q onset ) everywhere between the source and destination lakes.After each successive iteration of Eqs.(2-4) (or Eqs.8-10), we recalculate Q RC and N RC (or Q CC and N CC ) along the flow path using a shooting method: beginning with an initial guess of Q RC (or Q CC ) at the outflow point, we use Eq. ( 7) (or Eq. 11) to calculate dN RC /dx (or dN CC /dx) locally on the staggered grid and the corresponding values for N RC and T RC (or N CC and T CC ) at the next point downstream on the regular grid.Using these newly calculated values for N RC and T RC (or N CC and T CC ), we calculate dQ RC /dx (or dQ CC /dx) on the regular grid and then Q RC (or Q CC ) at the next point downstream on the staggered grid.This process is repeated downstream until we calculate N RC (or N CC ) using Eq. ( 5) (or Eq. 9) at the destination lake, at which point it is a known quantity.We compare our calculated N RC (or N CC ) with the "known" N RC (or N CC ) to obtain a misfit.We then use a Newton's method iteration on N RC (or N CC ) at the downstream lake (treated as a function of Q RC (or Q CC ) at the source lake outlet).After a maximum of 12 iterations we have arrived at a value for Q RC (or Q CC ) at the source lake that results in a value for N RC (or N CC ) at the destination lake that is within 0.01 m w.e. of the known N RC (or N CC ).Once this value is obtained we then iterate Eqs. ( 2), (6), and (7) (or Eqs. 8, 10, and 11) and the corresponding sheet-flow evolution (Eqs.12-14).To improve model stability, we apply a modification of the Courant-Friedrichs-Lewy (CFL) (Courant et al., 1967) (19b) for the canal model, where "min" refers to the minimum of all expressions within the parentheses.

Real and idealized model domains
Our idealized model domain (Fig. 4a) was based on a simplified version of the flow path connecting the well-studied (e.g., Christianson et al., 2012;Horgan et al., 2012;Tulaczyk et al., 2014) SLW to the Ross Sea from Carter and Fricker (2012) (see Fig. 1b for location; Fig. 4b for comparison to the real domain).In this domain, the segment between the source lake and the seal, and the seal and the destination lake, is approximated with straight lines.We tested both the R-channel and canal models on this domain and compared the model outputs.y S and h i were held constant at 2500 and 500 m, respectively.A L was set at 100 km 2 .Based on water budgeting work by Carter et al. (2013), Q b was between 0.24 and 0.35 m 3 s −1 .Values for other constants in this model run can be found in Table 2.
To demonstrate the ability of the canal model to reproduce the timing and magnitude of actual observed lake-drainage events, we also applied it to several "real" domains for flow paths draining the lakes in lower Whillans and Mercer ice streams: see Fig. 4b (for SLW), c (for SLC and SLM), and d (for SLE).For the real-domain values, we obtained values for ice thickness, h i , ice base elevation, z b , and pathway width, y S , from RES-derived measurements of ice thickness and surface elevation made between 1971 and 1999 over multiple campaigns (see Carter et al. (2013) for a description of the interpolation strategy and Lythe et al. (2001) for descriptions of data used).For the idealized domain we used measured values of θ 0 and θ S at the source lake, seal, and destination lake and then fit a straight line between them, holding ice thickness, and flow path width constant.

Summary of experiments
Our experiments began with several tests comparing the simulated lake-volume time series output by the R-channel model against the volume change time series inferred from observations for SLW (Fricker and Scambos, 2009;Siegfried et al., 2014).After experimenting with the unaltered Rchannel model we then experimented with perturbations to parameters controlling the rate of channel grown (m RC ) and closure (C VRC ) following hypotheses put forth in Fowler (2009) and Evatt et al. (2006) (Sect. 4.1.1).
All experiments involving lakes on a real domain were run until the lake had undergone at least one complete fill-drain cycle and at least 10 years had elapsed.As the sensitivity studies typically involved a small lake, they were run for only 8 years.In cases where the channel failed to grow the model was allowed to run for 10 years.In cases where runaway growth in outflow rate led to unrealistically large ice surface drawdown, the model was stopped once the lake level had decreased by over 30 m.In order to facilitate easier comparison of modeled lake-volume change with that inferred from Ice, Cloud and land Elevation Satellite (ICESat) observations (which spanned 2003-2009) we referred to the time of observed volume change in years since 2000 and adjusted the timing of the first model year to coincide with the first observed filling cycle on the lake being simulated.
Following the suggestion that the channel closure rate inferred from lake-volume change observations would be more consistent with the rheology of soft sediment, we focused subsequent experimentation on the canal model.Lakevolume change as simulated by the canal model on a realistic domain was compared against lake-volume change as observed for lakes SLW, SLE, SLC, and SLM (Sect.4.1.2).In order to explore the possibility of implementing the lakedrainage model in areas where ice thickness measurements are more sparse, we compared the model's output on ideal-ized (straight lines) and realistic (interpolated at 1 km intervals from RES and satellite altimetry) flow-path geometries (Sect.4.1.3).The simplifications inherent in the idealized domain also proved useful for quantifying better how variations in the flow-path geometries between the source and destination lakes might affect the magnitude and recurrence interval of lake drainage (Sect.4.2.1).This idealized domain was also employed to explore how the time series output by the canal model varies in response to perturbations both to the input data, including inflow and lake area (Sect.4.2.2), as well as changes to parameters that are more difficult to measure directly, such as α T , d 15 , N ∞ , R 1 , and Q onset (Sect.4.2.3 and 4.2.4).

Overview
All time series for lake-volume change output by the models tested in this work were qualitatively similar to observed time series.Most models also reproduced key elements of the lake-drainage process including initiation of outflow before the lake level reaches the flotation height, and the beginning of refilling without complete drainage of the lake.We begin with the results of the traditional R-channel model (Sect.4.1.1)and then describe the results of the canal model (Sect.4.1.2).We compare results from the real and idealized model domains (Sect.4.1.3)and finally explore the model's sensitivity to changes in key input parameters on the idealized domain (Sect.4.2).

R-channel model vs. observed height anomalies
Our lake-drainage model based on outflow through an Rchannel on the realistic SLW domain (Fig. 4b) with realistic parameter choices (Tables 1, 2) produced a single, large drainage event over the 30-year model run (Fig. 5a).The long fill-drain cycle resulted in part from the presence of an adverse bed gradient along segments of the flow path.In these segments all of the turbulent heat generated by water traveling down flow went to maintaining the water at the pressure melting point rather than melting a channel (Fig. 5b, c).From the start of the model run, it took nearly 10 years for a significant channel to begin growing, by which time the stiffness of the ice was too large to halt the lake drainage once the lake drained back to its initial level (defined as 0 m).Only after draining for almost 10 years and losing almost 16 m of elevation from its highstand did Q out fall below Q in and lake volume begin increasing.
This modeled fill-drain cycle does not agree with observations from height anomalies at SLW (e.g., Fricker and Scambos, 2009;Siegfried et al., 2016), which show an ∼ 4-year fill-drain cycle (dashed line in Fig. 5a, d).In an effort to match the observed height-change time series, we adjusted the parameter choices that controlled channel opening and  closure (Fig. 5d).By lowering the value for the latent heat of fusion (L h ), we were able to increase the channel opening rate and produce a timing of highstand that better approximated observations.However, higher melt rates continued to exceed channel closure rates and outflow increased exponentially.We adjusted the rate of channel closure by increasing K RC , thereby increasing the rate of ice deformation into the channel.While softening the ice alone led to channels that never opened, when K RC was increased by a factor of 40 and L h was reduced by a factor of 400 (black line in Fig. 5d) we were able to achieve a recurrence interval and magnitude of volume change similar to previous observations.Given that these values were roughly consistent with those found by The Cryosphere, 11, 381-405, 2017 www.the-cryosphere.net/11/381/2017/Fowler (2009) (which suggested closure could occur if K RC was raised by 2 orders of magnitude), our finding supports the concept that most of the lake-drainage events inferred from currently available observations are unlikely to occur via an R-channel.

Canal model vs. observed height anomalies
For the real domains based upon SLW, SLE, SLC, and SLM, our canal model with realistic parameter choices (Tables 1, 2) produced a simulated time series of lake volume that appeared qualitatively similar to the observed lake-volume time series for each of these lakes with respect to the magnitude and recurrence interval of lake fill-drain cycles (Fig. 6).At low stand, distributed flow carried water toward the lake.As θ 0 at the lake increased with increasing z bL , the difference between it and θ 0 at the seal decreased until the seal was overtopped.At this point, a monotonically decreasing path from the lake over the seal could be maintained, allowing distributed flow out of the lake as lake level rose.
As the level of the lake increased further, the value for N S required to maintain a continuous path downstream decreased, and fluxes from the distributed system increased.Even after the onset of channelized flow (i.e., Q CC > 0.5 m 3 s −1 ), Q S initially remained closely correlated with lake volume, peaking slightly before the lake reached highstand.As channelized flow (Q CC ) increased, however, the channel soon evolved to be the more efficient outlet from the lake.With further growth in Q CC , lake level began to decrease more rapidly.Declining lake levels led to decreasing outflow via the distributed system as declining θ S upstream led to a reduction in water inflow and in water layer thickness at the seal.A decline in θ CC also caused a reduction in Q S , allowing the channel then to draw water away from what remained in the distributed system (Fig. 6).
With further lowering of the lake level, however, the gradient of θ CC decreased, limiting the erosive power of the channel.Simultaneously, the increase in N CC allowed for more sediment to creep into the channel, decreasing channel size.While in many model instances Q CC never ceased completely, it did eventually fall below the inflow level of the lake such that the lake could fill again.If a channel remained continuously open then a seal never developed in θ CC , only in θ S .With successive filling and draining cycles, growth in S CC lagged behind the level of the lake such that fluctuations in lake volume continued.In model runs where the channel remained open, the amplitude of volume change dampened through time (Fig. 6d), a phenomenon explored in more detail in our sensitivity studies (Sect.4.2.2).
We achieved a reasonable agreement between modeled and observed recurrence intervals and magnitude of lakevolume changes for SLW, SLC, and SLE using constant values for Q in at each lake (Table 2; Fig. 6a-d).However, our model could not reproduce observed periods of slower filling, such as was observed at SLE between 2006 and 2008 (Fig. 6b).The weakest fit between model and observations occurred for SLM, for which our modeled frequency of filling and drainage events was over 50 % higher than observed www.the-cryosphere.net/11/381/2017/The Cryosphere, 11, 381-405, 2017 and ultimately out of phase.Work by Carter et al. (2013) has inferred that the filling rate for SLM varied between 2.25 and over 50 m 3 s −1 and was controlled primarily by outflow from SLC, suggesting that the misfit could reflect the poor assumption of non-varying Q in .Achieving a simulated time series where the recurrence interval and magnitude of volume change were close to those observed required that we vary the input parameters substantially from lake to lake (Table 2).
We describe these sensitivities in more detail in Sect.4.2.

Canal model on the idealized vs. realistic domain
A comparison between the model output from an idealized domain and one that uses more realistic flow-path geometry showed almost no qualitative difference, with recurrence interval, flux, and volume ranges within 5% of one another (Fig. 7a, b).The primary difference between the two models was the length of time they took to run.While the linearly varying hydropotential of the idealized domain (Fig. 7c) resulted in a nearly uniform cross-sectional area of the distributed system (Fig. 7e), undulations in θ 0 between the seal and the destination lake in the realistic domain (Fig. 7d) resulted in large along-flow variations in distributed-system aperture (Fig. 7f) and corresponding effective pressure.In particular, the advection scheme employed resulted in minor oscillations in the amount of water stored along the flow path, with more water stored in areas where θ 0 was concave upward and less water stored in areas where θ 0 was concave downward.On these local maxima, small changes to water storage still comprised a large portion of the total and require a very short adaptive time step (see Eq. 17b).
In both domains, values for θ CC quickly converged upon a nearly constant gradient in θ CC between the seal and the des- tination lake (Fig. 7c, d), and thus rates of water flux, channel growth, and sediment deformation back into the channel were similar.As the canal transferred over 90 % of the outflow downstream, the evolution of this feature was the most significant component for determining the timing and magnitude of lake-volume change.

Sensitivity to flow-path geometry
Given that our lake-drainage model is a bottleneck-style model, in which the seal is stationary while the lake level rises and falls, we expect the rate of channel growth to correlate with hydropotential gradient, as well as the seal height that must be overcome initially.We tested the model's sensitivity to flow-path geometry by increasing (Fig. 8a-c) and decreasing (Fig. 8d-f) hydropotential gradients.We also explored both changes in total slope (Fig. 8a, d) as well as the slope between the source lake and the seal (Fig. 8b, e) and the seal and the destination lake (Fig. 8c, f).Overall we found that steeper slopes favored a more rapid, highervolume lake drainage with the slope downstream of the seal being the most important factor.Changes to the flow-path geometry upstream of the seal tended to be more important for determining the peak lake height before drainage would occur (Fig. 8b, e), with gentler slopes leading to higher highstands and only a minor effect on the total volume range.

Sensitivity to inflow and lake area
At lower Q in , there was an inverse linear relationship between recurrence cycle and input (Fig. 9a, b).At higher Q in , higher levels of inflow initially resulted in higher lake levels at highstand as the outflow channel took longer to grow large enough to drain the higher inflow of water.While the higher highstands and outflow levels initially provided a steeper hydropotential gradient resulting in a larger net volume loss, the higher inflow caused a net increase in volume when there was still a substantial quantity of water leaving the lake.This process resulted in a lower net filling rate and lower highstand volume during the next cycle, but a larger channel at Figure 9. Sensitivity of canal model to inflow and lake area.Variation in (a) surface elevation anomaly (z sL ) and (b) channelized outflow (Q CC ) due to changes in inflow (Q in ).Variation in (c) lake volume (V L ) and (d) mean elevation anomaly (z sL ) due to changes in lake area (A L ).The black lines denote the "control" run using our preferred parameter values (Table 2).
low stand.As a consequence the magnitude of the fill-drain cycle dampened more rapidly with each successive cycle for higher values of Q in .When we explored the effect of variations in A L (Fig. 9c, d), we found similar results with smaller lakes experiencing greater ranges in elevation initially, but then decreasing as outflow and inflow came into equilibrium.

Sensitivity to sediment and channel properties
Small changes to channel geometry, sediment effective pressure, and grain size (α T , N ∞ , d 15 ) significantly affected lake level at highstand, the magnitude of volume change, and the recurrence interval of drainage events (Fig. 10).Moreover, the variations in N ∞ to which our model is sensitive (Fig. 10b) are below the errors of current measurements (from borehole (e.g., Engelhardt, 2004) or seismic (e.g., Blankenship et al., 1987) experiments).Although we could find suitable values which resulted in a match between model runs and observation (Table 2), there is a significant element of non-uniqueness to these parameters, especially α T and d 15 , which are related through Eq. ( 8a) and (8b).

Lower sensitivity parameters
Any factor which affected the relationship between h S and N S (e.g., d(z s )/dx, R 1 , y S ) had an important influence on volume range because such factors control how easily water overcomes an apparent obstacle in the hydropotential.Any suite of parameter values that made it easy for water to clear apparent barriers had the effect of reducing the total amplitude of volume change, favoring steady drainage by sheet flow (Fig. 11a).Conversely, parameter values that made it difficult for water to overcome obstacles tended to increase the volume range and favor channelized drainage.
The range of values for Q onset was limited by the maximum inflow (e.g., ∼ 5 m 3 s −1 for SLW, inferred by Carter et al., 2013).Q onset had only minor effects when compared to the obstacle size (R 1 ) with respect to determining lake level at the onset of channelization and volume range.Given that differences between Q onset = 1 m 3 s −1 and Q onset = 3 m 3 s −1 were minimal, values for Q onset < 0.5 m 3 s −1 would likely be even less significant.Likewise, M c (inflow from other sources, such as basal melt and flow from outside the flow path) seemed to matter only when it was over an order of magnitude higher than values that would be reasonable for a location like Whillans or Mercer ice streams.Under higher values of M S , highstands were substantially lower but volume range was nearly unchanged (Fig. 11b).This low sensitivity may result from our low value for R kRC , which limited the transfer of water between the channelized and distributed systems.When values for M S exceed 0.1 m 3 s −1 km −1 it affected the rate at which the source lake filled. .The black lines denote the "control" run using our preferred parameter values (Table 2).

Discussion
Although our models suggest that an R-channel can theoretically grow and contract in the Antarctic subglacial environment, supporting previous work by Carter et al. (2009) and Peters et al. (2009), their rates of growth and shutdown do not match those inferred for most active lakes unless the modeled melt and closure rates are adjusted significantly.If drainage through a canal system is the dominant mechanism for "active" lakes, then it would suggest that such lakes are more likely to be found in areas where soft sediment is widespread.In such an environment many of the criteria used to locate lakes with airborne RES might fail.Here we explore in more detail (i) the validity and implications of the assumptions used in achieving our results; (ii) what might be inferred about Antarctic subglacial hydrology based on our results in the context of previous work; and (iii) the prospects for modeling subglacial lake drainage in ice-sheet models.

Reassessment of model simplifications and assumptions
We did not include a model in which erosion and deformation of both the ice and bed occurred simultaneously.Given the difference in recurrence intervals predicted for Siple Coast style lake-drainage events via the canal model vs. an Rchannel, we expect that the evolution of the canal will dominate short-term flow variability.However, as the canal appeared to remain open throughout the fill-drain cycle, it is possible that ice deformation, melting, and freezing would all play a role in longer-term water-flow evolution.In particular if the basal ice is softer than would be implied by our assumed K RC of 10 −24 Pa −3 s −1 (Budd and Jacka, 1989) and/or the flow path follows a substantial adverse bedrock slope and is subject to high rates of basal freeze-on such as observed in Recovery Glacier (Fricker et al., 2014), then deformation of the ice into a canal could significantly affect water flow and drainage system geometry.This process, which couples canals to ice deformation, should be considered in future iterations of our model.Additionally, we have not explored the sensitivity of the model to the initial values for Q RC and Q CC .Given that in most model runs the channel does not ever fall below Q shutdown , it seems that alterations to this initial value for Q RC and Q CC would affect only the timing of the first filling and drainage cycle and thus any issues can be resolved with careful model spinup.
One of the major challenges to our present model is the difficulty in measuring many of the properties to which our model is most sensitive.We have limited observations of grain size and sediment effective pressure (e.g., Tulaczyk et al., 1998;Engelhardt, 2004), limited channel geometry estimates from RES (e.g., Schroeder et al., 2013), and limited sediment effective pressure estimates from seismic surveys (e.g., Blankenship et al., 1987;Luthra et al., 2016).Even where present, estimates for these parameters are still not sufficiently accurate to constrain our model effectively.As a consequence, most of our simplifying assumptions were made due to data constraints, including those regarding the A L remaining constant over the filling-draining cycle, the till strength as a function of N ∞ remaining constant over time, and our assumption of a semi-circular channel geometry for the canal.For example, A L likely correlates positively with lake volume.This effect alone would only moderately affect the rate of draw down leading to an increase in surface lowering relative to what would be expected in a constant A L case.Lower lake levels, however, would lead to faster rates of channel closure, partially counteracting this effect; these issues remain to be explored.
If N ∞ is also increasing as the lake drains and shrinks, then we may expect the system eventually to come into equilibrium at low stand until flow from upstream increases.Other results exploring temporal changes in till strength (which indirectly relates to N ∞ ), however, suggest that N ∞ is likely to be changing significantly over the duration of our model  run as water is exchanged between the sediment pores and ice-bed interface (Christoffersen et al., 2014).When we consider our assumption that N ∞ remains constant over time in light of these results, we can speculate on how variations in N ∞ over time might affect a subglacial lake's filldrain cycle.Based on research investigating the exchange of water between the interfacial flow system and sediment (e.g., Christoffersen et al., 2014;Bougamont et al., 2015) we would expect N ∞ to co-vary with N S and N CC .In this situation, sediment strength would increase as lake level declined, N CC and N S increased, and water was removed from the surrounding sediment.As a result the channel might remain open longer than predicted by our model and low stands would last longer.A more detailed exploration of this process including in situ monitoring of sediment pressures may be key to improving the performance of our model.
The sensitivity of our model output to small variations in these parameters (as well as α T ) calls into question our assumption that they remain constant with time.Ideally our model would fully incorporate a more complex array of glaciofluvial processes with spatial variation in sediment effective pressure and cohesion, all of which are related to till water content.Erosion, rather than being uniform along the channel, would then be concentrated in specific areas.Additionally, our experiments on R-channels suggest that melting, refreezing, and deformation of the overlying ice may change the channel significantly, especially in light of our results where the channels rarely cease once initiated.For this work, however, the demonstration that reasonable rates for erosion and sediment deformation can produce changes in channel geometry sufficient to simulate volume change inferred from observations acts as a first step.

Comparison against other observations of lakes in Antarctica
Observations show that the highstand of active lakes is lower than the flotation heights of their dams (Christianson et al., 2012;Siegert et al., 2014).In all of our model runs, lake drainage initiated when θ 0 was still well below θ 0 at the seal (Fig. 3).Our model is the first to capture this important dynamic for Antarctic subglacial lakes.We note that onset of outflow before lake levels reach the flotation height has been previously observed and modeled for icemarginal lakes in alpine glaciers (e.g., Fowler, 1999) and the Grímsvötn caldera lake beneath the Vatnajökull ice cap in Iceland (Björnsson, 2003;Evatt and Fowler, 2007).
Our results are also consistent with a hypothesis put forth in Fowler (2009) that the change in lake level (and thus pressure) may be better explained if channels were incised into deformable sediments rather than ice.This mode of channelization has two major implications: (i) the mechanism that allows active lakes to be inferred from satellite observations might make them difficult to detect with more traditional RES technology; conversely (ii) lakes detected by RES might drain via a different mechanism.If the rates of volume change inferred from satellite observations are better accommodated by channels incised in sediment, then it is likely that, in addition to the outlet of such a lake existing in soft sediments, the lake itself exists in soft sediments.The type of sediments that would both erode easily and quickly deform closed would also tend to have a low angle of repose such that any lake formed within them would be shallow.This physical setting may explain why these lakes are not detectable using RES (Gorman and Siegert, 1999), as a shallow lake in the presence of saturated sediment would fail the specularity test (Carter et al., 2007).Additionally, a lake surrounded by saturated sediments might not be significantly brighter than its surroundings as the reflection coefficient between ice and saturated sediments is very close to that between ice and lake water (Schroeder et al., 2013).
Our model for the R-channel suggests that such a mechanism dominates Antarctica subglacial lake drainage only in very specific circumstances.While Evatt et al. (2006) explored the possibility of Antarctic subglacial lake drainage via such a mechanism, the recurrence interval they found (thousands of years between events vs. a few decades of observations) would make it possible that we have yet to observe a lake draining via an R-channel.However, the drainage of Lake Cook E2 (Smith et al., 2009;McMillan et al., 2013;Flament et al., 2014), which resulted in a surface drawdown of over 50 m (similar in vertical scale to observed subglacial lake drainages in Iceland), would be a promising candidate for drainage via an R-channel.Further study of lake-drainage events outside of the Siple Coast is necessary to understand fully the respective roles of channels carved into ice and sediment.

Lakes in an ice-sheet model
A number of recent models for ice flow have started to predict the formation of subglacial lakes in local hydropotential minima (e.g., Sergienko and Hulbe, 2011;Goeller et al., 2013;Livingstone et al., 2013;Fried et al., 2014).These models all assume that these lakes simply fill until the water level reaches the flotation height of the "static seal", at which point they drain steadily through a distributed network at the ice-bed interface.The results of both our R-channel and canal models indicate that channelized drainage is likely important for many of these lake systems.In particular, the fluctuations in volumes inferred for lakes in fast-flowing regions of the ice sheet are likely the result of channelization.As a consequence, ice flow in areas where subglacial lakes are common may have very different flow properties relative to a region where subglacial water flows in steady state (Stearns et al., 2008;Siegfried et al., 2016).
Regions of the ice sheet underlain by active subglacial lakes will likely exhibit more variability in ice flow rate with peak ice speed coinciding with peak distributed flow and peak lake volume (Fig. 6).In chains of lakes, the lubrication may also be related to the evolution of the region hydrological system rather than an individual lake (e.g., Siegfried et al., 2016).The exact degree to which lake drainage accelerates ice flow is, however, highly dependent on longitudinal stresses, which is outside the scope of this paper.
Spatial variations in basal traction in areas of fast flow have previously been proposed as a mechanism for forming lakes (Sergienko and Hulbe, 2011), with lakes forming in the lee of local traction highs.More recently inversions of basal traction have inferred bands of stiff sediment impounding water in areas of moderate to fast flow, but usually outside the regions where active lakes are found (Sergienko et al., 2014;Hulbe et al., 2016).It may be that these categories represent two different mechanisms of water storage, one stable and the other unstable, resulting from subtle changes in sediment properties.Along with findings of multiple other mechanisms of water storage beneath the Antarctic Ice Sheet (Ashmore and Bingham, 2014), these conclusions coupled with our results indicate that simply modeling the filling of an enclosed hydropotential depression, while an important first step, is not sufficient to simulate the full nature and impact of lake dynamics in an ice-sheet model.

Summary
We have developed a new model for subglacial lake drainage beneath Antarctic ice streams in which channels are mechanically eroded into deformable subglacial sediment and compared it to the lake-volume change predicted by the more conventional "R-channel" model, in which water incises into the base of the ice through melting.The "canal" model was better able to reproduce the timing and magnitude of subglacial lake-drainage events in the well-studied Whillans-Mercer ice-stream system.Due to effective-pressure differences between the lake center and the "seal", lake drainage begins well before lake levels reaches flotation height and accelerates once flow is sufficient to initiate erosion into the sediment.Peak distributed flow correlates with lake level, while channelization is dominant when the rate of volume loss is highest.This evolution in the drainage system and its impact on water pressures will likely influence local and regional ice dynamics.
The time series output by the model are highly sensitive to small changes in the properties of the subglacial sediments, in particular those related to sediment water content.Given recent work on the exchange of water between sediment pores and the interfacial flow system taking place over time frames of years to decades, it is likely that there are substantial changes in sediment strength over the filling or drainage cycle of a subglacial lake.It may even be the case that periodic filling and drainage of subglacial lakes requires very specific sediment properties.Intriguingly, recent analysis of sediment in lower Whillans Ice Stream has not found substantial evidence for large-scale fluvial deposition suggesting that the erosion necessary to create subglacial channels occurs only in limited places (Hodson et al., 2016) The requirement of sediments for the realistic drainage of active lakes may explain why such lakes often fail classic RES detection criteria as lakes in such environments are often neither deep (and therefore not specular; Gorman and Siegert, 1999) nor significantly brighter than their surroundings (Carter et al., 2007).Conversely lakes detected by RES are typically found in steep-sided bedrock troughs (Wright and Siegert, 2012) and their surroundings likely lack the deformable sediments that would facilitate episodic drainage on the time frames for which observations exist.However, occasional drainage of such lakes via R-channels remains a possibility.For all lakes, the mechanisms for storage and release of water respond more to internal variability of the ice sheet than they do to external forcing.Consequently, the mechanism by which they drain and the effect of this mechanism on ice flow remains a critical uncertainty in predicting future ice-sheet evolution.
The Supplement related to this article is available online at doi:10.5194/tc-11-381-2017-supplement.

Figure 2 .
Figure 2. Schematic diagram showing the various stages of lake activity (right panels) and the corresponding evolution of the channel (left panels).Red arrow denotes net erosion rate, white arrow denotes closure rate, brown arrow denotes canal effective pressure, and blue bar indicates hydropotential in the distributed system:(a) initially the lake is filling and no outflow takes place; (b) as lake level increases, effective pressure differences between lake and seal allow for some water to escape via a distributed system (a.k.a."water sheet"); (c) as the distributed system grows, a low pressure channel is eroded that begins to draw water from the surroundings; (d) while the low pressure channel allows for the lake to drain below the level at which outflow was initiated, as pressure lowers, sediment creeps inward; and (e) the channel contracts until its flow is negligible and the lake refills.
which describes creep closure as a function of effective pressure (N RC ) and aperture size (S RC ), where K RC and n are constants of Glen's flow law.Conservation of mass governs the change in flux along the flow path by

Figure 3 .
Figure 3. Along-flow hydropotential (left axis) in the distributed (solid) and channelized (dotted) drainage systems at four important time steps during a fill-drain cycle (initial state, highstand, peak outflow, low stand) for the four main lakes in the Whillans-Mercer subglacial system: (a) SLW, (b) SLE, (c) SLC, and (d) SLM.Thick black line (right axis) indicates overburden pressure calculated using Eq.(1a).

Figure 4 .
Figure 4. Ice surface elevation, hydropotential, and ice-bed elevation along the following flow paths (see Fig. 1b for location of flow paths): (a) an idealized lake flow path based on SLW, (b) the flow path draining SLW (from Carter and Fricker, 2012), (c) the flow path draining SLC and SLM, and (d) the flow path draining SLE.

Figure 5 .
Figure 5. Results of experiments with R-channel model on the SLW domain (see Fig. 1b for location and Fig. 4b for flow path geometry).(a) Modeled lake elevation anomaly and outflow obtained for the channelized and distributed drainage systems; observed lake elevations (Carter et al., 2013) are also shown.(b) Net melt rate along the flow path for five time steps shown in (a), with gray boxes denoting areas of net freezing; (c) viscous closure rate (C VRC ).(d) Modeled lake elevation anomalies resulting from varying parameters controlling channel closure (K RC ) and opening (L h ).

Figure 6 .
Figure 6.Results of experiments with canal model on the real domain (Fig. 1b), showing evolution of modeled lake volume and outflow via distributed (gray) and channelized (blue) systems for constant Q in for the four main lakes: (a) SLW, (b) SLE, (c) SLC, and (d) SLM.Also shown are the observed lake volumes from Carter et al. (2013).Colored vertical lines correspond to same time steps from Fig. 5.

Figure 7 .
Figure7.Comparison between canal model executed on the idealized and realistic domains.From top to bottom, pairs of panels show (a) modeled volume change and (b) modeled outflow; evolution of hydropotential in the canal (θ CC ) and in the distributed system (θ S ) compared to ice base hydropotential (θ 0 ) on the (c) idealized and (d) realistic domains; and evolution of cross-sectional area of the distributed system (S S ) on the (e) idealized and (f) realistic domains.

Figure 8 .
Figure 8. Sensitivity of canal model to flow path geometry in the idealized domain.Left panels (a, b, c) show sensitivity of model output to changes in vertical scaling: (a) vertical exaggeration applied to the entire domain; (b) vertical exaggeration applied upstream of seal; and (c) vertical exaggeration applied downstream of seal.Right panels (d, e, f) show sensitivity of model output to changes in horizontal scaling: (d) horizontal stretching applied to the entire domain; (e) horizontal stretching applied upstream of seal; and (f) horizontal stretching applied downstream of seal.

Figure 10 .
Figure10.Sensitivity of canal model to sediment and channel properties.Variation of elevation anomaly (z sL ) due to five different values of (a) channel geometry (α CC ), (b) sediment effective pressure (N ∞ ), and (c) mean sediment grain size (d 15 ).The black lines denote the "control" run using our preferred parameter values (Table 2).

Figure 11 .
Figure 11.Sensitivity of canal model to (a) mean obstacle size (R 1 ) and (b) mean basal melt rate (M S ).

Table 2 .
List of parameters used for each of the experiments."n/a" stands for not applicable.
criterion for the length of the time step.After calculating dS S /dt and dS RC /dt (or dS CC /dt) we then define the time step t by the following relation: