Interactive comment on “ Short term variations of tracer transit speed on Alpine glaciers ”

These are the authors’ final comments on the TCD paper “Short term variations of tracer transit speed on alpine glaciers” and contains the response to the comments by the two reviewers Doug Benn and Robert Bingham as well as Mauri Pelto’s interactive comment. We were very glad to see that our publication was received so well and we would like to thank all of them very much for their favourable opinion on our article and for their helpful comments which helped improved the manuscript further.


Introduction
The glacier drainage system governs how meltwater is routed through the glacier, which in turn influences ice flow dynamics by affecting basal sliding rates.Due to the virtual inaccessibility of the glacier interior and glacier bed, the investigation of the subglacial drainage system relies on observations of products and the subsequent inference of the underlying Correspondence to: M. A. Werder (mawerder@sfu.ca)processes (Clarke, 2005).Characteristics of the drainage system have been deduced from point measurements such as borehole water levels, slug tests (e.g.Iken et al., 1996), geophysical methods (e.g.Walter et al., 2008) or from bulk information such as discharge recession analysis and hydrograph separation (e.g.Collins, 1979).
Dye tracer experiments have proven to comprise a powerful tool for studying the sub-and englacial drainage system and have been performed to characterise the firn aquifer (e.g.Lang et al., 1979), to investigate the drainage system on seasonal (e.g.Nienow et al., 1998) and diurnal time scales (e.g.Nienow et al., 1996a;Schuler et al., 2004).The measured tracer transit speed is the quantity most readily compared to models of glacier hydraulics (e.g.Kohler, 1995), but other parameters (e.g.dispersion coefficient, returned tracer mass, double peaks) can also be used (e.g.Nienow et al., 1996b).The high frequency tracer experiments by Nienow et al. (1996a) and Schuler et al. (2004) revealed covariations between tracer transit speed and supraglacial discharge input into the injection moulin.Modulation of the flow in the tributary by the discharge conditions in the main channel has been suggested as a possible explanation for this.However, the partitioning of the total residence time into contributions from the tributary and main channel remained unresolved until recently (Schuler and Fischer, 2009).
This publication covers the following topics: the development of a two component hydraulic model, the presentation of results from a series of high frequency tracer injections on Gornergletscher (Switzerland), the interpretation of this series and two qualitatively different series (conducted on Unteraargletscher by Schuler et al., 2004) with the model, and design recommendations for future tracer experiments deduced from the model application.First we describe the field sites (Gornergletscher and Unteraargletscher, Switzerland), the field methods and the used terminology.We then introduce the two component hydraulic model which simulates tracer transport through a moulin and a subglacial channel.
Published by Copernicus Publications on behalf of the European Geosciences Union.(Sugiyama et al., 2009;Schuler et al., 2004).Marked are the moulins used for tracer injections (M4, moulin), the tracer transit distance l (dashed lines), the maximum extent of the lake -(a) grey shaded area, the borehole used for subglacial water pressure measurements -(a) BH1, and the automatic weather station -(a) AWS.The inset shows the location of map (a) and (b) in Switzerland.
We present the results of the high frequency series performed on Gornergletscher, consisting of ten tracer injections conducted over 24 h at intervals of about 3 h followed by two further injections on the subsequent day.These experiments used a moulin into which an ice marginal lake slowly and steadily drained.This situation provided a continuously high discharge input into the moulin with a relatively small diurnal amplitude and therefore prompted a comparison to similar experiments where the supraglacial discharge had a pronounced diurnal cyclicity (Schuler et al., 2004).The observations are interpreted by applying the two component model to investigate the hydraulic context leading to the observed variations of tracer transit speed, in both our observations and those of Schuler et al. (2004).Finally, the analysis of the model behaviour allows us to suggest a new experimental design for future tracer experiments.

Setting
Gornergletscher is a ∼60 km 2 valley glacier in the Valais Alps, Switzerland.It covers an altitude range from 2200 m a.s.l. to 4600 m a.s.l. and has a length of 14 km (Huss et al., 2007).Gornersee is an ice marginal lake located in the confluence area of the two main tributaries Gornerand Grenzgletscher (Fig. 1a).The lake has an elevation of 2530 m a.s.l. and lies 5.25 km upglacier from the terminus.
The greatest ice thickness of 450 m was measured 1 km downglacier of the lake (Sugiyama et al., 2009).During 2004During -2008, the yearly jökulhlaups (glacier lake outburst floods) of Gornersee were studied using an integrated approach that employed a range of different measurements (e.g.Huss et al., 2007;Sugiyama et al., 2007;Walter et al., 2008), including tracer experiments (Werder et al., 2009;Werder and Funk, 2009).The experimental data presented in this paper were collected during a field campaign in 2006 when the lake filled until its shore reached a small moulin (marked M4 on map in Fig. 1a) into which the lake spilled over.The moulin adjusted its capacity over one and a half days after the onset of the outburst (Werder and Funk, 2009); afterward the lake discharge stabilised and was limited by the rate of spillway incision (cf.Raymond and Nolan, 2000).
The average daily rate of spillway incision and of lake level lowering was slightly more than 1 m per day, translating into a discharge between 1 and 5 m 3 s −1 .It took about three weeks for the lake to empty.At the end of this time, the canyon was about 200 m long, 5 m wide and up to 50 m deep.Unteraargletscher in the Bernese Alps, Switzerland, was the site of two series of high frequency tracer injections conducted by Schuler et al. (2004).It is a ∼26 km 2 valley glacier fed by two tributaries (Fig. 1b).From the confluence zone at ∼2400 m a.s.l., it extends about 6 km eastward terminating at ∼1950 m a.s.l. and its maximal ice thickness is ∼400 m.

