Introduction

Recently discovered active (i.e. filling and draining) sub-glacial lakes under Antarctic ice streams (Reference Fricker, Scambos, Bindschadler and PadmanFricker and others, 2007; Reference Stearns, Smith and HamiltonStearns and others, 2008) and supraglacial lakes on the Greenland surface (e.g. Reference DasDas and others, 2008; Reference Hoffman, Catania, Neumann, Andrews and RumrillHoffman and others, 2011) have sparked interest in the processes controlling lake formation, dynamics and effects on ice flow. Among the additional 124 active subglacial lakes in Antarctica identified in NASA’s Ice, Cloud and land Elevation Satellite (ICESat) laser altimeter data (Reference Smith, Fricker, Joughin and TulaczykSmith and others, 2009), the majority are clustered in regions of fast-flowing ice streams and outlet glaciers. The majority of the ∼2000 supraglacial lakes in Greenland (Reference Selmes, Murray and JamesSelmes and others, 2011) are also located in areas where the ice flow is fast and primarily occurs by basal sliding. The mutual effects of fast ice flow driven by sliding and subglacial and supraglacial lakes is far from being fully explored. The goal of this study is to examine this interaction and to identify factors that control locations of supraglacial and subglacial lakes.

Numerous subglacial lakes detected under Antarctic ice streams are located near areas with enhanced basal traction (e.g. Reference Fricker, Scambos, Bindschadler and PadmanFricker and others, 2007; Reference Fricker and ScambosFricker and Scambos, 2009; Reference Walter, Brodsky, Tulaczyk, Schwartz and PetterssonWalter and others, 2011; Reference Christianson, Jacobel, Horgan, Anandakrishnan and AlleyChristianson and others, 2012). This correspondence is not a coincidence, because in addition to being a potential source of subglacial meltwater, areas of high basal traction (‘sticky spots’) modulate ice flow in a way that favors water ponding (Reference Sergienko and HulbeSergienko and Hulbe, 2011). In Greenland, seasonal supraglacial lakes form at the same location every year without being advected by ice flow (e.g. Reference Lüthje, Pedersen, Reeh and GreuellLüthje and others, 2006), and this suggests that they are tied to undulations of the subglacial bed (Reference Echelmeyer, Clarke and HarrisonEchelmeyer and others, 1991). The recent analysis by Reference Lampkin and VanderBergLampkin and VanderBerg (2011) confirms this suggestion. Using Fourier and wavelet analysis of the basal topography they find a strong correlation between supraglacial lake location and topography with wavelengths of 5.6 and 11 km.

In this study, I use numerical models of ice and water flow to examine how basal conditions (basal topography and basal traction) determine the locations of subglacial and supraglacial lakes. I further explore the hypothesis that subglacial and supraglacial lakes can potentially coexist, due to their common origin by interaction of ice flow with the bed. This hypothesis is plausible because basal conditions (both topography and stickiness) produce both the surface undulations that contain supraglacial lakes and the basal hydrologic potential field minima that contain subglacial lakes. The hypothesis, if true, is important for two reasons. First, it would allow observations of existing surface lake distribution and character to also inform on the conditions of more difficult-to-observe basal lakes. Second, it would imply that water transfer from surface lakes is, by virtue of the presence of coordinated basal lakes, draining preferentially into reservoirs at the bed that are able to collect and retain such water.

Geometry and Models

In order to gain insight into the physical mechanisms controlling lake locations, I consider idealized geometrical settings. For simplicity, I consider only idealized circumstances such as isolated Gaussian-shaped bumps and ‘sticky spots’ of circular footprint (Fig. 1). These two idealized configurations can be modified by reversing the sign of the perturbation to determine comparable effects of basal depressions and ‘slippery spots’ (Reference Sergienko, MacAyeal and BindschadlerSergienko and others, 2007), i.e. the effects caused by these basal perturbations are opposite to the effects of bumps and ‘sticky spots’. The choice of exploring the effects of isolated basal perturbations is motivated by the fact that subglacial and supraglacial lakes themselves appear in isolation, and are thus likely to be the response to a specific isolated feature of the bed. I consider a slab of ice that is 1000 m thick, 200 km long and 60 km wide, flowing on an inclined plane with an inclination angle of 1° (Fig. 1). A Gaussian-shaped bump and a patch with enhanced basal stickiness represent the perturbations at the ice bed. The bump elevation profile is

