The influence of a model subglacial lake on ice dynamics and internal layering

As ice flows over a subglacial lake, the drop in bed resistance leads to an increase in ice velocities and a draw down of isochrones and cold ice. The ice surface flattens as it adjusts to the lack of resisting forces at the base. The rapid transition in velocity induces changes in ice viscosity and releases deformation energy that can raise the temperature locally. Recent studies of Antarctic subglacial lakes indicate that many lakes experience very fast and possibly episodic drainage, during which the lake size is rapidly reduced as water flows out. Questions that arise are what effect this would have on internal layers within the ice and whether such past drainage events could be inferred from isochrone structures downstream. Here, we study the effect of a subglacial lake on ice dynamics as well as the influence that such short timescale drainage would have on the internal layers of the ice. To this end, we use a full Stokes, polythermal ice flow model. An enthalpy-gradient method is used to account for the evolution of temperature and water content within the ice. We find that a rapid transition between slow-moving ice outside the lake, and full sliding over the lake, can release considerable amounts of deformational energy, with the potential to form a temperate layer at depth in the transition zone. In addition, we provide an explanation for a characteristic surface feature commonly seen at the edges of subglacial lakes, a hummocky surface depression in the transition zone between little to full sliding. We also conclude that rapid changes in the horizontal extent of subglacial lakes and slippery patches, compared to the average ice column velocity, can create a traveling wave at depth within the isochrone structure that transfers downstream with the advection of ice, thus indicating the possibility of detecting past drainage events with ice penetrating radar.

locally. Recent studies of Antarctic subglacial lakes indicate that many lakes experience very fast and possibly episodic drainage, during which the lake size is rapidly reduced as water flows out. A question that arises is what effect this would have on internal layers within the ice, and whether such past drainage events could be inferred from isochrone structures downstream. 10 Here, we study the effect of a subglacial lake on ice dynamics as well as the influence that such short timescale drainage would have on the internal layers of the ice. To this end, we use a full-Stokes, polythermal ice flow model. An enthalpy gradient method is used to account for the evolution of temperature and water content within the ice.
We find that a rapid transition between slow-moving ice outside the lake, and full sliding 15 over the lake, can release considerable amounts of deformational energy, with the potential to form a temperate layer at depth in the transition zone. In addition, we provide an explanation for a characteristic surface feature commonly seen at the edges of subglacial lakes, a hummocky surface depression in the transition zone between little to full sliding. We also conclude that rapid changes in the horizontal extent of subglacial lakes and slip- 20 pery patches, compared to the average ice column velocity, can create a travelling wave at depth within the isochrone structure that transfers downstream with the advection of ice, thus indicating the possibility of detecting past drainage events with ice penetrating radar.

