Modelling the evolution of subglacial tunnels due to varying water input

Abstract

The time evolution of a subglacial tunnel cross-section is examined usine a two-dimensional finite-element ice-flow model coupled to an idealized drainage system. Simulations are driven by physically based calculations of surface water-input variations at Slorgiaciaren, Sweden. Highlights of the model are its ability to handle unsteady conditions and irregular tunnel shapes. Agreement between modelled water pressure and borehole water levels is good. The following conclusions are reached: (i) Tunnels adapt to fluctuating inflow on time-scales of days. Storms, during which effective pressure ranges from 0 to 0.9 MPa, cause significant adjustments but daily fluctuations due solely to melt-water inflow are minor, (ii) Open-channel flow may become commonplace late in the ablation season, (iii) Initial tunnel shape influences subsequent tunnel evolution and seasonal water-pressure variation. Over the course of a summer, tunnels retain some of their initial shape, though in all experiments the width-to-height ratio increased with time, (iv) Tunnel contraction forms broad low tunnels. However, (v) given two tunnels of equal initial area, the higher narrower one expands more rapidly. Thus, more semi-circular tunnels may capture How from broader neighbours early in the summer.

Introduction

Water in a subglacial drainage system may move in a thin water film at the ice-bed interface (Weertman, 1972), through a porous till layer (Paterson, 1994), through a network of interlinked cavities (Walder, 1986; Kamb, 1987) or in tunnels (Rôthlisberger, 1972: Walder and Fowler. 1994). The extent to which some or all of these exist under a particular glacier may vary over time. In the case of tunnels, tracer tests have revealed their existence under parts of many valley glaciers during the ablation season (Stenborg, 1969: Behrens and others, 1975; Collins, 1982; Burkimisher, 1983; Iken and Bindschadler, 1986; Willis and others, 1990; Fountain, 1993; Nienow and others. 1996). Tunnels are of particular interesi because they provide the most direct link to conditions at the ice surface and they convey most of the water through glaciers.

The model introduced in this paper represents a new tool for integrating the influences of many glaciological parameters on tunnel evolution and water-pressure regimes beneath glaciers. Water-pressure variations are strongly correlated with ice motion (Iken and Bindschadler, 1986; Kamb and others. 1994; Iverson and others, 1995; Jansson, 1995). Thus, the seasonal evolution of the drainage system has a direct influence on ice dynamics (Rôthlisberger and Lang, 1987). In general, systems tend to expand at the onset of the ablation season, as rates of tunnel melting due to dissipation of viscous energy exceed rates of closure by ice creep (Shreve, 1972). Superimposed on ibis broad trend are short-term changes due to fluctuating meteorological conditions. Our understanding of the response time of the network to these changes is limiied. Calculations by Spring (1980) suggest that circular englacial conduits adjust to discharge fluctuations over periods of a few days. Tracer experiments indicate similar time-scales for subglacial channels (Hock and Hooke, 1993; Kohler, 1995). The manner in which water pressure within the system responds to a given input is dependent on the size and shape of the channel and hence on cumulative effects of previous inputs. Until now these cumulative effects have been difficult to track.

Existing theory for subglacial tunnel hydraulics

Most theoretical examinations of channelized subglacial drainage have assumed a steady state (Rôthlisberger, 1972; Shreve, 1972; Hooke and others, 1990; Alley, 1992; Walder and Fowler, 1994). Steady-state conditions are- most likely under the inner parts of huge ice sheets or during the winter beneath valley glaciers. Fluctuations in discharge and, consequently, in water pressure, are common during summer on valley glaciers and at the margins of ice sheets. Rapid release of ice-dammed lakes is perhaps the ultimate unsteady condition. This has been examined theoretically by Nye (1976), Spring and Hotter (1981), Clarke (1982) and Fowler and Ng (1996).

To estimate the tunnel closure rare. r_t , most theoretical analyses use the relation:

(1)

(Nye, 1953) where r_t is the tunnel radius, n is a constant from Glen's flow law (Glen, 1955), B is an ice-viscosity parameter and P_e is the effective pressure, equal to P_i = P_w, where P_w is water pressure and P_i is ice pressure. This relation was derived for a cylindrical channel but can be applied to a semicircular subglacial tunnel if the ice-bed contact is considered to be frictionless. The principal limitation of Equation (1) is that the more a tunnel shape diverges from semi-circular form, the more tenuous the predicted closure rates (Hooke and others, 1990). As outlined further below, there is evidence that tunnels may become broad and low. Thus, a significant feature of the tunnel model introduced in this paper is its accommodation of non-semi-circular shapes.

Tunnel shape

