Introduction
Tracer techniques are powerful tools for investigating virtually inaccessible hydrological systems and processes. In such investigations, the transit velocity of a tracer through the drainage system and the shape of the tracer breakthrough curve (TBC) are the most commonly evaluated transport parameters. These parameters are governed by both water discharge through and the structure of the drainage system. A characteristic feature of discharge through a glacier is its pronounced variability on seasonal and diurnal time-scales. In addition, a glacial drainage system has the ability to adjust its geometry to prevailing hydraulic conditions. The structure of glacial drainage systems and their evolution with time can be characterized by performing tracer tests over a range of different discharges to investigate the variations of transit velocity and tracer dispersion with discharge (Reference CollinsCollins, 1982; Reference Seaberg, Seaberg, Le, Hooke and WibergSeaberg and others, 1988; Reference FountainFountain, 1993; Reference Hock, Le and HookeHock and Hooke, 1993; Reference KohlerKohler, 1995; Reference Nienow, Sharp and WillisNienow and others, 1996). In this paper, we focus on the investigation of the instantaneous structure of a subglacial flow path at Unteraargletscher, a temperate valley glacier in the Bernese Alps, Switzerland. We present results from series of tracer tests which were conducted in quick succession over two diurnal discharge cycles during the ablation season 2000 to assess the variability of transport parameters on a time-scale of 1 day.
Observations
Over two diurnal dischargecycles in August and September 2000, injections of sulforhodamine B into a moulin near the central flowline approximately 4.5 km up-glacier from the terminus were repeated at intervals of about 2 hours and were accompanied by simultaneous measurements of discharge in the proglacial stream. Figure 1 shows the concentration time series recorded in the proglacial stream which resulted from these injections. Both series yielded sets of easily separable and independent TBCs. The individual TBCs were analyzed using a simple advection–dispersion model (ADM). For an instantaneous input of a tracer mass m (kg), the one-dimensional analytical solution for the tracer concentration c (kgm–3) at a distance x (m) from the injection point (Reference Kreft and ZuberKreft and Zuber, 1978) is given by
Time series of tracer concentration in the proglacial stream of Unteraargletscher, resulting from injections of sulforhodamine B into a moulin ~4.5 km from the glacier terminus. Each experiment covers an entire diurnal cycle in August (a) and September 2000 (b).
where Q (m3 s–1) is the discharge, assumed to be steady, and t (s) denotes time. The transport parameters velocity v (m s–1) and the dispersion coefficient D (m2 s–1) are then determined by fitting the model to the rising limb of the TBC such that the peak concentration is conserved (Reference Sauty, Kinzelbach, Custodio, Gurgui and FerreiraSauty and Kinzel-bach, 1988). If the tracer transport is affected by storage retardation processes (e.g.Reference Willis, Sharp and RichardsWillis and others, 1990; Reference FountainFountain, 1993), the ADM can account only partly for the declining limb of the TBC. In this case, the difference of the area under the calculated and the measured curve is used to express storage retardation (Reference Hauns, Jeannin and AtteiaHauns and others, 2001).
Figures 2 and 3 show the derived transport parameters for the experiments in August and September, respectively. Also shown are the bulk discharges in the proglacial stream averaged over the duration of each individual tracer breakthrough. Note that due to the unequal time intervals between injections, the horizontal axes cannot be directly translated into time axes.
Results of repeated tracer injections over one diurnal discharge cycle in August 2000. (a) Bulk discharge in the proglacial stream averaged over the duration of tracer break-through. Transit velocities (b) and dispersion coefficients (c) are derived from fitting the ADM to individual TBCs. The horizontal axis on top refers to the time of tracer injections.
Same as Figure 2, but for September 2000.
The transit velocity during the experiment in August varies between 0.34 and 0.75 m s–1. In addition, the velocity variation precedes that of bulk discharge in phase such that the velocity at the rising limb of discharge is higher than that at the falling limb (Fig. 2). During the September experiment, the course of discharge is characterized by an increasing trend superimposed on a diurnal variation (Fig. 3). We surmise that this increase is due to an increase in meltwater input to the glacier with the return of sunny and warm-weather after a cold period. Despite the increasing trend in discharge, the transit velocity shows a diurnal variation similar to that in August.
The accuracy of the transport parameters v and D depends primarily on the degree to which the time of the tracer breakthrough can be determined. As a conservative estimate of this error, we assign an uncertainty of ±1min. With peak concentrations of the fastest TBCs occurring about 100 min after tracer injection, this uncertainty translates into inaccuracies ≤1% for both v and D. These inaccuracies are smaller than the size of the symbols plotted for each data point in Figures 2 and 3.
The observed fast transit velocities of the tracer through the glacier (Fig. 2 and 3) in conjunction with the reappearance of the tracer at the glacier portal in narrow TBCs (Fig. 1) suggest water flow through a channelized drainage system. For a given conduit with constant roughness and diameter, it has been found that D is approximately proportional to v (Reference TaylorTaylor, 1954). Investigations of the v–D relationship have been used to infer the dispersive behavior of subglacial flow paths from the slope of a linear regression of v on D (e.g. Reference Seaberg, Seaberg, Le, Hooke and WibergSeaberg and others, 1988).
Figure 4 illustrates the observed relationships between transit velocity and dispersion coefficients for the August and September experiments. The significant scattering in both v–D relationships (Fig. 4a and c) inhibits the determination of a linear functional relationship. In fact, if the chronological order of the data is taken into account (Fig.4bandd), the v–D relationships reveal hystereses in anticlockwise direction such that the dispersion coefficients are systematically higher during declining than during increasing velocity.
Velocity–dispersion relationships derived from the tracer tests in August (a, b) and September (c, d). The arrows indicate the chronological order.
Interpretation
Reference RöthlisbergerRöthlisberger (1972) postulated conduits with a dynamic geometry which is controlled by the counteracting processes of melt enlargement due to dissipation of potential energy of the flowing water vs creep closure of the viscous ice. While Reference RöthlisbergerRöthlisberger (1972) considered the steady state, Reference Spring and HutterSpring and Hutter (1982) investigated the time-dependent behavior of such conduits. It was found that in response to a diurnally varying discharge, the conduit diameter evolves such that it is larger during declining than during rising discharge. This in turn causes the velocity to be higher at the rising as compared to the declining limb of the discharge variation.
Such behavior is indicated by the observed variations of discharge and transit velocity as illustrated in Figure 2, whereas in response to a continuously increasing discharge, the conduit presumably grows steadily. In this case (Fig. 3), low discharge during the night would cause a pronounced velocity depression (Reference SchulerSchuler, 2002).
Based on this interdependence of discharge, conduit size and flow velocity for a dynamic Röthlisberger conduit, we suggest a positive relationship between conduit size and dispersion coefficient to explain the observed hysteresis in the v–D relationship (Fig.3).
Tracer transport simulation
To test the hypothesis stated above, numerical tracer tests were performed to study systematically the influence of changes of conduit geometry on tracer dispersion. Therefore, a model is adopted which was used by Reference Hauns and JeanninHauns and others (1998) and Reference HaunsHauns (1999) to investigate the retardation of water flow in karst conduits. The approach consists in describing the three-dimensional velocity field of turbulent discharge through a given geometry. Subsequently, tracer transit is simulated by a scalar volume which is propagated through this flow field by advection. Hauns and others (1998) have validated the model with laboratory experiments. In this study, we apply this model to several idealized conduits which differ in size, shape and roughness. For each conduit scenario, TBCs were generated for different hydraulic conditions.
Model description
The Navier–Stokes equations are solved with a k–ε–turbulence model (e.g. Reference Harlow and NakayamaHarlow and Nakayama, 1967; Reference Launder and SpaldingLaunder and Spalding, 1974) using the finite-volume method. We used the commercial package CFX-4.3 of AEA Technology, which provides a complete software environment to solve turbulent flow problems. In our study, the three-dimensional velocity field was computed for an idealized conduit which was represented by a 10m long rectangular geometry (Fig. 5). Along the conduit walls, a no-slip condition is applied. The amount and direction of water discharge is controlled by prescribing a constant velocity at the inlet and a no-pressure condition at the outlet plane.
The idealized conduit geometry through which numeri tracer tests were performed. Water flow is in x direction.
In the turbulent velocity field, transport of a conservative tracer is then simulated by advection of a scalar representing the tracer concentration. As such, the approach is valid only if the tracer’s influence on the fluid density is negligible. This is typically the case in the experiments at Unteraargletscher where maximum observed tracer concentrations were 530 ppm (Fig.1).
Tracer injection is simulated by applying a constant scalar flux over a 2 s interval at the conduit inlet. After transport of the tracer through the conduit, the scalar field is sampled at the outlet, and its time series is analogous to a TBC.
Explicit computation of the flow field and simulation of the tracer transport is restricted to idealized conduits at the 101m length scale because at larger length scales the computational effort becomes excessive. However, tracer tests in the field are typically performed over length scales of 102– 103 m. Therefore, to examine tracer transport in conduits at the scale of hundreds of meters, we used the transfer function approach (Reference Hauns, Jeannin and AtteiaHauns and others, 2001) to expand the model to the 102 m length scale. This approach entails that once the transfer function of a 101m length-scale conduit is obtained using the numerical model, simulations of tracer transport in conduits of any desired length can be generated by concatenation of such transfer functions (Reference Box, Jenkins and ReinselBox and others, 1994).
In this study, the transfer function approach was used for upscaling a 10m long conduit (Fig. 5) to a length of 240 m. Figure 6 illustrates how the shape of a return curve evolves with increasing distance from the injection point at the conduit inlet. Tracer concentration was logged at several distances along the conduit. The first breakthrough curve displays an asymmetrical shape which is characterized by a fast rise and a slow decline, indicative of pronounced retardation. With increasing distance, the asymmetry of the return curves decreases gradually, peak concentrations decline and curves become broader. This evolution was quantified using the ADM as described above. Figure 7 shows how the dispersion coefficient increases with distance to the injection-point length concurrent with a decrease of retardation. This observation is in agreement with findings of Reference Hauns, Jeannin and AtteiaHauns and others (2001) and demonstrates how tracer retardation at small scales contributes to hydraulic dispersion at larger spatial scales.
The evolution of TBCs along the conduit with increasing distance from the injection point (6, 60, 120, 180 and 240 m). Tracer concentration is shown relative to its input value.
The evolution of dispersion coefficient (solid line) and retardation (dashed line) of the TBCs shown in Figure 4 with increasing distance from the injection point. Retardation is quantified by the portion of the tracer mass, which is not accounted for by the ADM.
Numerical experiments
The influence on tracer dispersion was systematically studied by altering the size, shape and roughness of the conduit. To characterize the dispersive behavior, we derived velocity–dispersion relationships for each of these scenarios by using a number of different flow velocities.
In a first set of experiments, the height H of the conduit geometry was varied to generate conduits with different cross-sectional areas A. Table 1 lists details of individual experiments.
Height H, width W and cross-sectional area A of different conduits used to study the effect of conduit size on tracer dispersion. vin is the mean flow velocity
Since the flow velocity is zero at the conduit walls, the flow field and therefore also the tracer transport is affected by the fraction of water volume that is in contact with the wall. This influence can be examined by systematically varying the hydraulic radius Rh = A/U, where A is the cross-sectional area that is occupied by the water and U is the wetted perimeter of the conduit. In a second set of experiments, we studied this effect using four conduits with different cross-sectional shapes while keeping A = 0.25 m2 constant (Table 2). In each case, a mean flow velocity vin = 0.3 ms–1 was applied.
Height H, width W and hydraulic radius Rh of different conduits used to study the effects of conduit shape on tracer dispersion
In a final set of experiments, TBCs were generated to investigate the effect of roughness on tracer dispersion. In general, roughness enhances turbulence in the flow field and thereby affects tracer transport. Rectangular obstacles (0.2×0.2×0.1m) were introduced at the bottom of the conduit to represent roughness of a typical glacial stream bed (Fig. 8a). In another scenario, the influence of the bottom roughness was increased by widening the conduit. To conserve the cross-sectional area, the ceiling was shaped asymmetrically (Fig. 8b). These situations were simulated using again vin = 0.3 m s–1.
Conduit geometries containing roughness elements. Apart from the obstacles, the conduit in (a) is identical with that shown in Figure 5. The asymmetrical shape of the conduit in (b) was introduced to enhance the influence of bottom roughness. Main water trajectories are indicated by arrows.
Model results
Figure 9 shows v–D relationships resulting from the 11 experiments performed to study the interdependence of tracer dispersion and cross-sectional area. The curves feature a positive relationship between v and D, with the slope of the curves increasing towards higher velocities. The dependence of dispersion coefficients on cross-sectional area is displayed by several v–D relationships obtained for different conduit sizes. We find that a larger conduit yields systematically higher dispersion coefficients than a smaller one.
Velocity–dispersion relationships for three conduits of different cross-sectional area. A = 0.0576 m2 for solid circles, A = 0.0677 m2 for open diamonds and A = 0.125 m2 for solid triangles.
Results from varying the hydraulic radius at constant flow velocity and constant cross-sectional area are shown in Figure 10. We can extrapolate the graph with two additional points. First, the trivial case Rh = 0 yields D = 0. Second, we adopt a value determined by Hauns and others (2001) for a perfectly homogeneous velocity field. This condition is realized for an infinitely high and wide conduit. Therefore Rh → ∞ and it was found that D∞ = 0.045 m2 s–1 for v= 0.3 ms–1. Figure 10 demonstrates that with increasing hydraulic radius, dispersion approaches asymptotically D∞.
Table 3 lists the results of the three experiments performed to investigate the effect of roughness on tracer dispersion. Comparing the conduits with and without obstacles we find that introduction of roughness causes a substantial increase of the dispersion coefficient. When the influence of roughness is further increased as considered in the third scenario, D increases even more, although the hydraulic radius of this conduit is smaller than those of the previous scenarios.
Concluding Discussion
The v–D relationships derived from our tracer experiments at Unteraargletscher reveal pronounced anticlockwise hystereses. We suggested that the dynamic geometry of a Röthlisberger channel might be responsible for this behavior. To explain the hysteresis, we supposed a positive relationship between the dispersion coefficient D and the cross-sectional area A of the conduit. While the results of the numerical experiments shown in Figure 9 indicate such a relationship, the dependence is not very strong. In these experiments, the increase of A was accomplished by increasing the height of a rectangular conduit of constant width, which involves a change of the perimeter and therefore the influence of the conduit walls. The influence of the conduit walls was expressed in terms of hydraulic radius Rh. Figure 10 reveals that an increase of Rh at constant cross-sectional area also increases the dispersion coefficient. However, this effect is relevant only for hydraulic radii below the meter scale. Towards larger values of Rh, the increase of D is attenuated and D approaches a maximum value. This indicates that the increase of D in Figure 9 would also approach an upper limit since an increase of A is accompanied by an increase of Rh. The diameter of the subglacial conduit beneath Unteraargletscher is estimated to be on the order of meters (Reference SchulerSchuler, 2002). On this scale, the effect of changes in conduit size on tracer dispersion would be very weak.
Effect of hydraulic radius on dispersion. The flow velocity was kept constant in all scenarios. The dashed lines denote extrapolations of the relationship based on data from Reference Hauns, Jeannin and AtteiaHauns and others (2001). See text for details.
Based on these results, it is difficult to explain the v–D hysteresis by a change of the conduit size alone. However, the dynamics of a Röthlisberger channel can account for the observed behavior if it is accompanied by a roughness change. It is shown by our numerical experiments that D is strongly affected by changes in roughness such that increased roughness generates enhanced dispersion (Table 3). In our model we assumed constant size and roughness along a straight, rectangular conduit, which is certainly a strong simplification of the real flow path beneath Unteraargletscher. Instead, a sub-glacial conduit in nature might be curved and might exhibit variability in diameter and roughness along its length, thereby leading to significantly higher dispersion coefficients than those obtained from our numerical experiments. In addition, measurements of subglacial water pressure (e.g. Reference Hubbard, Sharp, Willis, Nielsen and SmartHubbard and others, 1995; Reference Murray and ClarkeMurray and Clarke, 1995) suggest that flow paths expand as water spreads out along the glacier bed in response to an increase in subglacial discharge.This in turn enhances the influence of the rough glacier bed. Thus, variations in discharge can directly lead to roughness changes as required by our interpretation. During the declining limb of a diurnal discharge variation, subglacial drainage is widespread and thereforemore dispersed, whereas it is more confined and less dispersed during the rising limb.
Comparison of dispersion coefficients for three conduits of different roughness. Here, roughness is expressed as the ratio of the obstacle volume Vo to the water volume Vw
Acknowledgements
The conduction of tracer tests over entire diurnal discharge cycles involved sleepless and cold nights on the glacier. The assistance of T. Khazaleh, J. Helbing and L. Mutter during our fieldwork is gratefully acknowledged. Moreover we thank M. Hauns for his model code and F. Hermann for support with the CFX software.Two anonymous reviewers made valuable suggestions to improve an earlier version of the manuscript. This study was partly funded by Eidgenossiche Technische Hochschule grant 0-20527-98. Swiss military helicopters provided logistic support in the field.
