Field site

We conducted our experiment in August 2020 on Rhonegletscher, a valley glacier in the central Swiss Alps covering an area of ~15 km² (as of 2019), elevations from 3600 to 2200 m a.s.l. and a length of ~8 km (GLAMOS, 1881–2020). The study site was located in the ablation zone between 2450 and 2500 m a.s.l. where a bedrock overdeepening is present. There the ice is ~200 m thick (Grab and others, Reference Grab2021) and the ice surface slope is 5° southwards (Fig. 1). In summer, local daily ice melt reaches 0.1 m d⁻¹, supplying meltwater to surface streams (Landmann and others, Reference Landmann2021). The streams relevant for our study site originate from several smaller streams ~700 m up-glacier. When crossing the region of the overdeepening, the streams merge into a few larger ones, with discharges of several hundred litres per second each (Fig. S1). In this region, ice flow is compressional and only a few small surface crevasses are present. Just downstream of the study site, these surface streams enter the glacier via moulins (Fig. 1).

Fig. 1.

(a) Location of Rhonegletscher (blue dot) within the borders of Switzerland (black); (b) outline of Rhonegletscher (blue) and indication of study site (black rectangle); (c) map of the study site indicated in (b), including positions of all boreholes and artificial/natural moulins, the closest surface streams, surface and bed topography as well as the glacier outline. The glacier outline in (b) and (c) refers to the state 2016 (Linsbauer and others, Reference Linsbauer2021). Coordinates are given in the CH1903+/LV95 coordinate system. Figure S1 provides an aerial image of the site.

Artificial moulins

We created artificial moulins by drilling boreholes of ~0.1 m diameter with a hot water drill (Iken and others, Reference Iken, Röthlisberger and Hutter1976). The boreholes were placed next to relatively large surface streams and were drilled down to the glacier bed. Boreholes which drained were interpreted as to have connected to the glacier's efficient drainage system and were turned into artificial moulins by diverting the neighbouring surface stream into them. Only one (BH15, Fig. 1) of the five drilled boreholes (BH11–BH15) drained directly after the drill reached the bed. The artificial moulin created from this borehole, named AM15, was used for the experiment during the following two weeks. Meanwhile, we kept monitoring the water level in the other boreholes and eventually turned the one with the largest diurnal water level fluctuations (BH13) into an additional artificial moulin (named AM13).

BH15 was drilled on 7 August 2020 to a depth of 187 m and turned into an artificial moulin on 8 August. For the first two days of the experiment (8 and 9 August) the flow in AM15 was pressurised, indicated by a completely water-filled borehole. During the following days, however, the water level dropped continuously. The limiting factor for supply of water to AM15 was the incision rate of the surface stream's spillway to the moulin. Only about half of the stream's water entered AM15. As our experiment relies on pressurised flow, we tried to maximise the water supply through frequent enlargement of the spillway with a pickaxe. However, starting from 10 August it was not possible anymore to maintain pressurised flow with this method.

BH13 was drilled on 6 August 2020, reaching the bed at a depth of 200 m. While there was no direct drainage, 20 h later the water level had dropped by 42 m. On 20 August, we diverted a stream into it (turning it into AM13, Fig. 2a) but measurable flow was only established during the following night. AM13 was completely water-filled for the entire duration of the experiment between 20 and 21 August.

Fig. 2.

(a) Artificial moulin AM13 (diameter ~0.2 m) and the feeding surface stream. The rope carrying the conductivity-temperature-pressure sensors (CTDs) is attached to stakes drilled into the ice. The green shovel on the right side of the picture is ~1 m long, for a measure of scale. (b) Schematic setup of the experiment. Two CTDs are installed at two different depths, z ₁ and z ₂, which span a test-section of varying length (between 30 and 100 m, see Table S1) and measure conductivity σ, water temperature T _w and water pressure p _w. We use these measurements to derive the hydraulic gradient ∂ϕ/∂z, flow speed v, mean discharge $\bar Q$ and cross-sectional area S as averages over the test-section.