Using Equation (1), in conjunction with an assumption of circular conduits, Rôthlisberger (1972) found that an unusually low value of B was required to explain observed water pressures beneath Gornergletscher. This problem may be alleviated by the choice of a more realistic tunnel shape (Lliboutry. 1983; Iken and Bindschadler, 1986). An obvious starting point is a semi-circular tunnel carved into the ice and underlain by bedrock. In the steady state, such a configuration should maintain its shape if the tunnel remains full and there is no basal drag inhibiting tunnel closure. If basal drag is important, a broader low geometry may arise (Hooke and others, 1990). Lliboutry (1983) and Hooke (1984) raised the possibility of frequent periods during which water may flow in partially tilled conduits at atmospheric pressure. This would focus melting low on the conduit walls. Such a process also favours transformation of semi-circular tunnels into broad low ones. Adapting Rothlisberger's theory for use with tunnels of this shape, Hooke and others (1990) found an improved agreement between observed and calculated water pressures. Hock and Hooke (1993) also found that, by assuming the existence of broad low tunnels in an arborescent drainage network, measured water velocities could be matched by calculations. Furthermore, Walder and Fowler (1994) demonstrated the feasibility of broad “canals” in deformable subglacial sediments and Fowler and Ng (1996) have subsequently argued that development of broad low tunnels during jokulhlaups from Grimsvotn, Iceland, provides a viable explanation for the rapid termination of lake drainage. The question remains: at what rate might broad low conduits evolve in settings where discharge varies hourly? One piece of information required to answer this is the temporal and spatial variation of melt rates on tunnel walls.

Energy available for melting tunnel walls

Potential energy is released by water as it flows down a potential gradient. Defining s as the streamwise direction, then the amount of energy released per unit length Δs is

(2)

(Nye, 1976), where Q(t) is time-dependent discharge and ?/s is the potential gradient. The latter is obtained from

(3)

(Shreve, 1972), where ρ^w is the density of water, g is the gravitational acceleration, z is the elevation of the bed above a datum (Fig. 1) and P_w is water pressure.

Potential energy is used both to melt ice and to maintain water temperature at the pressure-melting point as ice thins in the downstream direction, thus

(4)

(Shreve, 1972) where m_p is the rate at which tunnel walls melt, ρ_i is the density of ice, L_f is the latent heal of fusion, l_i is the wetted perimeter of the ice part of the channel cross-section, C_w is the specific heat of water, c_t , is the change in the melting point per unit pressure and H_s is the height of the ice surface above the chosen datum (Fig. 1). Advection of heat (Nye, 1976) is ignored in Equation (4). The present application of the model to Storglaciaren, Sweden, allows this simplification, as calculations based on water-temperature measurements (Hock and Hooke, 1993) have indicated a decrease of only 1 % in the energy available for melt after correcting for advection, and discharge is always small (<0.6 m³ s^-1) in the section of tunnel examined here. Nonetheless, Spring (1983) has argued that an upper limit to tunnel expansion exists as a result of loss of melt energy by advection out of the system during jokulhlaups, thus caution is necessary when applying Equation (4) to other settings where tunnels carry greater discharge. Advection into the system is also neglected here, though this factor may become important when modelling tunnels draining warm reservoirs (Spring and Hutter, 1981).

Fig. 1. Coordinate system used to define z and H_s at any point along the stream-wise (s) direction.

Inserting Equation (3) into Equation (2), defining

where h is pressure head (Rôthlisberger, 1972), combining the result with Equation (4), dividing by Δs and re-arranging, results in

(5)

Where

K₁ has dimensions of m^-1, whilst K₂ is dimensionless. (The latter differs by a factor of ρ_w / ρ_i from K in the analysis of Rôthlisberger (1972).) Taking ρ = 900 kg m^-3,ρ_w = 1000 kg m^-3, L_f = 334 000 J kg^-1, g = 9.81ms^-2,C_w = 4180 Jkg^-1 K-1 and ct = 0.074 x 10 K Pa ', then ?_γ = 3.26 x 10 ^:' m ', and K₂ = 0.278. The variables Q_: e h and ff_s introduce time-dependence into Equation (5). The last of these can be considered constant over time-scales of months or less.

The following assumptions are implicit in Equation (5):

• (i) All physical quantities are uniform across the channel cross-section.

• (ii) Effects of curvature on the velocity field in the water are neglected.

• (iii) In a partially full conduit only the submerged part of the wall is melting.

• (iv) No energy is lost to or gained from the bed of the channel.

• (v) Instantaneous transfer of energy occurs between water and ice: turbulent mixing is efficient throughout the channel.

• (vi) Both ice and water are at the pressure-melting point.

Uniform energy dissipation over the channel perimeter (embodied in assumption (i)) might not exist near the margins of broad low tunnels, where shallow flow is retarded by enhanced frictional resistance from the bed. This assumption may lead to overestimated melt rates close to the edge of the conduit.

A Transient Model of Subglacial Tunnel Evolution

Because Equation (1) is inapplicable to non-semi-circular conduits, it is necessary to use numerical techniques to study how such conduits evolve. In the following discussion, I describe a two-dimensional transient finite-element model designed to investigate subglacial tunnel evolution. This is an adaptation of a scheme developed by B. Hanson (Hanson, 1990; Hanson and Hooke, 1994; Grace, 1995) which solves the Stokes equations of conservation of momentum and mass for an incompressible medium (Hanson, 1990). The model is used to study (i) the time-scales over which tunnels respond to changes in water input, (ii)the influence of tunnel evolution on temporal variations in sub-glacial water pressure, and (iii) the potential for development of broad low conduits beneath glaciers. In particular, the drainage system below the lower half of the ablation area of Storglaciarcn, Sweden, is examined.

The model domain

