Skip to main content Accessibility help


  • Access
  • Cited by 6



      • Send article to Kindle

        To send this article to your Kindle, first ensure is added to your Approved Personal Document E-mail List under your Personal Document Settings on the Manage Your Content and Devices page of your Amazon account. Then enter the ‘name’ part of your Kindle email address below. Find out more about sending to your Kindle. Find out more about sending to your Kindle.

        Note you can select to send to either the or variations. ‘’ emails are free but can only be sent to your device when it is connected to wi-fi. ‘’ emails can be delivered even when you are not connected to wi-fi, but note that service fees apply.

        Find out more about the Kindle Personal Document Service.

        Numerical analysis of rapid water transfer beneath Antarctica
        Available formats

        Send article to Dropbox

        To send this article to your Dropbox account, please select one or more formats and confirm that you agree to abide by our usage policies. If this is the first time you use this feature, you will be asked to authorise Cambridge Core to connect with your <service> account. Find out more about sending content to Dropbox.

        Numerical analysis of rapid water transfer beneath Antarctica
        Available formats

        Send article to Google Drive

        To send this article to your Google Drive account, please select one or more formats and confirm that you agree to abide by our usage policies. If this is the first time you use this feature, you will be asked to authorise Cambridge Core to connect with your <service> account. Find out more about sending content to Google Drive.

        Numerical analysis of rapid water transfer beneath Antarctica
        Available formats
Export citation