Instrumentation

The instrumentation in each artificial moulin consisted of two Keller DCX-22-CTD loggers measuring electrical conductivity, temperature and pressure. We refer to these loggers as CTDs. The instrumental uncertainty of conductivity was ±5 μS cm⁻¹ (microsiemens per centimetre) and the one of temperature ± 0.05 °C. The pressure sensors of the two loggers measured with an uncertainty of ±490 and ±1470 Pa, proportional to their full scale range of up to 100 and 300 m water head, respecively. Each day, we attached the two CTDs at different positions on a static rope and installed them in the artificial moulin such that the lower sensor, which was attached to the end of the rope, was located at least 10 m above the glacier bed, i.e. at ~170 m depth for AM15 and at ~190 m depth for AM13. Once installed, we did not change the setup during one day of measurements. The weight of the sensors and the downward water flow ensured that the sensors’ positions were stable. The general setup is shown in Figure 2b and a list of the sensor depths on individual days in the Supplementary material (Table S1). We always used the same two loggers and never changed their relative position.

The rope itself served as a depth and sensor distance measure. Its length was measured under tension typical of that experienced during deployments and we estimate the uncertainty of the sensor distance to be ±1 m. The sensors were lowered by hand into the artificial moulin in the morning and retrieved in the evening on a daily basis, which was necessary to download the recorded data and reset the loggers for the next day. We used a 1 s sampling interval to allow for ~8 h of recording time per deployment.

Our individual tracer experiments characterise the section of the artificial moulin between the two sensors. We call this section the test-section. Under the constraint that both sensors had to be submerged, the length of this section was chosen as large as possible (up to 100 m) in order to minimise uncertainties per unit length. This constraint was the reason to reduce the length of the test-section in AM15: six days after activation (14 August), the water level had dropped to 130 m below the surface at mid-day, when the water supply reaches its daily maximum. As a consequence, the upper sensor (at 127 m depth) was left out of the water during the entire day. Eleven days after activation (19 August) the flow was only pressurised to less than 20 m above the lower sensor.

Experimental setup

In this study, artificial moulins represent pressurised englacial R-channels rather than moulins themselves. To investigate the properties of such an R-channel, we use the pressure and depth measurements to infer the hydraulic potential gradient and salt tracer experiments to obtain the discharge and the flow speed. Using these quantities we derive the rest of the hydraulic properties as listed in Table 1.

Table 1.

Quantities measured and derived through the experiment. The two uncertainties for the water pressure correspond to the upper and lower CTD sensor, respectively.

Water flow in an R-channel (and in general) is driven by the gradient of the hydraulic potential, or hydraulic gradient for short. Neglecting the kinetic potential, the hydraulic potential is given by the sum of elevation potential and water pressure

(1)$$\phi = \rho_{\rm w} g h + p_{\rm w},\; $$

where ρ _w is the density of water, g the gravitational acceleration, h the elevation relative to a reference point and p _w the water pressure. We neglect the kinetic potential as it is commonly done in glaciology (e.g. Röthlisberger, Reference Röthlisberger1972) and since its gradient is more than an order of magnitude smaller than that of the other components (see Supplementary material). A summary of all physical constants used in this study is given in Table 2 and the symbols and units for physical quantities are listed in Table 1. In our experiment, measuring water pressures at the two sensors with known depths enables us to calculate the hydraulic gradient along the flow direction of the test-section. Since we constructed the experiment to produce purely vertical flow, we can write the hydraulic gradient as

(2)$${\partial\phi\over \partial z} = {\partial p_{\rm w}\over \partial z} - \rho_{\rm w} g \approx {\Delta p_{\rm w}\over \ell} - \rho_{\rm w} g ,\; $$