Field methods
On Gornergletscher, we manually injected Uranine dye at 3 h intervals between 11:00 on 18 July and 14:00 on 19 July and on 20 July at 14:00 and 17:00.The dye was injected into the lake outlet stream, ∼200 m from its entry into the moulin M4.The flushing of the dye was always good as discharge was high.The detection of the dye at the gauging station was fully automated using a BackScat submersible fluorometer.A borehole (BH1, Fig. 1a) was hotwater drilled in 2005 and equipped with a vibrating-wire pressure transducer (Geokon 4500).Air temperature and precipitation were measured by an automatic weather station (AWS, Fig. 1a) at the northern margin of the glacier.The lake level was measured with a Keller pressure transducer (DCX-22) with 10 min logging rate.The elevation of the lake surface was determined daily by a differential GPS measurement to relate lake level to elevation.Aerial pictures were taken on 20 September 2006 which were used for photogrammetric determination of the lake basin hypsometry.Proglacial discharge was measured by a hydro-power company at a gauging station 1.25 km downstream of the glacier terminus.All times are given in Central European Summer Time (CEST, UTC + 2).More details of the experimental design used on Gornergletscher can be found in two related papers (Werder et al., 2009;Werder and Funk, 2009).
On Unteraargletscher, sulforhodamine B was used for injections into a moulin (marked in Fig. 1b).The dye concentration was continuously monitored close to the terminus (marked in Fig. 1b) using a Turner Model 10-AU field fluorometer.All injections were accompanied by simultaneous measurements of supraglacial discharge into the moulin and bulk runoff in the proglacial stream using a salt dilution method (e.g.Aastad and Sognen, 1954).

Terminology
In order to describe and discuss the tracer experiments and the accompanying model, a few concepts need to be elucidated and terms defined.The same definitions are used by Werder et al. (2009).We assume that the tracer and the water travel at the same velocity, thus the following definitions apply to both.
The time interval between the passage of the maximum concentration of the tracer cloud at two locations is the residence time ( t, sometimes called the dominant residence time).The shortest possible horizontal distance between those two locations is the transit distance ( l), and accordingly, the transit speed is the ratio of the transit distance and residence time ( v = l/ t).The actual distance travelled by the tracer is the flow path length (l).Of course, the time to traverse the flow path is also the residence time.Note that the flow path length will, in general, be longer than the transit distance due to the vertical distance covered, the geometry and the sinuosity of the flow path.The transit speed is therefore a lower bound on the channel cross-section averaged flow speed (v = l/ t).
Note that hydraulic models, including the one presented in this paper, use the flow path length and calculate flow speed, not transit distance and speed.For this reason, care must be taken when comparing experimental and model results.

