Existing theory for subglacial tunnel hydraulics

Most theoretical examinations of channelized subglacial drainage have assumed a steady state (Reference RôthlisbergerRôthlisberger, 1972; Reference ShreveShreve, 1972; Reference Hooke and KohlerHooke and others, 1990; Reference AlleyAlley, 1992; Reference Walder and FowlerWalder 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 Reference NyeNye (1976), Reference Spring and HutterSpring and Hotter (1981), Reference ClarkeClarke (1982) and Reference Fowler and NgFowler and Ng (1996).

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

(1)

(Reference NyeNye, 1953) where r_t is the tunnel radius, n is a constant from Glen's flow law (Reference GlenGlen, 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 (Reference Hooke and KohlerHooke 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, Reference RôthlisbergerRô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 (Reference LliboutryLliboutry. 1983; Reference Iken and BindschadlerIken 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 (Reference Hooke and KohlerHooke and others, 1990). Reference LliboutryLliboutry (1983) and Reference HookeHooke (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, Reference Hooke and KohlerHooke and others (1990) found an improved agreement between observed and calculated water pressures. Reference Hock and HookeHock 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, Reference Walder and FowlerWalder 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)

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

(3)

(Reference ShreveShreve, 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)

(Reference ShreveShreve, 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 (Reference NyeNye, 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 (Reference Hock and HookeHock 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, Reference SpringSpring (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 (Reference Spring and HutterSpring 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 (Reference RôthlisbergerRô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 Reference RôthlisbergerRô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.

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 Reference NyeNye'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 (Reference Walder and FowlerWalder and Fowler, 1994; Reference Fowler and NgFowler 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 (Reference HansonHanson, 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.

Reference HansonHanson (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 (Reference PatersonPaterson, 1994) and transience in the stress field (Reference Iken and BindschadlerIken and Bindschadler, 1986) are ignored. A constant n(equal to 3) is also adopted (Reference HookeHooke, 1981; Reference PatersonPaterson, 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. Reference LliboutryLliboutry, 1983; Reference HookeHooke, 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. Reference Walder and FowlerWalder and Fowler, 1994).) Q,_M1, can be related to properties of the channel using the Caucklcr-Manning- Stricklcr equation;

(6)

(Reference Roberson and CroweRoberson 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 (Reference Seaberg, Seaberg, R. LeB and WibergSeaberg and others, 1988; Reference Hock and HookeHock and Hooke, 1993). Although high, this is not unusual in a glacial setting (Reference NyeNye, 1976; Reference ClarkeClarke, 1982). For a given Q_out there are two unknowns: A and R which must be determined using Newtons iterative method (Reference Chow, Maidment and MaysChow 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 Reference Holmlund and HookeHolmlund 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)

(Reference Chow, Maidment and MaysChow 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.

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 (Reference Seaberg, Seaberg, R. LeB and WibergSeaberg and others, 1988; Reference Hock and HookeHock 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 Reference SpringSpring (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 Reference ShreveShreve (1972), and calculations by Reference SpringSpring (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 Reference SpringSpring (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 Reference Holmlund and HookeHolmlund and Hooke (1983). Early season water pressures may be around half overburden (Reference Hooke, Calla, Holmlund and StroevenHooke 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 (Reference LliboutryLliboutry, 1983). However, this may not be applicable to Storglaciaren, as Reference StenborgStenborg (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 (Reference CutlerCutler, 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 (Reference Behrens, Bergmann, Moser and JochutnBehrens and others, 1975; Reference CollinsCollins, 1982; Reference Hock and HookeHock 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. Reference KohlerKohler (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.