Introduction
During the 1988 field season, experiments were conducted to elucidate the subglacial drainage system in the lower ablation zone of Storglaciären, a small, temperate valley glacier in northern Sweden, and site of the University of Stockholm’s Tarfala Research Station. The principal objective was to determine the percentage of the drainage system in which conduits are pressurized rather than at atmospheric pressure.
Two approaches were employed. The first involves tracer experiments, the “traditional” method used to explore characteristics of an unknown drainage system. This work builds on the ideas presented by Reference Seaberg,, Seaberg,, Hooke, and Wiberg,Seaberg and others (1988), to the effect that the relationship between Velocity and discharge is distinctly different for pressurized and open-channel systems. The average velocities obtained from the tracer experiments are used with average discharges in a simple model of the drainage system to solve for the percentage of the total system length flowing as an open channel.
The second method involves a comparison of water inputs and outputs of the drainage system. Discharge signals were measured on the supraglacial melt stream emptying into one of the riegel moulins and on the proglacial stream to which the riegel moulins drain. These signals are used in a linear-systems analysis of discharge-wave propagation through the same model of the drainage system.
This paper is organized in two halves; in the first half, I present the motivation for this work, discuss the techniques and results of tracer experiments, develop a simple model for the drainage system and fit the free parameters of the model using the tracer results. In the second half, I extend the model to describe the passage of waves and use linear-systems analysis to fit the model to the measured input- and output-discharge signals.
Background
The exact nature of subglacial drainage systems is still poorly understood. Surface melt- and rainwater enter the system through crevasses or moulins that directly intercept the drainage system. Water leaves the system through streams emerging at or near the glacier terminus. Between entrance and exit the glacial drainage system is, with rare exceptions, inaccessible to direct observation; insight into the system must therefore be gained through a combination of theoretical considerations and indirect observations—remote sensing in the broadest sense.
During the summer melt season, water generated by melting of snow and ice within the glacial watershed, and from rainfall, flows over the glacier surface until it either enters the internal drainage network or flows off of the glacier. The most common input points fur the supraglacial water are moulins, deep vertical holes formed when a crevasse opens to sufficient depth to intersect the englacial drainage system. Water flows into the crevasse, enlarging it by melting the ice walls. When the crevasse eventually closes, the moulin remains and continues to expand, often reaching very large (and dangerous) proportions. Because the ice is in motion, the moulin is transported down-glacier while crevasses continue to open at the same relative position on the glacier surface. Eventually, a new moulin is formed upstream, cutting off drainage to the older moulin.
Moulins are found in areas of the glacier where extensional flow is sufficient to open crevasses and there a an englacial drainage system within reach of the crevasses. On Storglaciären, moulins form principally in one area, where ice flows over a large bedrock riegel, a transverse ridge formed in large part by an erodibility contrast in the underlying rocks Fig. 1. This moulin field is fairly static; Reference Stenborg,Stenborg (1969) reported moulin locations and surface drainage networks in the mid-1960s which are nearly identical to those today, more than later. There is abundant evidence that the drainage network, though, is anything but static. Average velocities, determined from travel times of tracers introduced into the moulins, have varied over the years (Reference Stenborg,Stenborg, 1969; Reference Hooke,, Miller, and Kohler,Hooke and others, 1988; Reference Seaberg,, Seaberg,, Hooke, and Wiberg,Seaberg and others, 1988), even under similar discharge conditions There are several places where subglacial streams exit the glacier; the relative amounts of discharge carried by each stream vary both during the melt season and from year to year. These streams may also shift location slightly, partly in response to retreat of thesnout, but also partly in response to changes in the subglacial drainage system. Finally, at certain times during the melt season, abnormally high levels of suspended sediment aredischarged from the glacier, an indication of drainage-system incursion into previously unexposed areas of the bed.
Fig. 1. (a) Ablation area of Storglaciären, showing bedrock contours (gray lines, contour interval 210 m), surface contours (heavy black lines, contour interval 210 m), moulins (black circles), tracer-injection sites (labels), proglacial streams, and Sydjåkk dam (X). Equilibrium line is at about 1400 m. (b) Cross-section along profile X-X'. The black and white bars show the system lengths of ihe study area flowing as pressurized conduit and open channel, for both the tracer (T) and discharge (Q) analysis. The gray bar denotes the range of reasonable values for the pressurized-conduit length, and the tick mark shows the best fit to the data
Changes in the conduit system occur as the conduit walls are melted by the viscous heat generated in turbulent shearing of the water as it moves through the glacier. Counteracting this melting is the plastic closure of the conduit due to the overburden pressure of the ice. Because water flux through a glacier is variable, conduit geometry is constantly adjusting to balance closure and melting; as a consequence, the most important control on the drainage system is the water flux. Some water is generated by geothermal heating at the bed and viscous-heat generation within the deforming ice mass, but this amount is negligible when compared to the surface-melt flux (Reference Östling, and Hooke,Östling and Hooke, 1986).
Reference Röthlisberger,Röthlisberger (1972) was the first to quantify the opposing effects of tunnel closure and viscous-heat generation in glacial conduits. He idealized englacial conduits as long straight tubes, either cylindrical or semicircular in cross-section, thus allowing him to use Reference Nye,Nye’s (1953) analytical solution for closure of such tubes. To balance closure, Röthlisberger posited that the energy loss in a section of conduit is converted directly into heat, part of which is used to maintain the water at the pressure-melting temperature, and the remainder of which is available to melt the conduit walls. Heat-transfer effects are neglected. A significant result of Röthlisberger’s analysis is that the pressure gradient, and thus the pressure at any arbitrary distance from the terminus, is inversely dependent on discharge. Thus, if two nearly identical tunnels lie side by side, one carrying more discharge than the other, then the resultant potential gradient between the two tunnels tends to drive flow into the larger tunnel. This has led to the idea that glacial drainage networks consist of arborescent systems of conduits that become progres-smaller and more widely distributed away from the glacier margin.
A fundamental question, however, is whether sub-glacial water actually flows in pressurized conduits or in open channels. This question was examined by Hooke who concluded that, throughout much of a glacier, conduits would be at atmospheric pressure. In the specific case of the lower part of the ablation zone of Storglaciären, Hooke found that pressurized conduits would be restricted to small patches of the two over-deepened areas above and below the riegel. Between the riegel moulins and the proglacial stream to which they connect, flow would occur in an open-channel drainage system.
While there is decidedly pressurized flow some distance downstream from the moulins, as evidenced by significant water levels in boreholes and moulins, just how far the region of pressurized flow extends is difficult to establish. During the 1988 field season, a borehole located only 50 m from the main moulin zone registered extremely low water levels, less than 20 m over bedrock, at a time when water inputs were relatively high. Either this might be confirmation of an extensive open-channel system, or it might simply reflect the existence of a pressurized-conduit system with extremely low energy gradients. There are ways of answering this question using other field data than isolated and possibly spurious borehole data; one possibility is to use information obtained from tracer experiments.
Tracer Studies
The most common technique used to evaluate characteristics of unknown flow systems is a tracer experiment; a suitable tracer is “injected” into the system and its return is monitored at the appropriate output locations. Tracer returns yield several types of information. In systems with multiple input or output points, general flow paths can he determined. By calculating the total amount of tracer passing a measurement site and comparing this quannty with the injected amount, it is possible to determine if there were significant tracer losses which might indicate diversion of the traced flow within the system. The return signal can be analyzed to determine the average velocity of the traced flow and the dispersion coefficient, a measure of the spread of the tracer within the system. The return signal can also be compared with theoretically derived returns, or with returns through systems with known geometries (Reference Davis,, Campbell,, Bentley, and Flynn,Davis and others, 1985).
A problem with tracer studies is that the results are particular to the pathway of the tracer. In a spatially homogeneous system, such as a laboratory column of uniformly sized sand, the results of a few tracer studies describe adequately the fundamental paramercrs of the system. Analysis of a spatially heterogeneous glacial drainage system requires many traces, however, with the extent and scale of the inhomogeneities determining the number of traces necessary to describe the system. A further problem concerns temporal heterogeneities, for glacial drainage systems not only are spatially heterogeneous, but are also highly dynamic.
With some exceptions (Reference Burkimsher,Burkimsher, 1983, Reference Seaberg,, Seaberg,, Hooke, and Wiberg,Seaberg and others, 1988; Reference Fountain,Fountain, 1993), most glacial tracing programs are conducted with fewer than 5–10 traces, often with long time periods between each trace (Reference Ambach,, Behrens,, Bergmann, and Moser,Ambach and others, 1972; Reference Ambach, and Jochum,Ambach and Jochum, 1973; Reference Krimmel,, Tangborn, and Meier,Krimmel and others, 1973; Reference Behrens,, Bergmann,, Moser,, Ambach, and Jochum,Behrens and others, 1975; Reference Lang,, Leibundgut, and Festel,Lang and others, 1979; Reference Brugman,Brugman, 1986; Reference Hooke,, Miller, and Kohler,Hooke and others, 1988; Reference Willis,, Sharp, and Richards,Willis and others, 1990; Reference Hooke, and Pohjola,Hooke and Pohjola, 1994). This is not due to slackness on the part of the investigators; in most of these studies it took several days before injected tracer appeared at the monitoring site and even longer periods for the tracer to be completely recovered. Starting another trace before the first has completely arrived can confuse interpretation of both returns.
Long recovery times are not a problem, however, in the lower part of the ablation zone of Storglaciären; here there is sufficiently rapid return in the subglacial drainage system to allow multiple tracer tests in a single day.
Tracing Studies on Storglaciären
Previous work
The earliest tracer work on Storglaciären was undertaken by Reference Stenborg,Stenborg (1969) in the mid-1960s, using salt as a tracer and a band-powered conductivity meter to measure the return. Based on the emergence of tracer at one of the two proglacial streams, Nordjåkk or Sydjåkk (Fig. 1), Stenborg concluded that a bipartite drainage system existed in the lower part of the ablation zone. Tracers injected in the moulins over the riegel and in crevasses on the south side of the glacier consistently came out in Sydjåkk, while those injected in crevasse on the north side ended up in Nordjåkk.
This pattern seems not to have changed in the nearly since Stenborg’s work. The peak arrival times that Stenborg reported are consistent with those found today, despite retreat of the glacier front, changes in the exact positions of moulins and changes in the relative discharges being carried in the two streams, which vary seasonally and from year to year (Reference Östling, and Hooke,Östling and Hooke, 1986).
The connection between the riegel moulins and Sydjäkk was reconfirmed by 15 traces conducted during i 984 and 1985 by Reference Seaberg,, Seaberg,, Hooke, and Wiberg,Seaberg and others (1988), by over 120 traces I performed in 1986, 1987 and 1988, and by 10 traces in 1989 (Reference Hock, and Hooke,Hock and Hooke, 1993). These traces generally yielded a single peak, with travel times of between 0.5 and 2 h.
Location of experiments
During summer 1988, I conducted over 70 tracer tests at sub-daily intervals. Tracer was injected in four large moulins, M1, M2, M3, and M4, all located roughly 50 m up-glacier of the riegel near the glacier center line Fig. 1. Note that previous studies have also referred to moulins in this area as M1, etc. (e.g. Reference Seaberg,, Seaberg,, Hooke, and Wiberg,Seaberg and others, 1988); no consistent labeling scheme has been adopted for this group of moulins, however, even for those active over several years, and thus the names in these various works may have no correspondence to one another.
The geometry of this group of moulins has been studied in detail. Reference Holmlund,Holmlund (1988) made several descents into active moulins, which typically plunge to a depth of before they widen out into a splash pool. He also spent several years mapping successively deeper horizontal cross-sections of older “fossil” moulins as they melted out at the surface. The moulins seem to drain at two levels; at the upper level, water empties through a system of lateral passaġeways that connect along the bottom of the crevasse to other moulins, while at the lower level water drains in a single conduit dipping 45° along the trend of the crevasse. The latter passageways are inferred to connect to other moulins at greater depth than observed.
The principal moulins in my experiments all lie within 100 m of each other along what is essentially the same crevasse. Based on Holmlund’s observations, I assume that the conduits leading from individual moulins along a crevasse system join at relatively shallow depths, and that the local plumbing between the four study moulins and their confluence exerts comparatively less influence on the tracer movement than the drainage system downstream of the confluence. With some exceptions, the tracer tests reported here support this assumption; calculated values of average velocity and dispersion are similar for traces conducted in different moulins at roughly the same time. Accordingly, in this paper I discuss the results of the traces conducted in the riegel moulins as if they were performed in a single moulin.
Tracer choice
The ideal tracer is conservative, i.e. it follows the exact path of the water molecules as they pass through the system and does not decay or chemically interact with any part of the system. Reviews of the tracers currently in general use can be found in Reference Davis,, Campbell,, Bentley, and Flynn,Davis and others (1985) and Reference Aley, and Fletcher,Aley and Fletcher (1976).
Two tracers are commonly used in experiments at Storglaciären: chloride ions, in the form of common table salt, and rhodamine WT, a fluorescent dye. Salt is inexpensive and is easily measured in the field with a conductivity meter. However, large quantities must he used in order to generate return signals that are detectable over background diurnal-conductivity variations without resorting to more elaborate laboratory ion-measurement procedures. For this reason, salt can he used only where the connection between input and output is well established and where flow through the system is sufficiently rapid. Rhodamine WΤ is more expensive, and data collection and analysis are more labor-intensive. However, concentrations as low as 0.1 ppb (parts per billion) can be detected, and it is an extremely conservative tracer in most environments (Reference Davis,, Campbell,, Bentley, and Flynn,Davis and others, 1985). This makes it useful in longer-term traces, where the pathways are not established and where dilutions can be up to eight orders of magnitude.
All of the traces discussed in this paper were conducted with salt, both because rhodamine traces were being conducted elsewhere on Storglaciären (Reference Hooke, and Pohjola,Hooke and Pohjola, 1994) and because the large number of traces could be monitored easily using a conductivity meter and a data logger.
Tracer injection
Chloride ions were injected by mixing 2–4 kg of ordinary table salt in a 101 jug of cold water (1–2°C) taken directly from the stream draining into the moulin. Complete dissolution was impossible at such high concentrations, so the best way to inject the mixture into the moulin was to agitate the jug for about 1 min to put the undissolved salt into suspension, and then pour as rapidly as possible. In all cases salt remained on the bottom of the jug, requiring a refill with water and an additional pour. The total time for the entire operation was typically 20–30s. A few salt traces were done using hot water (40°C) obtained from the drilling team when they were in the vicinity. With hot water, more salt can go into solution. Although there was usually some residue after the first pour, this amount was and a refill and second pour were not required. This reduced the total pour time to about 15–20s.
With either method, it was virtually impossible to put all of the salt into solution. Consequently, there was some concern that some undissolved salt might lag the dissolved salt, or that the dense saline solution might lag more conservative tracers, such as rhodamine, because of density effects. However, results from an experiment in which salt dissolved in cold water was injected simultaneously with two fluorescing dyes, rhodamine and lissamine, yielded comparable returns (Reference Kohler,Kohler, 1992). This suggests that density effects are insignificant.
Sampling
Conductivity was measured at a dam located on Sydjåkk, 300 m downstream from the glacier snout Fig. 1. This sampling site was farther downstream than desirable, but because the subglacial channels debouch at many points along the glacier front, data had to be collected downstream from the point where the channels meet.
A conductivity electrode was installed behind the Sydjåkk dam and the output was monitored with a Campbell Scientific 21X data logger. Two data sets were obtained from the data logger, one consisting of measurements taken at 3 min intervals, the other of measurements made whenever the signal changed by more than a certain amount. The latter record begins later in the field season, but the two records were combined as far as possible to give a more accurate registration of the signal on both the time and measurand axes. Figure 2 shows the 3 min records of stream conductivity and discharge for a typical trace.
Fig. 2. Discharge Q, conductivity C, and time-lo-peak tp for a typical trace
Data Analysis
Discharge versus velocity
The principal quantities of interest are the average discharge Q, which is calculated for the period between the injection and the arrival of the peak, and the minimum velocity U, which is estimated by dividing L, the straight-line distance between the pour sites and the Sydjåkk dam, by t p, the time elapsed between the injection and the tracer peak. With the combined conductivity record described above, precise identification of the arrival time of the peak is possible. In the 3 min record the peak is less accurately identified because it might lie anywhere within a 3 min period; however, the resultant error in estimating U would be no more than even 50%, in the worst case.
The relationship between average flow discharge and velocity has been used in the past to infer whether the subglacial drainage system occurs predominantly as a pressurized conduit or as an open channel (Reference Seaberg,, Seaberg,, Hooke, and Wiberg,Seaberg and others, 1988). Continuity requires that
where A is the average cross-sectional area of the drainage system. If the drainage system consists entirely of a pressurized conduit rather than an open channel, and if the system does not change significantly through time, then A is approximately constant and a linear relationship exists between Q and U.
This is decidedly not the case for flow in an open channel because both A and U vary with discharge Q, This is quantified by considering the Darcy–Weisbach equation for velocity in an open channel:
where the hydraulic radius of the flow R h = Α/Ρ, Ρ is the wetted perimeter, S is the bed slope, g is the acceleration due to gravity, and f 0 is an empirical friction factor. The friction factor for fully rough turbulent flow is:
where z 0 is the height of an “effective”-roughness element (Reference Henderson,Henderson, 1966, p.98). While this expression results from studies of flow in pipes, it is entirely appropriate, both from scaling arguments and from empirical verification, to use it in open-channel flows. Combining Equations (2) and (3) yields
Channel geometry influences the relationship between U and Q in an open channel. Consider two idealized channel types: V-shaped, representing a proglacial stream, and rectangular, representing a wide, shallow subglacial conduit with constraining side walls. It is straightforward to derive a relationship between A and R h for either channel type, which is used together with Equation (4)and Equation (1) to yield equations of the form
where the values of a and b depend on the cross-sectional geometry of the channel. Table 1 gives expressions for A, P, Rh, A = ƒ(R h). and U = f(Q) for the two open-channel types and for a fully pressurized semicircular conduit.
For the open-channel geometries considered, b ≤ 0.4, whereas in the case of the pressurized conduit, b = 1. It is this difference that provides the basis for distinguishing open channels from closed conduits.
Results
A piot of Q versus U for all of the riegel-moulin traces shows that most of the data fie aiong a finear trend (Fig. 3). While it appears that results of traces in individuai moulins group afong different trends, this is due to changes in the drainage system over time; at the beginning of the field season most traces were conducted in M4, while later traces were primarily in M1.
Fig. 3. Average discharge Q and velocity U for tracts in riegel moulins M1–M4. and in proglacial stream S2 (symbols). Dashed line is fitted to S2 data using Equation (5) for a V-notck channel
If the drainage system were flowing entirely as an open channel, one would expect discernible curvature in the Q versus U relationship. However, the observed distribution is, apart from the outlying points, quite linear. This implies that A is constant, and that the drainage system was flowing predominantly as a pressurized conduit. To highlight the difference in the relationship between discharge and velocity for the two different systems, Figure 3 also shows Q versus U for traces conducted in Sydjåkk between the glacier snout and the dam, and a best-fit line using Equation (5) for a V-shaped channel, with b = 0.25 and a adjusted to fit the data.
Table 1. Hydraulic parameter values for two open-channel models and closed conduit
Fig. 4. Upper graph: ratio of average discharge Q to average velocity U for the riegel moulins M1 (black squares), M2 (black triangles), M3 (crosses), and M4 (white triangles). Lower graph: smoothed discharge (dark gray line), hourly discharge (lighter gray line) and average discharge at the time of the trace for the riegel moulins (symbols). Dotted vertical lines demarcate times of tracer (T) and discharge (Q) analysis.
A conduit with time-varying geometry would appear in a Q versus U diagram as a curve whose exact position and shape would be a function of the instantaneous state of the system. Traces conducted quickly enough and over a sufficient range of discharges would provide a snapshot of the curve. If the state of the system changed more rapidly than the rate at which traces were conducted, then the snapshots would be able to record only a single point on the time-varying curve, and the results, when taken together, would show what appeared to be a non-systematic relationship of Q versus U. This would be particularly striking at higher discharges or at the beginning of the melt season when changes in the subglacial drainage system occur most rapidly.
Looking at system changes in a slightly different way, Figure 4 depicts Q/U = A as a function of time, together with the hourly mean discharge, the discharge smoothed with a 24 h moving-average filter, and the mean discharge Q during the time of the trace. Discharge varies on two time-scales; there is a diurnal variation, shown by the hourly discharge record, while a longer time-scale variation is illustrated by the smoothed discharge. The longer time-scale changes in discharge are echoed by changes in the cross-sectional area A; this suggests that system vanations occur predominantly at low frequencies, an explanation that is consistent with slow melting and closure in the subglacial conduits (e.g. Reference Röthlisberger,Röthlisberger, 1972). There may be sub-diurnal vanations as well, but the traces were performed too infrequently to support detailed analysis.
Figure 4 shows how the Q/U relationship changes over time, and highlights the differences between results from individual moulins. While I asserted earlier that the riegel moulins may be assumed to behave as a single effective moulin, this is clearly not the case for all traces; M1 and M4 vary inconsistently on 28 July, for example. If the conduits do join a short distance below the bottom of the moulins, this would imply that the early career of the travelling tracer cloud is important in determining what is ultimately detected at the snout. This can also be surmised from the apparent linearity of the Q versus U relationship. Because the tracer traverses an indisputable 300 m of open channel, or more than 25% of the straight-line system length, apparent Linearity implies that the tracer must have spent a disproportionate amount of time in the moulin-and-conduit part of the system. The shorter the inferred conduit section, the longer is the relative time the tracer must have resided there.
Modeling the Data
Model description
A simple model can be used to quantify the percentage of conduit flow necessary to produce apparent linearity in the Q versus U relationship. The drainage system is modeled as a single vertical moulin extending directly to the bed and connecting from there to a horizontal semicircular pressurized conduit draining into an inclined open channel (Fig. 5). Flow throughout the system is fully turbulent, the moulin and conduit have constant cross-sections, and the open-channel segment is a wide rectangular channel, with no roof constrictions. Conduit closure and melting are neglected; the system is assumed to be static.
Fig. 5. Simple model diagram. See text for explanation of symbols
The relevant parameters for the model are: the height of water m the moulin h, relative to the junction between the moulin and the conduit; the cross-sectional areas of the conduit and the moulin A c, and A m; the average velocities in the conduit and open-channel segments U c and U oc; the minimum velocity through the entire system U; the average discharge during the time of the trace Q, the straight-line horizontal system length L; and the system sinuosity v, defined as the ratio of the actual system length to L. If p is the fraction of L occupied by the conduit, then the length of the system occupied by an open channel l oc is given by
Model discharge versus velocity
I seek an expression that relates U to Q as a function of p and other geometric properties of the system. The average time, or time-to-peak t p, that it would take for a tracer to pass through the system can be partitioned among the three parts of the system, so that
where the subscripts m, c, and oc refer to the moulin, conduit and open channel, respectively. This assumes that tracer is transported instantly from the top of the moulin to the upper surface of the water backed up in the moulin. Dividing both sides of Equation (7) by L and inverting gives the minimum velocity
The individual travel times for each component of the system can also he written in terms of velocity to yield
From continuity,
and, from Table 1, for the case of a wide rectangular open channel,
where
S is the slope of the channel, W is the stream width, and z oc is the height of the effective-roughness element. These last three expressions are substituted into Equation (9) and rearranged to yield
The relationship between h and Q in the horizontal conduit segment is given by Bernoulli’s equation for pipe flow and the Darcy–Weisbach expression for frictional head loss (Reference Vennard, and Street,Vennard and Street, 1982, p.358):
where R c is the hydraulic radius and fc is the Darcy–Weisbach friction factor of the conduit. While Bernoulli’s equation is strictly valid only for streamline flow, the correction factor applied for turbulent flow is roughly unity, and is typically neglected (Reference Vennard, and Street,Vennard and Street, 1982, p.357). Substituting Equations (10) and (11) into Equation (15) results in
Because the first term within the brackets on the righthand side of Equation (16) dominates for all but the shortest of conduits, I drop the last two terms. For a semicircular closed conduit, R c = 0.244 √A c, while f c is given by Equation (3), with the appropriate change in subscripts. Combining these two expressions yields
where z c is the height of the effective-roughness elements in the conduit. Substituting Equation (17) into the truncated form of Equation (16) yields
which, when substituted into Equation (18), yields
where
The first, second and third terms in the denominator of Equation (19) are the contributions of the moulin, conduit and open-channel sections, respectively.
Taking each term separately, one sees that velocity in the conduit-flow part of the system is proportional to discharge and that velocity in the open-channel part is proportional to the two-fifths power of discharge, as noted previously. In the moulin term, U is inversely proportional to Q; this is because higher discharges require greater head to drive the water through the conduit section, resulting in increased distance that the tracer must travel. There is of course a physical limit; if the water in the moulin reaches the glacier surface, then this term becomes constant. For the proposes at hand, this is not a problem; water did not fill the moulins to the top during the period of these experiments.
Fitting data to the model
If either the open-channel or moulin terms in Equation (19) are of the same order as the conduit term, then there will be appreciable curvature in a plot of U and Q. This is clearly at odds with the observed Linear trend of the tracer results. The question is, what combination of parameters does explain the data, and what are the implications for the drainage system?
This can be answered by Siting Equation (19) to the tracer data to estimate possible values of K m, K c and K oc. In doing so, I use only those data taken after 2 August (Fig. 4), since during this period, with one exception, there are no significant excursions from the trend in the U versus Q relationship. Discharge during this period occurred as a steady diurnal signal, with no sudden increases due to rainfall or up-glacier releases of stored water. I assume that the conduit geometry had adjusted to an approximately steady-state configuration after a large rainfall that had occurred some days earlier. The one outlying tracer result, conducted in M4 on 5 August, was probably the result of a misrecorded pour time rather than a short-term change in the conduit system, and was removed from the reduced data set.
Fitting is performed by inverting Equation (19) to obtain
When written this way, any general linear least-squares fitting mutine might be used to solve for the coefficients K m, K c and K oc. However, this yields one negative and two positive values for K m, K c, and K oc, respectively. Negative coefficients are not physically plausible, and an alternative fitting scheme is employed. Fitting involves the application of a merit function to quantify the agreement between the data generated by a particular combination of model parameters and the data actually measured. One such merit function is
where yi are the observed data, and yf are the modeled data (Reference Press,, Flannery,, Teukolsky, and Vetterling,Press and others, 1989, p.551). In a general linear least-squares fit, the derivative of χ 2 with respect to the model parameters is set to zero and the resulting system of equations is solved exactly using linear algebra. The problem is that there may be other combinations of model parameters that give an equally reasonable fit. Various canonical methods exist for characterization of χ 2 in parameter space (Reference Press,, Flannery,, Teukolsky, and Vetterling,Press and others, 1989, p.551); as a “quick-and-dirty” alternative, I simply use Equations (23) and (24) to calculate χ 2 at discrete values of the coefficients K m, K c and Κ oc. This both reveals the range of satisfactory parameter values and permits the imposition of constraints on the parameter values.
As stated, negative values are physically implausible. A further constraint may be applied as well, namely that the proglacial stream has a measurable length with respect to the unknown total system length, so that parameter values indicating a pressurized section extending beyond the confines of the glacier are also physically implausible. Invoking this constraint is essentially equivalent to adjusting the velocity data to account for the time the tracer spends in the proglacial stream, as described by Reference Seaberg,, Seaberg,, Hooke, and Wiberg,Seaberg and others (1988). I rearrange Equation (6) and insert it into Equation (22) with L = 1300 and l oc = 300 to yield the condition that
Rather than estimating values for the parameters W, z 0 and S in Equation (12) to obtain a value for K oc I fit k oc to the U and Q values obtained from the traces conducred in the proglacial stream. While stream geometry may not be quite the same beneath the glacier, and while recognizing that the proglacial trace data comprise only four points, which do not occur over a very wide range of discharge, the value of k oc ohtained is unlikely to he significantly in error. This is hccausc W, S and z 0 contribute to k oc in fractional powers; for a rectangular channel, a two-fold error in estimating these parameters contributes to an error in k oc of only 25%, 23% or 7%, respectively. In estimating k oc from the proglacial-stream data, it does not matter whether the V-shaped or the rectangular channel model is used. In either case, k oc ≈ 0.7, which, when inserted into Equation (25), results in the Limitation that
Results
The range of constrained best-fit valuer lying within 68.3% confidence limits (one standard deviation) are K m = 0.0−0.2, K c = 1.4−1.7, and K oc = 0.33−0.7, with the minimum value of χ2 in the vicinity of K m = 0.05, K c = 1.6 and Κ oc = 0.4. Figure 6 shows the curve resulting from Equation (19) with these coefficient values, plotted together with the field data and the 68.3% confidence limits.
Fig. 6. Average discharge Q and velocity U for traces conducted 2–10 August (open circles), restricted best fit of data to Equation (24) (solid line), and 68.3% confidence limits (light gray lines)
To solve for open-channel length l oc Equations (22) and (6) are combined to yield:
Using the range of fitted values of Κ oc and k oc = 0.7 implies an open-channel segment of 300–640 m, with the best-fit value at 360 m. Subtracting the proglacial-stream segment of 300 m yields pressurized-conduit flow occurring to within 340–0 m of the snout, with the best fit at 60 m upstream of the snout (Fig. 1).
While there is no independent way of verifying this result directly, it is possible to check the validity of the model by comparing fitted values of K m, which reflect water levels in the moulin, with water-pressure data from the same area as the riegel moulins. This can be seen by recasting the expression for K m in terms of h and Q (cf. Equations (14) and (19), so that
Although A m was not measured during these experiments, observations from previous years can he used to make a crude estimate. In his descent to the bottom of one of the riegel moulins, Reference Holmlund,Holmlund (1988) mapped a cross-sectional area of roughly 5 m2 in the plunge pool, at a depth of roughly 40 m. He also reported that the conduits leading away from the splash pool of fossil moulins have cross-sectional areas of about 0.1 m2. Computing an “average” value for A m is difficult in this case, but to maximize the estimate of h, I will assume that A m is better represented by the small conduits draining the bottoms of the moulins. If a single moulin receiving roughly 5% of the observed downstream discharge has a cross-sectional area of 0.1 m2, then total moulin area A m will be about 2 m2. If K m = 0.05, and L = 1300 m, then Q in the observed rangr of 0.2–1.0 m3 s−1 leads to h varying from 5 to 100 m.
The low values determined for K m connote low water levels in the moulin, and thus in water pressures measured anywhere along the conduit. The water-level record in a borehole roughly 50 m down-glacier from M4 supports this conclusion; during 2–10 August water-level was at or below the level of the transducer, roughly 20 m over the bed.
The results of the tracer analysis suggest that a significant percentage of the drainage system flows within a pressurized conduit. I will now turn to the other experimental method used to determine this percentage.
Discharge Waves and Linear-Systems Analysis
Linear-systems analysis, better known to hydrologists as the “unit-hydrograph” method, has been used to analyze hydrologie processes in a variety of media (Reference Himmelbau, and Bischoff,Himmelbau and Bischoff, 1968; Reference Dooge,Dooge, 1973; Reference Neuman, and de Marsily,Neuman and de Marsily, 1976; Reference Dreis,Dreis, 1982). The drainage system is simply viewed as a “black box” with an input signal and output signal, in this case, discharge. For an unknown system, the two signals can he analyzed to determine empirically the characteristics of this black box, and the results can he compared with theory.
In linear systems, some process converts an input signal f i(t) into an output signal f o(t). The elemental-input signal is a single spike, and the system’s response to this elemental-input signal, i.e. the system’s output signal, is the transfer function (Reference Papoulis,Papoulis, 1962). The central tenet of linear-systems analysis is that it is possible to combine the output signals of any number of arbitrarily scaled elemental-input functions and sum them to obtain the total system output; this is what is expressed in the convolution integral. More relevantly, one can also deconvolve two signals, that is, solve for an unknown transfer function, using known input and output signals. Convolution and deconvolution are readily computed in the frequency domain using the Fast Fourier Transform (Reference Press,, Flannery,, Teukolsky, and Vetterling,Press and others, 1989).
The two types of drainage system, open-channel and pressurized-conduit, have characteristic ways of transforming an input discharge signal into an output discharge signal For example, if a sudden increase in discharge is introduced into a steady-state system, there is a characteristic delay in the arrival of this step increase at the downstream end of the system. This delay represents the time needed for the system to adjust its energy gradient to accommodate the new water flux. Because the two systems have different means of transmitting energy down-gradient, the delay depends primarily on the type of flow in the system. In the case of a pressurized conduit, this delay represents the time it takes for the moulin to fill up to the new level required to drive the additional discharge through the conduit. Because water is largely incompressible, transmission of the pressure wave is nearly instantaneous through the conduit, provided that the walls are relatively impermeable and rigid, and that further storage does not take place along the flow path. In an open-channel system, however, the storage takes place within the flow itself, and the storage “location” is advected with the flow. The characteristic speed at which the change in discharge travels is a function of the average velocity of the flow, which is itself a function of discharge and channel geometry.
To quantify the way in which discharge waves move through a subglacial drainage system consisting of a combination of the two possible end-member systems, I will use the same simple model described above (Fig. 5). input discharge Q i enters the moulin and backs up to the height necessary to drive the output flow through the length of the pressurized conduit. The flow then continues through the open-channel segment until it arrives as output discharge Q o at the glacier portal.
This development assumes that: (1) the effects of unsteady flow are negligible in comparison to the frictional losses along the conduit length, (2) frictional losses in the conduit are much larger than in the moulin, (3) all flow parameters can be cross-sectionally averaged, and (4) during the time of interest the conduit can be treated as a rigid pipe.
Discharge waves in pressurized conduits
Consider first the passage of discharge waves in a system flowing entirely as a pressurized conduit. Because mass is conserved and the conduit walls are assumed to be impermeable, the relationship between input and output discharge in the system is given by
The height of water in the moulin h required to drive the flow Q o through the conduit is given by the truncated form of Equation (16), with the change of notation that Q = Q o:
Equating Equation (29) with the time-derivative of Equation (30) and rearranging yield
where
The term β embodies all the geometrical parameters of the system, and is taken to be a constant for the simple system envisioned, fn terms of K m, the parameter describing the relative amount of time that tracers spend in the moulin section of the system, β is given by
Equation (31) is non-linear, with no known analytical solution in closed form (Reference Kamke,Kamke, 1942). One can explore the behavior of such non-linear equations using standard perturbation techniques (Reference Nayfeh,Nayfeh, 1973). The input discharge is written as a small perturbation about a reference discharge Ǭ:
where ε is a small constant that scales the perturbation and G is any function of time, scaled to unity. Similarly, Q o is written as a power-series expansion about Ǭ:
These two expressions are substituted into Equation (34) and terms are collected in ε. For sufficiently small perturbations, terms higher than 0(ε) are neglected, resulting in:
Equation (36) is essentially the “linear-reservoir” model Reference Nilsson, and Sundblad,Nilsson and Sundblad, 1975). Numerical integration of Equation (31) for various Q i gives qualitatively similar results to the linearized Equation (36), even at fairly large values of ε.
The initial condition used in solving Equation (36) is that as t → ∞,
where
and u(t − τ) is the unit-step (Heaviside) function. The integral in Equation (37) is a convolution of the input-discharge signal with an exponentially decaying function, H, known variously as the system-response function, kernel function, or transfer function, and hereafter called the transfer function (Reference Papoulis,Papoulis, 1962). This particular transfer function describes how a pressurized-conduit system adjusts to changes in the input signal.
As an example, consider a step increase in the input discharge by an amount Q a over the background level at time t = 0,
in this case, the explicit solution of Equation (37) is
(Fig. 7). The characteristic time tc at which the conduit system fully adjusts to the increase in discharge is controlled by the parameters β and Ǭ. If the criterion for full adjustment is taken as the time at which output discharge is 95% of input discharge, then this occurs at the time tc = 3βǬ.
Fig. 7. Input and output discharge Qi and Qo in a pressurized-conduit system for a step increase in Qi, and time tc at which Qo is 95% of Qi.
Recalling the definition of β one can see that a long, rough conduit with a large ratio between the moulin and conduit cross-sectional areas would take comparatively longer to adjust to a sudden increase in discharge than would a short, smooth conduit with a small ratio between the moulin and conduit cross-sectional areas.
Discharge waves in open channels
Modeling the passage of discharge waves in open is more complicated, and generally requires the use of numerical techniques. For small increases in discharge, however, an approximate solution exists. At the upstream end of a steady flow m a uniform open channel with velocity U oc, a step increase in discharge is introduced; this disturbance moves downstream as a monoclinal kinematic wave travelling at a constant velocity c. If the wave travels downstream with no frictional losses or change in shape, and the channel is uniform, then
(Reference Henderson,Henderson, 1966).
If a step increase in discharge is introduced at the top of the moulin, and is transmitted directly to an open channel conduit connected to the bottom of the moulin, then the wave speed in the open channel is given by:
where t oc is the characteristic time delay for the arrival of the wave at a distance X downstream. Equations (41) and (42) together give
The relationship between velocity and discharge for a wide rectangular channel
is combined with Equation (43) to yield
Because the forward motion of the wave is assumed to be steady and unaffected by bed geometry or by diffusion of the wave crest, I set the output discharge in the open-channel system as simply being identical to the input discharge offset by the characteristic time of Equation (45). Thus Q i and Q o in the open channel are related by
The dependence of t oc on Q in Equation (45) means that Equation (46) is non-linear. Since t oc is proportional to only a fractional power of Q, however, this non-linearity is relatively insignificant and I linearize Equation (46) by rewriting the expression for t oc using mean discharge Ǭ:
Discharge waves through a combined system
To describe passage of discharge waves through a system comprising both pressurized conduit and open channel sections, I combine Equations (37), (46) and (47), making the changes of notation that Q o for the pressurized conduit is Q i for the open channel, and that X in Equation (47) is (1 − p)vL. This yields
where
and H is given by Equation (38). The transfer function for the combined system is a decaying exponential delayed by the amount of time necessary to traverse the open-channel segment of the system (Fig. 8).
Fig. 8. Input and output discharge Qi and Qo in a combined presurized-conduit–open-channel flow system for a spike input, and the time toc it takes the pressure wave to traverse the open-channel segment of the system.
Storgiären as a Linear System
In the following sections, I describe the time series of input and output discharges measured on Storglaciären, deconvolve the time series to obtain an experimentally derived transfer function for the drainage system, and finally, fit the theoretical transfer function (Equation (46))to the calculated transfer function. There are three major assumptions in this development:
-
1) The single measured input-discharge signal is roughly equivalent to any input signal entering any moulin in the drainage system, with adjustments only in the amplitude of the signal reflecting the size of the particular drainage area. Although supraglacial streams are are constantly reorganizing themselves as the glacier surface melts downward, thus decreasing or increasing discharge to an individual moulin, this assumption is probably not too unreasonable; I monitored discharges at two neighboring moulins simultaneously for a short period, and found that the two signals were roughly synchronous, with gradual changes taking place over longer time-scales than the response time of the system.
-
2) The simple model a representative of the subglacial drainage system; a single input-discharge signal enters one group of moulins located at the upstream end of the system, with no other sources or sinks along the length of the system. Clearly this is not so; there are contributions from several smaller moulins downstream of the riegel, as well as from supraglacial runoff at the snout. However, as previously mentioned, the riegel moulins drain an area that is roughly 70% of the total glacier surface between the Sydjåkk dam and the upper ablation zone, while the lower moulins and supraglacial runoff areas are roughly 20% and 10%, respectively. In addition, it is reasonable to suppose that the lower moulins do not differ radically in their characteristics from the riegel moulins; although they are Farther downstream in the system, they are also smaller and thus might have similar transfer functions. Although the same argument cannot be made for the supraglacial runoff, the influence of a signal comprising roughly 10% of the total discharge, whatever the shape of its transfer function, should be minimal.
-
3) The geometry of the drainage system does not undergo significant change during the period of the experiment. Thus the parameters β and t oc are treated as constants and the entire time series is in one step.
Discharge Time Series
The discharge data consist of two records spanning the period 5–9 August. The output discharge was measured at the Sydjåkk dam. Water level over the dam was monitored through time, and periodic measurements of discharge, using the constant-injection fluorescent-dye method (Turner Designs, 1981), allowed the construction of a rating curve.
The input discharge was measured at M1 using a “moulin bag”. A large fabric funnel connected to a section of pipe with a small hole at the bottom was suspended within the mouth of a moulin; water flowing over the lip was channeled into the funnel and backed up in it to the level required to drive a particular discharge through the hole at the bottom of the pipe (Fig. 9). A pressure transducer in the pipe connected to a data logger was used to monitor the water level.
Fig. 9. Moulin bag: a fabric funnel is connected to a section of pipe with a hole at ike bottom and is suspended within themouth of a moulin. A pressure transducer in the pipe monitors the water level within the bag.
Two moulin bags were deployed during the 1988 field season. The first was unsuccessful because the length of pipe at the bottom of the funnel was insufficient to measure the full range of discharges; either the water filled all the the way to the top or it all drained out of the bottom. The second bag was modified by adding a 2 m length of pipe with two holes located at different heights from the bottom; low discharges would be measured by the lower hole, while higher discharges would drain through both holes. This arrangement proved satisfactory for measuring the range of water heights in the bag, although it made interpretation of the rating curve more difficult. Seventeen discharge measurements were made in the supraglacial stream to define the rating curve for the MI moulin bag. The discharge passing through the bag is, in theory, proportional to the square root of the height of water driving the flow. The constant of proportionality is a function of such geometrical factors as the cross-sectional area of the pipe and of the drain hole. Because there were two holes in the M1 moulin bag at different heights within the pipe, two curves were fitted to the data, each with the requisite half-power relationship, and a smooth interpolating curve was drawn at the junction where the curves meet (Fig. 10).
Fig. 10. Rating curve for the moulin bag in M1. Two curves art used because the moulin bag had two holes at different heights from the bottom
Analysis
The two time series consist of 2048 points of discretely sampled discharge data taken at 3 min intervals during 5–9 August. The input data are first normalized so that the total volume under the curve corresponds to that of the output discharge (Fig. 11a), in accordance with the assumption that the single measured input signal is representative of all signals entering the system. Low-pass filtering the time series reveals differencesin the low-frequency components of the two signals (Fig. 11b) which are most likely the result of gradual changes in the drainage area upstream from the moulin. Because the response time of the system is several hours, and because the changes occur over periods of days, I subtract these low-frequency components from both signals (Fig. 11c).
Next, the data are prepared for application of the fast Fourier transform (FFT). The first 20 and last 20 points are multiplied by a sinusoidal smoothing function to remove any sudden jumps in the time series at either end. The time series are then “padded” with 1024 zeros to eliminate “wraparound”, an effect whereby the operations performed at one end of a time series influence the other end, and a by-product of an assumption implicit in the FFT, namely that the time series is periodic. Finally, deconvolution is performed in the frequency domain to yield the transfer function (Reference Press,, Flannery,, Teukolsky, and Vetterling,Press and others, 1989).
Fig. 11. Input-discharge (gray line) and output-discharge (black line) signals after normalization of the input signal. Time is relative to midnight, 5 August 1988. (b) Low-frequency component of the input-discharge (gray line) and output-discharge (black line) signals, (c) Input-discharge (gray line) and output-discharge (black line) signals with the low-frequeney components subtracted.
The resultant transfer function has a large component of high-frequency noise that obscures the underlying structure (Fig. 12), which nevertheless does resemble the decaying exponential of the theoretical transfer function Equation (48). High-frequency noise is an inherent problem in deconvolutiun in the frequency domain, attributable mainly to random noise in the original time series (Reference Press,, Flannery,, Teukolsky, and Vetterling,Press and others, 1989). In dividing the output by the input in the frequency domain, equal weighting is given to all frequencies contained in the two signals. Although the transformed output and input signals may have low amplitudes at higher frequencies, when one is divided by the other a relatively high amplitude value results.
Noise is typically handled by filtering out the high-frequency component either before or after the transfer function is calculated. For the purposes at hand, however, this a unnecessary; instead, Equation (48) is fitted directly to the first 120 points of the calculated transfer function using the Levenberg-Marquardt algorithm (Reference Press,, Flannery,, Teukolsky, and Vetterling,Press and others, 1989 p.572ff.).
The range of acceptable fits defined by one standard deviation is βǬ = 3000–8000 and t oc = 6–15 min, with the best fit to the data occuring at roughly βǬ = 4500 and t oc = 9 min. Combining Equations (6) and (49) yields the expression for open-channel length:
The fitted values of t oc together with k oc = 0.7 and Ǭ = 0.5 indicate a total open-channel length of 285–715 m, with the best fit to the data at 430 m, or, upon subtracting the proglacial-stream length, a subglacial open-channel length of between −15 and 415 m from the snout, with the best fit occurring at 130m (Fig. 1). These values are consistent with the range estimated with the tracer analysis, negative values notwithstanding.
Fig. 12. Calculated transfer function, with close-up of the first 120 points (gray line) and the theoretical transfer function using best-fit values of the free parameters in Equation (48) (dashed line)
The range of values for βǬ is rather large, however. The best fit is equivalent to β = 9000 (using Ǭ = 0.5), which when inserted into Equation (33) yields K m = 3.5. This is well outside the highest value for K m determined from the tracer results and implies a physically implausible range of water height in the moulin. Inserting K m = 3.5 into Equation (28) with A m = 2 m leads to h varying from 450 to 2250 m over the observed discharge range. One possibility is that the estimate for A m is too low. The pressure record from the borehole in the vicinity of the riegel moulins indicates that h was roughly 50 m at a time when Q was 1.0 m3s−1; inserting these values into Equation (28) yields A m ≈ 35 m2. This value is closer to the cross-sectional dimensions of moulins above the splash pool reported by Reference Holmlund,Holmlund (1988). The problem remains that a cross-sectional area of such magnitude requires tracers to spend a significant amount of time in that segment. This would lead to a strong inverse relationship between U and Q, which is not supported by the tracer data.
The difficulty is that the calculated transfer function has too long a “memory”. A plausible explanation is that the calculated transfer function is polluted by “under-flow”, i.e. seepage of water from the drainage system farther up-glacier of the riegel moulins. When high discharge enters the riegel moulins, it backs up in the moulin to form a hydraulic “dam” which prevents the up-glacier water from draining into the conduit system. When discharge entering the riegel moulins decreases, the water level falls, the hydraulic dam is lowered and the up-glacier waters are free to drain again. This mechanism would lead to an extended tail in the transfer function, particularly if the drainage system up-glacier from the riegel is distributed in large numbers of small channels as is proposed in several studies (e.g. Reference Hooke,, Miller, and Kohler,Hooke and others, 1988; Reference Hooke, and Pohjola,Hooke and Pohjola, 1994).
For the purposes of this paper, however, an extended tail in the transfer function does not pose a serious problem. The simple model captures the essential shape of the transfer function of more complex multiple-input systems, and the important parameter for estimating the open channel length remains the time-to-peak in the transfer function.
Discussion and Conclusion
Using two different methods, I conclude that an appreciable length of the subglacial conduit system on Storglaciären occurs as a pressurized conduit. This is quite different from the completely open-channel system Reference Hooke,Hooke (1984) predicted using Reference Röthlisberger,Röthlisberger’s (1972) unmodified expression for a semicircular basal conduit.
This is by no means a new result, for it has long been recognized that the Röthlisberger model underestimates water pressures. Röthlisberger found that to model observed water pressures on Gornergletscher, he needed to use an abnormally low value for the ice-viscosity parameter in the flow law for ice (Reference Röthlisberger,Röthlisberger, 1972). Studies on other glaciers have yielded similar results, namely that basal ice needs to be in the range of 4–150 times softer than is normally assumed in order to successfulfy predict observed borehole water pressures with semicircular conduits (Reference Hooke,, Lauman, and Kohler,Hooke and others, 1990).
One possibility is that closure really is larger than assumed; closure rates in ice tunnels located beneath Bondhusbreen and Engabreen, Norway (160 and 200 m ice thickness, respectively) are observed to be 4–16 times larger than those predicted using Nye’s circular-tunnel closure model with customary values of the viscosity parameter in the flow law (Reference Hagen,, Liestøl,, Sollid,, Wold, and Østrem,Hagen and others, 1993).
Reference Hooke,, Lauman, and Kohler,Hooke and others (1990) summarize a number of proposed mechanisms for softer basal ice, such as higher water content, impurities, and tertiary creep effects, but conclude that these cannot adequately explain the low predicted water pressures. They propose instead a simple geometrical modification to the Röthlisberger expression that would result in more rapid closure rates. By assuming that subglacial conduits are broad and low rather than semicircular, they were able to match observed pressures.
The results of this work do not shed further light on which, if any, of these mechanisms may be responsible for enhanced closure, but do support previous findings using two new and different approaches.
Acknowledgements
This research was funded by grants from the Geological Society of America, Sigma Xi, the University of Minnesota, and the U.S. National Science Foundation. Many people at Tarfala helped with the fieldwork, more than I can thank in this space, but special recognition is extended to V. Pohjola and P. Jansson for carrying their daily salt quota and for the stimulating field conferences. Comments from R. LeB. Hooke, C. Paola, E. Isaksson, and two anonymous referees helped to improve the manuscript.
