1. Introduction

Glacial floods, often known by their Icelandic name of jökulhlaups, occur when large volumes of water drain catastrophically beneath an ice mass from either ice-marginal or subglacial ice-dammed lakes (Reference Tweed and RussellTweed and Russell, 1999). The best-studied jökulhlaups are those from the subglacial caldera lake Grímsvötn, beneath Vatnajökull, Iceland (Reference BjörnssonBjörnsson, 2003). The high geothermal heat flux causes enhanced subglacial melting, and subglacial hydraulic gradients direct water to Grímsvötn where it collects. Jökulhlaups are triggered when the lake level rises to within 60–70 m of ice overburden near the lake outlet (Reference BjörnssonBjörnsson, 1988, Reference Björnsson1992). Water is released from the lake and flows subglacially to the ice-cap margin (Reference Björnsson and GuðmundssonBjörnsson and Guðmundsson, 1993). A positive feedback results during the rising limb of the flood hydrograph, as the frictional heat of flowing water enlarges the ice-walled conduits by melting faster than they can close down by ice deformation. This is primarily controlled by the effective pressure, the difference between the ice overburden pressure and the conduit water pressure. After peak discharges have been reached, water pressure drops rapidly in the conduit which closes down and eventually terminates the flood before all the water in the lake has emptied. Consequently, the flood hydrographs are typically asymmetrical, with a gradual but accelerating rise to peak followed by a rapid fall to zero discharge. Jökulhlaups from Grímsvötn have occurred with a periodicity of about 5 years.

Many of the theories used to describe jökulhlaups have been developed and tested using data gathered at Vatnajökull. The starting point for most theories is that developed by Reference RöthlisbergerRöthlisberger (1972). Based on the Gauckler–Manning–Strickler formula of turbulent flow, it models steady-state water flow in a single conduit. The theory accounts for the equilibrium between frictional melting of the conduit walls and the closing of the conduit due to the effective pressure. Since jökulhlaups are only in steady state momentarily (at their peak), the theory has limited use.

Reference NyeNye’s (1976) theory introduces time dependence. The result is a series of five coupled partial differential equations which can be solved numerically for the key parameters of a jökulhlaup: conduit area, discharge, water pressure, melt rate and water temperature. Nye’s theory assumes the water in the conduit is at a uniform temperature, except for a thin boundary layer where heat is transferred to the surrounding ice. Heat transfer is assumed to be instantaneous; advection of energy along the conduit is not accounted for. Despite this simplification, Nye’s theory is extremely successful in many cases. In his own words: ‘the excellence of the fit between observation and theory up to the flood peak is rather astonishing’ (Reference NyeNye, 1976, p.194), although the theory does tend to overestimate peak and total discharge.

Reference Spring and HutterSpring and Hutter’s (1982) theory is more complex than Nye’s. Derived from first principles of continuum physics it is a general theory, which can be reduced to Nye’s with certain assumptions. The key advance over Nye’s theory is that it accounts for the advective heat transfer along the conduit; released frictional energy leads to an increase in water temperature, which leads to melting. In practice, it is found that many of the improvements over Nye’s theory are relatively small. The lack of experimental data associated with jökulhlaups means that the full generality of the Spring and Hutter theory can rarely be applied. However, it does provide a slightly better fit with experimental data from Grímsvötn compared to Nye’s theory.

Reference ClarkeClarke (1982) added to Nye’s theory by including the effects of reservoir geometry and lake temperature. A major advance was made by Reference ClarkeClarke (2003) who modified the advective treatment of heat transfer in the Spring and Hutter theory, allowing frictional energy release to be directly converted to conduit wall melting. The theory predicts smaller values for conduit roughness compared to earlier theories, values that are comparable to those for natural streams and rivers. It also predicts that the locations of flow constrictions, which control flood magnitude, can migrate along the conduit through time and are not restricted to the upper end where the conduit leaves the lake as in earlier theories.

Other theoretical developments have included work by Reference FowlerFowler (1999) who added to Nye’s theory, considered the flood release mechanism more explicitly and was able to reproduce the periodic behavior of the Grímsvötn jökulhlaups. Reference Fowler and NgFowler and Ng (1996) substantially modified Nye’s theory to deal with broad low-sediment floored canals rather than the circular or semicircular ice-walled conduit of earlier theories. Further complexity has been added by considering flow in a basal sheet coupled to a system of conduits (Reference Flowers, Björnsson, Pélsson and ClarkeFlowers and others, 2004). Flow begins as a wide sheet across the base of the glacier, which then rapidly feeds the growing conduit system. Compared to earlier theories, where conduits take days to weeks to enlarge, conduits in the coupled model enlarge over a few hours as the source water is spread along the length of the conduits. The model is able to reproduce the main features associated with the atypical 1996 Grímsvötn jökulhlaup.

A substantial amount of work has been carried out in applying the above theories to jökulhlaups from subglacial lakes such as Grímsvötn, where lake areas are 5–40 km², drainage distances are about 50 km, ice thicknesses are typically 300–600 m (e.g. see cross-section in Reference BjörnssonBjörnsson, 1998) and water exits the glacier margin where pressures are atmospheric. However, recent findings suggest it is possible for jökulhlaups to occur beneath Antarctica from lakes that are up to 600 km² in area, where water flow extends for hundreds of kilometers, ice thicknesses are up to 4 km and the system is completely closed with water being transferred between pressurized lakes (Reference Wingham, Siegert, Shepherd and MuirWingham and others, 2006; Reference Fricker, Scambos, Bindschadler and PadmanFricker and others, 2007).

Reference Wingham, Siegert, Shepherd and MuirWingham and others (2006) analyzed European Remote-sensing Satellite-2 ERS-2) satellite altimeter data from the Adventure subglacial trench region in East Antarctica between 1995 and 2003 (Fig. 1). In one region (inferred to be a ∼600 km² lake and known as lake L), surface elevation dropped by ∼3 m over a 16 month period between late 1996 and early 1998. Approximately 290 km away, three areas (known to overlie subglacial lakes from an independent ice radar survey and referred to as lakes U1, U2 and U3) rose by a comparable amount over the same time period. All lakes lie on a subglacial drainage pathway as determined from subglacial hydraulic potentials, with lake L lying upstream (i.e. up-potential, but at a lower elevation) from lakes U1, U2 and U3 (Fig. 1).

Fig. 1. Surface and bed topography of Adventure subglacial trench region of East Antarctica. Subglacial lakes (from Reference Siegert, Carter, Tabacco, Popov and BlankenshipSiegert and others, 2005) are denoted by triangles. Approximate locations of the source lake (L) and sink lakes (U1, U2 and U3) from Reference Wingham, Siegert, Shepherd and MuirWingham and others (2006) are shown as circles. The black lines are the predicted locations of the main drainage pathways across the subglacial hydraulic potential surface. The position of U3 is different to that shown in Reference Wingham, Siegert, Shepherd and MuirWingham and others (2006) since the location of our predicted subglacial pathway is different in this region. It is the same distance down-hydraulic potential as Reference Wingham, Siegert, Shepherd and MuirWingham and others’ (2006).

2. Methodology

In this paper, we use a numerical modeling strategy to estimate steady-state and time-dependent water flow between the subglacial lakes in the Adventure subglacial trench, as identified by Reference Wingham, Siegert, Shepherd and MuirWingham and others (2006). We begin with a simple energy-conservation model and a model based on Reference RöthlisbergerRöthlisberger’s (1972) theory, which allow us to estimate the steady-state behavior of the system. We then introduce time dependency with a model based on an implementation of Reference NyeNye’s (1976) theory for non-steady flow in a subglacial conduit.

One of the key assumptions throughout this paper is that the ice at the base of the ice sheet in the Adventure subglacial trench is at the pressure-melting point. Output from a three-dimensional (3-D) balanced ice-flux flow model (personal communication from R. Hindmarsh, 2007) suggests this is a valid assumption. For realistic boundary conditions of bed and surface geometry and geothermal heat flux, the model calculates spatial patterns in the basal shear strain rate by multiplying the balance velocity by the driving stress. It also accounts for horizontal advection of heat. It predicts temperate basal ice throughout the Adventure subglacial trench region and a basal temperature gradient of the order of 0.01 Km ^{− 1}. Application of our simple energy-conservation model (section 3) and Röthlisberger’s theory (section 4) shows that the flood between the lakes is accommodated in a semicircular conduit with a steady-state cross-sectional area of 36 m², equivalent to a radius of <5 m. Application of Nye’s theory (section 5) shows that the flood is likely to be accommodated in a semicircular conduit with a maximum cross-sectional area (at maximum discharge) of 250 m², equivalent to a radius of 13 m. Even at maximum discharge, the conduit is therefore likely to extend no more than 13 m above the bed. With a temperature gradient of 0.01 Km ^{− 1}, the ice temperature at the top of the conduit will be just −0.14°C. This will increase the viscosity at the top of the conduit by less than 3% compared to its value in our implementation of the Röthlisberger and Nye theories, where we assume ice is temperate throughout. Given other uncertainties in the boundary conditions for the Adventure subglacial trench, we feel we are justified in assuming ice is temperate throughout the depth of the conduit.

3. Calculations Using a Simple Energy-Conservation Model

Before launching into the full time-dependent theory of Nye, it is instructive to analyze water transfer using a more simplistic approach. This provides an estimate of the results we might expect to obtain from Nye’s theory, and also places some constraints on the input parameters. Figure 2 depicts a system schematic. To simplify the problem, we assume water transfer occurs between two lakes. Furthermore, we assume the source lake L and the sink lake U are joined by a single conduit, an assumption often made when dealing with jökulhlaups. It is justified by considering the pressures in a multi-conduit system. If one conduit in the system is enlarged, then the water pressure in that conduit will decrease. The resulting pressure gradient will drive water from the small conduits into the large conduit. Increased water flux will result in further enlargement of the conduit through frictional melting. The pressure gradient is therefore maintained and the large conduit grows at the expense of the smaller ones. As a consequence, drainage is likely to be dominated by one conduit.

Fig. 2. Simplified schematic diagram of the two subglacial lakes, L and U, joined by a single conduit. Water flow is from the lower lake to the upper lake as the ice surface slope and subglacial hydraulic potential gradient both trend in that direction.

A simple energy-conservation model can be used to estimate the size of the conduit produced during the flood. Treating the ice-bounded lakes as pistons, the net energy associated with the flood (E _T) is the change in gravitational potential energy of the system minus the energy required to keep the water at the pressure-melting point as it flows from an area of thicker ice to one of thinner ice. This can be related to the energy required to melt a conduit through the ice (E _M) to obtain an expression for conduit size S. The detailed derivation of the simple energy-conservation model can be found in Appendix A.

To solve the energy-conservation equations, values are required for the total volume of water discharged V, the length of the conduit l and the change in surface elevation of the source lake Δh _L. Reference Wingham, Siegert, Shepherd and MuirWingham and others (2006) were able to constrain values for these parameters based on observational evidence to V = 1.8 km³, l = 290 km and Δh _L = 3 m. The values of ΔH (difference in bed height) and Δz (z _L − z _U, i.e. change in ice thickness) are less well constrained by available data. Errors in excess of ±30 m are possible (Reference Lythe and VaughanLythe and others, 2001).

Figure 3 shows how the predicted conduit cross-sectional area S varies with ΔH and Δz, using values for the other parameters as specified by Reference Wingham, Siegert, Shepherd and MuirWingham and others (2006). For ΔH = 262 m and Δz = 449 m, a cross-sectional area of 36 m² is predicted. This is in fairly good agreement with the 28 m² calculated by Reference Wingham, Siegert, Shepherd and MuirWingham and others (2006) using a different energy-conservation model with the same values. Figure 3 shows that there is a cut-off beyond which a conduit cannot exist. Within the energy-conservation model, this occurs when the total energy released ET becomes negative. We can therefore constrain conduit creation to a region in the ΔH−Δz plane (Fig. 4). Given the simplicity of the energy-conservation model, particularly its exclusion of dissipative effects such as friction, it is expected to overestimate the cross-sectional area of the conduit.

Fig. 3. Predicted values of conduit cross-sectional area in the energy-conservation model for different values of ΔH and Δz.

Fig. 4. Areas of ΔH and Δz space, in which a semicircular conduit can exist for the energy-conservation model, and the Rthlisberger model where k = 12.5 m^1/3 s⁻¹ and k = 30 m^1/3 s⁻¹. A conduit can only exist for values above the plotted lines.

4. Calculations Using Röthlisberger’s Theory

Based on the Gauckler–Manning–Strickler formula for the mean velocity of turbulent flow, Röthlisberger’s theory describes steady-state flow through an ice-walled conduit. Steady-state flow occurs when the rate of conduit enlargement due to flowing water equals the rate of closure due to ice deformation and discharge is constant. During a jökulhlaup, this condition is met only at peak discharge. Röthlisberger’s theory therefore provides an alternative to the energy model for predicting the maximum cross-sectional area of the conduit. For a conduit of circular cross-section, we have

(1)

and for a conduit of semicircular cross-section, we have

(2)

where Q is the discharge through the conduit, k is the inverse of the Manning roughness parameter M, Φ is the hydraulic potential in the conduit, x is the distance along the conduit (x = 0 at lake L) and ρ _w is the density of water. The detailed derivation of Equations (1) and (2) is provided in Appendix B. It should be noted that Equation (2) differs from the equations used by Reference Wingham, Siegert, Shepherd and MuirWingham and others (2006; supplementary methods 2, equation (4)) due to what we believe are two mistakes in their equation. First, their relationship seems to be based on the hydraulic radius of a circular (rather than semicircular) conduit and, second, they have made a numerical error in the factoring out of the hydraulic radius in the Gauckler–Manning–Strickler equation. These errors result in estimates of conduit cross-sectional area that are a factor of approximately 4 too small.

Reference Wingham, Siegert, Shepherd and MuirWingham and others (2006) predicted a cross-sectional area of 26 m² for a semicircular conduit. This was based on ΔH = 262 m, Δz = 449 m, Q = 50 m³ s ^{− 1} and k = 12.5 m^1/3 s ^{− 1} (i.e. M = 0.08 m^−1/3 s, a value appropriate for a relatively rough-walled conduit). Using the same values in the corrected equation gives a cross-sectional area of 98 m² for a semicircular conduit and 91 m² for a circular conduit. Figure 5 shows how the predicted value of S changes with ΔH and Δz. The values are much larger than those predicted using the energy-conservation model for the same values of ΔH and Δz (cf. Fig. 3). Since it was expected that the energy model would overestimate the size of the conduit, it seems unlikely that these results are correct. However, by choosing a higher value for k (corresponding to a smoother conduit), the Röthlisberger model can be made to produce lower estimates of conduit area than the energy model. A value of k = 47 m^1/3 s ^{− 1} (M = 0.02 m^−1/3 s) produces a conduit area of 36 m², the value calculated by the energy-conservation model.

Fig. 5. Predicted values of conduit cross-sectional area in the Röthlisberger model for different values of ΔH and Δz.

The values of ΔH and Δz can be better constrained using the energy-conservation model in conjunction with the Röthlisberger model. Since the energy model should overestimate the cross-section of the conduit, we can stipulate that allowed values of ΔH and Δz must give smaller areas for the Röthlisberger model than for the energy model. Figure 4 depicts the constraints on ΔH and Δz for k = 12.5 m^1/3 s ^{− 1} and the much larger region of allowed values for k = 30 m^1/3 s ^{− 1}.

5. Calculations Using Nye’s Theory

While Röthlisberger’s theory is useful for estimating some of the flow parameters for a steady-state discharge, it cannot be used to study the time evolution of the jökulhlaup. Furthermore, it can only be used to predict the cross-sectional area of the conduit, which cannot be measured directly. To address these limitations, we now turn to Nye’s theory. We begin by defining some notation. A[x, t] refers to the value of a parameter (which we arbitrarily refer to as A) at position x and time t, where x is the distance from lake L; t = 0 corresponds to the initialization of the parameters. x = L and x = U are the positions of lakes L and U, respectively. A particular value of a parameter which has only spatial dependence (which we arbitrarily call B) is referred to as B[x]. For the purposes of numerical analysis, the conduit is divided into a finite number of segments of length δx. The parameters are calculated at each gridpoint (integer numbers of δx). Each gridpoint is referred to by an index n, where x = L + nδx. A similar scaling is used for the time coordinate.

5.1. Solving Nye’s equations

Nye derives five coupled partial differential equations to describe non-steady flow in a conduit, which account for the geometry and flow of ice, continuity of water, flow of water and energy and heat transfer (Reference NyeNye, 1976, equations (16–20)).

For reasons outlined in section 2 above, the temperatures of the lake water and the ice are likely to be similar (at or close to the pressure-melting point) and so we neglect the heat transfer equation (Reference NyeNye, 1976, equation (20)). We are also able to simplify the energy equation (Reference NyeNye, 1976, equation (19)) by neglecting the final term. As the specific heat capacity of water is ∼4000 J K ^{− 1} kg ^{− 1} and the latent heat of melting for ice is 3.3 × 10⁵ J kg ^{− 1}, this is a reasonable assumption. In order to remove all temperature dependencies, we also assume that the internal energy is constant. Under these conditions, Nye’s equation (19) reduces to:

(3)

where m is the melting rate (kg m ^{− 1} s ^{− 1}). Reference Spring and HutterSpring and Hutter (1982) found that neglecting the changes in internal energy leads to the overestimation of peak discharge and total discharge volume. We must keep this in mind when analyzing and interpreting our results.

Solving the four remaining differential equations (Nye’s equations (16–18) and Equation (3)) allows the time evolution of S, Q and m to be modeled given initial values and the pressure. All that remains is to determine how the pressure behaves. Ordinarily, jökulhlaups drain to a sink of constant atmospheric pressure. In this case, however, the pressure of the sink lake changes with time. Nye’s equations solve the discharge parameters iteratively from source to sink. We overcome the problem of variable sink pressure by re-scaling the pressure of the source lake at each time-step, based on the pressure in the sink lake. The same is done with the ice overburden pressure, i.e.

(4)

(5)

where P _{i s}[x, t] is the ice overburden pressure at distance x from lake L and at time t, scaled to equal zero at lake U. P _s[L, t] is the scaled water pressure at lake L at time t scaled to equal zero at lake U.

The pressures at the lakes P can be calculated by considering the change in hydraulic head as the lake fills or empties:

(6)

(7)

where A _L and A _U are the areas of lake L and lake U, respectively. Using Equations (5–7), it is possible to calculate the scaled pressure at lake L at any time. The scaled pressure along the conduit can then be found by rearranging Reference NyeNye’s (1976) equation (18):

(8)

where g_s is the component of gravitational acceleration parallel to the conduit, and N is a constant for any conduit shape; N = (iR^{− 2})^2/3 ρ _wgk^{− 2}, where i is the conduit hydraulic gradient and R is the conduit hydraulic radius, defined as the cross-sectional area divided by the wetted perimeter.

Evaluating Equation (8) midway between gridpoints allows the scaled pressure at each gridpoint to be calculated by linear interpolation. The remaining difficulty is that Nye’s model cannot be used to calculate water flow along sections of conduit where the potential hydraulic gradient is positive, since this would require taking the square root of a negative number.

Analysis of the surface, bed and derived hydraulic potential gradients along the Adventure subglacial trench shows that sections of positive potential gradient exist as water flows out of bedrock overdeepenings (Fig. 6). However, the overall hydraulic potential gradient between the two lakes is negative and we therefore model flow along a conduit of constant bed elevation gradient (ensuring potential hydraulic gradient is always negative) between the lakes (Fig. 6). The time evolution of S, Q, m and P along the conduit can then be found using a finite-difference time domain (FDTD) method.

Fig. 6. Surface and bed profiles and subglacial hydraulic potential for the proposed path of the jökulhlaup along the Adventure subglacial trench shown in Figure 1. Also shown is the smoothed bed used in the Nye model.

5.2. Sensitivity analysis using Nye’s theory

Specifying values for the input parameters is difficult since they are poorly constrained by observational data. It is crucial to understand how this uncertainty may affect the results. Here we examine the sensitivity of the model to four key parameters: 1. k, the inverse of the Manning roughness parameter; 2. S ₀, the initial conduit cross-sectional area; 3. the conduit geometry; and 4. the ratio of the source- to sink-lake radii.

For the purposes of this sensitivity analysis, we assumed that the discharge occur between lakes L1 and U3, a distance of 267 km. (Note that this is slightly shorter than the distance of 290 km calculated by Reference Wingham, Siegert, Shepherd and MuirWingham and others (2006) and used in our analysis so far. This figure is based on our own analysis of the ice-sheet surface, bed and hydraulic potential; see Fig. 6.) A single conduit of constant elevation gradient is assumed to link the two lakes. A time-step size of 48 hours and a segment size of 500 m are used. Unless otherwise stated, k is assumed to be 30 m^1/3 s⁻¹, the initial cross-sectional area is set to 1 m², the default conduit geometry is semicircular and both lakes are assumed to have a radius of 14 km and a cylindrical geometry. A radius of 14 km is consistent with an area of 600 km², the area of lake L estimated by Reference Wingham, Siegert, Shepherd and MuirWingham and others (2006) from their observational data.

Reference Wingham, Siegert, Shepherd and MuirWingham and others (2006) estimated that k could lie anywhere in the range 10–100 m^1/3 s⁻¹. Figure 7 shows how the predicted discharge hydrograph and change in lake depth vary with k. As k increases from 20 to 80 m^1/3 s⁻¹, peak discharge increases from ∼100 to 1000 m³ s ^{− 1} and the time to peak drops from >20 to ∼5 years. Similarly, the water-level drop (rise) in the source (sink) lake alters from ∼25 to >40 m. The increase in peak discharges is associated with a doubling of conduit cross-sectional area from ∼150 to ∼300 m² (not shown in Fig. 7). Thus, k affects not only the magnitude of the peak discharge, but also the time taken for the flood to initiate and the lifetime of the flood. This is to be expected since k represents the smoothness of the conduit (a higher value of k is indicative of a smoother conduit). A smoother conduit will provide less resistance to the flow, allowing the discharge to increase more rapidly and allowing a conduit of smaller cross-section to support a larger discharge.

Fig. 7. (a) Discharge and (b) change in lake depth for different values of k using Nye’s model.

Since Nye’s theory does not model the triggering of the jökulhlaup, it is necessary to specify an initial conduit area, S ₀. Figure 8a shows how the choice of this value affects the predicted discharge. Peak discharges are not significantly affected, but the peak discharge timing falls from ∼15 to ∼7 years as S ₀ increases from 1 to 7 m². The initial cross-sectional area of the conduit affects the time taken for the flood to initiate, as would be expected. It affects the magnitude of the peak discharge only slightly; a longer lag time before the peak discharge allows a greater amount of water to have transferred between the lakes, lowering the peak discharge.

Fig. 8. Discharge in Nye’s model: effect of (a) initial conduit area S ₀;(b) conduit geometry, and (c) ratio of lake L radius to lake U radius.

The geometry of the conduit affects the value of the hydraulic radius R. Figure 8b shows the difference in flood hydrographs for a semicircular conduit (which we assume is most likely at the bed) and a circular conduit. The latter would be relevant if the conduit were surrounded by ice, which may be applicable if it did not everywhere follow the bed. A circular conduit has a smaller value of N, which results in a faster onset of the jökulhlaup and a higher peak discharge.

Our model assumes that the source lake has an infinite capacity. The jökulhlaup therefore terminates when the hydraulic potentials of the lakes are equal, rather than when the source lake empties. The flood termination requires the cumulative change in depth of the two lakes to reach a certain value. For a given discharge, a lake with a smaller radius will change depth more rapidly, requiring a smaller total discharge to reach the critical value for termination. The ratio of the lake radii therefore has an important effect on the flood hydrograph (Fig. 8c). The timing of peak discharge drops by ∼1 year and its magnitude declines eight-fold as the ratio of source-lake to sink-lake radii changes from 1 : 1.5 to 1 : 0.5. A larger sink lake will change depth less for a given discharge, therefore requiring a greater depth change of the source lake to reach termination. This results in a larger total flood discharge.

5.3. Optimizing the parameter set

It is also possible to begin to optimize values of the main parameters discussed above by comparing model calculations with observations. The key (indeed, only) observational data are the ice-sheet elevation changes observed using ERS-2 altimetry data by Reference Wingham, Siegert, Shepherd and MuirWingham and others (2006); these data can be compared with modeled change in lake depth through time.

In the modeling work discussed so far, we found the change in lake depth with time by integrating the discharge and dividing by the lake area. An assumption made by Reference Wingham, Siegert, Shepherd and MuirWingham and others (2006) (and also in our simple energy-conservation model and our steady-state model discussed above) is that the observed ice-surface deformation corresponds directly to the change in lake depth. However, our time-dependent modeling work shows that the changes in lake depths predicted by Nye’s theory (Fig. 7b) are approximately ten times larger than the observed surface deformations (Fig. 9).

Fig. 9. (a,b) Optimized modeled and measured ice-surface deformation at (a) lake L and (b) lake U3. Error bars on the observed data were determined empirically by analyzing errors in nearby stationary ice (Reference Wingham, Siegert, Shepherd and MuirWingham and others, 2006). (c) Discharge for the optimized model.

Little is known about how vertical deformations at the base of an ice sheet manifest themselves at the surface. Three factors that must influence this are:

1. 1. the lateral extent of the basal deformation anomaly compared to the ice thickness;

2. 2. the time over which the vertical deformation occurs (the deformation rate); and

3. 3. the ice viscosity, represented by the rate factor in the flow law.

These will control the importance of horizontal stress gradients and the ability of lateral ice flow to offset the effect of the initial vertical basal deformation anomaly (Reference Jarosch and GuðmundssonJarosch and Guðmundsson, 2007). To allow for these effects, the surface deformation is assumed to scale linearly with the basal deformation. An additional parameter (the ice-thickness scaling factor A) is introduced to model this relationship, i.e.

(9)

where z is the ice thickness.

It should be noted that the value of A will be relevant to our particular case study, but will not necessarily be universally applicable.

Linear scaling is supported on the basis of empirical observations of water level and ice-surface elevation changes around Grímsvötn (Reference BjörnssonBjörnsson, 1988, Reference Björnsson1992), although we recognize that the Adventure subglacial trench system is an order of magnitude greater than the Grímsvötn system. Theoretical considerations also suggest that linear scaling is appropriate; Reference GuðmundssonGuðmundsson (2008) has derived analytical solutions for the ice-surface response to perturbations in basal conditions using a shallow ice-stream approximation. The magnitude of surface vertical motion is dependent on the size and shape of the basal perturbation, the surface slope and the slip ratio, but the important point is that surface vertical movement scales linearly with the basal perturbation. This theory is applicable for basal perturbations of more than four ice thicknesses, which is the scale of the lakes beneath the Adventure subglacial trench. However, we recognize the theory might be more relevant to ice streams (with their high slip ratios) than the inner part of the East Antarctic ice sheet.

We can now optimize the main parameters in Nye’s model together with the ice-thickness scaling factor by matching modeled to observed ice surface deformations. A rigorous approach to determining the parameters would use Bayes’ theorem to find the posterior probability function for each parameter set (e.g. Reference MacKayMacKay, 2003). However, the large number of parameters makes calculation of the likelihood function problematic. A Monte Carlo technique is required, and this demands large computational resources. In this investigation we use a simpler, less computer-intensive approach. The mean square difference between the observed altimetric data and the surface deformations predicted by Nye’s theory (with the ice-thickness scaling factor) is used as a measure of model success. It is assumed that the parameter set that minimizes the difference between observations and predictions represents the most likely values of the parameters. Hence the inference problem is reduced to one of minimization.

The minimization is achieved using a Simplex-type approach (Reference VanderbeiVanderbei, 2001). Starting from an initial set of model parameters, each parameter is varied in turn by a fixed amount and the fit is calculated. If varying the parameter improves the fit, then the parameter’s value is updated and the iteration continues. If varying any of the parameters does not improve the fit, then a minimum difference has been reached. The resolution of the method is determined by the step size for each parameter. While this method requires significantly less computation than the Bayesian method, it cannot distinguish between local and global minima. Realistic starting values for the parameters are therefore critical to the success of the method. The ice-thickness scaling factor was not included as a parameter in the Simplex optimization. Instead, for each Simplex iteration, its optimal value is found by considering all of its possible values and finding the best fit. This decreases the number of Simplex iterations required, reducing run time.

Using the same system of lakes and step sizes as in section 5.2 and assuming a semicircular tunnel with an initial cross-sectional area of 1 m², the minimization was carried out over five different parameters: k; the ratio of the source- to sink-lake radii r _L /r _U; error in the bed height for both the source and sink lakes; and the ice-thickness scaling factor A. Each run was carried out for a total time of 30 years. The parameters were able to vary by the following resolutions: ±1 m^1/3 s ^{− 1} for k; ±0.01 for r _L /r _U; ±1 m for each of the lake bed height errors; and 2 x 10 ^{− 5} for A. The first four parameters were allowed to vary over any range necessary to achieve a good fit; A was allowed to vary between 1 × 10 ^{− 4} and 1 × 10 ^{− 2}. For the lake bed height errors, each of the lake heights was allowed to vary separately above and below the heights taken from the BEDMAP data (Reference Lythe and VaughanLythe and others, 2001).

Eight different sets of initial values were tested, leading to eight different minima being located with a range of 0.161–0.912 (Table 1). The smallest of these minima (a difference of 0.161) corresponds to a value of k of 55 m^1/3 s⁻¹, a r _L /r _U of 1 : 1.44, an error of +4 m in the height of lake U and 0 m for lake L and an ice-thickness scaling factor of 4.16 × 10 ^{− 3}. Figure 9 depicts the associated fit with the altimetric data for the two lakes and the associated flood hydrograph.

Table 1. Initial and final values of the parameters used in the optimization scheme, together with the associated degree of fit (minima between observed and predicted ice surface deformation). Lake offset is the difference between the bed height of the lake used in the model and that predicted by the BEDMAP data

Nye’s theory combined with the bed-to-surface scaling function tends to underestimate the maximum surface deformation (Fig. 9). This is because the lake U data show a decrease in ground deformation after the peak at 3 years (Fig. 9b). Nye’s theory provides no mechanism for this to occur, and the optimization instead fits a line through the middle of the points. Since it is the cumulative fit of both lakes that is used to optimize the parameters, this error has carried through into the lake L data (Fig. 9a). While it is not absolutely certain that these values correspond to a global minimum in the fit, many of the parameters took similar values at the other minima that were located (Table 1). All runs that gave a fit of less than 0.2 (five of the eight) had a value of k in the range 55–68 m^1/3 s⁻¹, r _L /r _U in the range 1.38–1.44, an error in the height of lake U in the range 0–4 m and of 0 m for lake L and an ice-thickness scaling factor of (3.62 × 10 ^{− 3}) − (4.16 × 10 ^{− 3}). The absolute size of the lakes determines the total volume of water released in the jökulhlaup (as discussed in section 5.2). This is independent of the surface deformation, so it is not possible to infer the lake sizes from the available data.

5.4. A multi-lake model

So far, the investigation has modeled the discharge through a single conduit between two lakes. The altimetric data show the existence of four connected lakes: one source lake and three sink lakes (Fig. 1). The deformation above the three sink lakes is very similar, which is why we have so far simply modeled them as one lake. However, we now turn our attention to modeling this four-lake system using Nye’s theory. The extension of the theory to a multi-lake system is relatively straightforward. The depth of each lake is found for a particular time-step by solving the equations for the conduits feeding and draining the lake. The entire system of lakes and conduits is solved for each time-step, and the results are used to solve for the next time-step. The same time-steps are used as in the two-lake model above. We assume all conduits are semicircular, with an initial cross-sectional area of 1 m² and k = 60 m^1/3 s⁻¹. Figure 10 shows the results for four lakes of equal radius (14 km) at distances of 0, 205, 246 and 267 km along the Adventure subglacial trench from the source lake.

Fig. 10. Discharge between lakes for the four-lake model: (a) lake radii = 1 km, S ₀ = 1 m²; (b) lake radii = 14 km, S ₀ = 40 m² for the conduit between U1 and U2 and 1 m² for the other two conduits; (c) lake U1 radius = 8 km, other lake radii = 14 km, S ₀ = 40 m² for the conduit between U1 and U2 and 1 m² for the other two conduits.

One of the most striking results of introducing multiple lakes is the second discharge peak between lake L and lake U1 (Fig. 10a). The first discharge from lake L triggers a discharge from lake U1. This then empties lake U1 sufficiently to allow a second discharge from lake L. The altimetric data show no evidence for multiple stages of surface deformation above lake L and therefore no evidence of multiple discharges (Fig. 9). This may be because the conduit closes completely after the first discharge (which is prevented in the implementation of Nye’s theory). If there was insufficient pressure in lake L to reopen a conduit, only a single discharge would be observed. If this happened, lakes L and U1 would remain in a meta-stable state with a potential difference between them until the lake L water level rose sufficiently to trigger another flood.

A second interesting feature of the multi-lake model is the time lag of a few years between the discharge from lake L and the discharge from lake U1 (Fig. 10a). Again, this is not supported by the altimetric data, which show no obvious evidence of surface lowering of lake U1 after the initial rise (Fig. 9). More importantly, the surface uplift at U2 occurs at the same time as that above U1 (Reference Wingham, Siegert, Shepherd and MuirWingham and others, 2006, fig. 2). This apparent time lag is most likely a failing of the model. The size of the conduit draining out of lake U1 does not respond rapidly enough to changes in the depth of the lake. This may be explained by considering the triggering of the jökulhlaup. On a superficial level, the triggering is caused by the flotation of the ice dam. The initial size of the conduit is therefore related to the amount by which the ice dam is lifted. If the flotation occurs on a rapid timescale, then flotation effects will dominate frictional melting effects. This is not accounted for in Nye’s theory, where the initial conduit size is taken to have a fixed value. The likely effect of incorporating flotation effects can be seen by setting the initial size of the conduit leaving lake U1 to a much higher value than that of the conduit between lakes L and U1. Compared to Figure 10a, Figure 10b shows how the lag time is reduced when the initial size of the second conduit is set to 40 m² instead of 1 m².