Introduction
Nearly 400 subglacial lakes have been identified in Antarctica based on radar data and 25 satellite measurements (Wright and Siegert, 2012), as well as at least two subglacial lakes Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | in Greenland, all of which satisfy established criteria for identification of subglacial lakes (Palmer et al., 2013). Subglacial lakes are typically located either in the ice-sheet interior, close to ice divides where surface slopes are low and ice is thick, or in areas of higher ice velocity where internal ice deformation and sliding add to geothermal energy to produce melting at the base (Siegert et al., 1996;Dowdeswell and Siegert, 1999). 5 Subglacial lakes have been hypothesized to play a role in the initiation of fast flow by (1) providing a steady stream of water downstream of their location, lubricating the base (Gray et al., 2005;Langley et al., 2014), (2) substantially influencing the thermal regime of the ice, gradually warming the basal ice from below, or (3) providing ice with enough thermal momentum to resist freezing on downstream of the lake (Bell et al., 2007). 10 The reason why so few lakes have been discovered in Greenland, despite considerable efforts and relatively dense surveying compared to Antarctica, is believed to be because of the generally warmer ice and higher surface slopes in Greenland, which favour rapid and more efficient drainage as well as increased vulnerability to drainage instabilities (Pattyn, 2008). In addition, the large supply of surface meltwater to the base of the Greenland 15 Ice Sheet through hydrofracturing means that the drainage networks in place are probably already highly efficient, and thus capable of draining subglacial water effectively and preventing or limiting subglacial lake development (Palmer et al., 2013).
Apart from Lake Vostok, the largest subglacial lake in the world, as well as a few other ones, the typical subglacial lake is around 10 km in diameter (Siegert, 2000). Recent studies 20 have presented compelling evidence of rapid transport of water stored in subglacial lakes, indicating that lakes can either drain episodically or transiently on relatively short time scales (Gray et al., 2005;Wingham et al., 2006;Fricker et al., 2007). Lakes appear to form a part of a connected hydrological network, with upstream lakes draining into downstream ones as subglacial water moves down the hydrological potential . Both ice 25 and water transport a substantial amount of sediment over time which is deposited in the lakes, thus reducing their size over long time scales. Sedimentation rates can vary from close to zero to several millimeters per year, with sediment layers estimated to be up to several hundreds of meters thick in some lakes (Christoffersen et al., 2008;Bentley et al., Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | 2013). The dominant mechanism in transporting sediment to subglacial lakes is thought to be influx of sediment-laden water for open-system, or active, lakes, and melt-out from the overlying ice for closed-system lakes, where water exchange happens purely through melting and freezing (Bentley et al., 2013). Subglacial lakes thus change size on a variety of time scales with different mechanisms, from fast drainage to slow sedimentation. 5 The presence of a subglacial lake is often signalled by a flattening of the ice surface, given that the lake is in hydrostatic equilibrium, and as a local speed-up of ice velocities. This draws down cold ice and deflects isochrone layers within the ice (Weertman, 1976).
Numerous studies have investigated the effect of spatially varying basal (or surface) conditions on ice dynamics and the internal layering of ice, such as melting or basal resis-10 tance. Previously, Leysinger Vieli et al. (2007) studied the effect of areas with basal sliding or melting on internal layer architecture, but not explicitly a subglacial lake. The effect of a subglacial lake on ice dynamics was investigated by Sergienko et al. (2007), who used a 2-D vertically-integrated flow equation to study ice-sheet response to transient changes in lake geometry and basal resistance, and by Pattyn (2008) who investigated the stability 15 of a subglacial lake to drainage events. The interaction between lake circulation and ice dynamics and its effect on the basal mass balance has also been studied by Thoma et al. (2010).
The aim of this study is to investigate the influence of a subglacial lake on ice dynamics and the effect it has on the isochrone structure within the ice. As all stress components 20 are important for such an interaction, we use a fully 3-D thermo-mechanically coupled full-Stokes ice-sheet model, implemented in the commercial finite element software COMSOL. We employ an enthalpy-gradient method to account for the softening effect of ice temperature and water content on ice viscosity (Aschwanden et al., 2012) and we show how temporally varying basal conditions can lead to the appearance of flow bands, or arches 25 and troughs, within the internal layering downstream of the original flow disturbance.
Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | 2 Model description

Ice flow
Ice is treated as an incompressible fluid with constant density, obeying conservation laws for mass and momentum: (1) 5 and where u is the velocity vector, g the gravitational acceleration, ρ the ice density and σ the Cauchy stress tensor. The Cauchy stress tensor is given by: 10 where τ and pI are the deviatoric and the isotropic parts, p is pressure, and I is the identity matrix. Inertial forces are assumed negligible and only body forces arising from gravity are taken into account.
Ice is assumed to follow Glen's flow law (Steinemann, 1954;Glen, 1955), in which deviatoric stresses are related to strain rates (ε) by: where η is the effective viscosity,ε e is the effective strain rate, (A t ) is a rate factor that depends on the homologous temperature (T ) and the water content of the ice (W ) and 20 Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | n is the power-law exponent in Glen's flow law. The homologous temperature T = T + γp corrects for the dependence of the pressure-melting point (T m ) on pressure, where γ is the Clausius-Clapeyron constant. The formulation of the rate factor follows (Duval, 1977): and A 0 is a pre-exponential constant, Q is the activation energy and R is the universal gas constant, with numerical values listed in Table 1 (Greve and Blatter, 2009).