Data processing
Processing of the experimental data follows Schuler et al. (2004) and Werder et al. (2009).For Gornergletscher, we correct the transit speed, as defined above, for the time the tracer spent in the proglacial stream (for details see Werder and Funk, 2009.Thus, the presented transit speeds are as if the fluorometer were positioned right at the glacier snout.This correction reduces the transit speed by about 0.05 and 0.1 m s −1 at low and high proglacial discharge, respectively.Furthermore, the tracer was not injected directly into the moulin but ∼200 m away, into the lake outlet stream.Flow in the stream was rapid, we estimate ∼1-2 ms −1 , thus the residence time in the stream was ∼1-3.5 min, compared to ∼2-3 h englacial residence time.The error in transit speed, including the correction and the effects of injecting into the lake outlet stream, is about 6% (Werder and Funk, 2009).For data from Unteraargletscher, no corrections were applied since tracer was injected only a few meters away from the moulin and dye detection was close to the glacier snout.
The lake discharge Q m into the moulin was calculated from where Q melt is the meltwater input into the lake, h lake is the lake level and A lake (h) is the hypsometric function of the lake.Q melt is calculated from a distributed temperature index model (Hock, 1999) coupled to a linear-reservoir model as applied to Gornersee by Huss et al. (2007).This model is driven by the measured air temperature from the automatic weather station and was calibrated during the filling period of the lake by matching measured and calculated lake level.The hypsometric function A lake (h) expresses the relationship between lake level and lake area.It was determined from the hypsometry of the empty lake basin (from photogrammetry) together with the known lake level-elevation relationship.The error in the lake discharge Q m can be estimated by comparing the measured and modelled filling rate of the lake before the drainage initiated.The absolute error is up to 1 m 3 s −1 and varies diurnally due to shortcomings of the linear-reservoir model.
The time series and additional data of the tracer experiments, air temperature, borehole water level and lake discharge are provided in the supplement.
p Fig. 2. The lumped element model consisting of a moulin element (left) and a channel element (right).It is driven by the measured discharge into the moulin Q m and the proglacial discharge Q p .Note that the additional discharge Q a , required for conservation of water mass, is not of interest for our purposes and is not calculated.

Model
Our aim is to simulate the passage of the tracer through the glacier drainage system for a given moulin input and proglacial discharge which is achieved using a lumped element model (c.f.Clarke, 1996) combined with a water flow speed calculation.We envisage that the traced water enters a moulin which connects to a Röthlisberger-type of channel (R channel; Röthlisberger, 1972), through which the bulk of proglacial discharge is routed (Fig. 2).
The moulin element has a cone-like geometry and is fed by Q m .The water level in the moulin is equal to the hydraulic head h at the upper end of the R channel, which depends on the discharge conditions and resistance to flow in the R channel.A turbulent flow resistor of resistance R is used for the R channel element, carrying a given proglacial discharge Q p .This model is governed by the equations which are solved for the hydraulic head h, the discharge exiting the moulin Q and the cross-sectional area of the R channel are constants (cf.Table 1), h max is the maximal possible filling height of the moulin, A(h) is the cross-sectional area of the moulin as a function of height, h p is the hydraulic head at the glacier terminus, h = h − h p is the head drop in the R channel, l is the R channel flow path length, h ob = ρ ice h ice /ρ w is the hydraulic head corresponding to flotation pressure of the ice above the channel and h = 1 2 h + h p is the mean hydraulic head in the R channel.Equation (2) assumes that the input discharge into the moulin reaches the water level in the moulin without delay.Table 2 summarises the model parameters and constants.The resistance R for a channel of circular cross-section is given by where n man is the friction factor used in the Gauckler-Manning-Strickler formulation (Chow et al., 1998).Using Eq. ( 3), the hydraulic head h can be calculated that is required to drive the given discharge Q p through the R channel element.This hydraulic head defines the water level in the moulin.Thus in this model, the R channel is independent of the moulin but not vice versa.
Exploiting our finding that model performance is not compromised by assuming a static R channel (thus now referred to as channel), Eqs.(2-4) can be simplified further.The steady-state value of the channel cross section (S c ) is related to R, l and the mean proglacial discharge Qp by setting dS/dt = 0 in Eq. (4).Furthermore, the case distinction in Eq. ( 2) is removed and h p is set to atmospheric pressure (≈0 m).The system thus simplifies to Equation ( 6) can be solved for Q where dh dt is determined by differentiating Eq. (7).

Transit speed calculation
To compute the transit speed, for comparison with the observations, we first determine the residence time of a parcel of water.Tracer diffusion and storage-release processes are neglected.The exit time of the tracer t out j from a lumped element j obeys where Q j is the discharge into the element, V j is the water volume in the element and t in j is the tracer entry time.Thus, the volume of water flowing into the element between entry and exit of the tracer is equal to the volume of water in the element at the exit time.The total tracer residence time t for an injection at time t in 1 is then given by the sum of the residence times of each element where t j = t out j − t in j .In the moulin, pressurised flow conditions prevail only below its filling height h, and the model (Eqs.2-5) assumes that injected water reaches h instantaneously.Thus, the volume of water in the moulin at the exit time (t m ) of the water parcel can be obtained by integrating the moulin cross-sectional area A from the bottom to h(t m ).We prescribe A as a linear function of height above the bed z (Fig. 3) where A t is the cross-sectional area at the glacier surface and A b is the cross-sectional area at the glacier bed.Integrating Eq. ( 11) from 0 to h(t m ) and substituting it into Eq.( 9) gives where t inj is the injection time.This equation can be solved for t m once Q m is specified and then gives the moulin residence time t m = t m − t inj .
The channel cross-sectional area S and its flow path length l are constant, hence the water volume in the channel is equal to Sl, assuming fully pressurised flow.Equation ( 9) becomes www.the-cryosphere.net/4/381/2010/The Cryosphere, 4, 381-396, 2010 and can be solved for the channel exit time t c once Q p is specified.The (static) channel volume Sl is obtained by solving Eq. ( 4) with dS dt = 0: Thus by specifying R and Qp , the channel volume is fixed without prescribing either S or l.It follows that the equations for channel cross-sectional area S, channel flow path length l and channel roughness n man comprise an under-determined system given by Eqs. ( 5) and ( 14).With t c = t c − t m , the (total) residence time t is given by which is a function of t inj .The solutions of Eqs. ( 12) for t m and t c and (13) will have to be determined numerically for general Q m and Q p .Finally, the transit speed is given by where l is the transit distance (Fig. 1).

Model configuration
The lumped element model is driven by the discharge into the moulin Q m and the proglacial discharge Q p , and uses atmospheric pressure at the terminus as the lower boundary condition.Defined in this way, three free parameters remain: A t and A b describing the geometry of the moulin and the channel resistance R. As the channel parameters S, l and n man are non-unique, we will present the range of n man and S corresponding to a range of sinuosities σ between 1 and 2.
To compare the model results to the measurements from Gornergletscher and Unteraargletscher (Schuler et al., 2004), the model was run with the measured proglacial discharge, the derived (Eq. 1) discharge into the moulin for Gornergletscher and the measured discharge into the moulin for Unteraargletscher.The measured data were interpolated using a piecewise cubic Hermite polynomial to obtain a continuous discharge function.These piecewise polynomial representations of the discharge were then integrated numerically for use in Eqs. ( 12) and ( 13).On Gornergletscher, the transit distance was l = 5250 m, the mean proglacial discharge Qp was 25.3 m 3 s −1 and we ran the model from 10:00 on 18 July to 23:59 on 20 July 2006.On Unteraargletscher, l was 4450 m, the two model run time intervals were 27 h long starting at 09:00 on 2 August 2000 and 8 September 2000 and the corresponding Qp were equal to 10.9 and 14.0 m 3 s −1 , respectively.In the Appendix A, further results are presented exploring the model behaviour with synthetic data.

Fitting procedure and error estimates
The tuning parameters (A t , A b , R) were fitted by minimising the least square differences between modelled and measured transit speed.For the experiments on Unteraargletscher we found that using different A t and A b did not improve the fit and thus we set A t = A b .On Gornergletscher however, it was essential to have different A t and A b .For Unteraargletscher, model runs denoted U1 and U2 correspond to the experiment performed in August and September, respectively.For Gornergletscher, we employed two fitting strategies.First, we fitted the model using all available tracer experiments (model run G1) and second, we used only the tracer experiments during the first 24 h and also excluding the experiment which was conducted at 17:00 18 July, just before an iceberg was blocking the spillway (model run G2).
The errors in the measured transit speed are small (1% for Unteraargletscher and 4% for Gornergletscher experiments) and were thus ignored in our estimation of the errors of the fitting parameters.However, the errors in the discharge data are larger.For the discharge measurements by salt dilution at Unteraargletscher, Schuler et al. (2004) estimate an error of 5%.The proglacial discharge on Gornergletscher was measured by the hydroelectric power company for which we estimate an error of 10%.The discharge into the moulin was derived from lake level measurements and a modelled meltwater inflow (Eq. 1) and has estimated errors of up to 20%, which are probably systematic.In addition, the interpolation gives rise to further errors.
A simple scheme was employed to estimate the influence of these errors in the model estimate of v.The moulin input and proglacial discharge data were modified by increasing and decreasing their mean and amplitude by the errors given above.For the moulin discharge on Gornergletscher, we also used a ± 2.5 h phase shift of the modelled water input into the lake Q melt (Eq. 1) to account for inaccuracies in the linear-reservoir model.The model was then fitted to the transit speed using the modified discharges.Thus, for Unteraargletscher the model was fitted to 3 4 different combinations of discharges and for Gornergletscher to 3 5 different combinations of discharges.The range of v obtained by this procedure we take as the error bounds and are presented as bands in the figures.

Results
We first present the tracer experiments and related measurements conducted on Gornergletscher and we give a short summary of the tracer experiments conducted by Schuler et al. (2004), then the results of applying the model to the experiment series are given.

Observations
The results of the tracer experiments conducted on Gornergletscher and other related measurements are presented in Fig. 4.During the observation period (18 to 20 July 2006), the weather conditions were stable without precipitation The Cryosphere, 4, 381-396, 2010 www.the-cryosphere.net/4/381/2010/and air temperature in the range 6-15 • C (Fig. 4d) and the drainage system was well established in its summer configuration (Werder, 2009).At this time, the lake had been draining into the moulin for two weeks and hence the moulin had adjusted to the greatly enhanced input.Due to the stable weather, the hydraulic conditions in the glacier drainage system were also stable, as can be inferred from the regular diurnal variations of proglacial discharge (Fig. 4a) and englacial water pressure (Fig. 4c).The water pressure head in the borehole fluctuated between ∼310 m in the morning and ∼350 m in the late afternoon.Proglacial discharge varied between 15 and 35 m 3 s −1 and was in phase with the water pressure fluctuations in the borehole.The lake discharge varied between 1.8 and 4 m 3 s −1 , except on 18 July in the afternoon, when it suddenly dropped to 1 m 3 s −1 followed by a subsequent rise to 5 m 3 s −1 .This erratic fluctuation was caused by an iceberg collapsing (17:00) at the lake outlet and blocking the spillway.The iceberg obstructed the discharge out of the lake for about four hours, and therefore was also responsible for the enhanced discharge once the lake outlet was cleared.
Figure 4b shows that the tracer transit speed varied between 0.49 and 0.76 m s −1 .The lowest tracer transit speed was measured in the 17:00 experiment on the first day when the tracer was injected two minutes before the blockage of the lake spillway.This unanticipated blockage led us to conduct more experiments on the following days to fill the gap left by this unrepresentative experiment.The replacement injection was performed on the third day at 17:00 and yielded a transit speed of 0.65 m s −1 compared to 0.49 m s −1 just after the blockage.Thus, the highest transit speed of 0.75 m s −1 was attained at 11:00 when the subglacial water pressure started to rise; the speed then dropped to 0.65 m s −1 in the afternoon (neglecting the transit speed during the blockage), rose again to 0.75 m s −1 during the night and lowered again to 0.68 m s −1 in the morning.

Summary of observations on Unteraargletscher
Schuler et al. (2004) conducted two series of tracer injections (to which we will apply our model too), each over a diurnal discharge cycle in August and September 2000 on Unteraargletscher, a temperate valley glacier in the Bernese Alps, Switzerland.Their main finding was that transit speed covaried with supraglacial discharge input into the moulin but not with proglacial discharge.The measured transit speed varied between 0.75 m s −1 in the afternoon (12:00-16:00) and 0.15 m s −1 in the early morning (00:00-04:00).The discharge input into the moulin varied over an order of magnitude between 0.3 m 3 s −1 in the afternoon and 0.01 m 3 s −1 in the early morning.

Model applied to Gornergletscher
Figure 5 shows input data, measurements and the results of fitting the model to all measured transit speeds (model run time t m varies between ∼5 min and ∼105 min, its minimum is at 06:00 and its maximum at around 16:00 (Fig. 5f solid and dashed line, respectively), except during the iceberg blockage event when it reaches 200 min.The channel residence time t c varies between ∼60 and ∼130 min and is in antiphase with t m .The total residence time t displays two maxima, one at each maximum of its components t m and t c .The range of t is between 115 and 160 min, and is smaller than the ranges of t m and t c .Table 3 summarises the fitting parameters and derived quantities for both model runs G1 and G2.The confidence intervals for R are smaller than 10% and the estimates of R from the two model runs lie within 2% of each other.The moulin cross-sectional areas show a large confidence interval for G1 and G2 of ± 30 m 2 .A t is between 35 and 110 m 2 and A b between −40 and 50 m 2 .Note that the negative value of A b does not cause unphysical effects because the total volume of water contained in the moulin is always larger than zero (i.e.h > h min at all times, cf.Fig. 3).
The root mean squared error (RMSE) of v in G1 is 0.1 m s −1 and 0.07 m s −1 in G2.Values for channel crosssectional area are 22 > S > 11 m 2 and those for roughness 0.24 > n man > 0.062 m −1/3 s for a given channel sinuosity 1 < σ < 2. Table 4 summarises the ranges and means of the fitting parameters when the uncertainties related to Q m and Q p are taken into account.It shows that the parameters are constrained better for G1 than for G2.However, even in G1, the range of values of the moulin cross-sectional area is large.R is similar for G1 and G2 and better constrained than A t and A b .

Model applied to Unteraargletscher
Figure 6 shows the input data, measurements and the results of the model fitted to the data from tracer experiments performed at Unteraargletscher in August 2000 (U1, Fig. 6a-g) and in September 2000 (U2, Fig. 6h-n).In U1, the model input Q p displays a maximum at 00:00 and a minimum at 08:00.Q m has its maximum at 14:00 and minimum between 03:00 and 06:00.The measured variation of transit speed has two maxima (0.75 m s −1 ) at 12:00 and 16:00 and its global minimum at 03:00 (0.34 m s −1 , Fig. 6b).For U1, our model reproduces the measured transit speeds well using the fitting parameters presented in Table 3.The largest discrepancies are for the injections conducted at 00:00 and at 10:00 on the second day, where the model underestimates the observed v.The error bounds on v are generally below ± 0.025 m s −1 except for the time of low moulin discharge when they reach ± 0.05 m s −1 .Apart from the outliers mentioned above, all measured transit speeds are within or very close to the error bounds.The hydraulic conditions change considerably during the passage of the tracer (Fig. 6c-e).t c stays fairly constant between 70 and 100 min, whereas t m varies between 20 and 100 min (Fig. 6f).Thus, most of the variation of t stems from t m (Fig. 6g), although the small maximum (and corresponding dip in v) at 14:00 is caused by changes in t c .
In U2, the model input Q m has its maximum at 14:00 and minimum between 22:00 and 08:00, and its amplitude is larger than in U1 (Fig. 6h).Q p increases between 09:00 and 18:00, stays at an almost constant level overnight and continues to rise again after 10:00 on the next day which is quite different to Q p of U1.The measured transit speed has its maximum at 12:00 (0.55 m s −1 ) and its minimum at 00:00 (0.15 m s −1 ), thus they are both lower than the respective extrema in August but their amplitude is the same (Fig. 6i).U2 overestimates the transit speed of the first tracer experiment and underestimates the transit speed of those conducted during the night as well as that of the last one.The error bounds on v are generally smaller than ± 0.025 m s −1 except at 20:00 and 22:00.U2 fits the measured transit speed slightly better than U1 (Table 3).The difference in hydraulic conditions at t inj and t m between U1 and U2 is most apparent in Q (Fig. 6e,l).For U2, only tracer injected between 20:00 and 23:00 actually exits the moulin during the time of low Table 3. Parameters from fitting the model to measurements from Gornergletscher (G1 and G2) and Unteraargletscher (U1 and U2).The fitting parameters and their 95% confidence intervals: A t and A b are moulin cross-sectional area at top and bottom, respectively, and R is channel resistance.The root mean squared error (RMSE) of the fit of v. Derived parameters for a sinuosity σ of the channel between 1 and 2: manning roughness n man and channel cross-sectional area S.  Q m .Tracer injected later does not exit the moulin before 08:00 after Q m has increased again.t c stays fairly constant between 150 and 100 min, whereas t m varies between 20 and 500 min, thus the variation t is dominated by t m (Fig. 6m,n).Table 3 lists the fitting parameters and their ranges for U1 and U2.The estimated moulin cross-sectional areas are 1.2 and 1.5 m 2 for U1 and U2 respectively, with an error of about ± 0.35 m 2 for both.The values for resistance are 2.1 and 1.4 s 2 m −5 for U1 and U2, respectively, both with an error of about ± 0.05 s 2 m −5 .Thus, even though the moulin size is similar for U1 and U2, the channel resistance is quite different.The RMSE for both is smaller than 0.04 m s −1 .Assuming a channel sinuosity 1 < σ < 2, the Manning roughness for U1 is 0.27 > n man > 0.076 m −1/3 s and the cross-sectional area of the channel is 12 > S > 6 m 2 .For U2, the corresponding ranges are 0.48 > n man > 0.13 m −1/3 s and 22 > S > 11 m 2 .Thus the channel properties differ significantly between August and September.Table 4 summarises the ranges and mean of the fitting parameters when uncertainties related to Q m and Q p are taken into account.The fitting parameters and RMSE for both U1 and U2 are stable within the error estimate.

Implications for experiment design
To infer the flow conditions and evolution of the channelised component of the drainage system with tracer experiments, the effects of all the system's constituents must be accounted for.Here we investigate how many tracer experiments and discharge measurements are actually needed to estimate the model parameters (A t , R) accurately.We fitted the model to all combinations of three or more tracer experiments chosen from one of the two injection series of Unteraargletscher (with 12 and 11 injections, for a total of 6142 combinations), (a) using discharge data collected only at the injection time of the chosen experiments and (b) using all the available discharge data.For (a) we find a wide range of RMSE for the predicted v when using less than ten experiments, accompanied by a wide distribution of estimated model parameters.For (b), the range of RMSE for v is much reduced compared to (a) and the distribution of estimated parameters is much narrower.With eight tracer experiments in (b), the fit becomes almost as good as with all 11 or 12 experiments: the standard deviation of the distribution of A t is 0.1 m 2 , of R is 0.01 s 2 m −5 and of the RMSE for v is 0.002 m s −1 (cf.Table 3).for three suitably chosen experiments (e.g., at 10:00, 16:00 and 20:00) the model parameters are well estimated, but only when using all discharge data available (b).

Model
An observation we made during the model development is that for a given channel resistance R of an R channel in steady state, the sinuosity σ , the channel cross section S and the Manning roughness n man are not independent.This follows from Eq. ( 14), which shows that for a given R and Qp the channel volume (needed to calculate the residence time, Eq. 9) is fixed, irrespective of σ , S and n man .This arises because channel enlargement is proportional to l −1 (Eq.4, first term on right hand side) and channel closure is proportional to S (second term on right hand side).All other variables only depend on R or Qp .Channel residence times thus only depend on R and Qp , hence we fitted R and not σ and n man .But this also means that σ and n man cannot be unambiguously inferred from tracer experiments.Even in time dependent situations, this distinction is not clear as shown by Werder and Funk (2009) and by Schuler and Fischer (2009).The latter found similar responses of modelled transit speeds to perturbations applied to roughness, sinuosity and parameters controlling the channel geometry.
The model presented here is similar to that of Kohler (1995) which was used to interpret tracer experiments conducted on Storglaciären.Kohler (1995) also uses a static channel and a cylindrical moulin, but includes an open channel flow section and assumes a constant discharge along the tracer flow path.Kohler's (1995) primary aim was to determine the extent of open channel flow within the glacier.

Observations
The transit speeds measured on Gornergletscher and Unteraargletscher show both quantitative and qualitative differences.For the Unteraargletscher experiments, the transit speed correlates well with supraglacial discharge into the moulin, exhibiting one maximum in the afternoon and one minimum in the early morning.Thus, we hypothesise that the variation in transit speed are due to changing residence times in the moulin and not in the channel.In contrast, the Gornergletscher transit speeds display two maxima and two minima, the latter coinciding with both maximal and minimal lake and proglacial discharge.We interpret these observations as follows.High discharge in the main drainage system causes high englacial water pressures, and correspondingly, a high filling level of the moulin.Therefore, the injected tracer has a longer residence time in the moulin, that more than compensates the faster flow in the main drainage system.At low discharge the situation is reversed: the passage through the moulin is faster as its filling level is lower, however the flow speed in the main drainage system is low, which more than compensates the fast flow in the moulin.Maxima in transit speed occur when flow is moderately fast in both components.
The pronounced decrease in moulin discharge, caused by the blockage event, led to lower transit speeds as the passage through the moulin was prolonged.The range of measured transit speeds on Unteraargletscher was 0.1-0.75m s −1 compared to only 0.5-0.75m s −1 on Gornergletscher.The smaller observed range of transit speeds on Gornergletscher is caused by the relatively smaller range of moulin discharge and is further aggravated by the opposing processes described above.

Model applied to Gornergletscher
The two daily maxima and minima of the transit speed are reproduced by the model applied to Gornergletscher.It displays the complex interplay between different parts of the drainage system and their forcings as described in the last section (cf.Fig. 5).Further insights into these processes can be gained from the model run S1, presented in Appendix A, using artificial input data.S1 shows that two daily maxima and minima in transit speed (Fig. 7) are possible with aboveexplained processes given a sufficiently constant moulin input discharge combined with a matching moulin size.
The model reproduces the observed transit speeds within the error bounds (Fig. 5b).The characteristic occurrence of two diurnal maxima and minima in transit speed v is robust against the large errors in Q m .This is due to a combination of two factors.First, this type of behaviour is already possible with a constant Q m as is seen in the model S1 (cf.Appendix A, Fig. 7a-g).Second, the conical shape of the moulin makes the maximum of t m narrower and the minimum broader (Fig. 5f).Without this the two maxima and minima per day in t would vanish due to negative interference of t m and t c .A conical shaped moulin also makes sense from a physical perspective as creep closure effects should make the shaft narrower further down.The variation of subglacial water pressure head 100 < h < 400 m (Fig. 5d) is larger than that observed in BH1: 310 < h < 350 m (Fig. 4c).However, subglacial water pressure measurements from the same borehole in the previous year produced a range of 280 < h < 350 m (Werder and Funk, 2009), suggesting that BH1 was not well connected during the course of the tracer experiments presented here.Nevertheless, the calculated range of h is still too large, demonstrating that the predicted channel resistance R is too high.The iceberg blockage event produced a visible reduction in both measured and modelled v, clearly demonstrating that the water flux into the moulin is a key determinant of v.
The fitted moulin cross-sectional areas (A t , A b ) are large but A t matches the observed cross-sectional area of ∼80 m 2 at the glacier surface.The water draining into the moulin from the lake had a temperature of ∼1 • C (Werder and Funk, 2009); thus a large moulin cross-sectional area at depth seems plausible.For a sinuosity σ = 1, the Manning roughness n man = 0.24 m −1/3 s is improbably high and thus the channel is likely to be sinuous or to have a low and broad cross section.For a channel with a semi-circular cross section, n man is reduced by a factor of 1.1 compared to a channel with a circular cross section, and more for a low and broad geometry (Hooke et al., 1990;Werder and Funk, 2009).
The error ranges and the confidence intervals of A t and A b are large (Table 4).This poorly constrained moulin geometry in turn affects the estimates of R. Therefore, the estimated range of R is much larger than for U1 and U2.However, even though the errors in the input data are large, we submit that the model captures the major processes determining the variation of v because it reproduces many of the observed features.Nevertheless, it is possible that there are other important processes influencing the passage of the tracer.Furthermore, this large error shows that the discharge into the moulin is an important parameter to measure accurately when performing tracer experiments.

Model applied to Unteraargletscher
The model is successful in explaining the variation of transit speed v obtained from tracer experiments at Unteraargletscher.It supports our earlier hypothesis (Sect.5.2) that the variations of v are due to changing moulin residence times ( t m ) and not due to changing channel residence ( t c ), as can be seen in Fig. 6f,g.Thus, most of the change of t m is due to the changing input discharge and not due to the changing moulin filling height.However, it is important to note that the mean of t c is comparable to the mean of t m but that the variation of t c is not as large.Thus, the variation in v arises during the passage through the moulin but its mean value is determined during both the passage through the channel and the moulin.
The model produces reasonable values for the fitted parameters and hydraulic variables, and even succeeds in reproducing minor features in the transit speed variations.For example, in U1 (Fig. 6a-g), we modelled two maxima of v in quick succession, interrupted by a small local minimum at 14:00, in agreement with observations.This minimum is one of the few qualitative features caused during the passage of the tracer through the channel: at the time the tracer of the 14:00 injection passes through the channel (Fig. 6c, dash-dotted line), Q p drops from 11 to 10 m 3 s −1 causing channel residence time to increase.This shows that although the variation of v is largely dominated by effects that can be attributed to the tracer passage through the moulin, the model also captures small features caused by the channel.
The lower error bound on v in U2 (Fig. 6i) has a large jump at 19:00 due to the same process that causes the discontinuity in S2 model presented in the Appendix A (Fig. 7h-n).The lower error bound is produced by lowering Q m by 5%, which is low enough to allow upwelling of subglacial water into the moulin.The variation of t m is much larger for U2 than U1 due to the larger variation in Q m and, in particular, due to the very low discharges during the night.
The fitted values of A t are reasonable (Table 3) and their 95% confidence intervals are small.In September, A t is larger which can be attributed to the seasonal evolution of the moulin.However, R is quite different for U1 and U2 which suggests that the channel system either had a different morphology or that different subglacial discharge conditions prevailed.The weather before the September experiments was cold and discharge was low (Schuler et al., 2004); the temperature and associated meltwater production then increased just prior to the experiments and discharge steadily increased (Fig. 6h).The channel was therefore not in steady state but increasing in size as suggested by Schuler and Fischer (2009).This limits the applicability of our static channel element and leads to exceedingly high modelled water pressure head h at the end of the period considered in U2 (Fig. 6k).A large h leads to increased retardation in the moulin and therefore an underestimated v for the two last experiments.Schuler and Fischer (2009) also modelled the transit speeds of the measurements from Unteraargletscher with a model consisting of a moulin element as used in our model but coupled to an R channel with time dependent cross section.They found different values for the moulin crosssectional area: 4 and 5 m 2 for the August and September experiments, respectively, compared to 1.2 and 1.5 m 2 found here.On Unteraargletscher, the moulin diameter was never measured, its estimate has a range of 1-2.5 m and thus cannot be used to discriminate between the conflicting model results.This mismatch arises for two reasons: first, their modelled subglacial water pressure is lower, leading to a lower filling level of the moulin and shorter moulin residence time; second, their channel residence time is shorter than ours.These two effects are compensated, in order to fit the measured total transit time, by a larger moulin cross-sectional area.
The subglacial water pressure produced by our model is too high, especially for the September experiments where it reaches 1.5 times the overburden pressure.In this respect, the model of Schuler and Fischer (2009) is more accurate and, for that reason, our model underestimates the moulin crosssectional area.Conversely, there is evidence that the channel residence time in our model is more accurate: the double maxima of transit speed observed at 12:00 and 16:00 in the August experiments are produced during the flow through the channel (Fig. 6b).Our model run U1 reproduces this feature whereas Schuler and Fischer's (2009) model does not.One cause for the discrepancy between the model results could be that Schuler and Fischer (2009) tuned their model manually and, consequently, their results might not be close to the best fit.Another factor could be that our model does not capture all the relevant physical processes occurring in the drainage system.Hence, these differences merit further investigation.

Implications for experiment design
The fitting of the model to all possible combinations of three or more tracer injections chosen from one of the Unteraargletscher datasets showed that near-continuous records of both the moulin and proglacial stream discharge are a prerequisite for running the model, whereas tracer experiments do not need to be as frequent once it has been shown that the presented model can be applied.These insights allow the formulation of a measurement strategy tailored to probe the evolution of the drainage system over several days to weeks: first, determine if the model is applicable to a specific experimental situation with a series of at least eight injections over a diurnal discharge cycle which the model is fitted; second, once the applicability of the model is established, three injections per day at rising, high and low discharge are then sufficient to monitor the evolution of the drainage pathway.The model should be validated again by conducting another high frequency series injections after a time span of days to weeks or after an event with a high impact on the drainage system.For longer study periods, it would be necessary either to use this model with an evolving R channel, and possibly also with a dynamic moulin cross section, or to fit the model independently for each day.

Conclusions
Tracer experiments were conducted to investigate the diurnal variability of the glacier drainage system under different conditions at Gornergletscher and Unteraargletscher.The supraglacial input to the injection moulin at Unteraargletscher displayed a pronounced diurnal cyclicity.In contrast, the experiments at Gornergletscher were performed using a moulin into which an ice-dammed lake drained, providing input with a much smaller relative diurnal variability.The transit speed of water flow through a moulin-channel system is affected by discharge variations in each of the subsystems.Hence, on Gornergletscher, having a reduced variability in input discharge, we expected that the experiments would provide more direct information about the subglacial component.Instead, the observed transit speeds vary in a complex fashion, displaying two maxima and minima over one diurnal discharge cycle.The influence of the moulin and channel were found to be equally important with the flow speed variations in the two subsystems in antiphase, leading to the complex observed behaviour.Conversely, on Unteraargletscher, the transit speed correlated with discharge into the moulin and displayed only one daily maximum and minimum.
Our simple two-component model of the glacier drainage system simulates water flow through a moulin and a channel.Both the moulin and the channel element share the characteristic that water flow speed is proportional to the discharge and inversely proportional to the water volume they contain.The water volume in the channel element is constant whereas, in the moulin element, it changes with the filling level of the moulin governed by the subglacial water pressure, which is, in turn, proportional to the squared discharge of the channel element (Q p ).Thus, with increasing proglacial discharge, flow speed increases in the channel element but decreases in the moulin element.The interplay of these two opposing processes leads to the complex variation of observed transit speed in tracer experiments.The qualitative differences in the Gornergletscher and Unteraargletscher experiments are captured by our model, suggesting that the interaction between inflow modulation and channel flow are indeed responsible for the observed behaviour.
Our model simulates the dominant water flow process, which is captured by the presented tracer transit speed observations.However, tracer experiments yield additional information encoded in the shape of the breakthrough curve, e.g., tracer dispersion, retention and alternate flow paths.In the supplement we included the full tracer concentration time series as well as parameter describing dispersion, retention and recovered tracer mass (they are also presented in Werder, 2009).Such additional information could be used to infer further details the drainage system (e.g.Nienow et al., 1996b) and could also be compared to simulations (e.g.Werder et al., 2009;Schuler and Fischer, 2003).
The input discharge into the moulin is a key parameter for modelling the tracer transit speed.This is demonstrated by two findings: (a) the uncertainty related to the input discharge on Gornergletscher substantially reduces the performance of the model, and (b) fitting the model to only a selection of observations from Unteraargletscher showed that frequent discharge measurements are more important than frequent tracer injections.These findings emphasise the importance to record the moulin discharge at a high frequency during the tracer experiments and we propose a measurement strategy that permits to infer the evolution of the drainage system over time scales of days to weeks.Furthermore, our results demonstrate that the sinuosity and roughness of the channel are not independently constrained by measurements of transit speed and thus other experiments are required to discriminate between them.R = 0.25 s 2 m −5 which corresponds to, for example, l = 5000 m, S = 6.9 m 2 and n man = 0.04 m −1/3 s and for the moulin, we set A t = A b = 1 m 2 .To explore the model response, we performed model runs using different constant input discharges q m in the range between 0.005 and 1 m 3 s −1 , including two detailed model runs with input discharge q m = 0.2 m 3 s −1 and q m = 0.008 m 3 s −1 to which we refer to as S1 and S2, respectively.In S1, the variation of v displays two minima of similar size, coinciding with maximal and minimal Q p (Fig. 7b).The hydraulic variables Q p , h and Q at t inj , t m and t c are almost identical (Fig. 7c-e), thus the passage of the tracer is fast compared to the change of the hydraulic variables.Note that even though Q m is constant, the discharge out of the moulin Q varies with time (Eq.8, Fig. 7e).The moulin residence time t m varies between 5 and 25 min and is in phase with Q p (Fig. 7f).The channel residence time t c varies between 17 and 35 min but is in antiphase with Q p (Fig. 7f).Since the amplitudes of t m and t c are similar but vary in antiphase, the resulting variation of their sum ( t) and, consequently, v display two maxima and minima each.It should be noted that it is essential that t m and t c are not quite The Cryosphere, 4, 381-396, 2010 www.the-cryosphere.net/4/381/2010/sinusoidal, because otherwise their superposition would produce a constant function.The results of model run S2 are displayed in Fig. 7 (right).The resulting variation of transit speed v shows only one maximum and one minimum, both occurring in the morning (Fig. 7i).There is a discontinuity in v at 06:30.The mean v is about one order of magnitude smaller in S2 than in S1; however the amplitudes in both model runs are of similar absolute size.Since the total residence time t is up to 10 h long, the hydraulic variables Q p , h and Q differ substantially during the passage of the tracer (Fig. 7j-l).Q p at t c (dashed-dotted line in Fig. 7j) is similar to Q p at t m (dashed line), as in S1, since the passage through the channel is not affected by the choice for Q m .Conversely, the hydraulic variables change greatly between t inj (solid line) and t m (dashed line) because the moulin residence time ( t m ) is long.Q varies between −1 × 10 −3 and 17×10 −3 m 3 s −1 (Fig. 7l) and Q at t inj and at t m differ qualitatively as the latter has a kink at 06:30.Both the mean and amplitude of t m in S2 are much larger than those of t c , thus t is almost exclusively determined by t m (Fig. 7n).
Varying the input discharge Q m affects only t m but not t c , and accordingly, the influence of the moulin on the total residence time t changes.Figure A11a shows t inj against v for different values of Q m .The transit speed v is fastest for the largest Q m = 1 m 3 s −1 and has its maximum and minimum at 18:00 and 06:00, respectively, coinciding with the maximum and minimum of Q p .For Q m = 0.5 m 3 s −1 a second minimum appears at 19:00.Both maxima and minima become comparable in size for Q m = 0.2 m 3 s −1 (cf.S1, Fig. 7a-g).The minimum at 06:00 disappears for Q m = 0.08 m 3 s −1 and v is in antiphase with respect to Q p .For even smaller Q m = 0.008 m 3 s −1 , a discontinuity in v (cf.S2, Fig. 7h-n) appears.Figure A1b shows a phase portrait of the timings of the extrema of v in the t inj -Q m space.At low moulin discharge Q m , there is one maxima and minima occurring in the morning.As Q m increases both extrema migrate towards later times in the day.At Q m = 0.08 m 3 s −1 a new pair of extrema arises at 04:00.At Q m = 0.5 m 3 s −1 the original t ↓ inj branch annihilates with the new t ↑ inj branch (the graph is periodic in t inj ).
Even using simple and idealised input data, the model exhibits complex behaviour, thereby providing insights into several observed characteristics of the drainage system.The qualitative behaviour of the modelled transit speed (Fig. A11) can be understood in terms of the moulin t m and channel t c residence times whose superposition gives the total residence time t = t m + t c .They both vary nearly sinusoidally, at least for Q m > 0.01 m 3 s −1 , and are in antiphase.At low Q m ∼0.03 m 3 s −1 , the total residence time t is dominated by t m and the transit speed v has its maximum in the morning when t m is short.At intermediate Q m ∼0.2 m 3 s −1 (Fig. 7a-g), t m and t c are of similar magnitude and v has two maxima and minima, the latter coinciding with the maxima of t m and t c .The resulting variation of v is comparable to the observations on Gornergletscher with the timing of transit speed extrema reproduced correctly.At even higher Q m , t m becomes negligible, as if tracer was injected directly into the main channel and v has a single maximum in the afternoon.The model behaves qualitatively differently at very low Q m < 0.01 m 3 s −1 when a discontinuity appears in most of the model variables (e.g.S2, Fig. 7h-n).It can be seen in the results of S2 that this discontinuity is related to the negative Q between 12:00 and 15:00, i.e., water is flowing from the subglacial drainage system into the moulin which is caused by the quickly increasing subglacial water pressure head h.During such periods of upwelling, the tracer cannot exit the moulin and a discontinuity is produced in t m (Fig. 7m) which is propagated to the other variables.If such a situation was encountered during a tracer experiment, it would produce a double peaked return curve; part of the tracer cloud exits the moulin but the rest is pushed back up and exits later.However, only one tracer experiment, non-overlapping with others, could yield double peaks by this mechanism.

Fig. 1 .
Fig. 1.Map of the tongue of Gornergletscher (a) and Unteraargletscher (b) with solid contours of surface elevation and dotted contours of bed elevation(Sugiyama et al., 2009;Schuler et al., 2004).Marked are the moulins used for tracer injections (M4, moulin), the tracer transit distance l (dashed lines), the maximum extent of the lake -(a) grey shaded area, the borehole used for subglacial water pressure measurements -(a) BH1, and the automatic weather station -(a) AWS.The inset shows the location of map (a) and (b) in Switzerland.

Fig. 3 .
Fig. 3.The geometry of the moulin element, with water filling height h, top (A t ) and bottom (A b ) cross-sectional area.(a) A b > 0; (b) A b < 0, in which case the apex of the cone has height h apex and the volume below the apex is negative; (c) equivalent moulin geometry to (b) provided h does not fall below h min = 2h apex .

Fig. 6 .
Fig. 6.Input, results and comparison to transit speed of the model fitted to the data from Unteraargletscher experiments in August (left, U1) and September 2000 (right, U2).The layout is identical to Fig. 5, except in (a) where crosses mark the measurements of Q p and Q m and the line is the interpolation.

Fig. 7 .
Fig. 7.The model input and output variables against injection time t inj from running the model with synthetic data: left panels (a)-(g), model run S1 with Q m = 0.2 m 3 s −1 and right panels (h)-(n) model run S2 with Q m = 0.008 m 3 s −1 .The layout is identical to Fig. 5.

Table 2 .
Model constants and fitting parameters.

Table 4 .
Parameters from the error estimate taking uncertainties of Q m and Q p into account: range of A t , mean Āt with mean 95% confidence interval, similarly for A b , Āb , R, R. Range of RMSE and mean RMSE.