The model domain is a rectangular ice block, oriented transverse to the main ice-flow direction, with a half-tunnel located in the lower left corner (Fig. 2). The x direction is considered as down-glacier, y is transverse to ice motion and z is vertical: thus, the model simulates changes in the y plane. Node spacing in the y and z directions (Δy and Δz, respectively) is smallest in the immediate vicinity of the tunnel. Typically, 0.02 < Δy < 0.04 m, and 0.04 < Δz < 0.08 m. Maximum node spacings, Δy_max and Δz_max, are between 1.0 and 2.0 m. Choices of Δy, Δz, Δy_lwl* and Az_lunx arc governed by stability considerations, model accuracy and potential run times. All simulations discussed here used grids with approximately 1000 nodes. Run time fora 90 day simulation at 1 hour time steps was 2 days on a Silicon Graphics Indigo work station.

The model was tested against Nye's (1953) analytical solution (Equation (1)) for the case of a contracting semicircular tunnel. Using ice thicknesses ranging from 25 to 400 m, values of Δy and Δz from 0.01 to 0.10m, and values of Δy_max and Δz_max from 0.5 to 2.0 m, predictions of the model and Equation (1) agree to within 2%. Slight variability is introduced by the fact that a perfect semi-circular form is vertical at the ice-bed contact, whereas the discretized form must always have a slope less than vertical at this margin, because basal nodes cannot possess the same y coordinate. The disagreement can be reduced by decreasing node spacing.

Fig. 2. Example of a nodal grid configuration, with an enlargement of the region immediately surrounding the half tunnel. Only the basal nodes are shown.

Assumptions

The following assumptions are used in the model:

• (i) The ice is isothermal at the melting point.

• (ii) The ice is incompressible and therefore has a constant density.

• (iii) The ice is isotropic (constant B and n).

• (iv) One channel exists within the model domain and this channel has uniform geometry in the x direction.

• (v) Water pressure is distributed evenly over the entire ice water interface.

• (vi) Partially full conduits arc at atmospheric pressure.

• (vii) The ice lies directly over a flat undeformable bed.

• (viii) No frictional resistance to transverse ice motion is imposed by the bedrock.

• (ix) The normal stress is isotropic on the top of the model domain.

• (x) The stress field around the tunnel is symmetrical with the unseen half of the tunnel.

• (xi) Far-field transverse and longitudinal strain rates are zero.

• (xii) No exchange of water occurs between the tunnel and the ice-rock interlace,

• (xiii) The ice-flow direction is parallel to water flow.

• (xiv) No melt occurs above the water line from splashing during periods of open-channel flow.

As noted earlier, one common argument for the development of broad low conduits is that friction between ice and the bed limits the transverse How along the bed. In the absence of such friction, the calculated rate at which tunnels widen may be too low. On the other hand, the assumption of uniform energy dissipation may result in an overestimate of tunnel widening. Thus, confidence is placed only in qualitative trends in predicted tunnel evolution and resultant hydraulic characteristics. Fortunately, any influence of neglect of friction on tunnel evolution is reduced by the fact that calculated melt rates far exceed closure rates for at least the first half of all simulations.

The assumption that ice lies directly over bedrock ignores the possibility of an intervening till layer. Broad low-channel morphology may result from till creep and piping failure, and subsequent removal of material as bed load or suspended load (Walder and Fowler, 1994; Fowler and Ng, 1996). Inclusion of a till layer would therefore be desirable. Nonetheless, a number of fundamental questions can be explored with the present scheme. Formation of broad low tunnels in the absence of a till layer would only strengthen arguments for their existence if till were present. The same argument holds for predictions in the absence of frictional resistance to transverse ice flow.

Boundary and initial conditions

Hydrostatic pressure acts on the upper boundary of the model domain (Fig. 2) and water pressure within the tunnel can vary between zero and ice overburden. No vertical ice motion is allowed on the lower boundary outside the channel and no horizontal movement can occur on the left margin. Lastly, zero vertical shear stress is imposed on both lateral boundaries. Longitudinal and transverse strain rates ?_xx and ?_yy , respectively) can be applied as boundary conditions to the model. ?_xx acts parallel to the tunnel but it influences the effective strain rate within the y plane, and is therefore included (Hanson, 1990). During all runs discussed here ?_xx and ?_xy were zero, as earlier results turned out to be fairly insensitive to reasonable non-zero values (Culler, 1996). The influence ?f ?_xz on effective strain rate in the y— z plane is also neglected. No significant impact on model results is anticipated over the time-scales discussed here.

Hanson (1995) found that the ice-viscosity parameter B= 0.20 MPa a1/3 yielded the closest match between measured and modelled three-dimensional velocity distributions on Storglaciâren. This value is adopted as the standard value in model runs. In so doing, the potential influences on B from impurity content, non-random c-axis fabric, water content (Paterson, 1994) and transience in the stress field (Iken and Bindschadler, 1986) are ignored. A constant n(equal to 3) is also adopted (Hooke, 1981; Paterson, 1994).

Each run requires initial values for a number of parameters. Typical values are summarized in Table 1. The initial dimensions of the tunnel can be varied from a semi-circular case to any broad and low case. Alternatively, tunnel shape at the end of one simulation can be used as the initial condition for a later run. Lateral migration of the tunnel boundary occurs as follows; a node on the bed (Fig. 2) is considered part of the tunnel (and therefore susceptible to creep closure and melt) if a straight line joining it to the second node to its left lies below the first node to its left. Thus, the zero vertical-velocity boundary condition on basal nodes outside the tunnel is removed once a node satisfies this criterion and vice versa.