where z is the depth coordinate pointing downwards (in our case along flow), Δp _w the pressure difference and ℓ the length of the test-section. When deriving other quantities from the so-computed hydraulic gradient, we use the pressure averaged over a time period which starts when the tracer is first visible in the readings of the upper sensors and ends when it ceases to be visible in the readings of the lower sensor. This duration is mostly ~2–3 min in AM15 and between 20 and 30 min in AM13 where the velocities and discharges are smaller and the test-section is longer.

Table 2.

Physical constants used for deriving relevant quantities (listed in Table 1, Eqns (1–8)), running the size evolution models (Eqns (9) and (10)) and calculating the heat transfer parameters (Eqns (11–16))

The discharge and the water flow speed are obtained by salt dilution gauging using measurements of electrical conductivity. In the field, we performed discrete injections of NaCl with a known mass, typically 25–100 g, at time intervals spanning from 20 min to a few hours. The mass was measured with an analogue suspended scale, and we estimated the uncertainty to 2 g by repeated comparisons with a digital kitchen scale. Before injecting the salt, we dissolved it in ~1 L of water. For a well-mixed salt tracer in a confined channel, the injection mass m is equal to the time integral of discharge Q times the salt concentration c. Assuming that the discharge is constant over the course of a tracer signal, one can write

(3)$$m = Q \int_{ t_0}^{t_1} c( t) \, {\rm d}t,\; $$

with the integration interval [t ₀,  t ₁] over the time span during which the concentration is above the background level. The salt concentration c(t) is proportional to electrical conductivity σ(t), and the factor of proportionality is determined from four calibrations conducted in the field. In each calibration, we measured the conductivity at different known salt concentrations, which is used to fit a linear function c = a σ with a proportionality factor a for each sensor. The time series of c(t) in the individual tracer experiments further allow the calculation of the average flow speed v, for which we use

(4)$$v = {\ell\over \Delta t_{c_{\max}}},\; $$

with $\Delta t_{c_{\max }}$ being the time difference between the peak concentrations at the two sensors. Finally, the test-section's average cross-sectional area S is obtained with

(5)$$S = \overline{Q}/v,\; $$

where $\overline {Q}$ is the mean value between the discharges at the two sensors. Note that we only conduct this calculation when the discharges at the two sensors agree within their uncertainties. Finally, the Reynolds number, a measure of turbulence, is given by

(6)$${\rm Re} = 2\vert v\vert \sqrt{{S\over \pi}} {\rho_{\rm w}\over \mu_{\rm w}},\; $$

with μ _w being the dynamic viscosity of water.

Propagation of uncertainties

For the forward propagation of measurement uncertainties we use a Monte Carlo approach: quantities with uncertainties are represented with ensembles of 20 000 samples each from Gaussian distributions with standard deviations equal to the given or estimated uncertainties (Table 1). All calculations are then carried out with these ensembles, in turn yielding ensembles of results with the corresponding statistical distributions. These Monte Carlo uncertainty calculations are carried out using the software package MonteCarloMeasurements.jl (Bagge Carlson, Reference Bagge Carlson2020). Note that we treat the stated instrumental uncertainties (Table 1) as systematic and thus the mean of repeated measurements does not have a reduced uncertainty. The uncertainties of the derived quantities are stated as the standard deviations of the resulting distributions and are plotted with error bars in Figures 3–5.

Fig. 3.

Quantities derived for the test-section of the investigated artificial moulins (AM) (corresponding to Eqns (2–8)). The following quantities are shown: (a) hydraulic gradient ∂ϕ/∂z, smoothed with a moving average filter with window of size ±1 min (AM15) and ±15 min (AM13), which is the typical duration during which the tracer is visible; (b) discharge Q; (c) flow speed v; (d) cross-sectional area S; (e) Reynolds number Re; (f) Darcy–Weisbach friction factor f and (g) Manning roughness n′. Note that the scale of the y-axis differs between AM15 (left) and AM13 (right) and that measurements of the hydraulic gradient are logged continuously while the other quantities rely on manual tracer injections. For the hydraulic gradient, we only show measurements where the flow between the two sensors is pressurised. Uncertainties are given as one standard deviation and are represented by either a grey area (hydraulic gradient) or vertical lines (other quantities).

