Determining the evolution of an alpine glacier drainage system by solving inverse problems

Abstract

Our understanding of the subglacial drainage system has improved markedly over the last decades due to field observations and numerical modelling. However, integrating data into increasingly complex numerical models remain challenging. Here we infer two-dimensional subglacial channel networks and hydraulic parameters for Gorner Glacier, Switzerland, based on available field data at five specific times (snapshots) across the melt season of 2005. The field dataset is one of the most complete available, including borehole water pressure, tracer experiments and meteorological variables. Yet, these observations are still too sparse to fully characterize the drainage system and thus, a unique solution is neither expected nor desirable. We use a geostatistical generator and a steady-state water flow model to produce a set of subglacial channel networks that are consistent with measured water pressure and tracer-transit times. Field data are used to infer hydraulic and morphological parameters of the channels under the assumption that the location of channels persists during the melt season. Results indicate that it is possible to identify locations where subglacial channels are more likely. In addition, we show that different network structures can equally satisfy the field data, which support the use of a stochastic approach to infer unobserved subglacial features.

1. Introduction

Water flow at the ice/bedrock interface exerts controls on ice flow (e.g., Iken, 1981), glacier erosion (e.g., Herman and others, 2011), sediment transport (e.g., Beaud and others, 2018), as well as catchment hydrology (e.g., Verbunt and others, 2003). Hence, a significant scientific effort has been dedicated to characterizing the spatial distribution of subglacial systems and their temporal evolution by both field data acquisition and numerical modelling. Pioneer work conceptualized the subglacial drainage system with two main components: a network of subglacial channels with fast water flow (Röthlisberger, 1972; Nye, 1976) and a distributed system with the slow flow (Lliboutry, 1968; Weertman, 1972). In addition, it has been observed that areas of a subglacial drainage network can become temporally or permanently hydraulically disconnected, which can be described as the third component of subglacial systems (e.g., Iken and Truffer, 1997; Hoffman and others, 2016; Rada and Schoof, 2018).

The study of subglacial systems has focused on both field observations and numerical modelling. Field techniques that have been employed include water level measurements in boreholes drilled into the glacier (e.g., Hubbard and others, 1995; Rada and Schoof, 2018), dye-tracer tests in moulins (e.g., Nienow and others, 1998; Schuler and others, 2004; Werder and others, 2009; Chandler and others, 2013) and glacial speleology expeditions to survey portions of channel networks (e.g., Gulley and others, 2014). More recently, seismic and other geophysical methods have also provided insights into the evolution and structure of the subglacial systems (e.g., Gimbert and others, 2016; Church and others, 2019; Zhan, 2019; Lindner and others, 2020; Nanni and others, 2020).