Table. 1. Standard parameter values used during testing of the model

Calculations of variable water level in the tunnel

The following information is specified or calculated at the beginning of each time step and, in the case of open-channel flow (e.g. Lliboutry, 1983; Hooke, 1984), it is used to determine water level: (i) inflow from the glacier surface, Q-_m; (ii) cross-sectional area of the tunnel, A; (iii) potential gradient, δφ/ds; (iv) hydraulic radius, R; (v) channel roughness, ».,; and (vi) outflow, Q₀nt- If no water is backed up in the system and Q_m is less than the maximum possible discharge as open-channel flow, Q_maKil, then Q_msl will equal Q\_u. (This makes the implicit assumption that the contribution of tunnel melt to Q_mi is insignificant (cf. Walder and Fowler, 1994).) Q,_M1, can be related to properties of the channel using the Caucklcr-Manning- Stricklcr equation;

(6)

(Roberson and Crowe, 1985). Channel roughness, n_m is a composite of the roughnesses of the bed and the overlying ice. For the specific case of Storglaciâren, a value of n_m = 0.20 m-1/3 s yields a mean velocity that agrees with tracer results (Seaberg and others, 1988; Hock and Hooke, 1993). Although high, this is not unusual in a glacial setting (Nye, 1976; Clarke, 1982). For a given Q_out there are two unknowns: A and R which must be determined using Newtons iterative method (Chow and others, 1988) if open-channel flow occurs.

Calculation of tunnel water pressure

Water pressure must be calculated for cases where discharge exceeds the open-channel capacity of the tunnel. The open-channel capacity is the maximum discharge with a free surface for a given bed slope, bed roughness and tunnel shape, and is calculated using Equation (6). It is attained when the flow depth is approximately 95% of tunnel height (American Society of Civil Engineers, 1982). The water-pressure calculation requires some assumptions about the configuration of the drainage network. The specific ease of water entering Storglaciâren at moulins Ml—M4 (Fig. 3) is examined. These moulins are formed due to crevassing over a subglacial bedrock high or riegel. As a first approximation of the system, it is assumed that a single crevasse receives the combined inflow to Ml-M4. Holmlund's (1988b) observations in the area containing Ml- M4 indicate that subsurface connections exist between moulins, usually along the line of the crevasse in which the moulins originated. Hence, it is not unreasonable to combine the inflow of Ml-M4 in the present scheme.

For now, it is assumed that the crevasse penetrates to the bed. This is unrealistic, based on evidence from Holmlund (1988b); however, more realism will be attempted later. The single crevasse is a reservoir that captures the calculated water influx to M1-M4 (Fig. 4). Its surface area, A_c , was estimated using information from Holmlund and Hooke (1983), Holmlund (1988b) and field observations in 1992 and 1993. A value of A_c =100m² was adopted. From the crevasse base, a single straight channel drains to the glacier terminus. This channel has a uniform gradient, dzfdb. over its complete length, s_tot.The tunnel model simulates channel evolution at a distance Δ_s from the base of the crevasse. lce thickness, (H_s — z), is constant in the vicinity of the tunnel model. It is assumed that the potential gradient, ?φ/ds, is constant along the tunnel and that water emerges at atmospheric pressure at the terminus. The influence on tunnel evolution of the assumption that hydraulic head may decline to zero at some point before the terminus is tested in the next section. Internal head losses within the reservoir, primarily at the outlet, are negligible compared to losses along the tunnel itself and are therefore ignored.

Fig. 3. The ablation area of Storglaciâren, with one possible path of water from the riegel moulins (M1-M4) to the terminus indicated by the thick dashed line. Two potential options for the last section of this route, marked “A” and “B”, are shown. The lower overdeepemng is located to the east of moulins M1-M4. Location of a borehole from which a water-level record was obtained in 1993 (Iverson and others, 1995) is indicated by “BH”.

Fig. 4. A simple configuration for the drainage system between M1-M4 and the terminus. The proportions of individual components are not to scale. Symbols are defined in the text.

