Characterising englacial R-channels using artificial moulins

Abstract

The englacial and subglacial drainage systems exert key controls on glacier dynamics. However, due to their inaccessibility, they are still only poorly understood and more detailed observations are important, particularly to validate and tune physical models describing their dynamics. By creating artificial glacier moulins – boreholes connected to the subglacial drainage system and supplied with water from surface streams – we present a novel method to monitor the evolution of an englacial hydrological system with high temporal resolution. Here, we use artificial moulins as representations for vertical, pressurised, englacial R-channels. From tracer and pressure measurements, we derive time series of the hydraulic gradient, discharge, flow speed and channel cross-sectional area. Using these, we compute the Darcy–Weisbach friction factor, obtaining values which increase from 0.1 to 13 within five days of channel evolution. Furthermore, we simulate the growth of the cross-sectional area using different temperature gradients. The comparison to our measurements largely supports the common assumption that the temperature follows the pressure melting point. The deviations from this behaviour are analysed using various heat transfer parameterisations to assess their applicability. Finally, we discuss how artificial moulins could be combined with glacier-wide tracer experiments to constrain parameters of subglacial drainage more precisely.

Introduction

The experimental study of englacial and subglacial hydrological systems is inherently challenging due to the difficulty of accessing such systems. Various field techniques have been applied, including (1) water-level measurements, as a proxy for subglacial water pressure, in boreholes (e.g. Fountain, 1994; Hubbard and others, 1995; Harper and others, 2002) and moulins (e.g. Iken, 1972; Andrews and others, 2014; Covington and others, 2020), (2) tracer experiments between multiple boreholes, moulins or the glacier outlet to infer, for instance, discharges, transit velocities and dispersion coefficients (e.g. Sharp and others, 1993; Nienow and others, 1998; Hock and others, 1999; Schuler and others, 2004; Cowton and others, 2013), (3) direct exploration of moulins (e.g. Holmlund, 1988; Gulley and others, 2014; Mankoff and others, 2017) and (4) imaging through borehole cameras, to determine the structure of englacial drainage systems (e.g. Fountain and others, 2005; Church and others, 2019). Current developments include measurement techniques such as the Cryoegg (Prior-Jones and others, 2021), the Glacsweb wireless probe (Hart and others, 2019) or drifters (Alexander and others, 2020) that can be deployed into glacier streams to take in situ measurements. In addition to these in situ measurement techniques, remote geophysical methods have been used to infer properties of glacier drainage systems by radar and seismic imaging (e.g. Harper and others, 2010; Schroeder and others, 2013; Church and others, 2019, 2021; Egli and others, 2021) or by passively measuring seismic noise (e.g. Dalban Canassy and others, 2016; Gimbert and others, 2016; Röösli and others, 2016; Nanni and others, 2020).

However, with the current understanding of englacial and subglacial drainage systems, physical models that aim at representing these systems still require a number of assumptions and involve poorly constrained model parameters, which motivates field experiments such as the one presented in this study. The validation of such models is difficult since measurements of subglacial water pressure, temperature, discharge and flow speed are usually sparse in space and time. Moreover, these measurements often represent bulk values over large portions of the glacier (e.g. Werder and others, 2010). Irarrazaval and others (2021) used one of the most comprehensive sets of observations of a glacier drainage system (Huss and others, 2007) to invert a subglacial drainage model for channel locations and discharge conditions. They found that the available data are insufficient to uniquely constrain the model parameters and that more frequent observations at higher resolution are needed for appropriate constraint. Brinkerhoff and others (2021) arrived at largely the same conclusion as they coupled a subglacial drainage model (Werder and others, 2013) to an ice dynamics model thus enabling the use of surface velocity measurements to invert for subglacial drainage properties.

An example for the lack of experimental model validation is the concept of R-channels (Röthlisberger, 1972), (semi)circular, ice-incised, pressurised conduits that are opened by melt and closed by viscous ice deformation. For such features, melt is commonly calculated under the assumption that the water temperature follows the pressure melting point and that the frictional and sensible heat results in instantaneous melting of the ice walls. Since their first mathematical description in Röthlisberger (1972), R-channels have been incorporated into the majority of subglacial drainage models to represent the efficient drainage system (e.g. Flowers, 2015; De Fleurian and others, 2018). However, even after decades of applying the theory, the aforementioned assumptions are yet to be validated and the parameters still need better constraints. The latter is particularly important in terms of the friction factor, which controls frictional energy loss and thus the energy available for melting the surrounding ice by parameterising the hydraulic roughness of the channel.