We use a simple energy-conservation model and a model based on Röthlisberger’s theory for steady-state water flow in a subglacial conduit to model water movement between lakes in the Adventure subglacial trench region of East Antarctica during a 1996–98 jökulhlaup. Using available field evidence to constrain the models suggests that water flow would likely be accommodated in a tunnel with a cross-sectional area of 36 m2 and a value for k (the reciprocal of Manning’s roughness parameter) larger than the 12.5 m1/3 s−1 previously calculated. We also use Nye’s theory for time-dependent conduit water flow to model the temporal evolution of conduit discharge, cross-sectional area, water pressure and lake draining and filling during the flood. We initially assume one source and one sink lake. We perform sensitivity tests on the input parameter set, matching modeled source- and sink-lake depth changes with measured surface elevation data. Using a simple function for vertical ice deformation in which surface deformation scales linearly to the lake depth change, we find the scaling factor is of the order 4 × 10−3 of the ice thickness. The most likely value of k lies in the range 55–68 m1/3 s−1, and the ratio of source to sink-lake radii is approximately 1 : 1.4. Finally, we experiment using Nye’s theory to model water movement between one source and three sink lakes. The model fails to produce the observed patterns of water movement as indicated by the surface deformation data.

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 (Tweed and Russell, 1999). The best-studied jökulhlaups are those from the subglacial caldera lake Grímsvötn, beneath Vatnajökull, Iceland (Bjö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 (Björnsson, 1988, 1992). Water is released from the lake and flows subglacially to the ice-cap margin (Bjö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 Rö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.

Nye’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’ (Nye, 1976, p.194), although the theory does tend to overestimate peak and total discharge.

Spring 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.

Clarke (1982) added to Nye’s theory by including the effects of reservoir geometry and lake temperature. A major advance was made by Clarke (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 Fowler (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. Fowler 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 (Flowers 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 km2, drainage distances are about 50 km, ice thicknesses are typically 300–600 m (e.g. see cross-section in Bjö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 km2 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 (Wingham and others, 2006; Fricker and others, 2007).

Wingham 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 km2 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 Siegert and others, 2005) are denoted by triangles. Approximate locations of the source lake (L) and sink lakes (U1, U2 and U3) from Wingham 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 Wingham 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 Wingham 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 Wingham and others (2006). We begin with a simple energy-conservation model and a model based on Rö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 Nye’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 m2, 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 m2, 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. Wingham and others (2006) were able to constrain values for these parameters based on observational evidence to V = 1.8 km3, l = 290 km and Δh L = 3 m. The values of ΔH (difference in bed height) and Δz (z Lz U, i.e. change in ice thickness) are less well constrained by available data. Errors in excess of ±30 m are possible (Lythe 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 Wingham and others (2006). For ΔH = 262 m and Δz = 449 m, a cross-sectional area of 36 m2 is predicted. This is in fairly good agreement with the 28 m2 calculated by Wingham 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 m1/3 s−1 and k = 30 m1/3 s−1. 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


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


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 Wingham 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.

Wingham and others (2006) predicted a cross-sectional area of 26 m2 for a semicircular conduit. This was based on ΔH = 262 m, Δz = 449 m, Q = 50 m3 s 1 and k = 12.5 m1/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 m2 for a semicircular conduit and 91 m2 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 m1/3 s 1 (M = 0.02 m−1/3 s) produces a conduit area of 36 m2, 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 m1/3 s 1 and the much larger region of allowed values for k = 30 m1/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 (Nye, 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 (Nye, 1976, equation (20)). We are also able to simplify the energy equation (Nye, 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 × 105 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:


where m is the melting rate (kg m 1 s 1). Spring 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.



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:



where A L and A U are the areas of lake L and lake U, respectively. Using Equations (57), 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 Nye’s (1976) equation (18):


where gs 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 0, 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 Wingham 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 m1/3 s−1, the initial cross-sectional area is set to 1 m2, 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 km2, the area of lake L estimated by Wingham and others (2006) from their observational data.

Wingham and others (2006) estimated that k could lie anywhere in the range 10–100 m1/3 s−1. Figure 7 shows how the predicted discharge hydrograph and change in lake depth vary with k. As k increases from 20 to 80 m1/3 s−1, peak discharge increases from 100 to 1000 m3 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 m2 (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 0. 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 0 increases from 1 to 7 m2. 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 0;(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 Wingham 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 Wingham 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 (Wingham 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 (Jarosch 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.


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 (Björnsson, 1988, 1992), 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; Guð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. MacKay, 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 (Vanderbei, 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 m2, 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 m1/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 (Lythe 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 m1/3 s−1, 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 m1/3 s−1, 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 m2 and k = 60 m1/3 s−1. 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 0 = 1 m2; (b) lake radii = 14 km, S 0 = 40 m2 for the conduit between U1 and U2 and 1 m2 for the other two conduits; (c) lake U1 radius = 8 km, other lake radii = 14 km, S 0 = 40 m2 for the conduit between U1 and U2 and 1 m2 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 (Wingham 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 m2 instead of 1 m2.

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 (Wingham 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 Rö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 m2 , comparable to that calculated by Wingham 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 m1/3 s 1 calculated by Wingham and others (2006).

We have also used Nye’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 m1/3 s 1, peak discharges increase from ∼100 to 1000 m3 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 m2.

  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 m2.

  3. 3. A circular conduit reduces the timing of and increases the magnitude of peak discharge by ∼2 years and ∼50 m3 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 m2 are obtained for values of k of 55–68 m1/3 s−1, 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 Wingham and others (2006) of 12.5 m1/3 s 1. This range represents smooth to medium roughness ice-walled tunnels (Rö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 m1/3 s−1 ; Arnold and others, 1998) or jökulhlaups from Grímsvötn (k = 31 m1/3 s−1; Clarke, 2003; Flowers 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.


We thank R. Hindmarsh, T. Schuler, M. Shinwell and S. Taraskin for advice about various numerical and computational aspects of the work. We are grateful to R. Hindmarsh and G. Leysinger Vieli for providing us with unpublished basal ice temperature calculations from their balanced ice-flux flow model, and for discussing the implications of these for modeling water flow beneath the Adventure subglacial trench. Useful discussions with H. Guðmundsson helped to clarify how best to model the relationship between basal and surface vertical ice deformation simply. We thank M. Siegert for supplying the altimetry data reproduced in Figure 9. G. Flowers and R. Bell reviewed the paper; their detailed and constructive comments helped improve it significantly.


Arnold, N., Richards, K., Willis, I. and Sharp, M.. 1998. Initial results from a distributed, physically-based model of glacier hydrology. Hydrol. Process., 12(2), 191219.
Björnsson, H. 1988. Hydrology of ice caps in volcanic regions. Vísindafélag Ísl. Rit 45.
Björnsson, H. 1992. Jökulhlaups in Iceland: prediction, characteristics and simulation. Ann. Glaciol., 16, 95106.
Björnsson, H. 1998. Hydrological characteristics of the drainage system beneath a surging glacier. Nature, 395(6704), 771774.
Björnsson, H. 2003. Subglacial lakes and jökulhlaups in Iceland. Global Planet. Change, 35(3–4), 255271.
Björnsson, H. and Guðmundsson, M.T.. 1993. Variations in the thermal output of the subglacial Grímsvötn caldera, Iceland. Geophys. Res. Lett., 20(19), 21272130.
Clarke, G.K.C. 1982. Glacier outburst floods from ‘Hazard Lake’, Yukon Territory, and the problem of flood magnitude prediction. J. Glaciol., 28(98), 321.
Clarke, G.K.C. 2003. Hydraulics of subglacial outburst floods: new insights from the Spring–Hutter formulation. J. Glaciol., 49(165), 299313.
Flowers, G.E., Björnsson, H., Pélsson, R. and Clarke, G.K.C.. 2004. A coupled sheet–conduit mechanism for jökulhlaup propagation. Geophys. Res. Lett., 31(5), L05401. (10.1029/2003GL019088.)
Fowler, A.C. 1999. Breaking the seal at Grímsvötn, Iceland. J. Glaciol., 45(151), 506516.
Fowler, A.C. and Ng, F.S.L.. 1996. The role of sediment transport in the mechanics of jökulhlaups. Ann. Glaciol., 22, 255259.
Fricker, H.A., Scambos, T., Bindschadler, R. and Padman, L.. 2007. An active subglacial water system in West Antarctica mapped from space. Science, 315(5818), 15441548.
Guðmundsson, G.H. 2008. Analytical solutions for the surface response to small amplitude perturbations in boundary data in the shallow-ice-stream approximation. Cryosphere, 2(2), 7793.
Jarosch, A.H. and Guðmundsson, M.T.. 2007. Numerical studies of ice flow over subglacial geothermal heat sources at Grímsvötn, Iceland, using Full Stokes equations. J. Geophys. Res., 112(F2), F02008. (10.1029/2006JF000540.)
Lythe, M.B., Vaughan, D.G. and BEDMAP consortium. 2001. BEDMAP: a new ice thickness and subglacial topographic model of Antarctica. J. Geophys. Res., 106(B6), 11,33511,351.
MacKay, D.J.C. 2003. Information theory, inference and learning algorithms. Cambridge, etc., Cambridge University Press.
Nye, J.F. 1976. Water flow in glaciers: jökulhlaups, tunnels and veins. J. Glaciol., 17(76), 181207.
Röthlisberger, H. 1972. Water pressure in intra- and subglacial channels. J. Glaciol., 11(62), 177203.
Röthlisberger, H. and Lang, H.. 1987. Glacial hydrology. In Gurnell, A.M. and Clark, M.J., eds. Glacio-fluvial sediment transfer: an alpine perspective. Chichester, etc., Wiley, 207284.
Siegert, M.J., Carter, S., Tabacco, I., Popov, S. and Blankenship, D.D.. 2005. A revised inventory of Antarctic subglacial lakes. Antarct. Sci., 17(3), 453460.
Spring, U. and Hutter, K.. 1982. Conduit flow of a fluid through its solid phase and its application to intraglacial channel flow. Int. J. Eng. Sci., 20(2), 327363.
Tweed, F.S. and Russell, A.J.. 1999. Controls on the formation and sudden drainage of glacier-impounded lakes: implications for jökulhlaup characteristics. Progr. Phys. Geogr., 23(1), 79110.
Vanderbei, R.J. 2001. Linear programming: foundations and extensions. Second edition. New York, etc., Springer.
Wingham, D.J., Siegert, M.J., Shepherd, A. and Muir, A.S.. 2006. Rapid discharge connects Antarctic subglacial lakes. Nature, 440(7087), 10331036.

Appendix A. Energy-Conservation Derivation

Assuming a cylindrical geometry for both the lakes and the ice that moves above them, the energy released by the system shown in Figure 2 can be written as the sum of the following four terms.

  1. 1. Energy released by the lowering of lake L:


  2. 2. Energy used to raise lake U:


    Assuming that the volume of water is conserved, it follows that:


    and therefore


  3. 3. Energy used to raise the water by ΔH:


  4. 4. Energy used to maintain water at the pressure-melting point:


Subscripts L and U refer to the lower and upper lakes, respectively, r is lake radius, t is ice thickness, Δh is surface elevation change, ρi = 917 kg m 3 is the density of ice, ρ w = 1000 kg m 3 is the density of water, g = 9.81 m s−2 is gravitational acceleration, ΔH is bed height difference, cv = 4.2 kJ kg 1 is the specific heat capacity of water, T = 273 K is the melting point, L = 3.3 × 105 J kg 1 is the latent heat of melting for ice and V is the total volume of water discharged, estimated to be 1.8 km3 by Wingham and others (2006). V neglects the contribution from the melted walls of the conduit, as the volume of a 290 km long conduit of cross-sectional area 36 m2 is 10.44 × 10 3 km3 (or 0.58% of the total discharge).

The total energy released is then:


For temperate ice, the energy required to melt a conduit of cross-sectional area S is


where l is the length of the conduit (290 km in this case). Rearranging gives:


Appendix B. Röthlisberger Derivation

The idealized geometry of the lakes and conduit is shown in Figure 2. Röthlisberger’s theory is based on the Gauckler–Manning–Strickler formula for the mean velocity of a turbulent flow:


where Q is the discharge (m3 s 1), k is the inverse of the Manning roughness parameter, Φ is the hydraulic potential and R is the hydraulic radius, defined as cross-sectional area over wetted perimeter.

Assuming a constant hydraulic potential gradient along the conduit,


The hydraulic potential in each lake can be expressed as:



where t is ice thickness above the lake, w is the head of water in the lake and z is the elevation of the lake bed.

Assuming the lakes are the same depth and substituting Equations (B3) and (B4) into Equation (B2), the gradient may be written:


where Δt is the difference in ice thickness between the lakes and ΔH is the difference in bed height between the lakes.

For a conduit of circular cross-section,



Substituting Equation (B7) into Equation (B1) then gives


For a conduit of semicircular cross-section,



Substituting Equation (B10) into Equation (B1) then gives


Equations (B8) and (B11) correspond to Equations (1) and (2) in the main part of the text.