(1)

where A is the bump amplitude, taken to be 400 m, x _c and y _c are the coordinates of its center and r _c is its horizontal extent, taken to be 6 km. The sticky spot is a circular patch with 5 km radius. The stickiness coefficient (described below) is five times larger within the circular patch than on the rest of the bed. The choice of geometric and basal parameters, although arbitrary, resembles conditions characteristic of both Antarctic ice streams and locations in Greenland where supraglacial lakes are formed. The geometric parameters of basal perturbations (the bump amplitude and the stickiness of the sticky spot) are chosen such that they result in similar ice surface perturbations.

Fig. 1. Model geometry: (a) bump; (b) sticky spot.

This modeling study uses three model components: an ice-flow model and subglacial and supraglacial water-flow models. The basis of the water-flow models is Hele–Shaw flow (e.g. Reference LambLamb, 1932) for subglacial water, and lubrication theory of thin-film flow for the surface flow. The main purpose of these models is to specify how water moves, down the hydraulic potential gradient in the subglacial water model, and down the surface slope in the supraglacial water model. These models are not intended to comprehensively simulate subglacial and supraglacial hydraulic systems. Thus, it is assumed that water flows as a sheet, similar to its treatment in the formulations by Reference NyeNye (1976), Reference ClarkeClarke (2005) and Reference Le Brocq, Payne, Siegert and AlleyLe Brocq and others (2009).

Ice-flow model

I use a two-dimensional, plan-view, vertically integrated ice-dynamics model that accounts for vertical shear stresses in ice. This is a so-called ‘hybrid’, or L1L2, model (Reference HindmarshHindmarsh, 2004). This hybrid approach has significant advantages in the simulation of ice flow over variable basal traction: it is computationally cheap compared to a full-Stokes model, thus allowing fine spatial resolution, yet it retains the capacity of a full-Stokes model to resolve vertical shear stresses (albeit in approximate form) associated with transitions of ice flow over the bed with variable basal traction (Reference Goldberg and SergienkoGoldberg and Sergienko, 2011). I follow the hybrid model formulation of Reference GoldbergGoldberg (2011). The friction law at the ice bed is a linear function of shear stress:

(2)

where is the basal shear, is the horizontal ice velocity at the bed, with components {u, v}, and β is the stickiness coefficient. The value of β is 5 × 10¹⁰ Pa s m⁻¹ in all experiments, except inside the sticky spot, where it is five times larger (2.5 × 10¹¹ Pa s m⁻¹). I use Glen’s flow law, with exponent n = 3, and a vertically averaged ice stiffness parameter equal to 1.68 × 10° Pa s¹/³. The ice dynamics model is complemented by the mass continuity equation,

(3)

where H is the ice thickness, ∂_t is the partial derivative with respect to time, is the two-dimensional divergence operator, is the surface ablation rate (positive for melting, and assumed to be 1 m a⁻¹) and is the basal melting/freezing rate, assumed positive if melting occurs. The basal melting/freezing rate, , is determined from the energy balance at the bed:

(4)

where ρ _i is the ice density, L is the latent heat of fusion for ice, k _i is the heat conductivity for ice, T_z is the vertical gradient of ice temperature (assumed to be constant, T_z = −0.045 K m⁻¹; Reference Joughin, Tulaczyk, MacAyeal and EngelhardtJoughin and others, 2004), G is the geothermal heat flux (G = 42 mW m⁻²) and is the frictional energy dissipation rate. It is assumed that no melting occurs when ice flows over a subglacial lake, due to the absence of frictional energy dissipation and geothermal flux (it is assumed that the water layer in the lake isolates ice from effects of the geothermal flux and that the water residence time in the lake is less than the heat conduction timescale; lake circulation and heat transport are not simulated in this study). Justification of these assumptions and an exploration of the thermodynamic aspects of ice flow over subglacial lakes are given by Reference Sergienko and HulbeSergienko and Hulbe (2011).

Performance of the hybrid ice-flow model used in this study was assessed in a separate study that investigated the effects of basal undulations on ice flow (Reference SergienkoSergienko, 2012), where results of this model were qualitatively compared to the simulations of a full-Stokes ice-flow model. Simulations of ice flow over the bump with the hybrid ice-flow model compare well with the full-Stokes simulations, and the differences between the two models decrease as basal sliding increases.

Subglacial water model

The basis of the subglacial water model is the thin-sheet treatment described by Reference ClarkeClarke (2005). The evolution of the thickness of the water sheet, w _b, is described by the mass conservation equation

(5)

where is the water velocity, which is proportional to the gradient of the hydraulic potential, Φ,

(6)

where p _w is the subglacial water pressure, ρ _w is the water density, g is the acceleration due to gravity and B is the bed elevation. The subglacial water pressure, p _w, is approximated by the ice overburden, i.e. p _w ≈ ρ _i gH. The water velocity formulation follows Reference ClarkeClarke (2005):

(7)

where K is the hydraulic conductivity of the water sheet (K = 10⁻⁶ m s⁻¹). Further details of this model are given by Reference Sergienko and HulbeSergienko and Hulbe (2011).

In contrast to Reference CarterCarter and others (2011), who investigate possible water sources of active Antarctic subglacial lakes (i.e. whether water formed locally in the vicinity of the subglacial lakes or was transported from upstream catchment regions), this study focuses on physical controls that determine specific lake locations. It is a priori assumed that the source of subglacial water is local, and water is formed due to energy dissipation associated with basal friction and due to geothermal flux. This assumption is adopted to eliminate uncertainties associated with poorly known large-scale fluxes of subglacial water, and to focus on regions where subglacial lakes are formed. The question addressed here can be formulated as follows: given the existence of subglacial water (either formed locally or brought from afar) what makes it flow to a particular location to form a subglacial lake? The assumption that subglacial water is sourced locally does not affect the generality of the results and conclusions obtained in this study, which are valid regardless of whether the water source is local or remote.

Surface water model

The model simulating surface water flow is also based on a thin-film treatment. The thickness of the water sheet, w _s, is described by the mass continuity equation that takes the form

(8)

where is the surface water flux determined from the thin-film lubrication theory (e.g. Reference Landau and LifshitzLandau and Lifshitz, 1988, p. 84)

(9)

where η _w is the water viscosity (η _w = 1.8 × 10⁻³ Pa s) and S is the ice surface elevation. This model does not simulate development and flow of surface streams and rivers, nor surface runoff through moulins. The purpose of this model is to simulate surface water distribution, and identify places where it can collect into lakes.

Simulation procedure

A suite of numerical simulations is performed with the finite-element package COMSOL^™ to obtain the following four types of solution:

1. (A) the steady-state configuration of ice flow, surface elevation and thickness in response to basal perturbations of both types in the absence of water at the surface or bed (Fig. 2);

2. (B) the configuration of standing surface water in response to the surface topography obtained in (A) (Fig. 3);

3. (C) the configuration of subglacial water and ice flow, surface elevation and thickness (fully coupled to subglacial water), but with surface water assumed to be absent (Fig. 4);

4. (D) the configuration of the ice and water at the surface and base in a fully coupled simulation (Fig. 5).

Fig. 2. Ice bed and simulated surface topography with a steady-state ice flow: (a) with a Gaussian bump on the bed, and (b) with a circular patch (indicated by black shading) of enhanced basal traction. The mean ice thickness is 1000 m; horizontal dimensions of the bed perturbations are ∼10 km. The vertical dimension is not to scale.

Fig. 3. Ice surface elevation (shape) and surface meltwater depth (m) (colors) after 5 weeks of ablation with a constant melt rate, assumed to be 1 m a⁻¹. (a) Bump, (b) sticky spot.

Fig. 4. Depth of subglacial lakes. Blue arrows are water velocity vectors. (a) Black lines are contour lines of the bump (100 m spacing). (b) Black line outlines the sticky-spot boundary.

Fig. 5. (a, c) Same as Figure 3 when a subglacial lake is present at the ice bed. (b, d) Top view of the surface, depth of the supraglacial lakes (color), location of the subglacial lakes shown by blue contours, locations of the bump and sticky spot shown by black contours.

In simulations of the steady-state ice flow over the bed with irregular basal conditions (simulations (A)) it is assumed that the ablation rate, , is 0. Simulations of the surface and subglacial water flow (simulations (B) and (C)) use the steady-state ice configurations obtained in simulations (A) as initial conditions. Simulations (B) are run as time-dependent for 60 days (the ablation period) with surface ablation rate, . In simulations (C) the ice-flow and subglacialwater-flow models are run together for 10 years to reach steady state (i.e. where typical changes in all parameters are <1% a⁻¹). In this simulation, it was assumed that the ice bed is frictionless, i.e. β = 0, wherever the water depth exceeds 50 cm. In simulation (D), the surface-flow model was run for the same period of time with the same ablation rate as in simulations (B), 60 days and 1 m a⁻¹, using the surface elevation resulting from ice flow in the presence of the subglacial lake, obtained in simulations (C).

Results and Discussion

The numerical model simulations show that the response of the surface of the ice sheet to both forms of basal perturbation

is a dipole-like departure from the large-scale surface slope (Fig. 2a and b). A bulge is formed on the upstream ends of the bump and the sticky spot, and a depression is formed on their downstream ends. Similar surface configurations have been obtained in previous analyses of ice flow over an irregular bed (e.g. Reference JóhannessonJóhannesson,1992; Reference GudmundssonGudmundsson, 2003); however, these studies did not explore the ramifications of these configurations on the dynamics of the lakes.

The hydraulic potential of surface meltwater is controlled by the ice surface slope. The surface depressions caused by ice flow over bed irregularities are thus natural locations for accumulating surface meltwater. Figure 3 shows surface water depth produced under assumptions of a relatively mild (1 m a⁻¹) surface ablation rate. These surface depressions are not advected with the ice flow, due to their bedrock origin. The depth of a supraglacial lake is primarily determined by the depth of the surface depression, and the rate of their filling is determined by the ice ablation rate, both surrounding the lake and within the lake itself, where the ablation rate is amplified by radiative effects of water, and the surface slope (Reference Lüthje, Pedersen, Reeh and GreuellLüthje and others, 2006; Reference TedescoTedesco and others, 2012). In the simulations presented here, the surface depression fills within several weeks of the beginning of the ablation season. The water depth in the supraglacial lakes does not change, however, since, as the ablation season progresses, as soon as the lakes are filled to capacity, they start overflowing (Fig. 3). Simulations with increased surface ablation rate (2 and 3 m a⁻¹) show that the characteristic patterns and general dynamics in all simulations remain the same (i.e. the depth of the lake and its overflow). The increased ablation rate only affects time during which the lake is filled to its full capacity.

A subglacial lake is formed in the location where the hydraulic potential is minimum, and its gradient is zero (e.g. Reference ClarkeClarke, 2005). The hydraulic potential of subglacial water depends on ice overburden, which is determined by the ice thickness and the bed slope. Ice flow over basal irregularities (either the bump or the sticky spot) reduces ice thickness and leads to a surface depression, which then leads to a minimum in the hydrologic potential at the bed. The other component contributing to the hydraulic potential is the basal slope. Antarctic ice streams where active subglacial lakes tend to form, have relatively flat basal topography with small basal slopes. Therefore, the variability of the ice thickness plays a dominant role in the variability of the hydraulic potential. Figure 4 shows water depth in the subglacial lake. The simulated depth values are not representative of all observed lakes, however, due to their strong dependence on the hydraulic conductivity.

There are several positive feedbacks in the configuration, where both the supraglacial and subglacial lakes are present, that suggest that the two types of lakes could exist in pairs. The presence of a subglacial lake modifies the ice surface configuration in a way that stabilizes both the subglacial and supraglacial (if surface melting occurs) lake. This is because the subglacial lake is equivalent to adding additional basal irregularity: a slippery spot. Its effect on ice flow is opposite to the effect of the sticky spot, and results in a bulge formation on the ice surface on the downstream side of the subglacial lake (Fig. 5). On the one hand, this bulge plays the role of a surface levee that inhibits overflow of the supraglacial lake. On the other hand, this bulge is formed by thicker ice, therefore increasing the gradient of the basal hydraulic potential, to make the subglacial lake more stable. This feedback is seen in Figure 5, where the surface depression deepens in response to the presence of the subglacial lake, and collects more surface meltwater than the configuration where the subglacial lake is absent (Fig. 3). A hypothesis of coexistence of supraglacial and subglacial lakes is supported by observations. The presence of subglacial water on the downstream side of a basal topographic feature was revealed by a seismic survey in the vicinity of a supraglacial lake on Russell Glacier, West Greenland (Reference BoothBooth and others, 2011). The data quality was not sufficient to estimate the water depth, however. It should be pointed out that not every supraglacial lake will be accompanied by a subglacial lake. There are other conditions that influence the formation and stability of subglacial lakes (Reference Sergienko and HulbeSergienko and Hulbe, 2011). Additional field observations are needed to establish basal conditions in the vicinity of supraglacial lakes.

Stress distributions around basal irregularities calculated in this study (Fig. 6) indicate that the extensional stresses are sufficiently large to cause development of surface crevasses (Reference Benn, Warren and MottramBenn and others, 2007) that advect with ice flow to locations of the surface depressions, where ice is thinner, and where they may result in moulin formation. Field observations of persistent, multi-year moulins confirm these suggestions, and indicate that moulins tend to form in thinner ice, and are associated with locations of strong basal melting (Reference Catania, Neumann and PriceCatania and others, 2008; Reference Catania and NeumannCatania and Neumann, 2010).

Fig. 6. Surface ice velocity (m a⁻¹) (color) and the principal components of the horizontal deviatoric stresses (kPa) for (a) bump and (b) sticky spot. White color corresponds to extensional stresses; black color corresponds to compressional stresses. (a) Gray lines are contour lines of the bump, spaced at 100 m. (b) Gray line outlines the sticky-spot boundary.

The size and shape of the subglacial and supraglacial lakes are determined by the size and shape of the basal irregularities and the ice surface response to the presence of these irregularities. As analysis by Reference GudmundssonGudmundsson (2003) shows, basal irregularities with spatial extents of several ice thicknesses affect ice flow the most. The ice surface topography created by ice flow over such irregularities will have dominant length scales that reflect characteristic length scales of the basal irregularities. Thus, subglacial and supraglacial lakes can be viewed as a surface proxy that reflects conditions at the ice bed. Analysis of the spatial and size distribution of the active subglacial and supraglacial lakes, similar to one performed by Reference Lampkin and VanderBergLampkin and VanderBerg (2011), can provide valuable information about characteristic length scales of ice basal conditions that cannot be observed directly. It should be emphasized, however, that the surface undulations caused by variable basal conditions do not mirror the magnitudes and spatial distributions of the basal conditions.

Conclusions

The analysis of ice flow over an idealized irregular bed topography and stickiness reveals unifying principles that determine the location of subglacial and supraglacial lakes. These lakes cannot occur randomly across an ice sheet, but must be coordinated with the features that control surface topography and basal hydrologic potential. In other words, subglacial and supraglacial lakes are manifestations of ice-flow response to heterogeneous basal conditions. The results of this analysis also suggest that supraglacial and subglacial lakes could appear in pairs, originating from the same basal feature, and that they are mutually reinforcing. Of course, the presence of supraglacial lakes is determined by the availability of surface meltwater, and this may explain why twinned lakes are not commonly seen in Antarctica, where surface warming is not yet able to support extensive melt seasons. However, they may become common features in areas experiencing strong warming, such as the Antarctic Peninsula.