In this paper, we present a new approach to experimentally assess the properties of an englacial R-channel. By using boreholes connected to the subglacial drainage system and supplying them with water from a glacier surface stream, we create englacial features that route the water from the surface to the glacier bed; we use the term artificial moulin to refer to such a feature. The part of an artificial moulin which is below the water level experiences pressurised water flow and is thus an artificial, englacial R-channel, which is the focus of this study. Although the concept of artificial moulins as a field technique is not completely new, it has only been used for piezometric observations of the water level (e.g. Röthlisberger and others, 1979). Here, we combine pressure and temperature measurements with salt dilution gauging to derive the temporal evolution of the channel's hydraulic gradient, discharge, water flow speed and cross-sectional area. This way we are able to examine the hydraulic and thermodynamic properties of a single R-channel in detail. The highly resolved time series of both the hydraulic friction factor and of the cross-sectional area enable us to directly test the basic theory of Röthlisberger (1972). In contrast to a natural moulin with its complex geometry, an artificial moulin enables safe deployment of instruments and thus precise measurements, which in turn allows the calculation of R-channel properties.

Methods

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, 2021) 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, 2021). 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, 2021). 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, 1976). 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, 1972) 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, 2020). 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, 1845). 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, 2003)

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

Model for channel size evolution

According to Röthlisberger (1972), 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, 2015). 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, 2003). 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, 1975). 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, 1958) 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, 2010) with 10⁶ iterations via the software package KissMCMC.jl (Werder, 2021). 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, 2005; Sommers and Rajaram, 2020). 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, 1975)

(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, 1976; Spring and Hutter, 1982; Clarke, 2003) even though they were derived for engineering applications with experiments of water flowing through heated pipes. The coefficients of Lunardini and others (1986) 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, 2010; Ogier and others, 2021), 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 (2020) 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, 2019). 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.

Results

Measurements

In the time period 8–21 August 2020, we performed a total of 70 tracer injections in the two artificial moulins AM15 and AM13. After calculation of the discharge, we select 40 tracer experiments suitable for further considerations. The selection is based on two criteria: we require that (1) both sensors are at least 5 m below the water table, recognisable from the pressure record and (2) the discharge estimates at the two sensors are consistent, i.e. overlapping within their uncertainty ranges. In AM15, which we used for the experiment from 8 to 19 August, the limiting factor is the continuously dropping water level (see Fig. S2 and ‘Data/code availability’ section). In AM13, where the experiment lasted only for two days (20 and 21 August), the selection of only the second day is due to insufficient water flow on the first day.

Selected measured and derived quantities are summarised in Figure 3 (see also Table 1). In AM15 (left panels), discharge and flow speed measurements of 9 and 11 August capture daily variations, with peak values in the early afternoon. The flow speed decreases over the six day measurement period, whereas the discharge and the hydraulic gradient do not show a recognisable trend. Consequently, the cross-sectional area increases over time, with a clear linear trend on 9 August. Note that the inflow into AM15, and therefore the melt enlargement, was stopped during the nights from 8 to 9 and 9 to 10 August. In AM13 (right panels of Fig. 3), the absolute hydraulic gradient is low throughout the single day of measurements. Accordingly, the discharges, velocities and the opening rate of the cross-sectional area are much smaller than in AM15. The friction factor f fluctuates between 0.1 and 0.4 on 8, 9 August (AM15) and 21 August (AM13), corresponding to n′ between 0.02 and 0.04 s m^−1/3. In AM15, the friction increases from 10 August onwards to reach a maximum of f = 13 on 13 August.

Observed and modelled size evolution

To compare the modelled opening rates to the measurements, we focus on the two days with the most dense measurements, i.e. 9 and 21 August 2020 (using AM15 and AM13, respectively); we will refer to them as 9-Aug/AM15 and 21-Aug/AM13, respectively. On both days, the size evolution inferred from the c _t-gradient model overlaps with the measurements (Fig. 4) meaning that the c _t-gradient model is in principle successful at predicting the evolution of the cross-sectional area. However, the prediction of the free-gradient model follows the measurements more accurately.