Enthalpy balance
An enthalpy-gradient method (Aschwanden et al., 2012) is employed, as opposed to the 10 typically used cold-ice formulation, which is incapable of reproducing correctly the rheology of temperate layers within ice sheets. The enthalpy formulation allows for the possibility of including liquid water content within temperate ice, based on mixture theory, without explicitly tracking the cold/temperate transition surface (Aschwanden et al., 2012;Greve, 1997).
In the enthalpy-gradient method, enthalpy replaces temperature as the thermodynamical 15 state variable, such that: where H = H(T, W, p) is the temperature and water-content-dependent specific enthalpy and Q = 4ηε 2 e is the heat dissipation due to internal deformation. The conduction term in Eq. (8) depends on whether the ice is cold (H < H Tm (p)) or temperate (H ≥ H Tm (p)). The 20 conduction coefficient for cold is ice defined as K i = k/c, where k is the thermal conductivity and c is the heat capacity, both assumed to be constant (Table 1). H Tm (p) is the Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | specific enthalpy of the pressure-dependent melting point of ice. The diffusivity for temperate ice is poorly constrained as little is known about the transport of microscopic water within temperate ice (Aschwanden et al., 2012;Hutter, 1982). In practice, we use the value K 0 = 10 −3 K i , shown by Kleiner et al. (2015) to be sufficiently low to suppress transport of water by diffusion through the ice matrix, while still numerically stable. 5

Age equation
In addition to the balance equations above, we solve a separate equation for the ice age (χ) to determine the influence of the lake on the isochrone structure (Hindmarsh et al., 2009;Parrenin et al., 2006;Parrenin and Hindmarsh, 2007;Leysinger Vieli et al., 2007): 10 where χ is the age of ice and the second term on the right represents a diffusivity term needed for numerical stability, in which d χ is the numerical diffusivity.

Boundary conditions
At the surface, stresses arising from atmospheric pressure and wind can be neglected as they are very small compared to the typical stresses in the ice sheet (Greve and Blatter,15 2009), resulting in a traction-free boundary condition (BC), where n is the normal vector pointing away from the ice. Accumulation and ablation (a s ) at the surface are assumed to be zero, giving the kinematic surface BC as: 20 Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | where z s is the surface elevation. We employ an inverse Weertman-type sliding law (Eq. 12), where the basal drag (τ b ) is expressed as a function of the velocity of the ice (u b ) immediately above the ice/base interface, except over the lake surface where basal traction is set to zero (full slip). With basal sliding exponents (p, q) = (1, 0) appropriate for ice-streaming conditions, the sliding relationship simplifies to a linear relationship between basal sliding 5 and basal traction. The ice is assumed to be in hydrostatic equilibrium everywhere and the basal normal pressure (τ n ) taken as the ice overburden pressure. Ice accretion and melt at the base are assumed to be zero and along with the stress BC at the surface, a nopenetration condition is used to close the system: 10 and u · n = 0, where C is the sliding coefficient (Table 1). Periodic BCs for the inlet/outlet are used, such that the velocity, pressure and specific enthalpy are the same at the upstream and downstream extremes. On the side boundaries of the domain, symmetry for velocity and 15 thermal insulation is imposed: At the surface, a value is set for specific enthalpy (H s ) corresponding to a surface tempera-20 ture of T s = −30 • C (W s = 0, p s = 0) such that: Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | At the base, the geothermal flux (q geo ) is used for cold ice and a zero flux for temperate ice: where H i is the specific enthalpy of pure ice at the melting temperature. To determine correctly the basal boundary condition for the enthalpy field equation, a switching between a Dirichlet and a Neumann condition is necessary (Aschwanden et al., 2012;Kleiner et al., 5 2015), which depends on basal temperature, water availability at the base and whether a temperate layer exists immediately above it. Here we opt for a simpler, more computationally efficient boundary condition where the geothermal flux is gradually decreased in the , with a smoothed Heaviside function with continuous derivatives. 10 For the age-depth relation (Eq. 9), periodic boundary conditions are used for the inlet and outlet. At the surface, the age is set to χ s = 0. Zero normal fluxes are used for the side boundaries and the lower boundary.
− n · (−∇χ) = 0 15 We define the computational domain as a rectangular box 350 km long and 100 km wide with a fixed bed slope (α). Three different lakes sizes are used, based on typical sizes of subglacial lakes found in Antarctica (Siegert et al., 1996;Siegert and Kwok, 2000). The lakes are all elliptical in shape, with major and minor axes defined in Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | layers, which become thinner towards the base, where the thinnest layer has a thickness of ∼ 35 m.