A third feature of the multi-lake model is the broad width of the hydrograph from lake U1 lasting around 5 years (Fig. 10a and b). As mentioned above, this is not evident in the altimetric data, which show no surface lowering above U1 after the initial uplift (Fig. 9). Furthermore, the uplift at U2 lasts for around 2 years, not 5 years (Reference Wingham, Siegert, Shepherd and MuirWingham and others, 2006, fig. 2). The hydrograph width is strongly affected by the radii of the lakes, as can be seen in Figure 10c, showing the results when the radius of lake U1 is set to 8 km rather than 14 km. A fourth feature of the multi-lake model is that no discharge occurs between lakes U2 and U3 (Fig. 10). The reasons for this are not entirely clear, although it could be due to errors in the lake-height or ice-thickness values (i.e. errors in the BEDMAP data), which do not produce a large enough hydraulic potential gradient between U2 (when full) and U3 (when empty) to trigger a discharge. Further investigation is needed to confirm this.

6. Conclusions

We have used a simple energy-conservation model and a model based on Reference RöthlisbergerRöthlisberger’s (1972) theory to model steady-state water flow between two subglacial lakes in the Adventure subglacial trench, East Antarctica (sections 3 and 4). Constraining the models with available field evidence suggests that water flow could be accommodated in a tunnel with a cross-sectional area of ∼36 m² , comparable to that calculated by Reference Wingham, Siegert, Shepherd and MuirWingham and others (2006). Our results suggest, however, that the value for k (reciprocal of Manning’s roughness parameter) must be larger than the 12.5 m^1/3 s ^{− 1} calculated by Reference Wingham, Siegert, Shepherd and MuirWingham and others (2006).

We have also used Reference NyeNye’s (1976) theory to analyze the time-dependent transfer of water between the subglacial lakes (section 5). Our primary focus was to: (1) model water transfer between two lakes; (2) analyze the sensitivity of Nye’s theory to input parameters; and (3) use the theory in combination with a simple linear bed-to-surface ice deformation function to predict the values of these parameters, given observed surface elevation change data for the region. The main findings are summarized as follows.

1. 1. The results are highly sensitive to assumed values of k (Fig. 7). As k increases from 20 to 80 m^1/3 s ^{− 1}, peak discharges increase from ∼100 to 1000 m³ s ^{− 1} and the time to peak drops from >20 years to ∼5 years. Similarly, the water-level drop (rise) in the source (sink) lake changes from ∼25 m to >40 m and the peak conduit area doubles from ∼150 to ∼300 m².

2. 2. Results are also sensitive to assumed values of the initial conduit area (Fig. 8a). Peak discharges are not significantly affected, but the peak discharge timing drops from ∼15 to ∼7 years as initial conduit area increases from 1 to 7 m².

3. 3. A circular conduit reduces the timing of and increases the magnitude of peak discharge by ∼2 years and ∼50 m³ s ^{− 1}, respectively, compared to a semicircular conduit (Fig. 8b).

4. 4. Similarly, the timing of peak discharge drops by ∼1 year and its magnitude declines eight-fold as the ratio of source-lake to sink-lake radii changes from 1 : 1.5 to 1 : 0.5 (Fig. 8c).

5. 5. A Simplex-type approach is used to optimize the main parameters in Nye’s theory coupled to a simple vertical ice deformation function by matching model output to observed surface deformation data (Fig. 9). An ice-thickness scaling parameter is introduced to convert lake-level change to ice-surface elevation change. The best matches between modeled and observed data for a semicircular conduit of initial size 1 m² are obtained for values of k of 55–68 m^1/3 s⁻¹, an error of 0–4 m in the height of the sink lake and 0 m for the source lake, a ratio of lake radii (source/sink) of between 1 : 1.38 and 1 : 1.44 and an ice-thickness scaling factor of between 3.62 × 10 ^{− 3} and 4.16 × 10 ^{− 3}.

6. 6. The range of values for k supports the findings of the energy-conservation and Röthlisberger models, which predicted a value greater than that found by Reference Wingham, Siegert, Shepherd and MuirWingham and others (2006) of 12.5 m^1/3 s ^{− 1}. This range represents smooth to medium roughness ice-walled tunnels (Reference Röthlisberger, Lang, Gurnell and ClarkRöthlisberger and Lang, 1987) and compares with slightly lower values (representing slightly rougher conduits) required to model typical summer water flow beneath Haut Glacier d’Arolla, Switzerland, (k = 20 m^1/3 s⁻¹ _; Reference Arnold, Richards, Willis and SharpArnold and others, 1998) or jökulhlaups from Grímsvötn (k = 31 m^1/3 s⁻¹; Reference ClarkeClarke, 2003; Reference Flowers, Björnsson, Pélsson and ClarkeFlowers and others, 2004).

7. 7. Systems of more than two lakes are poorly modeled by Nye’s theory compared with observational evidence. Discharge between the first sink lake and the second sink lake occurs a long time after the initial flood from the source lake to the first sink lake. This suggests that the conduit from the second sink lake opens more rapidly than Nye’s theory predicts, which may be due to mechanical ice-flotation effects rather than thermal enlargement processes.