Fig. 4. Size evolution of the cross-sectional area S of artificial moulin (a) AM15 on 9 August 2020 and (b) AM13 on 21 August 2020 as calculated from measurements (black markers with error bars), the c _t-gradient model (blue area) and the free-gradient model (grey lines resulting from some of the individual MCMC iterations fitting ∂T _w/∂z and the initial S). The lower panels show probability densities for the temperature gradient for 9-Aug/AM15 (c) and 21-Aug/AM13 (d): temporal mean of measured temperature gradient $\overline {\Delta T_{\rm w}}/\ell$ (black line); temporal mean of melting point gradient as predicted by the c _t-gradient model (blue); posterior distribution of ∂T _w/∂z inferred by the free-gradient model (grey) with 95% confidence interval (CI) highlighted explicitly (dark grey).

Figures 4c, d show the temperature gradients corresponding to the c _t-gradient model, the free-gradient model and the measurements. We compare the posterior distribution of ∂T _w/∂z of the free-gradient model to the temperature gradient distributions both calculated from temperature measurements and predicted by the c _t-gradient model. The observed temperature gradients agree with both models’ results, but the measurement uncertainties are too large to make any further statements about the question whether the temperature gradient follows the depression of the pressure melting point. For 9-Aug/AM15, the uncertainty is notably larger since the length of the test-section ℓ is only 50.6 m (compared to 101.5 m for 21-Aug/AM13). With a test-section of ~100 m, achieving an observational accuracy comparable to the model results would require measurement uncertainties of absolute temperature smaller than ~0.002 °C.

The temperature gradient of the c _t-gradient model (the melting point gradient) overlaps with the 95% confidence interval of the posterior from the free-gradient model inversions on both days, meaning that the c _t-gradient model is able to explain the derived S time series for each day individually. However, the overlaps are only small and on opposite ends of the c _t range for the two days: for 9-Aug/AM15, only the models with c _t values close to the lower limit c _t (pure) (smaller absolute gradient) lie within the confidence interval, whereas for 21-Aug/AM13 only the values close to the upper limit at c _t (air) lie therein. Since there is no indication for the c _t values being different between the two days, it follows that the c _t-gradient model is not able to explain the time series of both days simultaneously. This means that the actual water temperature gradients are different from the ones predicted by the c _t-gradient model.

The models can also be used to calculate the relative contributions of frictional dissipation and sensible heat change to ice melt. The c _t-gradient model predicts that averaged over the day, sensible heat provides $66\pm 7\percnt$ of the total energy for 9-Aug/AM15 and $90\pm 11\percnt$ for 21-Aug/AM13, when the discharges and hydraulic gradients are smaller. In the free-gradient model, the contribution of sensible heat is $57\pm 4\percnt$ for 9-Aug/AM15 (smaller than in the c _t-gradient model) and $94\pm 6\percnt$ for 21-Aug/AM13 (larger than in the c _t-gradient model). Compared to the free-gradient model, the c _t-gradient model predicts total melt rates that are on average $26\pm 9\percnt$ larger for 9-Aug/AM15 and $40\pm 6\percnt$ smaller for 21-Aug/AM13. The opening rates of both models and their partitioning are plotted and listed in the Supplementary material (Fig. S3, Table S3).

Heat transfer

Since the heat transfer calculations depend on the results from the size evolution models, we focus on the same days here, 9 August (AM15) and 21 August 2020 (AM13). Figure 5 shows the daily mean values for the thermodynamic variables z _eq, τ _eq and τ _w, as well as τ _w − τ _eq for the different Nusselt number parameterisations. Note that we treat the different coefficients for the Dittus–Boelter equation as individual parameterisations.

Fig. 5. Daily means of thermodynamic properties, derived from six different ${\rm {Nu}}$ parameterisations (Eqns (15) and (16), Table 3): (a) Nusselt number $\hbox {Nu}$, (b) equilibrating length scale z _eq, (c) equilibrium offset-temperature τ _eq, (d) water offset-temperature τ _w and (e) difference between water offset-temperature and equilibrium offset-temperature. Panel (d) also shows the 95% confidence interval of the τ _w measurements. The figure only shows daily mean values because the daily variations are much smaller than the difference between the daily mean values of the $\hbox {Nu}$ parameterisations (see Fig. S6).