Q_out varies with ∂φ/∂s (Equation (6)), which is time-dependent. The contribution of ∂z/∂s to ∂;φ/∂χ can be prescribed bin (he remaining part of ibis gradient, ∂h_c/∂s, is a function of the hydraulic head in the reservoir, h_c. To determine It,. during each time step, water is routed through the system in Figure 4 using a third-order Runge-Kutta numerical approximation of the continuity equation:

(7)

(Chow and others, 1988), where Q-_in is reservoir inflow and dS/df is the time rate of change of water storage in the reservoir. In a glacial system, the time-dependent nature of the outlet geometry results in a variable head-discharge relationship.

Testing The Model

The response of the model to systematic inflow variations (Fig. 5a) is examined first. The chosen 120day water-input regime is loosely based on a typical regime for the riegel moulins on Storglaciâren. Four large “storms”, each lasting 24 hours and spaced 7 days apart, are prescribed. A steady sinusoidal regime is maintained for 85 days after the final storm in order to investigate whether an equilibrium pressure regime and shape can be achieved. For clarity, only the first 60 days of model output are shown in Figure 5a-d. Predicted water-pressure variations for days 117-120 are shown in Figure 5e. Three runs were performed to test model sensitivity to choice of initial tunnel shape and maximum potential gradient. In two runs (cases 1 and 3), the tunnel began with a ratio of height (b) to half-width (w) (henceforth termed the b:w ratio) of 1:2 (Fig. 5d). The hydraulic head was assumed to be zero at the terminus in case 1, whilst in case 3 the head declines to zero halfway down the tunnel so that open-channel flow occurs in the last half. Cases 1 and 3 therefore bracket the extent of open-channel flow determined by Kohler (1995). The b:w ratio for the other run (case 2) was initially 1:4, though the initial cross-sectional area of the conduit, A, was identical in all three cases (0.023 m²). The choice of both A and the b:w ratio is arbitrary at this point. Possible constraints on these will be examined later on.

Figure 5b illustrates the calculated variation of A over the first 60 days of the three runs. During the first 9-21 days of this period -(δφ/θ,^/ρ^^ is steady at its maximum value (0.15 in cases 1 and 2, and 0.25 in case 3), and the tunnel expands because the melt rate exceeds the closure rate. Overflow from the crevasse is indicated by the flat tops of some peaks in Figure 5c, particularly early in the season. This reflects the fact that Q_in > Q_out even with the maximum head provided. Due to overflow, the first of the four prescribed storms does not influence tunnel evolution in cases 1 and 2. The water level, and hence water pressure.

P_W, eventually declines as the tunnel continues to expand (Fig. 5c). The drop occurs after 13 days in case 1, 16 days in case 2 and 6 days in ease 3. The simulated temporal variation of P_W is qualitatively similar to the regime under many valley glaciers (Rothlisberger and Lang, 1987), although the early season overflow occurs for what is probably an unreasonably long time. The cause of the lower maximum P_w in case 3 is the assumption of a more steeply negative potential gradient in conjunction with the fact that tunnel c\-olution is simulated at a point 100 m from the base of the reservoir.

Clearly, the effect of a steeper potential gradient on tunnel evolution outweighs the impact of initial shape, at least for the combinations of parameters selected for the runs in Figure 5. The tunnel in case 3 expands more rapidly than the others, owing to the larger melt rates from the sleeper (more negative) ?φ/ds and achieves a larger maximum area (Fig. 5b). Tunnel expansion in case 1 is more rapid than for the broader case 2, because of the smaller initial wetted perimeter along which melt energy is distributed. This suggests that tunnels with higher initial b:w ratios early in a melt season would tend to expand at the expense of those with lower initial b:w ratios, other factors being equal. The resulting lower pressures in these tunnels (Fig. 5c) would allow them to capture flow from initially broader tunnels, further enhancing their growth rate.

From Figure 5d, it is clear that choice of initial tunnel shape strongly influences its form after 60 days. However, shapes after 60 days are strikingly similar between the three cases. This leads to increasingly similar pressure regimes through time. In fact, all three cases are completely pressurized after 120days (Fig. 5e). At this time, only minor adjustments in pressure regime and cross-sectional area (not shown) occur.

Application of The Tunnel Model to Storglaciaren

A model of water input to moulins M1-M4 (Fig. 3) during 1993 is now used to drive tunnel evolution. Water inputs Q_in, was calculated using a distributed model in which melt and precipitation are routed through supraglacial drainage pathways to the moulins (Cutler, 1996). The accuracy of Q_in depends on uncertainties in three factors: melt, precipitation and influences on routing. In 1993, the mean différence between calculated melt and ablaiograph measurements was 4±20% for 13 1 day periods at a single site in the ablation area of Storglaciâren. The calculated glacier-wide summer balance was within 3±2% of measurements (Cutler, 1996). A tipping-bucket rain gauge was used to monitor rainfall intensity in the centre of the area draining to Ml-M4. The measurement uncertainty was estimated at±7% (Cutler, 1996), though this value may double when extrapolating results over the catchment of Ml-Ml. The sensitivity of Q_in to poorly constrained influences on water routing-such as extent of channelized flow at the base of the snow-pack, and vertical and lateral snow permeability, becomes less problematic as contributing area increases (Cutler, 1996), On the scale of the area draining to M1-M4, errors are primarily expected in the magnitude of estimated peak inflow rather than in timing. With the addition of uncertainties in calculated melt and precipitation, the error in Q_in is estimated at ±20%.

Equipotential surfaces in the ablation area of Storglaeiaren (Holmlund, 1988a) suggest that one major sub-glacial corridor for pressurized water should exist parallel to and close to the southern margin of the glacier. Existence of a second route through the centre of the lower overdeepening (Holmlund, 1988a) is not supported by borehole evidence (Hock and Hooke, 1993), which indicates that most subglacial water is probably usually diverted around the southern margin of the overdeepening. A possible route for water from Ml-M4 is plotted in Figure 3. The location of its upstream end is poorly constrained. However, Holmlund (1988b) noted that, from the bottom of moulins, water initially seems to flow transverse to the ice-flow direction along old crevasse lines. Assuming zero sinuosity and that the sub-glacial route follows option “A” (Fig. 3) near the terminus (indicated by tracer tests in 1992), the total length of the route is approximately 1000 m. ?ζ/?κ is significantly sleeper (more negative) during the last 250 m of this route than for the first 750 m. Hence, the assumption of uniform dz/Os is violated. The combination of increased release of potential energy and thinner ice here (Fig. 3) probably results in tunnels which are pressurized far less of the time than is the case further upstream (Kohler, 1995).

Fig. 5. The influence of initial tunnel shape and head gradient on tunnel evolution. a. Inflow regime. b. Temporal variation of A for three channel configurations. Case 1 has a b: W ratio of 1:2, where b = 0.15 m and a maximum head gradient, (dh_c/ds)_max, of 0.1. For case 2. h;w= 1:4 and (uhJUs)_lu;^ = -0.1. For case 3. b ur = 1:2 and (Oh_r/<h)_mn). = - 0.2. All cases have the same initial area. c. Temporal variation of P_w.d. Tunnel shapes at the start and after 60 days for each case. e. Variation of P_w after 117 days.

The values summarized in Table 1 are retained as attributes of an idealized tunnel beneath the lower part of the ablation area of Storglaciaren. A reservoir with straight sides (Fig. 4) and a surface area of 100m² is also retained, as halving or doubling this area has little impact on tunnel evolution because overflow occurs in all cases. An identical initial tunnel shape to that of case 1 in Figure 5 (with a b: w ratio of 1:2) was selected for the run. The significance of this, and the assumption of a straight-sided reservoir, is examined later. In the context of an arborescent drainage system such as that proposed by Hock and Hooke (1993), the tunnel simulated here would probably tarry about a One-quarter of the discharge of Sydjokk at the terminus (Fig. 3).

Fig. 6. a. Calculated water-input variations to M1-M4 during 1993 and seasonal variation of Q_max0 due to tunnel evolution in response to this input. b. Calculated water pressures during 1993. c. Measured borehole water levels at location BH (Fig. 3). Initial b: w ratio was 1 ; 2 with b=0.1m. The tick marks on the time axis in each panel are 1 week apart.

A water-level record exists from a borehole that intercepted the bed at a depth of 118m at location BH (Fig. 3) during 1993 (Iverson and others, 1995). This is used for comparison with model results. No information exists regarding the character of the bed at the base of the borehole, though the fact that the borehole drained immediately on reaching the bed, and that subsequent water-level fluctuations occurred on a sub-daily time-scale, suggests that the borehole maintained a strong hydraulic connection to a nearby subglacial tunnel.

Results

Figure 6a shows the calculated water input to moulins Ml-M4 during 1993, along with the seasonal variation of calculated maximum open-channel discharge through the sub-glacial tunnel, Q_max . Figure 6b illustrates the concurrent seasonal variation of calculated P_W in the tunnel. This can be compared with water levels measured at location BH (Fig. 6c). The amplitude of measured water-level fluctuations is lower than the calculated amplitude of P_w . The former may have been damped by diffusion as the pressure wave migrated through a porous subglacial till layer or smaller conduits. Hooke and others (1989, Fig. 4a) have presented a record of water-level fluctuations in a borehole (83-2) that directly intercepted a subglacial channel in the vicinity of BH 10 years earlier. Their data are comparable in character to P_w in Figure 6b.This supports the suggestion that the record from BH is damped.

The following similarities between the measured water Level and P_w in Figure 6 are notable, (i) The timing and duration of pressure peaks during storms C through G match measurements closely. For example, the model succeeds in calculating broader peaks in P_w during storms D and F than during storms E and G. Additionally, simulated overflow during storm F is supported by observations of overflow from crevasses near BH. (ii) Lower peaks in ?_w during storms G and H than in storms D, E and F mimic the muted response in the measurements. (iii) The gradual climb of P_K up to storm G, with sub-daily oscillations super imposed on the signal matches the trend in water level. Oscillations in measurements and calculations are almost perfectly synchronous on an hourly time-scale, (iv) P_w is generally higher before storm G than after storms D and E, as are measured water levels. Furthermore, Q_max0 peaks at about the same time that the borehole water level reaches its seasonal minimum (Fig. 6a). The system apparently achieved its maximum efficiency at this time. In support of this prediction, the transition from early to late-season drainage configuration usually occurs in late July/early August (Seaberg and others, 1988; Hock and Hooke, 1993). This, and the preceding agreement between observations and calculations, allows us to draw a pair of solid conclusions. First, tunnel adjustment in response to storm inflow occurs over 2 or 3 days once the system has undergone its early season growth (Fig. 6a). This time-scale is in agreement with the calculations of Spring (1980). Secondly, open-channel flow is apparently possible for about 50% of the time from mid-July to mid-August, and then for most of the time until early September, as tunnel closure did not keep paee with declining water inputs.

The results in Figure 6 support suggestions by Rothlisberger (1972) and Shreve (1972), and calculations by Spring (1980), that discharge and P_w vary in phase with each other under transient conditions of diurnally varying discharge. Further broad agreement exists with Rothlisberger's scenario for phase differences between P_w and discharge over longer time-scales. Specifically, he proposed that low discharge is associated with high P_w through the spring and early summer whilst tunnels continuously expand. As discharge and tunnel size peak in mid- or late-summer, so P_w declines, on average, according to Rothlisberger. Finally, P_K falls more rapidly than discharge towards the end of the summer, as tunnel closure lags behind falling water levels. Predictions by Spring (1980) are also borne out in Figure 6. He calculated that time-scales for tunnel expansion are shorter than for closure and that diurnal fluctuations in P_w tend to be larger during periods of increasing mean discharge compared to periods of decreasing mean discharge. The latter pattern is seen in both calculated P_w and measured water level during the last 3 weeks of July.

Discussion of disagreements between modelled P_w and measured water level is warranted. Values of P_w before storms A and B (Fig. 6b) remain close to or at their maximum value, whereas the measured water level rose from approximately 40 m prior to storm A, peaked at 110 m (approximately 5 m above flotation) and then dropped back to about 40 m before storm B. Model results suggest that water overflowed from the reservoir for about 3 weeks before storm A. Such prolonged overflow is unrealistic, according to the observations by Holmlund and Hooke (1983). Early season water pressures may be around half overburden (Hooke and others, 1989) rather than the selected initia] value of zero. Adopting ibis as a more realistic initial condition enhances early season tunnel expansion and reduces by a few days the duration of overflow. Initial channel shape and size also influence the persistence of early season overflow (Fig. 5). These controls are further explored in the next section. A number of other factors could also reduce the duration of overflow, including an increased reservoir capacity, temporary storage of some water at the ice-bed interface and a steeper potential gradient in the tunnel. Additionally, the early season hydraulic system may consist of a number of smaller channels which allow better reservoir drainage than a single tunnel. The system may become organized into fewer, larger, tunnels as die summer progresses.

Constraints on the choice of initial tunnel shape and cross-sectional area

One approach to constraining initial tunnel shape and size is to estimate how contraction during the winter alters these parameters by the time the next ablation season begins. As noted earlier, the initial tunnel shape used for the run in Figure 6 had a b:w ratio of 1: 2. Figure 7a illustrates die evolution of this tunnel through periods of summer expansion and subsequent winter contraction, Figure 7b shows the evolution of an initially semi-circular tunnel with the same initial cross-sectional area as in Figure 7a, through which the same inflow was routed. After 88 days, both tunnels retained some evidence for their initial form. Thereafter, Q_in dropped to zero. Trends in borehole water levels reported by Hooke and others (1989, Fig. 3a, borehole 83-6) show that P_w gradually increases to approximately half overburden pressure during the winter. Irregularities in the tunnel system, which might lead to blockages and the subsequent build-up of P_w , are not addressed by the model. Thus, P_w was artificaily raised from zero to half overburden pressure over the next 240 days. After 157 and 188 days, in the cases of Figure 7a and b, respectively, there is little distinction between the two tunnels. Both develop much higher aspect ratios but unfortunately close completely before the end of the year. Thus, a reasonable initial tunnel size or shape cannot be determined from these experiments. More rapidly increasing water pressure in the autumn could help a tunnel sur-vive through the winter. This may occur in the natural system due to enhanced tunnel closure on the upstream side of bumps, particularly as water inputs decline. Alternatively, a constant inflow from ground-water sources might allow survival of some tunnels (Lliboutry, 1983). However, this may not be applicable to Storglaciaren, as Stenborg (1965) observed no water flowing out of the southern drainage stream, Sydjokk, in winter. Nevertheless, if some sections of the tunnel do survive, model results indicate that these may be broad and low. This conclusion is apparently applicable even to initially semi-circular tunnels (Fig. 7b).

Fig. 7. Comparison of the seasonal evolution of two tunnels both receiving the same inflow (Fig. 6b) for the first 88 days and then closing with no further inflow passing through them. a. A tunnel with a b: w ratio of 1:2. The shape in panel (a) was used in the run displayed in Figure 6b. b. An initially semi-circular tunnel. Initial cross-sectional area were 0.027 and0.023m² in (a) and (b), respectively.

Another approach to the problem of determining initial tunnel shape might be to find a shape that results in a better match between P_K and measured water levels at BH. Figure 8 illustrates such an experiment for tunnels with initial cross-sectional areas equal to those in Figure 6 but with initial b : w ratios of 1:1 and 1:4. The validity of this approach depends greatly on the suitability of the approximation of the drainage system, but we will proceed. It is apparent that the case of b:w = 1:1 does a better job of matching water-level fluctuations in the period containing storms A through D. This is also an improvement over the case in Figure 6b. However, nearly all of the inflow during storms F and G is accommodated as open-channel flow (Fig. 8c), whereas measured water levels suggest pressurized flow at this time. At the other extreme, the case of b :w= 1:4 remains at maximum P_w for too long (until storm D) but matches measured water levels better than b:w = 1:1 in the later half of the summer. The most appropriate initial tunnel shape may therefore lie within the range tested here. It is noteworthy that both cases in Figure 8 attain a similar cross-sec-tional area by the time storm G occurs, so différences in calculated P_w at this time are alone due to differences in tunnel shape. The lower hydraulic radius of the taller, narrower tunnel allows it to evacuate water more efficiently than the broader, lower tunnel.

Comparison of tunnel growth beneath Storglaciaren in 1992 and 1993