In addition to field methods, a range of subglacial numerical models have been built, which we divide here into two categories: (1) subglacial routing models and (2) subglacial evolution models. Routing models need to make an assumption for the hydraulic potential to compute the water flow direction. Most commonly, the water pressure is assumed equal to the ice overburden pressure to which the elevation potential is added to get the hydraulic potential, as was first proposed by Shreve (1972). Then, a routing algorithm (e.g., O'Callaghan and Mark, 1984; Tarboton, 1997) calculates a map of the upstream area or accumulated flow, where locations with higher upstream areas are conceptualized as channels. Routing models have been widely used to study subglacial systems, with applications including channel network topology determination (e.g., Richards and others, 1996), the occurrence of subglacial lakes (e.g., Livingstone and others, 2013) and delineation of subglacial basins (e.g., Chu and others, 2016). Furthermore, Banwell and others (2013) and Arnold and others (1998), constructed a channel network structure based on routing methods on which they then calculated the water pressure and channel radii evolution. One of the limitations of such routing models is that only one single subglacial channel network is obtained and then further used for all numerical experiments. Moreover, when combining the surface elevation and bedrock elevation, different sources and spatial scales must be either downscaled or upscaled to reach a common resolution, which in turn adds uncertainties. Efforts to explore the variability of the subglacial network structures have considered the effects of varying a spatially uniform flotation factor in Shreve's equation. They have found that channel topology and subglacial basin delineation are sensitive to the choice of a flotation factor (e.g., Chu and others, 2016). A first effort to consider a spatially two-dimensional (2-D) nonuniform term which adds variability to the channel network was presented in Irarrazaval and others (2019). However, it was tested in a synthetic setup and used to infer the channel structure and hydraulic properties rather than the actual channel location.

The second category of models represents, besides water flow, the process of the evolution of the subglacial drainage system. First, 1-D models represent the evolution of one subglacial channel, for example, for describing glacial lake outburst floods (e.g., Flowers and others, 2004). More recently, 2-D subglacial evolution models became capable of the spontaneous formation of channel network structures as the outcome of simulating subglacial drainage system equations. They successfully incorporated many of the known physical processes, including channel opening (by heat dissipation) and closing (by viscous-ice creep) allowing the formation and evolution of channels and distributed systems (e.g., Schoof, 2010; Hewitt, 2011; Werder and others, 2013; Rada and Schoof, 2018). However, in real-world applications, such approaches are difficult to be implemented and often cannot produce models that match observations. Exceptions include the probabilistic inference of a lumped glacier hydrology model by Brinkerhoff and others (2016), in addition to Koziol and Arnold (2017) with uses an optimization approach to infer certain model parameters from data.

However, the inversion is applied to a lumped (zero-dimension) glacial hydrology model, and, hence, the channel network is not explicitly modelled. In any case, it is not possible to validate the produced channel networks due to the restricted access to the subglacial drainage systems because these models often require a high number of parameters whose uncertainty cannot be fully constrained by observations (de Fleurian and others, 2018).

The aim of the present paper is to infer, based on field observations and prior knowledge, the location of subglacial channel networks and their hydraulic parameters at different snapshots during the melt season. As information is scarce and contains uncertainties, a unique solution is neither possible nor desirable, and the model has to be formulated as an inverse problem using a probabilistic approach based on Bayesian inference. The proposed framework is based on a modified version of the method introduced by Irarrazaval and others (2019). It allows inferring the main characteristics and hydraulic parameters of subglacial channels from real data, in turn enabling the assessment of the channel-network development across the melt season. The modelling approach includes a parsimonious, statistics- and a physics-based model which is computationally inexpensive, allowing us to perform a large number of model runs and, thus, enabling uncertainty quantification via Markov chain Monte Carlo methods. The framework explores a range of different network structures by including a spatially nonuniform perturbation of the hydraulic potential used in the subglacial routing method. Comparison with field data are achieved by using a water flow model that provides water pressure and tracer-transit time.

The framework is applied to the trunk of the Gorner Glacier, Switzerland, during the 2005 melt season. We consider five snapshots during which the hydraulic setting is considered constant in terms of hydraulic conditions, meltwater input and topology of the subglacial system. Our goal is to determine the subglacial drainage system for each of these snapshots. For each of them, water pressure from boreholes is used to condition the subglacial system through an inversion procedure. In addition, tracer-transit times from dye-tracer tests carried out from moulins are used to constrain the water velocities in the subglacial channels to a physically reasonable range.

The remainder of this paper is organized as follows. Section 2 describes the model framework, which has its roots in Irarrazaval and others (2019), however, incorporating new adaptations to fit real-world applications. Section 3 describes the field site, datasets, modelling choices and assumptions to apply the framework to the Gorner Glacier. Section 4 shows a set of subglacial systems conditioned to data across the melt season and parameter uncertainties provided by the inversion. The last section discusses the contribution and limitations of a physical-geostatistical subglacial model and trade-offs in integrating data in complex numerical models.

2. Methods

The methodology consists of four steps. First, we use a geostatistical subglacial channel generator to produce a set of unconditioned subglacial systems states for a reference snapshot. The subglacial channel generator approximates the hydraulic potential by Shreve's equation solely to produce channel networks based on a flow routing approach (see Section 2.1). Second, we compute the hydraulic potential and water discharge inside the unconditioned subglacial systems, consisting of the above-generated network and a homogeneous distributed system, by a water flow model. The water flow model solves laminar/turbulent flow equations in steady-state conditions to recover water pressure and discharge in the channel network (see Section 2.2) and the distributed system. Third, a formal parameter inversion is performed, which enables assessing the misfit of the modelled water pressure and tracer-transit times with field observations by a likelihood function. The outcome is a set of subglacial systems conditioned to data for the reference snapshot (see Section 2.3). Lastly, for the subsequent snapshots across the melt season, the temporal persistence of the channel network throughout the season is guaranteed by allowing only small perturbations to the maximum likelihood subglacial system of the reference snapshot. The outcome is a set of subglacial systems honouring the observations during each snapshot, which allow assessing the changes in the system across the melt season (see Section 2.4). The core of the methodology is based on Irarrazaval and others (2019), with two significant improvements: (1) the subglacial channel generator has been expanded to account for channel location. (2) The parameter inference strategy has been adapted and improved to suit the application to a real-world dataset and to assess the development of the channel network across the melt season.

2.1 Subglacial channel generator

The subglacial channel generator is a statistics- and physics-based tool that creates a channel network embedded in a 2-D distributed system (Irarrazaval and others, 2019). It depends on seven parameters whose values will be inferred in our inversion procedure. The subglacial channel generator presented in this study incorporates two novelties: (1) channel location variability is incorporated by extending the parametrization and adding a geostatistical parameter. (2) The flotation factor f in Shreve's equation (Shreve, 1972) is incorporated to add more variability to the channel networks.

Note that the subglacial channel generator produces a set of channel networks embedded in a distributed system. Water pressure and flow speed in the subglacial systems is then computed by a separate model via flow equations and recharge forcing (Section 2.2).

The subglacial systems are set up in a 2-D domain considering pressurized water flow, where water follows the hydraulic potential defined by

(1)$$\phi = \phi _z + p_{\rm w}.$$

Here p _w is the water pressure and the ϕ _z = ρ _wgB is the elevation potential with water density ρ _w, acceleration of gravity g and bedrock elevation B(x,y).

The first step is to generate the channel network topology or channel architecture. Similar to previous studies based on routing algorithms (e.g., Arnold and others, 1998; Chu and others, 2016), Shreve's approximation to the hydraulic potential (Shreve, 1972) is obtained by setting p _w = fp _i in Eqn (1). Where p _i = ρ _igH is the ice overburden pressure considering an ice density ρ _i and ice thickness H(x,y). The flotation factor f is usually set to 1 but values ranging from 0.6 to 1.1 have been employed previously (e.g., Chu and others, 2016). Here, we include a nonuniform noise or perturbation field ϕ _R(x, y) which is added to Shreve's hydraulic potential to obtain an approximated hydraulic potential field ϕ*,

(2)$$\phi ^\ast{ = } \rho _{\rm w}gB + f\rho _{\rm i}gH + \phi _{\rm R}.$$

The role of ϕ _R is key in determining channel location. Flow routing can be highly sensitive to small variations in the potential surface and a small perturbation could induce large downstream deviations of hydraulic paths. In this study, ϕ _R is modelled as a Gaussian random field (GRF). The GRF is a spatially correlated random field generated from random white noise and parametrized in terms of its variance and integral scales. This is useful, as varying the variance, integral scales and modifying the white noise allows for perturbations of different magnitude, correlation in space and location, which influence the channel topology. However, one of the challenges is to gradually modify the white noise to provide a continuous transition between channel networks. To overcome this issue, we create a larger GRF from which we crop ϕ _R by varying the shift s of the cropping-region along the y dimension (Fig. 1). Tests showed that the overall channel location probability is not affected by the dimension along which s is varied. This is expected as the GRF is isotropic and large enough to accommodate several independent cropping-regions or realization of ϕ _R. The inclusion of s allows for a smoother deformation of the network, useful when exploring the model space (see in Supplementary material: Subglacial channel generator video 1). As an analogy, the variance, integral scales and shift can be seen as changing the amplitude, wavelength and phase shift of a periodic wave, which results in a simple parametrization of the GRF. Note that Irarrazaval and others (2019) found that a sufficient channel network variability can be achieved with a fixed GRF variance. Here we chose a GRF variance of 0.0058 (MPa)² equivalent to a standard deviation of ~7.7 m water column. The selected value is similar to the small-scale variations of ϕ, such that it influences ϕ* but does not disrupt the large-scale trend. Then, the subglacial channel generator considers three parameters that modify ϕ*: flotation factor f, integral scales l_xy and shift s, hereafter referred to as the structural parameters of the subglacial channel generator.

Fig. 1. Synthetic example of subglacial channel generator (water flow from right- to left-hand side). On the right-hand side colour fields correspond to the GRF (ϕ _R), and on the left-hand side to unconditioned subglacial channels. The influence of the shift s parameter is shown in (a), (b) and (c), by shifting the GRF (f and l_xy are fixed) and computing the channels. (a) and (b) show minor differences when the shift is small. (c) shows different structures as it considers a completely different part of the GRF. In (d), shift is the same as in (c), but the integral scale (l_xy) of the GRF is smaller.

Once the hydraulic potential is obtained, the D8 flow routing algorithm (O'Callaghan and Mark, 1984) is computed over the hydraulic potential ϕ* and the upstream contribution flow accumulation is computed for each cell. In this step, the spatial variations in water recharge (i.e., moulins and water recharge by tributary glaciers) are considered to compute the upstream flow contribution. The flow paths with an accumulation higher than threshold parameter c correspond to the channel network. Parameter c is set to ensure that all moulins are connected to the subglacial channels. Flow routing and channel structure are computed using TopoToolbox 2 (Schwanghart and Kuhn, 2010; Schwanghart and Scherler, 2014).

The last step is to assign radii to each of the channel network segments r _i. This is done computing the Shreve's stream order u _i (Shreve, 1966) and transforming it to a radius using the following equation (Borghi and others, 2016),

(3)$$r\lpar {u_i} \rpar = a{\rm e}^{u_ib}\comma \;$$

where a and b are parameters determined in the inversion procedure. Then, the channel network is embedded into the distributed system modelled as a homogeneous 2-D layer with transmissivity T _d. Lastly, in order to model dye-tracer transit speeds, we account for the moulins' inflow modulation time. Inflow modulation time corresponds to the time water spends in the englacial drainage system until it reaches the subglacial channels. Moreover, inflow modulation time, which has been shown to significantly affect the overall tracer transit times (e.g., Werder and others, 2010), could be estimated if the recharge hydrograph, subglacial water pressure and moulin geometry are available. As these data are unavailable for the study site, we include parameter τ_i which accounts for the inflow modulation time where the subscript i enumerates the moulins with tracer data.

In summary, the channel generator uses seven parameters. Parameters a and b determine the channel radii and T _d the transmissivity of the distributed system. The structural parameters l_xy, s and the flotation factor f control the channel location and/or network topology. The first two modify the GRF, whereas the flotation factor controls an overall trend on Shreve's equation. Finally, τ_i accounts for the inflow modulation time. All parameters are inferred in the inversion procedure.

2.2 Water flow model

A water flow model is used to compute water pressure p _w (Eqn (1)) and water discharge inside previously generated subglacial systems. Water flow is solved in steady-state conditions using the finite element GROUNDWATER code (Cornaton, 2007) based on the framework presented in Irarrazaval and others (2019). The subglacial drainage system is modelled as two coupled components: distributed and channelized. Both components of the subglacial system are coupled using a finite-element mesh with shared nodes, assuming continuity of the pressure field. This allows water and mass exchanges between the distributed system and channels that are driven by the pressure gradient.

The distributed system is modelled as a 2-D layer, discretized in 25 m edge quadrangles and with homogeneous transmissivity T _d. Water mass is conserved assuming incompressibility and pressurized flow, then Darcy-laminar flow is considered under the assumptions of a nondeformable porous medium and vertically-integrated flow equations:

(4)$$q = {-}T_{\rm d}\nabla \phi \comma \;$$

with q the flux, T _d the transmissivity of the distributed system and $\nabla \phi$ the gradient of the hydraulic potential.

The channelized system is modelled as a network of 1-D cylindrical elements of radius r. Note that the radii are fixed as we are considering a snapshot model (no time derivatives). Discharge in the channel network is computed using the nonlinear Manning-Strickler empirical formulation for turbulent flow assuming mass conservation, incompressibility and pressurized flow. The water flow inside channels Q is computed by

(5)$$Q = {-}K\nabla \phi \comma \;$$

with

(6)$$K = \displaystyle{{\alpha {\lpar {r_i/2} \rpar }^{2/3}} \over {n_m\sqrt {\vert {\nabla \phi } \vert } }}\comma \;$$

where K corresponds to the channel hydraulic conductivity, with a circular cross-section α = πr_i ², and n_m = 0.04 m^−1/3s is the Manning friction coefficient for subglacial channels (e.g., Gulley and others, 2014).

The model is set up with prescribed water recharge from the tributaries and moulins that enters via subglacial channels. In addition, discharge at the glacier's outlet is modelled as atmospheric pressure Dirichlet boundary condition. The total residence times of tracers are calculated in two steps. First, tracer-transit times are computed by integrating the advective velocity along the subglacial channels and secondly, by adding the inflow modulation time (τ_i), which is a parameter to be inferred.

2.3 Parameter inversion strategy

The inversion aims to find the subglacial networks and hydraulic parameters that enable the model to fit the observations. As data are scarce and uncertain, we expect the solution to be highly nonunique, meaning that many networks can honour the observations. Therefore, we use a maximum likelihood approach where the goal is to explore the model space and find multiple conditioned subglacial (CS) systems. Here, we define a likelihood function, which is a measure of the misfit of the field data with the model outputs, to evaluate model parameters. A general formulation of a likelihood function, assuming independent Gaussian uncertainties for the data (e.g., Mosegaard and Tarantola, 1995) is given by

(7)$$L\lpar { {\bf m} {\rm \vert }\bf d} \rpar \propto {\rm exp}\lpar {-S} \rpar \comma \;$$

where

(8)$$S\lpar {\bf m} \rpar = \displaystyle{1 \over 2}\mathop \sum \limits_{i = 1}^n \left({\displaystyle{{F_i\lpar {\bf m} \rpar -{\bf d}_i} \over {\sigma_i}}} \right)^2\comma \;$$

with S the misfit function, m = [a, b, l_xy, s, f, T _d, τ_i] the model parameters and F the forward model (water flow model and subglacial channel generator). The forward model outputs F(m) (water pressure and tracer transit times) are based on a simulation run using m as parameters. Field data d = [d₁,…d_n] correspond to the observations, with n the total number of data. The standard deviation of the errors σ = [σ ₁,…σ_n] have to be estimated and should reflect not only instrumental errors, which are relatively small, but also the effect of averaging a point measurement over a model cell and model errors. The choice of the uncertainties is a subject of debate; for example, Brinkerhoff and others (2016) chose uncertainties similar to the magnitude of daily variations. However, the setup here differs as the inversion is performed for a snapshot of time. We run several inversions to test the effect of σ in the likelihood function and select a value that avoids overfitting but represents the uncertainties that each data source is subjected to, such as field measurement errors, model errors and limitations (e.g., mesh resolution).

Each data source is subjected to field measurement errors, which are added to model errors and limitations (e.g., mesh resolution). Note that Eqn (9) is a general formulation that is subsequently adapted for each kind of data, as is commonly done for joint inversion. Assuming independence, the total likelihood is the product of the individual likelihoods, which can be stated as

(9)$$L\lpar {\bf m{\rm \vert }\bf d}\rpar \propto {\rm exp}\left({-\mathop \sum \limits_j S_j} \right)\comma \;j\in \lcub {\rm p_{\rm w}\comma \;{\rm ts}\comma \;z} \rcub .$$

Here each misfit function S_j is designed to account for the particularities of the observational data type, where S _pw, S _ts and S _z correspond to water pressure observations, tracer-transit speeds and a correction term to account for piezometric level exceeding the surface elevation (see below).

Borehole water pressure data are subject to different sources of uncertainty. Furthermore, it has been observed that neighbouring boreholes can show significant differences in water level and behaviour (e.g., Rada and Schoof, 2018). Therefore, we have low confidence in the spatial representativeness of a measurement. The latter is implemented in the likelihood function by adding a tolerance to the borehole location, i.e., we select the modelled value within a radius r _bh which is closest to the measurement and use this in the likelihood. Then, the misfit function is stated as:

(10)$$\eqalignb{S_{\rm p_{\rm w}} &= \displaystyle{1 \over 2}\mathop \sum \limits_{i = 1}^n \left({\displaystyle{{{\rm min}\lpar {\vert {F_i\lpar {\bf m} \rpar -\;{\bf d}_{i\comma {\rm \gamma }}^{\rm p_{\rm w}} } \vert } \rpar } \over {\sigma_i^{{\rm p}_{\rm w}} }}} \right)^2\;\comma \;\;\; \cr &\quad {\rm \gamma }\in \lcub {{\rm \;}\forall \;x_{\rm \gamma }\comma \;y_{\rm \gamma }{\rm \;with\;}{\lpar {x_i-x_{\rm \gamma }} \rpar }^2 + {\lpar y_i-y_{\rm \gamma }\rpar }^2 \lt r_{{\rm bh}}^2 \;} \rcub \comma \;$$

where subscript γ refers to the neighbouring cells (coordinates x_γ, y_γ) within a distance r _bh of a borehole location (x_i,y_i). The standard deviation of borehole data are considered as $\sigma _i^{p_{\rm w}}$.

Tracer-transit times provide an estimate of subglacial water flow conditions. However, the total tracer resident times are strongly influenced by the inflow modulation time (time spent inside the moulin), which is estimated from water recharge into the moulin and its geometry. Within the data used in this study, there are no observations about recharge conditions into moulins. Consequently, observed transit times are only used to constrain or regularize the model into a range of physically plausible transit speeds (Schuler and others, 2004; Werder and others, 2010). We define the binary variable I^ts as an indicator of the tracer-transit speeds within a physically reasonable range defined by ${\bf d}^{{\rm ts}}_{{\rm min}}$ and ${\bf d}^{{\rm ts}}_{{\rm max}}$. Transit speeds outside the physically reasonable range are then penalized through a likelihood function with standard deviation $\sigma _i^{{\rm ts}}$:

(11)$$\eqalign{S_{{\rm ts}} = & \displaystyle{1 \over 2}\mathop \sum \limits_{i = 1}^n \left({{\bi I}^{{\rm ts}}\displaystyle{{{\rm min}\lpar {\vert {F_{\bi i}\lpar {\bf m} \rpar -\bf d^{{\rm ts}}_{{\rm min}} } \vert \comma \;\vert {F_{\bi i}\lpar {\bf m} \rpar -\bf d^{{\rm ts}}_{{\rm max}} } \vert } \rpar } \over {\sigma_i^{{\rm ts}} }}} \right)^2\;\comma \;\cr \;{\bi I}^{{\rm ts}} = & \left\{{\matrix{{0\comma \;\;\;{\bf d}^{{\rm ts}}_{{\rm min}} \le F_i\lpar {\bf m} \rpar \le {\bf d}^{{\rm ts}}_{{\rm max}} \comma \;}\hfill \cr {1\comma \;\;\;{\bf d}^{{\rm ts}}_{{\rm min}} \gt F_i\lpar {\bf m} \rpar \;{\rm or}\;{\bf d}^{{\rm ts}}_{{\rm max}} \lt F_i\lpar {\bf m} \rpar .\;\;} \cr } } \right.} $$

Additionally, we assume that it is unlikely that the hydraulic head exceeds the surface elevation. To this end, we implement an equation that is activated only if water pressure exceeds surface elevation and penalizes the likelihood when the water level is higher than the surface elevation. This is implemented as a truncated Gaussian given by

(12)$$S_{\rm z} = \displaystyle{1 \over 2}\mathop \sum \limits_{i = 1}^n \left({{\bi I}^{\rm z}\displaystyle{{F_{\bi i}\lpar {\bf m} \rpar -{\bi z}_i} \over {\sigma_i^{\rm z} }}} \right)^2\;\comma \;\;\;{\bi I}^{\rm z} = \left\{{\matrix{ {0\comma \;\;\;F_i\lpar {\bf m} \rpar \le z_i\;} \cr {1\comma \;\;\;F_i\lpar {\bf m} \rpar \gt z_i\comma \;} \cr } } \right.$$

where I^z is an indicator corresponding to the index of nodes where F_i(m) has exceeded the surface elevation z_i with a standard deviation $\sigma _i^{\rm z}$.

One of the challenges in the model inversion is that we expect different channel networks with the correct hydraulic parameters to have a similar hydraulic response (i.e. to honour the observations). In other words, we expect to find several local maxima, which can be distant in the model space. In addition, there are nonlinearities and discontinuities in the likelihood function. For example, gradually modifying structural parameters (l_xy, s, f) can modify the channel network until a large reorganization is introduced. This is, for instance, the case when two parallel channels are being gradually pushed together until they merge, culminating in a larger restructuring of the network.

Consequently, we choose an inversion algorithm that explores the space thoroughly and avoids getting trapped in one local maximum. This is commonly done by Markov chain Monte-Carlo algorithms such as the Metropolis-Hasting algorithm (Metropolis and others, 1953). Here, we use a more advanced algorithm, the Markov chain Monte-Carlo DREAM_(ZS) algorithm (Laloy and Vrugt, 2012). DREAM_(ZS) explores the model space similarly to a Metropolis-Hasting algorithm but includes multiple chains an adaptive sampling which facilitates a robust exploration of the model space and therefore it is able to find several local maxima more readily. In addition, we configure DREAM_(ZS) with different starting points to ensure an extensive model space exploration. As a result, DREAM_(ZS) captures the multimodal distribution of our problem, where each mode corresponds to a local maximum and therefore to a CS network. Moreover, it allows uncertainty quantifications of the model parameters for each explored local maximum. Note that this involves running the inversion procedure numerous times, which is feasible as our model is relatively computationally inexpensive.

2.4 Channel network structure persistence across the melt season

The proposed methodology finds sets of subglacial systems (or model parameters) for one snapshot. We refer to the set of inferred subglacial systems as CS systems, as they honour the observations. While this is acceptable for a single snapshot, if independent inversions are carried out for different snapshots, the spatial and topological dependence of the channel network structure across the melt season is not guaranteed. Similar to previous studies (Arnold and others, 1998; Banwell and others, 2013), we assume that the channel network structure does not change radically in the course of the melt season. While, changes in recharge during the course of the season can affect the hydraulic parameters and channel radii.

To ensure channel-architecture dependence between snapshots across the melt season, we force the parameters that define structural parameters of the subglacial channel generator (l_xy, s and f) to remain within a vicinity. For this, we implement the inversion in two steps: First, we infer a set of subglacial systems conditioned to a reference snapshot (labelled with subscript R: a _R, b _R, l_{xy, R}, s _R, f _R, T _d,R, τ_{i, R}), select manually one model (e.g., maximum likelihood mode) and extract the structural parameters (labelled with prime: l_xy’, s’, f’). The second step consists in inferring independently the subglacial systems of all other snapshots sε{1, …, 5}, limiting the exploration of the structural parameters to the vicinity of the selected model (l_xy’, s’, f’). The bounds are set in order to allow perturbations such as channel translation and local deformation to avoid large network reorganizations. This is done by sensitivity tests that constrain the bounds of l_xy’, s’ and f’ such that no large reorganizations are observed in the network. As a result, the inferred structural parameters for each snapshot (l_,xy,s’, s_s’, f_s’) are near the selected model. The inferred structural parameters for all snapshots are therefore correlated because they are forced to be in the surroundings of the same local maximum of the likelihood function, imposing a dependence across the melt season. Note that the channel radius parameters (a and b), transmissivity T _d of the distributed system and τ_i are not constrained to remain in the vicinity of the selected model (Fig. 2).

Fig. 2. Parameter inversion strategy and snapshot dependence. First, the parameters for a reference snapshot are inferred and one model is selected (e.g., maximum likelihood model). For the following snapshots, the structural parameters are constrained to the vicinity of the structural parameters from the selected model.

3. Application to the trunk of the Gorner Glacier

3.1. Field site and dataset

The study site and model area correspond to the trunk of Gorner Glacier, Switzerland (Fig. 3). Gorner Glacier is located in the Valais Alps at the foot of Monte Rosa and has experienced substantial mass loss since the 1980s (Huss and others, 2012). The glacier system consists of five tributaries covering an area of ~50 km² (in 2003). The model extent comprises the ablation zone of the main glacier tongue, with an area of ~6.5 km² and a surface elevation range between 2200 and 2700 m a.s.l.. Recent observations indicate surface velocities of up to 200 m a⁻¹ near the confluence of the tributaries Grenz and Zwillings Glacier, 35 m a⁻¹ near bh1 (Fig. 3) but much slower surface movement over much of the lower part of the glacier (Benoit and others, 2019). In addition, bedrock elevation maps show an extended over-deepening where ice thickness exceeds 400 m (Sugiyama and others, 2008). A distinct feature of the Gorner Glacier is the polythermal composition in the ablation zone. A cold ice body extends between the Grenz and Gorner confluence (near boreholes bh1–7 expecting bh5 in Fig. 3a), occupying 80–90% of the ice thickness (Ryser and others, 2013). The cold ice is impermeable to meltwater, except for fracturing or water flowing permanently through cracks and moulins, leading to supraglacial streams and ponds.

Fig. 3. (a) Model configuration, the location of boreholes (bh), and moulins (M) with tracer test data. (b) Overview of the Gorner Glacier system in 2005. The investigated area is shown in grey. Abbreviations correspond to the tributary glaciers Unterer Theodul (The), Breithorn (Bre), Schwärze (Sch), Zwillings (Zw), Grenz (Gre), and Gorner (Gor) Glacier.

Gorner Glacier was extensively studied over the years 2004–08 with the main focus in glacial lake outburst floods, which occurred periodically from the drainage of the Gornersee. Extensive fieldwork and data collection were done providing one of the most complete of such datasets (e.g., Huss and others, 2007; Sugiyama and others, 2008; Walter and others, 2008; Werder and others, 2009; Riesen and others, 2010). For this study, we use datasets acquired in 2005 which include borehole water pressure data (Huss and others, 2007; Riesen and others, 2010), dye-tracer experiments (Werder, 2009; Werder and others, 2009) and discharge of the catchment. In addition, hourly surface melt on a 25 m grid is available from a melt model calibrated to in situ ablation measurements and catchment discharge (Huss and others, 2008). Estimates of ice thickness distribution are available for the entire glacier system based on ground-based and airborne ground-penetrating radar surveys and spatial interpolation of the radar profiles (Farinotti and others, 2009; Rutishauser and others, 2016). Furthermore, repeated DEMs providing information on ice surface topography are available (Bauder and others, 2007). Key aspects of the datasets are summarized below and illustrated in Figures 3, 4.

Fig. 4. (a) Observed borehole water levels, measured outlet discharge and modelled discharge. Tracer-transit times correspond to vertical green bands defined by the injection time and the time of exit at the outlet (for reference, thicker line is ~13 h and thinner range from 2 to 4 h). Each modelled snapshot is denoted by black dashed line. (b) Zoom for each snapshot.

Water levels are available from seven boreholes drilled down to the glacier bed and instrumented with pressure transducers (Riesen and others, 2010). Boreholes bh1, bh3, bh4 show strong diurnal amplitudes of water head (sometimes in excess of 100 m) almost in phase with the discharge measured at the gauging station (Fig. 3, Fig. 4). Conversely, boreholes bh5, bh6 and bh7 show minor or no diurnal oscillations across the season, with a relatively higher water level. Note that at the beginning of August 2005, mean water levels rose and oscillations reduced significantly in boreholes bh1 and bh4. Borehole bh2 shows an anticorrelated behaviour with an amplitude of 3 m. Such an anticorrelated signal has been attributed to stress re-distribution at the glacier bed (e.g., Murray and Clarke, 1995; Lefeuvre and others, 2018; Rada and Schoof, 2018). Overall, boreholes with a steady and relatively high-water level can be interpreted as situated in partly or fully disconnected areas, whereas boreholes with strong diurnal-oscillating water level are situated in hydraulically well-connected areas. In addition, the rise in the water level of bh1 and bh4 in August 2005 can be interpreted as a transition from a channelized to a distributed or less efficient system in their vicinity.

Discharge for the entire catchment is available from the hydropower gauging and water catchment station, located ~1 km downstream the glacier outlet by the time of the study (2005). The discharge time series shows strong daily fluctuations and seasonal variations (Fig. 4). Note that proglacial discharge is in phase with the subglacial water pressure and thus we take this as the recharge rate, in which our model directly enters the subglacial system. However, the recharge needs to be spatially partitioned between the moulins on the glacier trunk (which is the model domain) and between the tributaries (input via inflow boundary conditions). We do this partitioning according to the melt rates calculated by a calibrated, distributed glacier surface melt model (Huss and others, 2008). According to the results, the trunk of Gorner Glacier (i.e. the area investigated in this study) contributes ~22% of the total discharge during the summer months, whereas the remaining catchment runoff stems from higher regions of the catchment and the tributary glaciers.

Dye-tracer tests were carried out throughout the melt season from four different moulins (m1–m4, Fig. 3). Tracer-transit times were extracted from the breakthrough curves of dye-tracer tests and corrected for proglacial residence time following Werder and others (2009). Overall, tracer-transit times for most of the injection points across the season ranged between 2 and 4 hour. The only exception corresponds to moulin m2, the most distant to the glacier outlet, which showed a transit time of ~13 hour in June and September 2005. Tracer-transit times are presented in Figure 4 as green vertical bands, defined by the injection time and the time of exit at the glacier outlet. Water flux entering the moulins was not measured during the tracer experiments. In addition, a mapping of moulins location was carried out in the field in 2005. We expect that several moulins may have been missed by the survey, but it is likely that their spatial density is well represented.

In addition, the seasonal drainage of a large ice-marginal lake impacts on the evolution of the subglacial system (Fig. 3). Note that lake drainage was observed between 13 and 15 June 2005. It is suggested that the lake drained via an R-channel which could have started opening a week before the event (Huss and others, 2007). However, it is estimated that before the outburst flood (13 June), the lake contribution prior to the main phase of the drainage was stored sub- or englacially and is therefore not accounted for in our model.

From the dataset, we select five snapshots across the melt season. We chose the timing of each snapshot such that it is at the end of a relatively stationary period, in order to make it representative of a stable situation of the subglacial system (Fig. 4 and Supplementary Table S1).

3.2 Modelling choices and settings

Total water recharge in the subglacial system is set to the discharge measured at the outlet for each snapshot (Fig. 4). The relative contribution of each tributary is extracted from the distributed surface melt model results and input from tributary glaciers as Neumann boundary conditions. In addition, basal melt, which is negligible compared to the total recharge, is assumed to correspond to a uniform 5 mm a⁻¹, following the same order of magnitude as presented in Herman and others (2011) and Cuffey and Paterson (2010). Water exits the subglacial system at the glacier outlet, which is modelled as a Dirichlet boundary condition at atmospheric pressure.

As precise moulin locations and water recharge are unknown, sensitivity tests are performed to evaluate their influence in the subglacial channel generator. Here, we carry out a sensitivity test evaluating the unconditioned networks as it allows inspecting the channel networks before running the more computationally expensive water flow model.

First, the subglacial channel generator is used to produce multiple realizations of unconditioned subglacial systems by randomly sampling structural parameters (l_xy, s, f). In addition, we explore the influence of the variance of the GRF on the network structure (Eqn (2)) and set it similar to small-scale changes in the hydraulic potential because the variance should not modify the general trend given by the other terms in Eqn (2). This step allows defining the structural parameter ranges where networks are realistic.

Second, we construct a moulin probability density map from a survey carried out in 2005 (Loye, 2006). Third, we randomly sample moulin locations from the density map (Supplementary Fig. S1) and assign a water recharge contribution proportionally to the Voronoi cell area of each moulin. Several random realizations including different number of moulins were used to determine the consequences of such assumptions (Supplementary Figs S2–S5). To compare the different channel network realizations (e.g., Fig. 5a–d), we compute the likelihood for each cell to be occupied by a channel weighted by the accumulated water flow (e.g., Fig. 5e, f), obtained by routing the hydraulic potential (Eqn (2)). Note that the normalized accumulated flow differs from the number of upstream-cells as is commonly found in the literature (e.g., Arnold and others, 1998). Here, the accumulated flow is weighted by the pointwise recharge from the moulins and tributary glaciers. Overall, sensitivity tests found that moulin number, location and relative recharge does not greatly impact the channel network structure (see Supplementary material). This can be explained by the fact that in this case study, water recharge into moulins within the model domain only accounts for 22% of the total subglacial discharge. Moreover, unconditioned channel locations for all configurations are well constrained by ϕ* and result in one or two main channels receiving most of the accumulated flow. For the following analysis and inversion, we thus run the model using a single set of 60 moulins whose positions are randomly sampled from the moulin probability map. A selection of unconditioned channel network realizations is presented in Figure 5.

Fig. 5. (a–d) Examples of unconditioned channel networks. (e) Colour scale shows the normalized sum of the accumulated flow for 500 realizations. (f) Colour scale shows the logarithm base 10 of the normalized sum of accumulated flow for 500 realizations.

Next, the parameter inversion strategy requires defining the misfit function of each component (Eqns (10–12)). Here, two sources of errors related to borehole observations are considered: borehole location (or spatial representativeness of a punctual measurement in a grid) and water pressure. As model errors are unknown and are usually larger than the instrumental error, we considered borehole water level standard deviations ${\sigma _i}^{{\rm p_{\rm{w}}}} = 10\,{\rm{m}}$, and spatial tolerance of r _bh = 100 m (Eqn (10)).

Tracer-transit times are largely dependent on water recharge into moulins and their geometry (e.g., Werder and others, 2010). As these data are not available, tracer-transit times are only considered as a measure of physically reasonable velocities in the subglacial system. In general, tracer experiments carried out in Alpine glaciers have shown that transit speed (straight path from the injection point to the glacier's snout divided by transit time) rarely exceeds 1 m s⁻¹ (e.g., Nienow and others, 1996; Schuler and others, 2004; Werder and others, 2010; Chandler and others, 2013). Indeed, tracer tests carried out in 2005 from four different moulins (Fig. 3) resulted in transit speeds in the range from 0.1 to 1 ms⁻¹ (Werder and others, 2010). Thus, we chose to bound tracer-transit speeds at the moulins with observations to a range of ${\bf d}^{{\rm ts}}_{{\rm min}}$ and ${\bf d}^{{\rm ts}}_{{\rm max}}$ (Eqn (11)) corresponding to transit speeds of 0.1 and 1.2 m s⁻¹. Note that transit speed is computed from the total residence time, including inflow modulation time. Water recharge into the moulin as well as moulin volume is unknown, hence inflow modulation time is unknown and thus we incorporate a fitting parameter that accounts for inflow modulation time (τ_i), with a range from 0 to 20 h. This is critical in order to avoid model overfitting. For transit speeds that fall outside the physically reasonable range, we consider a standard deviation of $\sigma _i^{{\rm ts}}$ = 0.25 m s⁻¹ (Eqn (11)).

Water pressure exceeding surface elevation has only rarely been observed on glaciers and has not been reported for Gorner Glacier. Here we decrease the likelihood of a model when the hydraulic potential exceeds the surface elevation. This is done by evaluating Eqn (12) at 108 points dispersed at regular 250 m intervals over the entire surface of the glacier and using a standard deviation of $\sigma _i^{\rm z}$ = 10 m. The value of ${\sigma _i}^{\rm{z}}$ has been chosen to rapidly decrease the likelihood of water pressures exceeding this bound. However, we note that this term plays a secondary role by driving the model towards a solution closer to the actual borehole observations.

Lastly, we select snapshot dates across the season as well as a reference snapshot. The snapshots are defined as the average condition at midday within a 15 minutes window. The main criteria to select snapshot dates is to find periods where water pressure, discharge and meteorological variables are as stationary as possible. Then the snapshot is considered at the end of the stationary period. It is likely that the channel network is well developed in midsummer. We thus select snapshot 3 as the reference snapshot. However, tests (not shown) indicate that the choice of the reference snapshot does not greatly influence the results. As shown in Figure 2, each snapshot is an independent inversion where the structural parameters are bounded in the vicinity of the selected CS system structural parameters (l_xy, f, s). The bounds considered are s ± 50 m, l_xy ± 50 m and f ± 0.1. A summary of the data used for the inversion is presented in Table 1.

Table 1. Model coefficients and parameters

^a Corresponds to a variable of the subglacial channel generator as well.

^b Superscripts correspond to ^pw water pressure, ^ts transit speed and ^z hydraulic head exceeding glacier surface elevation.

4. Results

4.1 Reference snapshot inversion (snapshot 3)

A set of 15 independent inversions were run for the reference snapshot (snapshot 3). Each inversion is configured to start with a different set of parameters. From the 15 independent inversions, 11 different local maxima or CS systems were found (some independent inversions converged to the same local maxima). Note that for each of the inversions a likelihood distribution of the inferred parameters is calculated by our numerical scheme. The maximum likelihood CS systems are shown in Figure 6a (see Supplementary Fig. S7 for additional results). Note that all but one inversion had a similar maximum likelihood value, whereas the remaining inversion converged to significantly lower likelihood value (CS 9 in Fig. S7). Overall, it was found that water flow occurs predominantly in channels that run-in depressions of the hydraulic potential (e.g., CS 1–5) and areas between channels are subject to a relatively higher hydraulic potential. CS system CS 6 does not fully capture this pattern, as its southern channel is not located in a depression of the hydraulic potential (Fig. 6a). Note that CS 6 transmissivity is relatively large compared to other CS, implying that more flow occurs in the distributed system and consequently the hydraulic potential field is smoother in comparison to other CS (Fig. 6c). All CS systems exhibit the presence of high discharge channels passing near boreholes showing diurnal oscillations in water pressure (bh1, bh3 and bh4). Boreholes with a high and steady hydraulic head are located in areas distant from predominant channels. One of the most significant differences between the inversions is the point where the main channels merge: for example, CS 1 has two main parallel channels that join near the glacier outlet whereas CS 2 has only one predominant channel already 3.5 km in front of the glacier snout (a). Note that for CS 3 parameters a and b converged to a range that differs from the rest of the models. The a and b parameter values indicate that CS 3 presents smaller radii for channel starting points and larger radii near the outlet in comparison to other CS. However, note that the sum of the volume of all subglacial channels (TCV) in CS 3 is similar to CS 6, indicating that the channel networks of CS 3 and CS 6 share similar overall volumes. This difference indicates that the channel structures influence the hydraulic efficiency of the network and therefore particular channels radii are needed for conditioning the model to observations.

Fig. 6. Conditioned subglacial systems (CS) for six independent inversions of the reference snapshot (snapshot 3). CS 1 is highlighted in grey as it is propagated throughout the rest of the melt season. (a) Maximum likelihood of CS hydraulic potential field and channel discharge. Note that CS 1, 3 and 6 present two main parallel channels, whereas CS 2 and 5 present one dominant channel downstream of the boreholes. (b) Borehole water pressure residuals (boreholes in x-axis). Observations are displayed as blue vertical bands and black dot indicates the maximum likelihood model displayed in (a). Boxplots indicate the uncertainty with boxes corresponding to the 25^th and 75^th percentiles, respectively. The whiskers extend to the most extreme data points (~±2.7σ for normal distribution). Outliers are plotted as red cross. (c) Inferred parameter distributions (CS in x-axis). Black dots indicate the maximum likelihood model displayed in (a), boxplots the parameter residuals. The total volume of subglacial channels (TCV) is shown in the right-hand side in m³.

Borehole water pressure residuals show that most observations lie within the uncertainty bounds (Fig. 6b), except for bh3, which is difficult to match because it implies a sharp hydraulic gradient. Indeed, boreholes 2 and 3 are at a distance of 220 m, but show a dissimilar water pressure level and behaviour. The hydraulic parameters (a, b and T _d) of the CS systems are well constrained. However, these parameters differ from one conditioned model to another. The latter implies that each channel network architecture has a particular geostatistical radii relation (Eqn (4)) that best fits the observations. The transmissivity (T _d) of the maximum likelihood models ranges between 10^−7.9 and 10^−7.2 m² s⁻¹. Lastly, the structural parameters (l_xy, s, f) converged to a different value for each CS, showing that each network is different by construction (c).

4.2 Inversion of subglacial systems across the melt season

We selected CS 1, as the maximum likelihood CS system model to carry on the inversion across the melt season (snapshots 1–5). This choice is supported by its high likelihood and the fact that two other independent runs (CS 3 and Supplementary Fig. S8) converged to a similar network structure. Additionally, we performed the inversion across the melt season assuming channel network persistence of CS 2 (Supplementary Fig. S9).

Inversion results show that the maximum likelihood CS system of the five snapshots presented, as expected, minimum variations in the network architecture throughout the melt season (Fig. 7a, c). Similarly, the CS systems present channels in areas with lower hydraulic potential, separated by high-pressure areas of the distributed system.

Fig. 7. Conditioned subglacial system (CS 1) across the melt season (five snapshots). Grey bands indicate reference snapshot (snapshot 3). (a) Maximum likelihood hydraulic potential field and channel discharge. Discharge at the outlet is provided in parenthesis (b) Borehole water pressure residuals (boreholes in x-axis). In snapshot 1, boreholes 4 and 5 do not have data. Observations are displayed as blue vertical bands and black dot indicates the maximum likelihood model displayed in (a). Boxplots indicate the uncertainty with boxes corresponding to the 25^th and 75^th percentiles, respectively. The whiskers extend to the most extreme data points (~±2.7σ for normal distribution). Outliers are plotted as red cross. (c) Parameter distributions (CS in x-axis). Black dots indicate the maximum likelihood model displayed in (a), boxplots the parameter residuals. The total volume of subglacial channels (TCV) is shown in the right-hand side in m³.

Across all snapshots, a majority of borehole pressure observations are within the uncertainty bounds of the model (Fig. 7b). Note that snapshots 2, 3 and 4 present a similar magnitude in water recharge, however, the observed borehole water pressure differs for each snapshot (Fig. 7b). The channel radius parameters (a, b) are subject to relatively small differences, which allow honouring the observations of each snapshot (Fig. 7c). Snapshots 4 and 5 indicate similar pressure observations, however, the recharge forcing of snapshot 4 is more than triple that of snapshot 5. Consequently, parameter a (radii linear scaling) is smaller in snapshot 5. This is clearly indicated by the total volume of subglacial channel networks (TCV), where the TCV of snapshots 2, 3 and 4 is more than double of that for snapshot 5. The transmissivity of the distributed system does not vary significantly across the inferred snapshots. As the structural parameters have been bounded, l_xy, s and f do not show any major difference between the snapshots.

5. Discussion

In this study, we infer properties of the subglacial drainage system of Gorner Glacier based on fitting our model to observations at different snapshots across the melt season. Our approach focuses on building a parsimonious model that allows the exploration of multiple solutions to explain observations. In this section we explore the limitations, advantages and further work associated with our framework.

The proposed model makes the following necessary simplifications: (1) it assumes a spatially uniform transmissivity to represent the distributed system, (2) it produces the channel network structure by a flow routing method, (3) it uses a statistical relation for the channel radii, (4) it assumes constant recharge for the duration of a snapshot period, which is necessary to resolve water flow in steady state. Conservatively, we selected the snapshots after a relatively stable recharge period in order to be relatively close to steady-state conditions. These assumptions make the model computationally inexpensive, and as a result, its parameters can be inferred from data, along with uncertainty quantification.

An important challenge when applying our model regards the uncertainties in recharge conditions. We use a distributed melt model calibrated with measured discharge and in situ ablation measurements, to estimate how recharge is partitioned between different locations of the model domain. Tributaries account for 76% of the total recharge and moulins account for 22% of the recharge. The remaining 2% of the discharge is incorporated outside the model domain (between the glacier snout and the gauging station). To assess the impact of the subglacial channel generator on the number of moulins, location and relative contribution, we carried out a sensitivity analysis. It was found that these factors do not greatly influence the network structure (see Supplementary material). This result, however, cannot be generalized, as it is particular to this site where channels are well constrained by Shreve's hydraulic potential and moulin recharge corresponds to only about a fifth of the total recharge. One particularity of Gorner Glacier is its polythermal nature in the ablation zone (Ryser and others, 2013). The cold ice body acts as an impermeable barrier for meltwater except when water flows continuously into pre-existing moulins or cracks. Such moulins in the cold ice area are rare, hence, they drain a larger surface area providing sufficient meltwater to maintain them. This feature could imply a more stable supra- and englacial system, and consequently that recharge locations into the subglacial systems are persistent across the melt season.

Results and previous studies have stressed that flow routing models are sensitive to small-scale perturbations (e.g., Chu and others, 2016; Irarrazaval and others, 2019). This is particularly relevant when downscaling (or interpolating) bedrock elevation data to match the higher spatial resolution of a surface DEM, as flow routing becomes biased towards small-scale features from the high-resolution DEM. Moreover, modelled bedrock elevation can significantly differ depending on the model and available data used to determine the ice thickness (e.g., Farinotti and others, 2017). The subglacial system generator adds variability, which allows generating sets of unconditioned subglacial systems. The core of the channel network generation is based on the addition of a nonuniform spatial perturbation term, modelled as a GRF, to the hydraulic potential field. The GRF is parametrized by structural parameters l_xy,, f and s, allowing for a low dimensionality of the inverse problem. As a result, modifying structural parameters l_xy, f and s results in global changes in the network. A drawback of this approach is that fitting independently individual water pressure observations requires locally modifying a channel, which is not possible with only three global parameters. To avoid overfitting, we add a spatial tolerance to the borehole locations (Eqn (8)), allowing for a broader range of networks to match observations. Although the choice of adding a GRF to Shreve's equation satisfied the objectives of this case study, further work should investigate the variability incorporated with different parametrizations or directly using a nonspatially uniform flotation factor, for example. Another drawback is that the exploration of the model space is challenged by discontinuities in the likelihood function as a consequence of modifying parameters (l_xy, s, f). This points out the importance of having a fast forward model which enables a large number of runs necessary for uncertainty quantification.

We designed a two-step inversion strategy that first allows inferring multiple subglacial systems models of a reference snapshot. Here we selected snapshot 3 (10 August 2005), i.e. a midsummer date, as it is more likely that channel networks are well developed. For this date, we obtained different CS systems which honour observations, demonstrating that data are insufficient to constrain a well-defined or unique channel network. Secondly, we selected two CS systems (CS 1 maximum likelihood model, and CS 2 in Supplementary material) to force the results of following snapshots across the melt season to resemble their conditioned system. The assumption that the channel network's main structure persists across the melt season is practical, as it avoids inversion parameters that introduce discontinuities in the likelihood function. As a consequence, convergence is much more easily reached in the subsequent snapshots. Moreover, the hydraulic parameter variability found within the reference snapshot is large (see Fig. 6) and thus considering one channel network allows for a better understanding of changes on the hydraulic parameters (enlargement and shirking of channel radii) across the melt season. However, this strategy certainly underestimates uncertainty as it only explores few models across the melt season. Further work, with more computational resources, should implement inversion frameworks to infer independently each snapshot (or all snapshot simultaneously), as well as not enforcing channel network persistence across the melt season. In that case, channel network persistence across the melt season will rely on data and its uncertainties

Subglacial transit times from four moulins were used to constrain the channel network parameters to a physically reasonable range. A number of assumptions, which were made to avoid model overfitting, needed to be made as the available tracer observations only loosely constrain the transit times: using a large range for the inflow modulation time (0–20 hour), as well as providing extended bounds for the subglacial transit speeds through channels (0.1–1.2 m s⁻¹). Nonetheless, including even these broad constraints on subglacial transit and inflow modulation facilitated model space exploration and convergence. Experiments performed without these constraints (not shown) converged to subglacial models composed of small channels dominated by large transmissivities in the distributed system. This configuration resulted in unrealistically fast water transit speeds and a smooth hydraulic potential field that did not capture sharp hydraulic potential gradients observed in nearby boreholes. Clearly, if a more accurate solution (which we found) exists in the transit-time constrained space, it also exists in the unconstrained space, but it appears that the capability of the DREAM_(zs) algorithm to find it is reduced. Thus, in our case, limiting transit speeds facilitated convergence to local maxima and forced discharge into channels. However, other possible strategies could have been considered such as adding a lower limit for channel radii, lowering the upper bound of T _d or utilizing another algorithm to explore the model space.

The high sensitivity of the inversion to subglacial transit time constraints also highlights that having more accurate observations available would constrain the inversion significantly better. This finding is coherent with previous work of Irarrazaval and others (2019) on a synthetic inversion, where tracer tests considerably reduced the uncertainty of the model parameters. This means that future tracer experiments need to separate total transit time into inflow modulation time and subglacial transit time, which necessitates observations of moulin recharge and obtaining constraints of moulin geometry (e.g., Nienow and others, 1996; Schuler and others, 2004; Werder and others, 2010).

Conditional subglacial systems show that water is concentrated in subglacial channels with low water pressure surrounded by areas of poor connectivity and high-water pressure. The subglacial channels were enlarged at the beginning of the melt season (from snapshot 1 to snapshot 2), maintaining their total volume over the summer (from snapshot 2 to snapshot 4) then substantially shrink at the end of the season (from snapshot 4 to snapshot 5) as shown in Figure 7c. This is compatible with current understanding of the seasonal evolution of the glacier drainage system (e.g., Fountain and Walder, 1998).

The result of our framework is a set of solutions for subglacial systems that can match the data with a high likelihood (Fig. 6). The set of inferred subglacial systems depends on the definition of the likelihood function, which includes errors from different sources (e.g., observation and model errors) that must be estimated. In our experiments, it occurred that the likelihood function accepted subglacial systems presenting a good agreement with data but a poor representation of the expected physics, such as channels not capturing water from the surroundings (e.g., CS 6). Further restrictions, for example, limiting the minimum radii (i.e. a parameter) and regularization may overcome this issue. Note that the variability within the set of CS system of snapshot 3 (Fig. 6) is larger compared to the variability of CS 1 across the season (Fig. 7). This implies that not considering a multi-solution approach could strongly underestimate the possible diversity of network configurations.

6. Conclusions

We presented a framework which enables the inference of subglacial channels from water pressure and tracer-transit times at different snapshots across the melt season. First, we introduced an updated subglacial channel generator, which produces a broad spectrum of likely channel networks. This is relevant as flow routing algorithms are sensitive to small-scale features and uncertainties in bed elevation datasets can be significant. Previous studies found that varying the spatially homogeneous flotation factor in Shreve's equation can divert subglacial flow routing paths (e.g., Banwell and others, 2013; Chu and others, 2016). Here we added a GRF that accounts for a spatially nonuniform term. The choice of a GRF is practical in an inversion framework, as it reduces the dimensionality of an unknown 2-D spatial field. Consequently, the incorporation of a GRF adds variability to the simulated subglacial systems, which are typically missing in flow routing models solely based on Shreve's equation.

The inversion framework allowed conditioning subglacial systems to observed water pressure and to a range of plausible tracer-transit speeds. Our main finding is that multiple channel networks can match observations within their uncertainties. Gorner Glacier presents one of the most comprehensive datasets, but data are still too sparse to result in one single solution. This highlights the importance of working in a multi-solution framework, particularly when using results for prediction as distinct channel networks could influence, for instance, basal sliding or respond differently under a glacial lake outburst flood scenario.

Benefit of the use of tracer experiment data for uncertainty reduction was limited. To avoid overfitting, tracer experiments from four moulins were used to constrain modelled water transit speeds to a physically reasonable range. For further uncertainty reduction, as pointed out in previous studies, it is recommended to measure the moulin recharge hydrograph during the tracer experiment. Last, our approach favoured a physically simple model for subglacial hydrology. Further frameworks should attempt finding a balance between adding more complex physics (and more parameters) versus physically simpler models capable of multi-solution, uncertainty quantification and data integration.