The results can be divided into two groups: a low heat-transfer group, comprising the Dittus–Boelter equation with coefficients either from the standard formula (e.g. Clarke, 2003), from Sommers and Rajaram (2020) or from Lunardini and others (1986), and a high heat-transfer group comprising the Dittus–Boelter equation with coefficients from Vincent and others (2010) or Ogier and others (2021), or the Gnielinski correlation (Gnielinski, 1975). The low-heat transfer group shows lower$\hbox {Nu}$, longer equilibrating length scales z _eq, higher equilibrium offset-temperature τ _eq and higher water offset-temperature τ _w than the latter group, and vice versa. The discrepancy between the two groups is large: for example, τ _w is ~0.12 °C for the former and ~0.02 °C for the latter. Note that the close correspondence of the results within each of the two groups are particular to our setting and that the grouping would break apart at much higher or lower Reynolds numbers (the parameterisations depend on this number with an exponent varying from 0.58 to 1).

The Reynolds numbers that we obtain are between 6 × 10⁴ and 3 × 10⁵ (Fig. 3) and thus mostly in the ranges over which the Nu correlations are stated to be valid (Table 3). The only exception is the correlation of Ogier and others (2021), but their Re range only represents the conditions under which the experiments were conducted and therefore the validity is not necessarily exclusive to this range.

Between 9-Aug/AM15 and 21-Aug/AM13, all parameterisations show the same differences for all thermodynamic variables except for τ _w. For instance, all parameterisations predict a lower τ _eq for 21-Aug/AM13 than for 9-Aug/AM15. In terms of the difference between τ _w and τ _eq (Fig. 5e) all parameterisations predict a water offset-temperature above the equilibrium offset-temperature for 9-Aug/AM15 and below the equilibrium offset-temperature for 21-Aug/AM13. Comparing the measured τ _w with values from the different$\hbox {Nu}$ parameterisations (Fig. 5d) shows that all parameterisation are able to explain the measurements to within the 95% confidence interval except for the one by Sommers and Rajaram (2020).

Discussion

The results show that artificial moulins can provide insights into the hydraulics and evolution of an englacial R-channel. To our knowledge, the time series that we derive for hydraulic gradient, discharge, flow speed and cross-sectional area are the first of their kind.

Friction factor

Within the first two days after activation of the artificial moulins, the derived Darcy–Weisbach friction factor lies between 0.1 and 0.5. This is well within the range of values that are typically used in modelling studies, between 0.01 and 0.5 (e.g. Clarke, 2003; Werder and others, 2013). However, between 10 and 13 August the friction factor determined for AM15 increases considerably, reaching a value of 13. A possible explanation for this increase could be the development of scallops, periodic cavities eroded into the ice, which have been shown to form spontaneously when a flat ice surface is in contact with flowing water and which can reach a steady-state geometry (e.g. Bushuk and others, 2019). To further investigate the evolution of the friction factor, in particular to determine whether it eventually reaches a steady value, it would be necessary to carry out measurements over a longer time period. This would require a site with higher subglacial water pressure to keep the flow in the artificial moulin pressurised such that it continues to operate as an englacial R-channel.

However, even if such a steady-state was reached, the artificial moulin would still differ from an R-channel that has formed by itself. The main difference would likely be the shape, as natural R-channels are not perfectly straight and can have wide and flat cross-sections (Hooke and others, 1990; Church and others, 2019). In particular the sinuosity can have a major impact on the friction factor, at least in the case of subglacial conduits (Chen and others, 2018). Both an elliptic cross-section and sinuosity would increase the friction, thus the values we derive for the friction factor likely represent lower bounds. Note that the inclination of the channel does not play a role as long as the flow is pressurised, and thus our results are transferable to horizontal channels.

The comparison to subglacial channels is difficult since the contact to the bed adds complexity and also a dependence on the type of the glacier bed. Chen and others (2018) derived friction factors of up to 2.34 for a pressurised channel with sediment bed and ice roof, whereas Gulley and others (2014) reported friction factors of up to 75 (for open channel flow conditions). The comparison to our results emphasises that a single friction factor, constant in space and time, cannot be effectively used to parameterise the dissipation of energy in both englacial and subglacial systems.

Comparison of size evolution models

The time series of the hydraulic variables in the test-sections allows us to test the most commonly used size evolution model, where the temperature gradient ∂T _w/∂z is fixed at the melting-point gradient c _t∂p _w/∂z (the c _t-gradient model). To this end, we compare the c _t-gradient model to the free-gradient model, which uses ∂T _w/∂z as a parameter to fit the observed channel opening rates (note that the potential energy dissipation term is identical in both models).