Model experiments
The aim of the study is to show how the presence of a subglacial lake affects ice dynamics and thermal resolution, and to follow its effects on internal layering through simple tempo-5 rally dependent and steady state experiments.
All transient simulations are started from an initially steady-state configuration, where equations for mass, momentum and enthalpy are solved jointly with a direct solver along with equations for surface or grid evolution. All presented results are based on simulations with a fully nonlinear ice viscosity (n = 3) unless otherwise stated explicitly, and all numeri-10 cal values for model constants are defined in Table 1.
The lake itself is modelled as a "slippery spot" (Pattyn et al., 2004;Sergienko et al., 2007). The lake surface is assumed to be fixed to the bed plane. In reality, any changes to lake size or volume would lead to vertical movement of the lake surface, which in itself would induce an expression at the surface of the ice. Sergienko et al. (2007) found that vertical movement 15 of lake surfaces due to drainage resulted only in small changes to ice velocity and that the ice surface reached its initial stage after just a few years, indicating that it would not lead to a significant disturbance in isochrone structure compared to changes in the horizontal lake extent. Here we consider only planar changes in lake geometry, or changes in lake size fixed to the bed plane. 20 Water draining from subglacial lakes is likely to have an impact on ice dynamics and isochrone structures downstream of the lake, either through melting of ice or changes in basal water pressure and sliding, although no representation of subglacial hydrology is included here. The effect the draining water will have on isochrone layers will depend on the state of the drainage system and how the water is transported downstream (Stearns 25 et al., 2008) but will likely be limited by the swiftness of the ice response compared to the time needed to perturb isochrone layers considerably. Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | The lack of basal friction over the lake results in increased velocities, not just over the lake but also a considerable distance both upstream and downstream of the lake that depends on ice thickness, basal traction (Hindmarsh, 2006;Kamb and Echelmeyer, 1986) and the size of the lake. This distance is around 15 km for the L S lake size, 50 km for L M and 100 km for the L L lake size, defined as the distance it takes for the horizontal velocity to 5 drop below 1.05 times the background velocity. Figure 1a shows the surface velocity with black lines indicating streamlines and Fig. 1b shows the relative change in surface elevation caused by the presence of the lake, where the general slope has been subtracted from the surface. The surface responds to the change in basal conditions by becoming flatter and with more than a doubling of horizontal velocities over the lake. Streamlines contract and 10 ice is brought in towards the lake from the sides. Figure 2 shows different model output in vertical cross sections, through the center of the lake, in the direction of flow for the L M lake size. Over the lake, the ice basically moves as an ice shelf with more or less uniform horizontal velocity throughout the ice column (Fig. 2a). At the upstream end of the lake, a strong downward movement of ice causes cold ice from the 15 upper layers to be drawn down towards the bottom, steepening the temperature gradient close to the base. Conversely, at the downstream end a strong upward flow restores internal layers to their prior depths (Fig. 2b). The temperature and microscopic water content within the ice are shown in Fig. 2b. Black lines represent streamlines which coincide, or line up, with isochrone layers for steady-state 20 simulations .
The intense internal deformation close to the borders of the lake, where ice velocities change significantly over short distances, gives rise to a thin temperate layer of ice both at the upstream and downstream ends of the lake (Fig. 2b). The vertical resolution is somewhat limited in the model but simulations with higher resolution do indicate that a temperate 25 layer should form for the given velocity field, although its thickness might be slightly overestimated. This is strongly affected by the chosen sliding law and sliding coefficients and should not be considered representative for subglacial lakes in general. Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Figure 2c shows the internal deformation energy, with high values near the base where velocity gradients are high. Maximum values are reached near the edges of the lake, in the transition zone between low and full sliding. Increases in temperature, pressure and water content, and larger effective stress, all have the effect of decreasing the viscosity, which is why the lowest viscosity values are obtained at the base of the ice sheet, in particular 5 around the edges of the lake (Fig. 2d). Generally low viscosity is furthermore obtained throughout the ice column over the lake boundary, extending all the way up to the surface, effectively creating a vertical zone of softer ice. Figure 3a shows profiles of horizontal surface velocity for the three different lake sizes considered. As expected, velocity peaks over the lake surfaces and two fringe peaks at the 10 lake edges are discernible for the smallest lake size as well.
A characteristic surface dip feature is observed at the upstream end of many lakes in Antarctica, for instance the Recovery Lakes (Bell et al., 2007;Langley et al., 2011) and Lake Vostok (Studinger et al., 2003) as well as a small hump at the downstream end. Figure 3b shows surface profiles for the three different lake sizes, each with the characteristic 15 flattening of the surface as well as the dip and ridge features on each side of the lake. Figure 4a shows a comparison of surface profiles for the L M lake size, made with a nonlinear viscosity (n = 3 in Glen's flow law as used in our other experiments) versus a constant viscosity (η const = 10 14 Pa s) and a pressure-and temperature-dependent viscosity (n = 1) but otherwise using the same setup. The hummocky feature is noticeably absent 20 from simulations with a constant viscosity. Figures 4b and c show profiles of horizontal velocity and viscosity over the centre of the lake for the three different viscosity cases. The horizontal velocity has been scaled with the horizontal velocity at the surface (the uppermost point of the vertical profile). For both the n = 1 and n = 3 cases, the horizontal velocity at depth is larger than at the surface. 25 A weaker contrast in basal traction inside and outside the lake, or a less pronounced switch in flow mode such as one might expect beneath an ice stream (Leysinger Vieli et al., 2007;Parrenin and Hindmarsh, 2007), decreases the amplitude of the dip and ridge features as can be witnessed in Fig. 5a, where surface profiles for simulations with different Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | sliding coefficients and ice thickness are compared (n = 3). Isochrones relative to the ice base are shown in Fig. 5b for two of the simulations presented in Fig. 4, (n = 3 (blue line) and η = 10 14 (red line)), in addition to isochrones for a simulation with a lower contrast in basal traction and a lower bed inclination than the other two (cyan line). The structure of internal layers depends heavily both on the contrast in basal traction and the viscosity for-5 mulation used. Low contrast results in small deflections of internal layers whereas a fixed viscosity (compared to flow with nonlinear viscosity and a similar surface velocity) results in larger deflections of internal layers due to the generally smaller ice flux at depth outside the lake compared to over it.
Lakes drain and fill on different time scales. Several studies have documented relatively 10 rapid drainage events for subglacial lakes in Antarctica (Gray et al., 2005;Wingham et al., 2006;Bell et al., 2007). Typically, drainage occurs over the course of several years but refilling takes much longer. To simulate such an event and what effect it could have on the internal structure of the ice, we set up a model run where the lake diameter shrinks during a 10 year period from the maximum (L L ) to the smallest lake size used in the paper (L S ). 15 Figure 6 shows four different time slices of horizontal velocity, with black lines indicating isochrone layers. As the velocity field adjusts in time to the new boundary conditions, a travelling wave is created at depth within the isochrone structure that transfers downstream with the flow of ice.

20
The frictionless boundary levels the surface over the lake, changing surface gradients and causing the ice to speed up in the vicinity of the lake. The increase in velocity is further amplified by the effect of velocity gradients on ice viscosity. Outside the lake, the velocity field is affected primarily by changes in surface gradients and longitudinal stresses and, to a lesser degree, by changes in ice viscosity and basal velocities. Over the lake, the ice 25 moves more like an ice shelf, with almost uniform horizontal velocity throughout the ice column. The increase in velocity (Fig. 2a) is due to the lack of basal friction over the lake, Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | causing the highest velocities to appear there, but secondary velocity peaks (Fig. 3) can also be discerned at the lake edges, which result from the interaction of surface evolution and ice dynamics. The velocity increase for the two fringe peaks propagates from the surface and downward whereas the velocity peak over the lake is mostly caused by acceleration in the basal layers of the ice sheet.

Thermal regime
For the given model setup, the geothermal flux is sufficient to ensure that the basal temperature reaches the pressure melting point everywhere in the model. Internal deformation then adds to the available thermal energy with the potential to form a temperate layer with non-zero microscopic water content in the vicinity of the lake as seen in Fig. 2b, where a 10 thin temperate layer has formed (∼ 70 m, 2 vertical cells). Lakes in ice-streaming areas, where a considerable portion of the surface velocity arises due to sliding at the base, can be expected to have a much gentler transition in flow mode at the lake edges and much less deformational energy available. No temperate layer forms for simulations where the transition in sliding is considerably less sharp and more appropriate for streaming areas as 15 in the cases with C = 10 10 (Fig. 5).
The steeper temperature gradient over the lake efficiently leads away excess heat created by internal deformation at the upstream end, refreezing whatever microscopic water that was created upstream. Large quantities of microscopic water, in reality, would drain to the base ( 3 %), although drainage is not included here as water content is relatively low. 20 Two factors would contribute to facilitating freeze-on at the interface between ice and lake water, if it were included in the model. Firstly, draw-down of cold ice from above increases the temperature contrast in the lower part of the ice sheet, secondly, heat from internal deformation ceases with the removal of basal traction over the lake surface. Both result in a steepening of the temperature gradient close to the base, which removes heat more 25 efficiently from the base and shifts the balance in favour of freeze-on.
As lake size increases so too do horizontal velocities over the lake (Fig. 3), as well as the effective strain rates and the available deformational energy for internal heating. Accretion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | rates at the ice/lake interface would be limited by the amount of latent thermal energy that the ice above it can lead away efficiently, and any temperate layer formed by internal deformation at the lake boundary would, during its existence, block heat flow and thus accretion, until it completely refreezes. 5 The softening effect of local increases in effective stress at the lake edges effectively creates a vertical layer of soft ice in between higher viscosity ice (Fig. 2c). Strong vertical flow at the edges of the lake results from the localized lack of basal traction. A clear difference can be seen between simulations with constant viscosity and viscosity that depends on pressure and temperature (Fig. 4). The softening effect of increasing temperature and pressure with 10 depth not only causes velocity changes in the vertical to be concentrated in the lower layers of an ice sheet but also means that for areas with varying basal traction, such as subglacial lakes, the ice at depth will support less lateral shear and lower longitudinal stresses compared to the upper layers where viscosity is higher and the ice stiffer. This in return means that as the ice encounters a slippery spot or a spot with a sharp decrease in basal traction, 15 such as a subglacial lake, that the force balance will be different than in the isotropic case, where viscosity is constant. The imbalance in mass flux at depth must be compensated by a more localized increase in vertical flow and a subsequent drop in surface elevation at the upstream side and an upwelling at the downstream side of the lake due to the limited vertical extent of a typical ice sheet. 20 For large subglacial lakes, where the basal traction is zero, the horizontal velocity at depth is predicted to be slightly larger (here about 0.1 % for n = 3, Fig. 4b) than at the surface.

Transition zone
Lakes, such as Subglacial Lake Vostok, should therefore experience extrusion flow at the base, where the basal horizontal velocity exceeds that at the surface. Extrusion flow is not, however, a requirement for the formation of these dip and ridge features. To form them it is 25 sufficient to have a sharp transition in sliding along with a negative downward gradient in ice viscosity. Our model further predicts that dips and ridges can form in situations where there is merely a strong decrease, not necessarily a complete disappearance, in basal traction. Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper |

Drainage experiment
Subglacial-lake drainage cycles (including both draining and filling) can have frequencies on decadal to centennial time scales or potentially even larger (Pattyn, 2008;Wingham et al., 2006). Although the drainage occurs over 10 years here, it can be seen as instantaneous given the generally slow flow of ice. The velocity field adjusts rather rapidly (∼ 100 a) to the 5 new basal boundary conditions relative to the time it takes for the isochrones to respond, as the flow of ice is relatively slow. For episodic drainage cycles, the strength of the response will depend partly on how long it takes for the lake to refill. After roughly 2000 years, the wave has moved far enough downstream to be more or less separated from its initial location (Fig. 6d). As both the upstream and downstream lake boundaries move during the drainage 10 event, a two-troughed wave is created. In draining lakes where one end is much deeper than the other a single wave would be expected, as only the shallower end is likely to move significantly. The velocity of the moving boundary also affects the amplitude of the resulting isochrone disturbance (Wolovick et al., 2014), where a slip boundary, moving with the ice, is capable of distorting isochrone layers to a much greater extent than stationary slip 15 boundaries. For maximum effect, the boundary should be moving at a velocity comparable to the averaged ice column velocity. Only slip boundaries moving with the ice are capable of distorting layers to a greater extent than stationary boundaries. Boundaries moving in opposite directions to ice flow, like our downstream lake boundary, will have a smaller impact on internal layers than a stationary boundary as the relative horizontal velocity increases. 20 Water in very active subglacial systems, such as recently discovered in West Antarctica (Gray et al., 2005;Fricker et al., 2007) has relatively short residence times and fast circulation, contrary to previous beliefs. The impact of such short drainage cycles, where both drainage and refilling happen on decadal time scales, is unlikely to have a strong effect on the internal isochrone layers as it takes a long time for the ice to respond. Regular drainage 25 events in subglacial lakes that have a much shorter cycle than the time it takes for a particle of ice to be brought up by vertical flow at the edge of the lake, will therefore probably not be easily detectable in the isochrone structure. If the frequency of drainage cycles is high, the Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | ice will have little time to respond and the amplitude of the resulting wave will be small. On the other hand, drainage of large lakes in areas with low basal melt rates and consequently long filling times could be expected to generate travelling waves, downstream of the lake, with a sufficiently large amplitude to be detectable within downstream isochrones. The amplitude of the travelling wave would be expected to be similar in magnitude as the steady 5 state isochrone disturbance over the lake itself (Fig. 6). In our simulated scenario of Fig. 6 the travelling-wave amplitudes of ∼ 100 m would be well within the bounds of detectability with modern radar systems. Drainage of subglacial lakes where the transition in flow mode is less abrupt would result in smaller amplitudes (Fig. 5b). As layer stratigraphy is often quite complex, a numerical model of ice age and velocities would be needed to separate the ef-10 fect of temporally changing lake size, or basal conditions, from layer deflections caused by varying basal topography, or rheology.
For our particular setup, the isochrone disturbance, or the travelling wave, should eventually overturn and create a fold as the stress situation downstream of the lake is essentially one of simple shear without any longitudinal extension (Waddington et al., 2001;Jacobson 15 and Waddington, 2005). Resolving this in the model would require an equally fine resolution downstream of the lake, as over the lake itself, which is not done here. For a subglacial lake situated at the onset of streaming flow, a fold might not be expected though, as overturning would be counteracted by longitudinal extension and vertical compression, which would tend to flatten all layer disturbances. In general, both horizontal shear and longitudinal ex-20 tension can be assumed to be present and thus, whether a layer disturbance develops into a fold or flattens out, eventually to disappear, will be decided by the balance between the two (Waddington et al., 2001).

Conclusions
Subglacial lakes represent areas of the ice-sheet base which are incapable of exerting any 25 horizontal stress on the overlying ice. A rapid transition between little to full sliding at the lake edges causes intense deformation and internal heating, leading both to an increase in Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | enthalpy and strain softening of ice. The decrease in viscosity with depth leads to preferential deformation of the lower layers of an ice sheet which, along with a sharp transition in sliding, results in a hummocky surface feature above the upstream and downstream edges of subglacial lakes. These dip and ridge features can be taken as evidence for a rapid transition in basal sliding, with strong vertical movement of ice as a result. Lakes without this 5 feature would be expected to experience a more gradual increase in sliding approaching the lake, a lower contrast in basal traction or continued fast flow downstream, as with subglacial lakes at the onset of or in ice-streaming areas. Over the lake itself, the ice moves more like an ice shelf, with zero basal traction and virtually no internal strain heating. As a result, temperature gradients over the lake steepen, making them more prone to freeze-on 10 at the ice/water interface. A rapid decrease in lake size (or basal friction) causes a travelling wave to be created at depth within isochrone layers, suggesting that certain aspects of lake history are preserved within them and could potentially be deciphered from radio echo sounding data downstream of lake locations when combined with output from numerical models. 15 Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper |