Introduction
More than 145 subglacial lakes have been identified underneath the Antarctic ice sheet, mainly by analyzing airborne radio-echo sounding profiles. Subglacial lakes are characterized by a strong basal reflector, a constant echo strength (corroborating a smooth surface) and a flat surface compared to the surroundings with a slope around ten times, and in the opposite direction to, the lake surface slope (Reference Siegert, Carter, Tabacco, Popov and BlankenshipSiegert and others, 2005). The ice column above a subglacial lake is therefore in hydrostatic equilibrium. Most lakes lie under a thick ice cover of >3500 m and are therefore situated close to ice divides. Although subglacial lakes are usually small (<20 km in length), the largest, Vostok Subglacial Lake, contains 5400 ± 1600 km3 of water (Reference Studinger, Bell and TikkuStudinger and others, 2004). Assuming an average water depth of 1000 m for large lakes and 100 m for shallow lakes, the known volume of water in Antarctic subglacial lakes is ∼22 000 km3 or approximately 25% of the global water in surface lakes (Subglacial Antarctic Lake Environments (SALE) Workshop report, http://salepo.tamu.edu/saleworkshop2006). This is equivalent to a uniform sheet of water ∼1 m thick if spread out beneath the entire Antarctic ice sheet.
Subglacial lakes are defined as relatively closed and stable environments with long residence times and slow circulations (Reference Kapitsa, Ridley, Robin, Siegert and ZotikovKapitsa and others, 1996; Reference SiegertSiegert and others, 2001; Reference Bell, Studinger, Tikku, Clarke, Gutner and MeertensBell and others, 2002). Such a huge amount of subglacial meltwater would pose no threat to the stability of the ice sheet if it were not moving through a hydrological network. This view has recently been challenged through the observation of rapid subglacial lake drainage events (Reference Gray, Joughin, Tulaczyk, Spikes, Bindschadler and JezekGray and others, 2005; Reference Wingham, Siegert, Shepherd and MuirWingham and others, 2006; Reference Fricker, Scambos, Bindschadler and PadmanFricker and others, 2007). The observations by Reference Wingham, Siegert, Shepherd and MuirWingham and others (2006) suggest the rapid transfer of subglacial lake water and the periodic flushing of subglacial lakes connected to other lakes that consequently fill through a hydrological network. For the draining lake (Adventure Trench Lake), they observe a surface lowering of 4 m, corresponding to a water discharge of ∼1.8 km3 over a period of 16 months at a rate of 50 m3 s-1. The analysis by Reference Fricker, Scambos, Bindschadler and PadmanFricker and others (2007) points to another event of similar magnitude, but underneath an ice stream (Whillans Ice Stream, West Antarctica, just upstream of the grounding line). In this case, a surface lowering of 9 m is observed over Lake Engelhardt, corresponding to a discharge of ∼1.25 km3 over a period of 18 months at a rate of 40 m3 s-1. The events are therefore similar in timing and magnitude.
The above-mentioned observations stem from rapid changes at the surface from which the drainage events are inferred. Direct observation of such subglacial events, and hence knowledge of their triggering mechanisms, is lacking. However, one way to investigate possible mechanisms leading to rapid subglacial drainage is through ice-sheet modeling in which ice flow across and interactions with a subglacial lake are taken into account. A full Stokes ice-sheet model, taking into account all stress components acting on the ice, is probably the most suitable tool for this (e.g. Martín and others, 2003; Reference Zwinger, Greve, Gagliardini, Shiraiwa and LylyZwinger and others, 2007). Therefore, such a model is developed and coupled to the subglacial lake environment to investigate the stability of subglacial lakes and their drainage through a sensitivity analysis. It is shown that the coupled ice-sheet/lake system, due to hydrostatic equilibrium, exhibits a feedback mechanism that leads to periodic lake flushing on decadal timescales.
Model Description
In the ice-flow model, the different surfaces are expressed with a right-handed (0, x, y, z) Cartesian coordinate system as depicted in Figure 1. The ice-surface and ice-bottom Cartesian representations are defined as z = S(x, y, t) and z = B ice(x, y, t), respectively, from which it follows that ice thickness equals H(x, y, t) = S(x, y, t) − B ice(x, y, t). The underlying bedrock is defined by z = B(x, y). For ice resting on the bedrock it follows that B ice(x, y, t) = B(x, y). When a subglacial lake is present, however, B ice(x, y, t) > B(x, y) and the water depth of the lake is defined by W(x, y, t) = B ice(x, y, t) − B(x, y).
Fig. 1. Ice-sheet surfaces expressed in a Cartesian coordinate system.
The model approach is based on continuum modeling and encompasses balance laws of mass and momentum, extended with a constitutive equation. Treating ice as an incompressible fluid with constant density, the equation for conservation of mass implies that
where (vx
, vy
, vz
) denote the respective x, y, z coordinates of the velocity vector v and 
are the normal strain rates. Ice incompressibility implies that the stress tensor must be split into a deviatoric part and an isotropic pressure P, i.e.
The constitutive equation for ice then links deviatoric stresses to strain rates:
where T′ and 
are the deviatoric stress and strain-rate tensor, respectively, and η is the effective viscosity.
Stokes’ equations
Neglecting acceleration terms, the linear momentum balance is written as:
where g is the gravitational acceleration (9.81 m s−2) and ρ ice is the ice density (910 kg m−3). Since only acceleration due to gravity in the vertical is considered,
By definition, the isotropic pressure is equal to P = −1/3 Σ i τii . Making use of Equation (2), we derive
An expression for the vertical normal stress τzz is obtained through vertical integration of Equation (7) from the surface S to a height z in the ice mass:
where the dynamic boundary condition| (17) is invoked to eliminate the contribution of the free surface. As such, P can be eliminated from Equation (7) and the force-balance equations further reduced to two equations in the horizontal direction (see below).
Rheological non-linearities are ignored, which is an acceptable approximation as long as the stress perturbations remain small enough for the flow law to be linearized around the mean stress field (Reference GudmundssonGudmundsson, 2003). A few experiments were performed with values of n = 3, i.e. the exponent in Glen’s flow law 
. Although the timing of drainage events was slightly different in these cases, due to the different parameter settings, the overall sensitivity of the system to parameter perturbations was the same.
For a linear flow law relating deviatoric stresses to strain rates, the effective viscosity in Equation (3) equals
The flow rate factor A is temperature-dependent (Reference PatersonPaterson, 1994). For this study, temperature–flow coupling is not taken into account and ice is considered to be at the pressure-melting point at the base all the way (a prerequisite for having subglacial lakes). Nevertheless, the flow parameter A is used as a major sensitivity parameter to investigate a large range of ice viscosities. In reality, however, thermomechanical coupling would allow for changes in viscosity over the vertical, hence influencing the vertical profile of the velocity field, with stiff ice on top of a more deformable ice layer. This would certainly result in different timings of subglacial lake drainage but not alter the trends in sensitivity of the system to parameter perturbations, explored later.
Velocity field
Combining the force-balance Equations (5–7) with the constitutive Equation (3) and replacing the isotropic pressure by Equation (8) leads to
and
Strain rates are related to velocity gradients by definition, so
Combining Equation (13) with Equations (11) and (12) leads to an expression for the horizontal velocity field, which is presented in Appendix A. An expression for the vertical velocity vz is obtained through vertical integration of the incompressibility condition (1).
Boundary conditions
Stress-free conditions apply at the upper surface. When neglecting atmospheric pressure, this leads to
Written in terms of deviatoric stress components, this becomes
Boundary conditions at the base of the ice mass are more complicated, as they involve either the presence of water (subglacial lake) or basal sliding over the bedrock. In the case of a subglacial lake, boundary conditions are
where nB is the outward-pointing unit normal vector to the bottom surface of the ice sheet (Reference Mayer and SiegertMayer and Siegert, 2000). The assumption is made that the pressure of the lake water at the bottom of the ice shelf is equal to the hydrostatic pressure necessary to float the ice mass on top, i.e. the ice sheet is floating in a medium with a hydrostatic pressure field. Written in terms of deviatoric stress components, this becomes
Basal sliding occurs outside the subglacial lake region, which is defined by
where F is a basal sliding coefficient and τ B is the basal drag or the sum of all basal resistance. Following the analysis by Reference Muszynski and BirchfieldMuszynski and Birchfield (1987) and by combining Equations (19) and (21), the sliding law can be written as:
For the series of experiments carried out in this paper, a viscous sliding law is assumed, implying m = 1. However, since basal sliding is not a dominant process here and inland ice flow is considered, the type of basal sliding law (viscous or plastic) does not really matter. A linear sliding law simplifies the computation.
Basal hydrology
Full basal hydraulics are not implemented in the model as such. A basal hydrological algorithm is used to check whether the lake ‘seal’ breaks or not, based on the hydraulic potential gradient described in Appendix B, but active water transport across the domain is not considered. Other hydrological features, such as tunnel formation and closure and associated water discharge (Reference NyeNye, 1976), are not taken into account either but are discussed at the end of this paper. Nevertheless, as the lake empties, the marginal ice should collapse or fracture to form a cauldron, a common feature of Jökulhlaups. The rate of subsidence is directly related to the rate of lake discharge, but requires an underpressure in the lake to draw down the overlying ice (Reference Evatt, Fowler, Clark and HultonEvatt and others, 2006; Reference Evatt and FowlerEvatt and Fowler, 2007). However, as small floods creating small surface depressions are considered here, only a flotation criterion for the ice over the lake is considered, as in Reference NyeNye (1976).
Ice-sheet evolution
Using a kinematic boundary condition at the upper and lower surface of the ice mass, the mass conservation Equation (1) is integrated along the vertical in order to obtain an expression for the change of local ice thickness in space:
Local ice-thickness variations are therefore in balance with the divergence of the ice flux. Surface and bottom accumulation, ablation and melting are not considered.
Numerical solution
The equations are solved on a regular finite-difference grid in the horizontal and an irregular grid in the vertical. First and second central-difference approximations on a regular and irregular grid are given in Reference PattynPattyn (2002, Reference Pattyn2003). At the horizontal domain boundaries, periodic boundary conditions are applied.
The numerical solution of the velocity field is obtained by following the method described in Reference PattynPattyn (2003). Equations (A1) and (A2) in Appendix A are solved as a sparse linear system using the conjugate gradient method, whereas the non-linear parts of Equations (A1) and (A2), represented by the righthand sides, were solved using the unstable manifold correction (Reference Hindmarsh and PayneHindmarsh and Payne, 1996; Reference PattynPattyn, 2002, Reference Pattyn2003). The continuity Equation (25) is reformulated as a diffusion equation for ice thickness H and solved using a semi-implicit method, also making use of sparse matrix algorithms.
Conservation of water volume in the subglacial lake implies the definition of a buoyancy level (for a floating ice shelf this is sea level, but for a subglacial lake this is defined as B ice + Hρ ice /ρ w). In an iterative procedure, the buoyancy level is determined for the local ice thickness to be in hydrostatic equilibrium on top of the lake for a given water volume. This procedure determines the position of both the upper and lower surfaces of the ice sheet across the subglacial lake in the Cartesian coordinate system. The contact surface between the ice sheet and the lake is then set to β 2 = 0.
Model set-up
The basic model set-up is an idealized subglacial lake underneath an idealized ice sheet (Fig. 2). The lake is defined as a Gaussian cavity. Initial upper and lower surfaces are defined by:
where α is the general slope, H 0 is the initial ice thickness, C 0 is the depth of the water cavity (taken as 400 m), σ is the approximate radius of the lake, and x 0, y 0 the coordinates of the center of the domain. Model parameters for the standard experiment and the sensitivity analysis are given in Table 1. Initially, and according to the standard experiment, the lake cavity was filled with 40 × 109 m3 of water. The horizontal domain is L × L, where L = 80 km, and a grid size of 2 km was used (order of magnitude of ice thickness). Model runs were performed with a time-step of 0.04 years. In order to determine the effect of domain size and the influence of the periodic boundary conditions on the model results, preliminary experiments were conducted with different sizes of model domain up to 160 km × 160 km. It appeared that the model behavior was not influenced by this choice.
Table 1. Standard model parameters and sensitivity values
Results
Subglacial lake stability
Beginning from the initial conditions, the model ran forward until a steady state was reached (Fig. 3). The general characteristics of the resulting ice-sheet geometry are typical for those of a slippery spot, as is a subglacial lake (Reference PattynPattyn, 2003, Reference Pattyn2004; Reference Pattyn, De Smedt and SouchezPattyn and others, 2004), i.e. a flattened surface of the ice–air interface across the lake and the tilted lake ceiling in the opposite direction to the surface slope, due to hydrostatic equilibrium. The tilt of this surface in the direction of the ice flow will determine the stability of the lake, since the hydraulic gradient is dominated by surface slopes and therefore the flatter this air–ice interface, the more easily water is kept inside the lake cavity. Changes in the surface flatness will therefore be crucial in subglacial lake stability. The purpose of the sensitivity analysis is to determine this stability as a function of model parameters (Table 1).
Fig. 3. Steady-state geometry calculated using the model. The two planes represent the upper and lower surface of the ice mass. A cross-section is shown to the right. Note the flattened surface of the ice above the lake and the tilted ice–water interface due to hydrostatic equilibrium. Ice thickness has been reduced to 300 m only for better rendering. Vertical exaggeration is a factor 100 in the three-dimensional plot, as the horizontally scaled distance corresponds to 80 km for the standard experiment.
Figure 4 summarizes the results of the sensitivity analysis, displaying the ice surface slope across the lake in the direction of the ice flow as a function of parameter settings. It follows that lake stability is prevalent for large lakes, shallow ice, small ice surface slope and fast ice flow. However, the most dominant parameter here is general surface slope. Neither Adventure Trench Lake nor Engelhardt Lake is the end member of stability or instability, albeit that Adventure Trench Lake might be regarded as a more stable situation. Small lakes tend to be less stable than large lakes, which implies that the lake region in the area around Dome Concordia, East Antarctica, should be more vulnerable in terms of drainage as it consists mainly of relatively small subglacial lakes. On the other hand, large lakes such as Vostok Sub-glacial Lake should be very stable and exhibit a pronounced flatness at the surface. In the next series of experiments, the system is perturbed to estimate its sensitivity to drainage.
Fig. 4. Steady-state modeled surface slope across the lake area for different parameter settings. ATL and EL correspond to Adventure Trench and Engelhardt Lake geometries, respectively. Best-fit lines demonstrate the tendency of the relationship. Solid circles and lines correspond to results of the standard experiments. Open circles and dotted lines refer to perturbation experiment P3, where 20% of the water volume has been removed from the lake.
Lake drainage
Our hypothesis is that (partial) lake drainage is a common feature of the Antarctic subglacial lakes and that only small perturbations are necessary to lead to a sudden outburst. Three perturbation experiments were therefore carried out: experiment P1 where the surface of the ice sheet is slightly perturbed through a perturbation in ice thickness; experiment P2 in which subglacial water is gradually added to the lake system; and P3 where 20% of the lake water volume is removed at once. The thickness perturbation according to P1 is defined by:
where L is the length of the domain, ΔH the amount of ice added to or subtracted from the surface every time-step and Δt the time-step (1 year). Perturbation P1 slightly changes the ice-flow pattern and can be regarded as a perturbation in accumulation pattern or gradual changes of the interior ice sheet as a response to the glacial/interglacial contrast, which might still occur today (Reference HuybrechtsHuybrechts, 2002). In P2, water is added to the lake at a rate of 0.1 m3 s−1. This corresponds to melting at the ice/lake ceiling at a rate of 2.5 mm a-1 for the standard experiment, which is of the same magnitude as basal melting occurring under the ice sheets. Thus, both P1 and P2 are regarded as common gradual changes in the glacial environment which are hardly noticeable on decadal timescales, hence forming part of the natural variability of the system.
Once the drainage condition is fulfilled (i.e. the hydrostatic seal is broken), the lake is drained at a rate of 50 m3 s−1 for a period of 16 months, a value given by Reference Wingham, Siegert, Shepherd and MuirWingham and others (2006) and which is of the same order of magnitude as the drainage rate of Lake Engelhardt (Reference Fricker, Scambos, Bindschadler and PadmanFricker and others, 2007). This essentially means that enough energy is released to keep the subglacial tunnel open for a while and a siphon effect occurs. It is not clear whether this effect is temporary, or ends when the lake is completely emptied as suggested by Reference Wingham, Siegert, Shepherd and MuirWingham and others (2006). However, the drainage event of Lake Engelhardt suggests that water is still present in the lake after the event, as the surface aspect (flatness) has not changed over this period of time. Furthermore, Reference Evatt and FowlerEvatt and Fowler (2007) suggest that the drainage event of Adventure Trench Lake is due to channelized flow over soft sediments at low effective pressure, rather than the emptying of a very shallow lake.
The time-dependent response of the lake system to such small perturbations for P1 is shown in Figure 5. Similar results are obtained for other sets of parameter values and also P2. Not all experiments exhibit lake drainage; certain parameter settings, such as small ice thickness, render lakes very stable (Fig. 5). When drainage occurs episodic events take place, even though the initial geometric conditions have not been met. These events occur at a higher frequency at the beginning, i.e. when more water is present in the lake, than later on and occur with frequencies of less than a decade. Most experiments show stable conditions after 200 years of integration, either because the drainage conditions are not fulfilled anymore, or because the cavity is devoid of water and a complete drainage has taken place. The episodic drainage results in variations in the surface elevation across the lake; specifically a rapid lowering followed by a gradual increase until the initial level is more or less reached. This surface rise is clearly non-linear, and the rate of uplift decreases with time. The frequency also decreases with time, as less water is present in the system, hence hampering drainage. It is interesting to note, however, that the surface rise is not caused by an influx of subglacial water entering from upstream. On the contrary, less water is subsequently present in the whole subglacial system. Surface increase is due to an increased ice flux, filling up the surface depression created by the sudden lake drainage. Evidently, the maximum speed of surface increase is significantly less than the surface lowering due to drainage. After the drainage event, surface slope slightly decreases before increasing again, when upstream ice flows in, as depicted in the phase space diagram (Fig. 6).
Fig. 5. Time series of the perturbation experiment according to Equation (28). Time evolution of (a) lake water volume; (b) ice surface elevation on top of the lake; and (c) ice surface slope across the lake for the ice thickness sensitivity. The sensitivity according to P2 provides similar results. No lake drainage is observed for H 0 =1000 m.
Fig. 6. Phase space representation of the evolution of surface slope vs ice-sheet elevation across the lake according to experiment P1 with H = 2500 m. Each star is plotted at equal time-step intervals of 1 year. Time integration starts in the upper left corner of the graph.
Figure 7 shows the lake volume after 200 years of model integration for different parameter settings and both perturbation experiments. Results for P1 and P2 are more or less the same, indicating that the perturbation is only a minor trigger to initiate episodic lake drainage. If one assumes the lake volume after 200 years is a measure for subglacial lake stability, then Figure 7 can be compared directly with Figure 4. The relationships are somehow different compared to the previous stability test: large lakes in slow-moving areas under a thick ice cover drain much more easily. The relationship is inverted for small lakes, which are more stable. Large lakes can drain more easily, despite their pronounced surface flatness.
Fig. 7. Lake water volume after 200 years of calculation for perturbation experiments P1 (solid circles and lines) and P2 (open circles and dotted lines). Lines represent the best fit to show the tendency of the relationship.
A third perturbation experiment, P3, involves the sudden removal of 20% of the lake volume (Fig. 4). The parameter/slope relations are similar to those of the standard experiment. However, the slope values are significantly higher due to a reduced amount of water in the lake cavity, hence a smaller contact area between the ice sheet and the lake. As seen from the initial sensitivity experiment, smaller lakes tend towards higher surface slopes. This also explains the episodic drainage events in phase space (Fig. 6). After the lake drains for the first time, the surface slope adjusts and reaches a critical state at a higher value before the hydrostatic seal breaks again, since less water is present in the system and the lake surface is smaller. As seen from the initial sensitivity test (Fig. 4), smaller lakes lead to steeper surfaces. This also accounts for ice thickness; the ice is thicker across the lake after a drainage event, since more ice flowed into the cavity due to the water release. Thicker ice also increases the surface slope (Fig. 4), leading to a positive feedback favoring subsequent lake drainage.
Discussion
Effect of lake shape on drainage
Only experiments on circular-shaped lakes are presented in this paper. In reality however, many lakes lying in sub-glacial trenches have an elongated form, with the tilting ice–water slope in the direction of the main ice flow. Moreover, the preconditioned Gaussian shape also influences the response to perturbations. The removal of a constant quantity of water will lead to changes in the bed slope by lowering water levels, which inevitably affects the slope of the ice–water interface and influences the trigger mechanism. Figure 8 presents time series for different lake forms, i.e. an elongated Gaussian form in the direction of as well as perpendicular to the ice flow, and a conical-shaped lake characterized by a constant bed slope, defined by:
within the limit range of 
Outside this range, B(x, y, 0) = S(x, y, 0) − H
0.
Fig. 8. Time series of the perturbation experiment depicted in Figure 5, for different types of lake cavities: (a) elongated Gaussian shape in the direction of ice flow; (b) elongated Gaussian shape perpendicular to ice flow; and (c) conic lake form.
All of these lake shapes exhibit episodic drainage, but the drainage behavior differs significantly. Elongated lakes in the ice-flow direction are relatively unstable compared to those perpendicular to the ice flow, and their behavior is similar to that of larger lakes. Stability is therefore governed by the lake dimension in the direction of the ice flow. The conic shape exhibits a more regular drainage pattern over time, as the bed slope remains constant and the only parameter that will account for drainage is the surface slope variability.
Origin of subglacial lakes
Reference Siegert, Le Brock, Payne, Hambrey, Christoffersen, Glasser and HubbardSiegert and others (2007) show that most subglacial lakes are lying close to ice divides and therefore in areas of low surface slopes and high hydraulic fluid potential. A fundamental question arises (Priscu and others, 2007): does the evolving ice sheet control the location of subglacial lakes or does the lithospheric character necessary for lake formation, such as geothermal heat flux and the presence of basal cavities and sub-ice aquifers, constrain the evolution of these catchments?
According to the model simulations, geometric effects such as ice thickness, mean surface slope, ice viscosity (hence ice temperature) and lake size play an important role in the stability of subglacial lakes. They favor the a priori existence of subglacial lakes near ice divides in areas of low viscosity and low surface slopes. Lake filling might therefore happen when the ice sheet is in full expansion during sustained cold periods and when surface slopes are very low, such as a glacial period. During interglacial periods, the smaller and warmer ice sheet characterized by slightly higher surface slopes in the interior facilitates lake drainage. This may explain why subglacial lakes are not observed below the Greenland ice sheet. Several reasons can be put forward, such as the lack of suitable subglacial basins or cavities to collect the subglacial water. Another factor might be given by geometry: the Greenland ice sheet is generally warmer than the Antarctic ice sheet (hence lower viscosity), thinner at the interior and characterized by higher surface slopes as the ice sheet is substantially smaller than its Antarctic counterpart. These factors lead to increased lake instability and favor rapid drainage, which may have occurred at the end of the last glacial period at the beginning of the Holocene.
Nye theory for lake discharge
Nye’s hydraulic theory offers an appropriate description of the drainage of subglacial lakes, especially where the lake is not connected to the ice surface and the ice over the lake responds dynamically to lake level fluctuations (Reference NyeNye, 1976). According to that model, the magnitude and duration of floods appear to be controlled by the channel hydraulics, but the peak discharge is essentially related to lake volume, i.e. large lakes produce large floods (Reference Evatt, Fowler, Clark and HultonEvatt and others, 2006). Therefore, small floods apparently observed by Reference Wingham, Siegert, Shepherd and MuirWingham and others (2006) are associated with drainage through channels at low effective pressure.
While the Nye theory offers the possibility to predict the flood discharge, the present study focuses on the trigger mechanism for lake discharge, i.e. when the hydrostatic seal is broken. It is therefore supposed that once the hydrostatic seal is broken, a subglacial channel is established in which the water pressure is lower than the ice pressure. However, in order to keep the channel from closing through viscous creep of ice, closure should be balanced by the melting of the channel walls. This energy is supplied by viscous dissipation of the turbulent water in the channel. As long as creep closure is balanced by the dissipated heat, the channel will remain open and discharge can take place. Reference Evatt, Fowler, Clark and HultonEvatt and others (2006) have shown that lake drainage events, as described by Reference Wingham, Siegert, Shepherd and MuirWingham and others (2006), can occur over sustained periods of tens of months without the complete emptying of the lake. Furthermore, they also corroborate the fact that lake drainage is a common mechanism underneath the Antarctic ice sheet.
Conclusions
An existing three-dimensional higher-order ice-sheet model (Reference PattynPattyn, 2003) has been modified and converted to a full Stokes model, further extended with a model for subglacial water storage and release that is in hydrostatic equilibrium with the overlying ice sheet. Although the model does not take into account thermomechanical effects nor the Nye theory for tunnel formation, it is capable of exhibiting the major characteristics of drainage due to changes in ice-sheet and/or subglacial lake geometry.
Results demonstrate that only small changes are necessary to provoke episodic lake drainage. These events only lead to a partial drainage of the lake, after which the ice sheet adapts to the changed conditions and ice flow ‘fills’ the surface depression after drainage. Lake drainage on a decadal timescale can therefore be regarded as a common feature of the subglacial hydrological system and may influence, to a large extent, the behavior of large ice sheets through such jökulhlaups (Reference Evatt, Fowler, Clark and HultonEvatt and others, 2006). This has potential implications for ice-stream dynamics, especially in areas that are underlain by sedimentary material. Unless this water escapes to the margin via well-defined channels, it will act to increase pore-water pressures, reducing the strength of sediments, hence increasing ice velocity (Reference Siegert, Le Brock, Payne, Hambrey, Christoffersen, Glasser and HubbardSiegert and others, 2007). Finally, the experiments show that large lakes under thick ice in slow-moving regions drain easily, which should be a common occurrence for Vostok Subglacial Lake.
Acknowledgements
This paper forms a contribution to the Belgian Research Program on the Antarctic (Belgian Federal Science Policy Office), project No. SD/CA/02 ‘Antarctic Subglacial Processes and Interactions: role of transition zones in ice sheet stability (ASPI)’. Thanks to M. Siegert, R. Bell and A. Salamatin for fruitful and illuminating discussions and to R. Greve, F. Saito and T. Zwinger for helpful reviews.
Appendix A
A complete derivation of the horizontal velocity field is obtained by combining Equation (13) with Equations (11) and (12). After some algebraic manipulation, this leads to
and
which in the case of a linear rheology reduces to
and
An expression for the boundary conditions at the upper surface (stress-free) is obtained by employing Equation (17) and expressing deviatoric stresses in terms of velocity gradients. The horizontal components Equations (15) and (16) therefore transform to:
and
An expression for the boundary conditions at the lower surface in the case of the presence of a subglacial lake is obtained by employing Equation (21) and expressing deviatoric stresses in terms of velocity gradients. The horizontal components Equations (19) and (20) therefore transform to:
and
When the ice does not flow across a subglacial lake, a sliding-law relationship holds, i.e. 
. A complete expression for the basal velocity field is obtained by replacing the deviatoric stresses by velocity gradients in Equations (23) and (24), so that both equations transform to:
and
where 
is a friction parameter. From Equations (A8) and (A9) it follows that when β
2 = 0, the base of the ice sheet experiences no friction and they reduce to Equations (A6) and (A7), respectively.
Appendix B
In the model, basal hydrology is introduced to check whether subglacial water captured in the lake cavity is transported outside this cavity or not, and therefore serves to detect whether the hydrostatic seal is intact or broken. It does not play an active part in the lake drainage. A continuity equation for basal water flow is written as
where w′ is the thickness of the subglacial water layer outside the lake (m), v
w the vertically integrated velocity of water in that layer (m a−1) and w the potential water release rate of the lake (m a−1). Since water that leaves the lake flows rapidly to the downstream edge of the domain, we can assume that the system is in steady state, or 
.
Subglacial water tends to move in the direction of decreasing hydraulic potential φ (Reference ShreveShreve, 1972). If we assume that basal water pressure is equal to the ice pressure, the hydraulic potential gradient (or water-pressure gradient) is written as ∇φ = ρ wg∇B ice − ∇P, where ρ w is the lake water density, assumed to be 1000 kg m−3, B ice is the lower surface elevation of the ice mass and P is the ice pressure. The steady-state basal water flux ψ w = (v w w′) is obtained by integrating w over the domain, starting at the hydraulic head in the direction of the negative of the hydraulic potential gradient ∇φ. This is done with the computer scheme for balance-flux distribution of ice sheets described in detail by Reference Budd and WarnerBudd and Warner (1996). The hydrostatic seal is thus broken when water originating from the lake potentially arrives at the downstream edge of the domain. Once this occurs, lake drainage begins at the rate of 50 m3 s−1 for a period of 16 months.