The opening rates of the c _t-gradient model differ by $\sim \! 35\percnt$ from the opening rates of the free-gradient model (cf. Table S3). When conducting a glacier-wide modelling study, employing a more advanced model including temperature as a field variable (e.g. Nye, 1976; Spring and Hutter, 1982; Clarke, 2003) would be physically more correct, but we suspect that the overall uncertainties would dominate any gains of modelling the temperature explicitly. In general, a model like the c _t-gradient model is accurate for situations similar to our field experiment, i.e. where the R-channel is at a roughly constant angle, the water temperature is already close to the equilibrium temperature when entering the channel and the equilibrating length scale z _eq is small compared to the domain of interest.

In other scenarios, for example for lake drainages involving relatively warm (say 2 °C) water, the temperature gradients will differ significantly from the melting point gradient over several equilibrating length scales (Eqn (13)) and thus provide a considerable amount of sensible heat that is not accounted for in the c _t-gradient model. In these situations, modelling the water temperature explicitly is indispensable (e.g. Clarke, 2003). Slightly different conditions arise in the case of channels along adverse slopes, where the pressure melting point increases along flow, reducing the amount of energy available for melting. As a consequence, τ _eq is smaller (see Eqn (12)), becomes zero when the adverse slope reaches the supercooling threshold (e.g. Hooke, 1991; Alley and others, 2003; Werder, 2016) and will be negative beyond that threshold. Indeed, beyond the threshold the only available energy for water warming and wall melting is the heat advected by the water. If that heat is not sufficient, refreezing and thus channel shutdown will occur. Modelling such a scenario will benefit from an explicit treatment of the temperature if the lengths of the adverse slopes are comparable to or smaller than the equilibrating length scale; this would, for instance, allow the model to capture channels existing on short, adverse slopes beyond the supercooling threshold.

Heat transfer

Using the thermodynamic variables z _eq, τ _eq and τ _w, we try to explain why the fitted temperature gradients are smaller than the melting point gradient for 9-Aug/AM15 and larger for 21-Aug/AM13. Theoretically, Eqn (14) predicts that the water offset-temperature τ _w approaches the equilibrium offset-temperature τ _eq exponentially. For 9-Aug/AM15, τ _w must be approaching τ _eq from lower temperatures, as c _t ∂p _w/∂z < ∂T _w/∂z (Fig. 4c) means that τ _w is increasing, and vice versa from higher temperatures for 21-Aug/AM13 (Fig. 4d). This is consistent with the result that irrespective of the$\hbox {Nu}$ parameterisation used, τ _w < τ _eq for 9-Aug/AM15 and τ _w > τ _eq for 21-Aug/AM13 (Fig. 5e). However, we cannot confirm this with direct observations due to the large uncertainties in τ _w (compare Figs 5c, d). The question is now whether the discrepancy is due to the two days and sites having different τ _w, τ _eq or both.

The heat transfer calculations (Figs 5c–e) show that, irrespective of the used$\hbox {Nu}$ parameterisation, the equilibrium offset-temperature τ _eq is higher for 9-Aug/AM15 than for 21-Aug/AM13, while τ _w is very similar on both days (for each$\hbox {Nu}$-parameterisation, τ _w is in-between the τ _eq values of 9-Aug/AM15 and 21-Aug/AM13). This suggests that the discrepancy in the fitted temperature gradients is due to a difference in τ _eq between the two days and sites, rather than due to a difference in temperature at the entry of the test-section. To answer why τ _eq is different between the two days and sites, we turn to Eqn (12) which shows that τ _eq is proportional to the ratio between the opening rate predicted by the c _t-gradient model (see Eqn (9)) and$\hbox {Nu}$. Both the c _t-gradient model opening rate and$\hbox {Nu}$ are larger for 9-Aug/AM15 than for 21-Aug/AM13. However, the opening rate is about six times larger for 9-Aug/AM15 (see Supplementary material, Table S3) while$\hbox {Nu}$ is only about three times larger (Fig. 5a). Thus, τ _eq is larger for 9-Aug/AM15 than it is for 21-Aug/AM13. The interpretation that the test-section entry temperature is similar on both days and sites is in line with the observation that the conditions in terms of meteorology (air temperature, cloudiness), topography (surface channel slope and morphology) and flow characteristics (fully pressurised conditions throughout the artificial moulin) were similar between the two days and sites, such that there is no strong indication that the surface water temperature nor the test-section entry temperature should be different between the two days and sites.

To interpret the temperature gradients, we can additionally take the temperature equilibrating length scale z _eq into account. Using Eqn (14), we can calculate the required water temperature at the surface entrance of the artificial moulin in order to reach the inferred τ _w at the mean depth of the test-section. Unfortunately, this back-calculation relies on a growing exponential function and thus the uncertainties also grow exponentially. The uncertainties for the high heat-transfer cases (small z _eq) are particularly high, precluding any inference on the surface temperatures. For the low heat-transfer cases, the temperature of the water entering the artificial moulin is estimated to be 0.05 ± 0.04 °C (9-Aug/AM15) and 0.35 ± 0.12 °C (21-Aug/AM13), which is plausible. Due to the large uncertainties of the surface temperatures in the high heat-transfer cases, it is not possible to discriminate between the high and low heat-transfer Nusselt number parameterisations.

The second method of analysing the plausibility of different ${{\rm Nu}}$ parameterisations is to compare their predicted offset-temperatures τ _w with values obtained from measurements (Fig. 5d). Here, the possible conclusions are also limited due to the large uncertainties of the water temperature measurements. However, we can at least exclude the Sommers and Rajaram (2020) parameterisation since the corresponding τ _w value lies outside of the 95% interval of the measurements. In our case, sensible heat advection dominates heat from frictional dissipation (Fig. S3), thus it is not surprising that this parameterisation – which is based on a case where frictional dissipation is the only heat source – is not applicable. For the other parameterisations, a clear statement is difficult albeit the parameterisations from the high heat-transfer group are closer to the mean of the measurements and thus likely more appropriate. They contain the parameterisations that are based on field studies (Vincent and others, 2010; Ogier and others, 2021), which could be a motivation for further field experiments specifically targeted to quantify heat transfer.

To summarise, we cannot ascertain the full picture as to why the inferred temperature gradients differ from the melting point gradient and why the one for 21-Aug/AM13 is below and the one for 9-Aug/AM15 above the melting point gradient. To gain more insight, future experiments in artificial moulins would need temperature measurements with higher accuracy, and they could involve several test-sections to get spatially resolved melt rates and temperature gradients. Such spatially resolved measurements would allow to constrain z _eq and thus, using Eqn (13), also the$\hbox {Nu}$ parameterisation. More accurate temperature measurements, at an accuracy of ~0.01 °C (compared to 0.05 °C in this study), would allow to discriminate between the two groups of heat-transfer parameterisations (Table 3) using Eqn (11). Hence, it would be possible to calculate$\hbox {Nu}$ and infer its parameterisation in two independent fashions, similar to the approach by Ogier and others (2021).

Summary and outlook

We found that the hydraulic friction factor of the investigated englacial R-channels is not well constrained and that it evolves substantially over the course of a few days. Modelling the evolution of the channel cross-sectional area showed that the simpler channel model by Röthlisberger (1972), which assumes that the water temperature follows the pressure melting point, cannot fully explain observations. However, it is sufficiently accurate to model the size evolution of R-channels in cases where the entry water temperature is close to the melting point.

The investigation of different parameterisations for the Nusselt number was inconclusive but showed that slightly more accurate temperature measurements or several test-sections would have allowed discriminating between them. However, since the uncertainties in the Nusselt number would probably remain considerable, using a range of the Nusselt number values in modelling studies would be prudent. The same consideration applies for the friction factor. In general, this study does not give rise to major concerns about the models that are used today in englacial and subglacial hydrology. However, we suggest that the modelling should follow a more stochastic approach (e.g. Brinkerhoff and others, 2021; Irarrazaval and others, 2021) instead of choosing a single, best-fitting value for each poorly constrained model parameter.

The potential of artificial moulins does not only lie within the investigation of the drainage in englacial R-channels or moulin-like features. We envision them to be useful also for the study of the subglacial drainage system when combined with glacier-wide tracer experiments as the latter have been shown to be significantly impacted by the flow processes within the injection moulin (Schuler and others, 2004; Werder and others, 2010). Characterising the discharge and flow speeds within the moulin using the methods described here would thus allow to partition the effects of subglacial and englacial flow in the corresponding glacier-wide tracer experiments. Moreover, accompanying such tracer experiments with novel autonomous sensing platforms such as the Cryoegg (Prior-Jones and others, 2021), the Glacsweb wireless probe (Hart and others, 2019) or drifters (Alexander and others, 2020) would greatly enhance the scientific yield. Their deployment using artificial moulins would also ensure that these sensing platforms readily reach the subglacial system, without getting stuck or being destroyed on impact when traversing step pool sequences of natural moulins.