Friction factor

Discharge and flow speed are related to the hydraulic gradient, which is commonly captured in terms of the empirical Darcy–Weisbach equation (Weisbach, Reference Weisbach1845). In the case of a circular, vertical conduit it reads

(7)$${\partial \phi\over \partial z} = f {\rho_{\rm w}\over 4} \sqrt{{\pi\, v^{5}\over Q}},\; $$

where f is the Darcy–Weisbach friction factor, which can be understood as an effective factor absorbing any effects of non-circular conduits. However, in our setup we do not expect major deviations from the circular cross-section since the conduits of interest are both pressurised and vertical, which prevents local incision, and since we conducted the experiment within one week after drilling, a timescale where viscous ice deformation should be of minor importance. Since we measure and derive all remaining quantities, we can compute f directly. To do so, we use the spatial discretisation from Eqn (2) for ∂ϕ/∂z and take $\overline {Q}$ for the discharge. Thus the derived friction factor is also an average value over the test-section. The calculation of f through Eqn (7) is only applicable when the flow is pressurised rather than a waterfall.

In order to facilitate comparisons with other studies, we also compute the Manning roughness n′, which is often used as an alternative way to relate the hydraulic potential to discharge. For a circular channel, the conversion between f and n′ is (e.g. Clarke, Reference Clarke2003)

(8)$$n' = \left({S\over 4 \pi}\right)^{1/12} \sqrt{{\,f\over 8g}}.$$

Model for channel size evolution

According to Röthlisberger (Reference Röthlisberger1972), R-channels are opened through melting the ice walls. In Röthlisberger's formulation, the amount of melt is obtained from considerations of energy balance, where the energy for melt is provided both by frictional dissipation and by sensible heat changes due to the water temperature changing along the channel. In the majority of subglacial drainage models, the temperature is assumed to be at a constant offset from the pressure-dependent melting point (Flowers, Reference Flowers2015). Under this assumption, the temperature gradient is equal to the gradient of the pressure melting point and obtained directly from multiplying the pressure gradient with the pressure melt coefficient c _t. The fact that this approach avoids solving for the water temperature as a field variable makes it numerically cheap and attractive to apply.

In order to test this assumption, we simulate the size evolution of our test-sections with two different models. In the simpler one, we follow the aforementioned assumption, therefore fixing the temperature gradient at the melting point gradient (referred to as the c _t-gradient model), and in the second model we perform a Bayesian inversion with the temperature gradient as a parameter that we fit to the experimental data of the cross-sectional area (the free-gradient model). The latter allows us to determine a confidence interval for the temperature gradient, and we can test whether the temperature gradients predicted by the c _t-gradient model are included within that interval. In both models, we neglect the closure of R-channels through viscous ice deformation since it is insignificant in our case: after calculation it turns out that the largest closure rates are three orders of magnitude smaller than the opening rates (Table S3, Figs S3 and S4).

In the c _t-gradient model, the evolution of the cross-sectional area is described through

(9)$${\partial S\over \partial t} = -{Q\over \rho_{\rm i} L}\left({\partial\phi\over \partial z} + c_{\rm w} \rho_{\rm w} c_{\rm t} {\partial p_{\rm w}\over \partial z}\right),\; $$

where ρ _i is the ice density, L the latent heat of fusion and c _w the heat capacity of water. The first term, proportional to ∂ϕ/∂z, describes frictional heat dissipation while the second one, proportional to the melting point gradient c _t∂p _w/∂z, is related to the sensible heat. Here, we take c _t as negative, noting that it is also common to define it as positive (e.g. Clarke, Reference Clarke2003). Since c _t depends on the air saturation of the water, its value is not accurately constrained. Two values are typically used: for pure water c _t (pure) = −7.4 × 10⁻⁸ K Pa⁻¹ and for air saturated water c _t (air) = −9.8 × 10⁻⁸ K Pa⁻¹ (Harrison, Reference Harrison1975). We forward-propagate uncertainties through the c _t-gradient model by assuming that the c _t values are uniformly distributed over the interval [c _t (air),  c _t (pure)] and by using the above-determined uncertainties for Q, ∂ϕ/∂z, ∂p _w/∂z and the initial S. Both models predict the evolution of S over the time period of tracer experiments on a single day, therefore the initial condition describes the state when the first tracer was injected on the corresponding day.

For the free-gradient model, we replace the melting point gradient in Eqn (9) by a general temperature gradient

(10)$${\partial S\over \partial t} = -{Q\over \rho_{\rm i} L} \left({\partial \phi\over \partial z} + c_{\rm w} \rho_{\rm w} {\partial T_{\rm w}\over \partial z}\right),\; $$

where T _w is the water temperature. Note that more negative water temperature gradients lead to larger opening rates. We fit the evolution of S modelled through Eqn (10) to the experimental data of S by a Bayesian inversion (Bayes, Reference Bayes1958) using ∂T _w/∂z and the initial S as the fitting parameters. For the likelihood function, which compares modelled and measured S, we impose a Gaussian distribution on the uncertainties. This is done by assuming independent measurements and by using the measurement uncertainty as the distribution's standard deviation. As a prior for the initial S, we choose a Gaussian distribution around the measurement with its uncertainty as standard deviation. For the prior of ∂T _w/∂z, we impose a uniform distribution which is non-zero only for ∂T _w/∂z ∈ [−10⁻², −10⁻⁶] K m⁻¹. (See the Supplementary material, in particular Eqns (S7–S8), for a detailed description of the mentioned priors and likelihoods.) Uncertainties of the input parameters $\overline {Q}$ and Δp _w/ℓ are neglected as otherwise the number of fitting parameters would increase from 2 to ~30. We run the Bayesian inversion by employing an affine invariant Markov chain Monte Carlo (MCMC) sampler (Goodman and Weare, Reference Goodman and Weare2010) with 10⁶ iterations via the software package KissMCMC.jl (Werder, Reference Werder2021). To discretise the spatial derivative in Eqns (9) and (10), we use Eqn (2). For the time stepping we use an explicit forward Euler scheme with the time step equal to the interval between tracer injections (i.e. ~20 min).

Heat transfer

When water flows in a long R-channel under steady conditions, the water temperature approaches an equilibrium temperature T _eq (Isenko and others, Reference Isenko, Naruse and Mavlyudov2005; Sommers and Rajaram, Reference Sommers and Rajaram2020). In our case, where the channel is vertical and pressurised, T _eq changes along the flow due to the depression of the pressure melting point. Here, it is rather the offset of T _eq from the melting point that approaches an equilibrium. In the following, we refer to this offset as the offset-temperature τ = T − T _i, where T _i is the ice temperature. Note that T can either refer to the water temperature T _w or to the equilibrium temperature T _eq with corresponding offset-temperatures τ _w and τ _eq, respectively. The ice temperature T _i is assumed to be at the pressure melting point since Rhonegletscher is a temperate glacier.

From the steady-state energy balance along the channel, a water temperature equation can be derived (Eqns (S10–S16)). When stated in terms of the offset-temperatures, it reads

(11)$$z_{\rm eq} {\partial \tau_{\rm w}\over \partial z} = \tau_{\rm eq} - \tau_{\rm w},\; $$

where

(12)$$\tau_{\rm eq} = -{Q\over \pi k {\rm{Nu}}} \Big({\partial \phi\over \partial z} + \rho_{\rm w} c_{\rm w} c_{\rm t} {\partial p_{\rm w}\over \partial z}\Big),\; $$

and z _eq the equilibrating length scale

(13)$$z_{\rm eq} = {\rho_{\rm w} c_{\rm w} Q\over \pi k {\rm{Nu}}}.$$

Here, k is the thermal conductivity of water and $\hbox {Nu}$ the Nusselt number quantifying turbulent heat transfer. Solving Eqn (11) gives

(14)$$\tau_{\rm w}( z) = \tau_{\rm eq} + ( \tau_{{\rm w}0} - \tau_{\rm eq}) \exp\Big({-z\over z_{\rm eq}}\Big),\; $$

where τ _w0 is the water offset-temperature at z = 0.

If we know$\hbox { Nu}$, we can gain insights into the thermodynamics of the system by using both our measurements and the free-gradient model outputs to derive the thermodynamic variables z _eq, τ _eq and τ _w. In particular, the free-gradient model allows us to constrain the temperature gradient ∂T _w/∂z more accurately than the measurements, thus we use ∂τ_w/∂z = ∂T _w/∂z − c _t∂p _w/∂z. The Nusselt number$\hbox {Nu}$ is not accurately constrained but can be obtained from hydraulic variables via empirical formulas. Most commonly used is the Dittus–Boelter equation

(15)$${\rm{Nu}} = A\, {\rm{Pr}}^{\alpha} {\rm{Re}}^{\beta},\; $$

where ${\rm {Pr}} = 13.5$ is the Prandtl number and $\hbox { {Re}}$ is given by Eqn (6). Another formula to obtain Nu is the Gnielinski correlation (Gnielinski, Reference Gnielinski1975)

(16)$${{\rm Nu}} = {( {\,f}/{8}) ( {\rm{Re}} - 1000) {\rm{Pr}}\over {1 + 12.7 ( {\,f}/{8}) ^{{1}/{2}}( {\rm{Pr}}^{{2}/{3}}-1) }},\; $$

which additionally depends on the hydraulic friction factor f. To obtain the thermodynamic variables described above, we calculate Nu with six different parameterisations, summarised in Table 3. Five of them correspond to different sets of coefficients for the Dittus–Boelter equation (Eqn (15)) and the sixth parameterisation is the Gnielinski correlation (Eqn (16)). The parameterisation of Nu depends on the flow regime and thus the different correlations are always associated with certain ranges of Re values, which are also given in Table 3.

Table 3.

Parameterisations used to calculate $\hbox { Nu}$ and the range of Re values for which they were determined. For the ones based on the Dittus–Boelter equation (Eqn (15)), dimensionless coefficients A, α and β are given. The different parameterisations are categorised into low and high heat-transfer, corresponding to low and high Nu values, respectively, for the conditions encountered in this study.

First, we use a set of coefficients that has been the standard in glaciology (e.g. Nye, Reference Nye1976; Spring and Hutter, Reference Spring and Hutter1982; Clarke, Reference Clarke2003) even though they were derived for engineering applications with experiments of water flowing through heated pipes. The coefficients of Lunardini and others (Reference Lunardini, Zisson and Yen1986) were also obtained from laboratory experiments but specifically in an ice-water related setup. In recent years, the Dittus–Boelter coefficients were fitted to glaciological field data (Vincent and others, Reference Vincent, Auclair and Meur2010; Ogier and others, Reference Ogier2021), albeit for surface streams with relatively warm water coming from a lake. Furthermore, we apply a set of coefficients that was proposed by the theoretical study of Sommers and Rajaram (Reference Sommers and Rajaram2020) for the case where all heat is generated by frictional dissipation. Finally, we apply the Gnielinski correlation, which is a refinement of the Dittus–Boelter equation with standard coefficients but which has had limited application in glaciology to date (Ancey and others, Reference Ancey2019). For each of these Nu parameterisations, we calculate τ _w using Eqn (11) and compare it to our temperature measurements in order to test whether the different parameterisations for$\hbox { Nu}$ are applicable to our test-section.