Seasonal tunnel evolution in 1992 and 1993 is now compared. A more realistic reservoir geometry is used in these experiments. This involves retaining a uniform cross-sectional area to a depth of 35 m (based on observations from Holmlund (1988b)) but then linearly decreasing the area from 35 m to the bed. Basal area is initially 0.5 m² but this is allowed to grow as tunnel width increases. The cross-sectional area of the top of the reservoir is set at 150 m ² to adjust for the decrease in storage capacity deeper in the glacier (to avoid even greater overflow problems early in the season). The reservoir is still a crude approximation of reality but its form (its better with Holmlund's (1988b) observations.

Figure 9 illustrates the calculated seasonal evolution of a tunnel in 1992 and 1993. Initial conditions were identical for both years, with a b:w ratio of 1:2. Calculated inflow in early and late 1992 (Fig. 9a) is on a daily time-scale due to lack of meteorological measurements on the glacier surface (Cutler, 1996). This does not seriously compromise results, because water pressure lends to respond lo input variations on time-scales of more than a day early and late in the season. No borehole water-level records were obtained in 1992. Thus, the simulated variation of P_w cannot be assessed. However, dye-tracer experiments were performed on live occasions during the summer liable 2). Though diurnal variations in dispersivity and velocity (Behrens and others, 1975; Collins, 1982; Hock and Hooke, 1993; Nienow and others, 1996) make interpretation of such a small dataset difficult, the injection times are similar enough that these factors do not alter the broad conclusion that the drainage system became increasingly efficient between 8 and 24 July. The model calculates a doubling of A to its peak value during this period (Fig. 9b.

Fig. 8. The influence of initial tunnel shape (either b:w =1:1 or b: w = 1:4) mi the agreement, between (a) measured borehole water levels and (b) calculated water pressures. (c) Seasonal variation of cross-sectional area of the two tunnels.

The pattern of calculated channel evolution is similar in 1992 and 1993. An initial period of high P_w and continuous expansion lasting 3-4 weeks (Fig. 9b and d) precedes a period of 4-6 weeks during which the channel expands further but primarily in response to repealed large storms. Higher melt inputs and more intense storms in the first half of the summer of 1993 cause A to surpass 1992 values by the end of July. On days without major contributions from rainfall, tunnels do not respond in a significant way to sub-daily variations in water input. No stable channel size is attained in either year. Instead, a single peak m A is Followed by gradual closure as water inputs decline towards the end of the summer. Closure is less rapid in 1992, because of occasional water inputs in October. Calculated tunnel shapes on 7 September in 1992 and 1993 are compared in Figure 9e, Despite different inflow regimes, the two tunnels are surprisingly similar in form, furthermore, while larger in size, neither is very different from its original shape. The 1992 tunnel is slightly broader and lower because open-channel flow first occurred approximately 1 month earlier in that year.

According to the model, pressurized flow dominates the early season regime in both years. As the summer progresses, there is increased potential for open-channel flow, even in tunnels close to the riegel. Kohler (1995) has sug-gested that water may flow under pressure over most of the distance from the riegel moulins to the terminus. His inferences do not conflict with model calculations, because his dataset extends only until early August, when water inputs arc still high. Only as inputs tail off are long periods of open-channel flow expected.

Fig. 9. A comparison of calculated tunnel evolution during 1992 and 1993. Panels (a) and (c) display modelled water-input variations into moulins M1-M4 in 1992 and 1993, respectively. Panels (h) and (d) contain simulated changes in water pressure and cross-sectional area during 1992 and 1993, respectively. Initial conditions were identical in both years (see Table 1 and text). The tick marks on the time axis are at 1 week intervals. Panel (e) contains a comparison of estimated tunnel shapes on 7 September 1992 and 1993.

Conclusions

A model has been developed to examine the seasonal evolution of a subglacial tunnel in response to fluctuating water inputs from the surface. The model can deal with (i) unsteady flow, and (ii) non-uniform tunnel shapes. Despite a crude drainage configuration composed of a reservoir feeding a single subglacial tunnel, the model is able to simulate many of the features in a record of borehole water levels from Storglaeiiiren.

The following model results are considered significant to our understanding of glacier hydrology. (i) Tunnels evolve on time-scales of days in response to fluctuating inflow. Storms cause the most significant change in the system, whilst typical daily fluctuations, due solely to meltwater inflow, play a minor role. Early in the ablation season, even the smaller inputs become backed-up in the inefficient drainage system, causing continuous, accelerating, tunnel expansion. During this period, die drainage system probably consists of many smaller channels which are replaced by fewer larger channels later in die summer. (ii) After the initial period of net tunnel expansion, lasting about 2 months, an abrupt change to contraction is indicated. However, (iii) the rate of closure is insufficient to keep pace with declining water inputs and open-channel flow may become common in the later half of the ablation season. This is suggested even for tunnels with initial width-to-height ratios as high as 8:1.(iv)Initial tunnel shape significantly influences subsequent tunnel evolution and, hence, daily and seasonal variation in water pressure. Over the course of a single ablation season, tunnels retain some semblance of their initial shape, though in all experiments the end result possessed a higher width-to-height ratio. However, a single summer is insufficient time for the development of channels with heights on the order of 0.1 in and widths of > 10 m as envisaged by Hock and Hocke (1993). (v) Contraction increases the width-to-height ratio of tunnels. These forms should be retained until the next ablation season if water pressure rises rapidly enough to prevent complete closure during the winter. However, (vi) given two tunnels of equal initial area, the one with a lower width-to-height ratio will expand more rapidly. Thus, more semi-circular tunnels may capture flow from broader neighbours early in the summer, thus limiting the latter's importance to subglacial drainage.

Table. 2. Summary of tracer experiments during 1992