A subglacial hydrological model dedicated to glacier sliding

Introduction Conclusions References


Introduction
The flow of glaciers is a combination of viscous ice deformation and subglacial phenomena such as basal sliding and/or deformation of the sediment layer if it exists.The sliding component is particularly important for temperate glaciers, and can account for up to 90 % of the total surface speed (Cuffey and Paterson, 2010).As shown Introduction

Conclusions References
Tables Figures

Back Close
Full by earlier theoretical considerations (Weertman, 1957;Lliboutry, 1968;Schoof, 2005), water pressure is the key variable explaining most of the modulation of basal sliding.
For high water pressure, the strength of the bed resistance is reduced, either by cavitation in the case of hard beds, or by dilatation of the sediment layer for soft beds, inducing an increase in sliding speed.Observations of surface velocity and basal water pressure at the same location on glaciers confirm the importance of water pressure in controlling the flow of glaciers (e.g., Walder, 1982;Mair et al., 2003;Anderson et al., 2004;Bartholomaus et al., 2008).The aim of our work is to develop a subglacial hydrological model which can simulate basal water pressure and couple it to an ice flow model.The general framework is therefore in line with work by e.g.Schoof (2010) and Pimentel and Flowers (2010), but with a different approach for the hydrological model.Nevertheless, as in these earlier works, supraglacial and englacial water systems are extremely simplified and our focus is on the subglacial system.The characterization of the subglacial system is rather ill-defined.Whatever the draining systems, they can be classified into two main groups: (i) inefficient draining systems, exhibiting high water pressure, such as water film (Weertman, 1972;Walder, 1982), linked cavities (Lliboutry, 1968;Walder, 1986) or diffusion in a sediment layer (Shoemaker and Leung, 1987) and (ii) efficient draining systems, exhibiting lower water pressure, such as ice-walled channels (R öthlisberger, 1972), channels opened in the bedrock (Nye, 1976) or at the interface between ice and sediment (Walder and Fowler, 1994).
Drainage systems under glaciers are a combination of these inefficient and efficient systems (Shoemaker and Leung, 1987;Boulton et al., 2007).These two types of system, as a consequence of their efficiency to drain water, have different impacts on water pressure and consequently on basal sliding.Inefficient drainage systems are highly pressurized which results in relatively fast sliding speeds, whereas efficient draining systems allow water to drain at lower pressures.These two types of system, where they coexist, are tightly coupled and the efficient drainage system will tend to drain Introduction

Conclusions References
Tables Figures

Back Close
Full water out of the inefficient one which in turn induces a decrease in the basal velocities (Bj örnsson, 1974;Magn ússon et al., 2010).
Recently published subglacial hydrological models take into account inefficient and efficient components for the drainage system (Pimentel and Flowers, 2010;Schoof et al., 2012;Hewitt et al., 2012;Werder et al., 2013).In our approach, a sediment layer is used to model the inefficient drainage system (IDS), and, rather than actually modelling a network of channels to represent the efficient drainage system, we use an equivalent porous layer (EPL).This approach, known as the dual continuum porous equivalent approach in hydrogeology has been developed for karstified aquifers (Teutsch and Sauter, 1991).Karstified and glaciological hydrological systems share some distinctive features that motivates this approach.Specifically, they both consist of systems with an inefficient drainage component and a more efficient one which is activated only under some water head conditions (Hubbard and Nienow, 1997;White, 2003).
The model is tuned and validated by performing a series of three experiments of increasing complexity using data obtained on the Haut Glacier d'Arolla.Data-sets containing both hydrological and ice flow observations are rare.The Haut Glacier d'Arolla data-set is one of the most complete, although it suffers from some non-synchronous measurements.
The present paper describes the double continuum approach and the numerical methods that are used for its treatment in Sect. 2. The ice flow model equations and boundary conditions are introduced in Sect.3. Section 4 presents the simulation settings and the results leading to the selection of hydrological parameters.Finally, the coupling between the hydrological and ice dynamic model is presented in Sect. 5.

Hydrological model
The basal drag of glacier is strongly modulated by the effective pressure (e.g.Schoof, 2005), i.e. the difference between the overburden ice normal stress σ nn and the water Introduction

Conclusions References
Tables Figures

Back Close
Full (1) High water pressure induces low effective pressure, and where N = 0 the ice is locally floating.The effective pressure is the key variable for the coupling between glacier sliding and the subglacial hydrological system.In Eq. ( 1), the basal water pressure p w is positive for compression whereas the normal Cauchy stress σ nn is negative for compression and defined as: with n the upward pointing normal vector to the bedrock and σ the Cauchy stress tensor.
Note that the definition of N in Eq. ( 1), using the normal Cauchy stress instead of the overburden hydrostatic ice pressure p, is more rigorous and fully accounts for the stress distribution at the base of the glacier.Moreover, solving the Stokes equations for the ice flow often results in a normal Cauchy stress at the base which differs slightly from the hydrostatic pressure, justifying the use of Eq. (1) to define the effective pressure.
Using the dual continuum porous equivalent approach, the inefficient and efficient drainage components are both modelled as sediment layers with the use of a specific activation scheme for the efficient drainage system.This approach, rather than pointing the position of individual channels at the glacier base will define, in a continuous manner, the location where the efficient drainage system is most likely to develop.This approach has the advantage of requiring a lower bedrock topography resolution than the discrete approaches describing each channel individually (Schoof et al., 2012;Hewitt et al., 2012;Werder et al., 2013).Furthermore, the use of a diffusion equation to compute the water head distribution in both systems reduces the computational cost.Introduction

Conclusions References
Tables Figures

Back Close
Full

Water distribution in a porous medium
Hereafter, the index j (subscript or superscript) may either refer to the IDS (j = i) or to the EPL (j = e) and the term porous medium is used to describe both systems.The two main assumptions of the model are that (i) the porous medium is always saturated with water and (ii) the aquifer is confined, assuming that the overlying glacier and underlying bedrock are impermeable.Considering these assumptions, mass conservation for both the porous medium and the water has to be considered (Bear, 1988).Mass conservation for the water is given as For the porous medium, mass conservation reads where U w and U j are the filtration velocities of the water and of the porous media, respectively.The filtration velocity corresponds to the velocity that the material would have should it use all the available section.This definition relies on the first experiments made by Darcy in 1856 where the velocity of the fluid and solid were computed as a flux divided by the surface of the sample.Densities of water and porous medium are ρ w and ρ j , respectively and ω j represents the porosity of the porous medium.Finally, q j represents an inflow (outflow if negative) by unit of volume that represents the transfer of water from the inefficient to the efficient drainage system and/or the recharge of the inefficient system (see Sect. 2.2 and following).Darcy's law in its classical form (Darcy, 1856) reads

Conclusions References
Tables Figures

Back Close
Full where k j is the tensor of intrinsic permeability of the porous media, µ w is the viscosity of water, g is the norm of the acceleration of gravity, p j is water pressure in the porous medium j and z j is the elevation at which the pressure is computed.In Eq. ( 5), U d is the mean macroscopic velocity of the fluid with respect to the mean macroscopic velocity of the solid.
It is a common assumption in hydrology to consider that water density shows very limited spatial variations and that the velocity of the solid is negligible with respect to that of the liquid (Bear, 1988).Introducing these assumptions, the combination of conservation Eqs. ( 3) and ( 4) with Darcy's law Eq.( 5) gives the diffusion equation for a confined aquifer as follows In Eq. ( 6), the water pressure p j is expressed in terms of water head h j , the altitude of the water free surface for an equivalent unconfined aquifer.On a glacier, h j would be the altitude of the water surface measured in a borehole connected to the subglacial draining system.With z j defined as the elevation of the observed point from a reference level (here the mean sea level), the relation between water head and water pressure is Equation ( 6) introduces the two main physical parameters for the porous media, namely the tensor of hydraulic conductivity K j and the specific storage coefficient S j s .These two parameters are defined as: and,

Conclusions References
Tables Figures

Back Close
Full where β s and α are the compressibility of the solid phase and the porous media, respectively, while β w is the compressibility of water.The compressibility β s defines the compressibility of the solid phase of the media (i.e., it can be assessed by a compression experiment on pure material), whereas the compressibility of the sediment, α, takes into account the compressibility due to the rearrangement of the grains.As it is usually done, the β s term is dropped from the expression of the specific storage coefficient since it is negligible relative to the water compressibility (β s ≈ 1/25 β w , Freeze and Cherry, 1979).Furthermore, Eq. ( 6) is vertically integrated so that the hydrological model is one dimension lower than the ice flow model and can be solved only over the bottom boundary of the ice flow model.Doing so, the problem simplifies to depend only upon the horizontal coordinates, which is consistent with the goal of simplicity of our approach.With the assumption of z-independent terms in Eq. ( 6), the integrated values reduce to where z tj and z bj are, respectively the altitudes of the top and base of the considered layer and e j its thickness.Using these expressions, Eq. ( 6) is rewritten in its vertically integrated form This last equation gives the water head at each point of the domain within a porous sediment layer under a given flux per unit of surface Q j , as a function of the layer transmitivity tensor T j and storage coefficient S j .The behaviour of both the inefficient and the efficient drainage systems are given by Eq. ( 11).In the following sections, the differences between the two systems are presented along with their coupling scheme.

Water routing through the inefficient drainage system
Darcy's law is commonly used in glaciology to express the draining of water through a sediment layer (e.g.Boulton and Dobbie, 1993;Fischer and Clarke, 2001).It describes inefficient drainage in connection with high water pressure (Walder and Fowler, 1994).
Consistent with this last assumption, the subglacial flow in a sediment layer has an important impact on glacier sliding.
In our approach, the distinctive feature of the IDS is that the water head h i is bounded by an upper limit h max such that the effective pressure at the bedrock stays larger than or equal to zero (N ≥ 0).Using the definition of the effective pressure, Eq. ( 1), and the definition of the water head, Eq. (7), the upper limit reads To conserve water volume, the volume of water contributing to a water head larger than h max is transferred to the efficient drainage system as a flux Q xs .The h max limitation on the IDS water head is imposed as a Dirichlet boundary condition on the system for the nodes where h i ≥ h max and the corresponding residual is used to compute Q xs (see Sect. 2.5 for numerical details).This leads to an iterative method similar to the one used by Zwinger et al. (2007) for the treatment of temperature fields in glaciers where the temperature is limited by the pressure melting point and the excess of energy used to melt ice into water.Introduction

Conclusions References
Tables Figures

Back Close
Full

Water drainage through the efficient drainage system
As stated before, the dual continuum porous equivalent approach adopted here relies on the use of equivalent porous medium to model the efficient drainage system.The modelling of an efficient drainage system by means of a physically inefficient system is not straightforward but we believe that this approach could lead to convincing results.
In this approach, the efficient drainage system is used as an incidental system whose goal is to drain the excess of water from the sediment layer.Keeping that in mind, places where the EPL is activated are more representative of zones where efficient drainage is likely to occur rather than of actual channel positions.The use of a diffusion equation, Eq. ( 11), to model the efficient drainage system requires the development of a special treatment to reproduce the specificities of this system.The physical parameters of the EPL (i.e.T e , S e ) are adjusted to account for the high hydraulic transmitivity and the low storage coefficient which characterize such efficient drainage systems (Hubbard and Nienow, 1997).Moreover, the resolution of the equation itself is subject to certain conditions.Indeed, activation of the EPL is not needed if the sediment layer alone can drain all the water produced.Therefore, the EPL is activated only where the water head in the sediment layer first reaches the maximum water head h max defined by Eq. ( 12), leading to two distinctive states for the EPL: (i) the EPL is not active where h i < h max (see Fig. 1a) or the (ii) the EPL is activated where h i = h max (see Fig. 1b and c).
The first of these two states could represent a winter configuration, when the amount of water driven to the base of the glacier is small enough to be solely drained by the sediment layer.Once h i reaches h max at some places the efficient drainage system is activated.However, the EPL passes through a transitional state before being able to drain water from the sediment layer leading to two sub-states for the activated EPL.The transitional state represents the time required for the efficient drainage system to extend enough to reach an infinite reservoir.In glaciology, an infinite reservoir can be a large subglacial lake, the snout of the glacier or the ocean.At these places, the Introduction

Conclusions References
Tables Figures

Back Close
Full water head is known, allowing a Dirichlet boundary condition to be imposed to the hydrological model.The transition between the two sub-states of the active EPL is illustrated in Fig. 2. Considering a glacier of horizontal domain Ω (the ice/bedrock interface on which the hydrological system is defined), the hydrological boundary condition are of Neumann type on the sides (Γ 1 ) and of Dirichlet type for any infinite reservoir (Γ 2 ), such as the snout.In its active transitional state, the EPL diffusion equation is solved on a domain ω with a zero flux boundary condition on all its boundaries (γ).The EPL becomes efficient when the boundary of its active domain ω reaches the Dirichlet boundary condition Γ 2 .Considering this, the transition between the two sub-states of the active EPL is defined on the domain as i. the EPL is active in a transitional state if γ Γ 2 = 0 (Fig. 2a), ii. the EPL is active in an efficient state if γ Γ 2 = 0 (Fig. 2b) The transitional state of the EPL represents a growing phase during which the water head in the EPL is maintained at a high level.These high water heads are due to the incoming water flux from the IDS which is not evacuated due to the zero flux boundary condition.The spreading of the EPL is controlled by the maximum water head h max .
Once the water head in the EPL h e reaches h max , the EPL is activated downstream and the volume of water is redistributed all over the active domain.
Once the EPL becomes active and efficient, its functioning is the same as that of the inefficient system.At this point, the only differences between the two systems are the source flux and the value of the physical parameters T j and S j .The method used for the estimation of the source flux Q t for the EPL is presented in the next section.So far, no closing procedure of the EPL has been implemented in the model.This feature would be required to perform simulations longer than one melt season, but can be ignored for the shorter simulations presented in this paper.Introduction

Conclusions References
Tables Figures

Back Close
Full

Coupling of the inefficient and efficient drainage systems
Once the EPL is activated, a transfer flux is established between the two different systems.This flux Q t is illustrated in Fig. 1b and c.The expression of the transfer term is a function of the water head in the two systems, of the characteristics of the inefficient drainage system (thickness e i and transmitivity T i ) and of the leakage factor ϕ, such that The leakage factor ϕ is a characteristic distance that the water has to cross to pass from one to the other drainage system.The introduction of the storage coefficient S j (j = i, e) is needed to convert the water head into volume.The value is dependent on the source of the water, thus if the water comes from the inefficient drainage system, S i is used, and if the water drains from the EPL to the sediment layer S e is used instead.

Numerical methods
The diffusion equation (Eq.11) is solved using the finite element method.The variational formulation is obtained by integrating over the domain Ω Eq. ( 11) and multiplying it by the test function φ, such that This equation is further transformed by applying Green's theorem to the second term, so that where Γ is the boundary surface of the domain Ω.Discretisation of this system finally leads to a formulation of the problem such that where H j is the vector solution, M a mass matrix, A j is the system matrix defined by the second term in Eq. ( 15) and B j is the force vector constituted by the two last terms.
A backward difference formula is then applied to discretizate the time derivative, with, the new matrix of the system and the new force vector.H p j and B p j are the solution and force vector at time step p and ∆t the time step.The treatment of this equation is the same for the two systems as long as h j remains lower than h max .The way the upper limit is imposed for each system requires two different treatments.
For the IDS, as stated in Sect.2.2, a Dirichlet method is applied to the water head to limit its height to the value h max .After the first iteration of the system, if any element h ip of the solution vector H i is greater than h max , then the system is manipulated such that: where the K i matrix is the same as the K i matrix except for the p th line which is fixed to zero apart from the value on the diagonal which is fixed to unity.Similarly, the force vector F i is equal to the F i vector except for its p th value which is fixed to h max .From this new system, a residual vector R i is computed such that: This is repeated until the relative change of H i falls below a given threshold.For the converged solution, the residual R i represents the necessary sink per node needed to keep the local water head below its maximum limit.Due to the assumption that all the water is drained into the effective layer, the residual R i is then added to the resolution of the EPL equation as follows; For the EPL, the volume of water above the given maximum limit is used to increase the size of the efficient drainage system.After the first iteration of the system, if any value h ep of the solution vector H e is greater than h max , the EPL domain (see Fig. 2) is increased in the direction of the lowest hydrological potential and the system is iterated until each element of H e satisfies the fixed upper limit.
The coupling of the two hydrological systems requires iteration between the two draining systems to achieve stability.This outer iteration loop is completed by two inner ones on each subsystem which ensures that the water head in each system does not overflow the maximum boundary h max .These inner loops are also needed to compute the water transfer between the two systems.A schematic view of this iterative process is presented in Fig. 3.

Ice flow model
The finite element code Elmer/Ice is used to solve both the hydrological and ice flow equations.The governing equations for the ice flow model are briefly summarised be-Introduction

Conclusions References
Tables Figures

Back Close
Full fore describing in detail the basal boundary condition which links the hydrological and ice flow models.Further information on the numerics and capabilities of Elmer/Ice can be found in Gagliardini et al. (2013)

Governing equations
The problem to be solved is the one of a gravity-driven flow of incompressible and non-linear viscous ice sliding over a rigid bedrock.
The ice rheology is given by Glen's flow law, defined as: where τ is the deviatoric stress tensor, εij = (u i, j + u j, i )/2 are the components of the strain-rate tensor, and u is the ice velocity vector.
The effective viscosity η in Eq. ( 23) is expressed as a function of the fluidity parameter A as where ε2 e = tr( ε2 )/2 is the square of the second invariant of the strain-rate tensor.Ice is assumed to be isothermal so that the fluidity parameter A is a uniform constant.
Moreover, the commonly used value for the exponent, n = 3, is adopted.In Eq. ( 26), ρ ice is the ice density, g the gravitational acceleration vector and the Cauchy stress tensor σ is linked to the deviatoric stress tensor such that σ = τ − pI, where p = −trσ/3 is the isotropic pressure.More details regarding the numerics of the ice flow model can be found in Gagliardini and Zwinger (2008).Solving for changes in the upper surface elevation z = z s (x, y, t) involves a local transport equation which reads where a is the accumulation/ablation function given as a vertical flux at the upper surface.Due to the duration of the simulation performed here (less than a year), we further assume a = 0.

Boundary conditions
Upper and lateral boundaries are treated as stress-free surfaces.The bedrock boundary is used to couple the ice dynamics with the subglacial hydrology, through the effective pressure N. In return, the mass redistribution derived by the ice flow model influences the hydrological model by modifying the Cauchy normal stress at the bedrock.The relation between the mean basal drag τ bi , basal velocity u bi and effective pressure N was first introduced by Lliboutry (1968).Recent studies from Schoof (2005) and Gagliardini et al. (2007) provide a friction law based on three parameters which depend only on the bedrock geometry.The proposed formulation fulfills the upper limit on the basal drag for a finite sliding velocity known as Iken's bound (Iken, 1981) and leads to a decrease in basal drag for low effective pressures and/or high velocities.This Coulomb-type friction law reads with In these relations, τ bi = t i • (σ n)| b is the basal shear stress in the tangential direction ) is the basal tangential velocity in direction i , with n being the upward-pointing normal to the bedrock surface.The parameters A s and n are the sliding parameter without cavities and the friction law exponent (n = 3 taken as in Glen's flow law), respectively.The exponent m controls the post-peak decrease in basal drag.Following Pimentel and Flowers (2010), we adopt the simplest form of the law setting m = 1, so that α = 1.Doing so, τ bi /N monotonously increases to its upper bound C which constitute with A s the last of the friction law parameters.

Field site and methods
Haut Glacier d'Arolla is an alpine glacier located in Switzerland (Fig. 4).This glacier is relatively small with a surface of 6.33 km 2 (Sharp et al., 1993) at altitudes ranging from 2560 to 3500 m a.s.l.The glacier is believed to be warm based and resting over unconsolidated sediments (Copland et al., 1996;Hubbard et al., 1995).The bed and surface of the glacier were mapped in 1989 by Sharp et al. (1993) and several surface DEMs have been created since then.In our study, we will use the 1989 bedrock DEM along with the surface elevation from 1993.
The main interest in this glacier for our study are the hydrological investigations that have been undertaken on it (e.g., Arnold et al., 1998;Gordon et al., 2001;Kulessa et al., 2003).These studies give a sound knowledge of the hydrological configuration in the area beneath the main tongue of the glacier, about 1.5 km from its snout and labeled borehole array in Fig. 4. A study by Hubbard et al. (1995)  1998) give a good insight to the evolution of the subglacial draining system during the melt season.

Strategy to estimate the hydrological parameters
The hydrological model has been designed to rely on a limited number of parameters.
As presented in Table 2, most of the hydrological parameters are well defined with the exception of the transmitivities of both layers T j (j = i, e), and the leakage factor ϕ.
We further assume that the transmitivity of both systems is isotropic and is therefore a scalar value T j .Estimating the parameter values for the hydrological model leads to two distinct problems.First, measurements of subglacial sediment transmitivity are rare and encompass a large range of values (Freeze and Cherry, 1979).Second, the use of an equivalent layer for the treatment of the efficient drainage system prevent us from directly using parameter values that would characterize a true discrete channel.
The EPL transmitivity and leakage factor have then to be estimated by comparing the model results directly to observations of the hydrological systems.We therefore adopted the following strategy.As first approximation, a broad range of values for the three unknown parameters is estimated from the available measurements.These ranges are large enough to produce very different model results.Then, these ranges of values are decreased by comparing the model results to large scale features of the hydrological system of Haut Glacier d'Arolla.
The sediment transmitivity is estimated using the hydraulic diffusivity at the bed of Haut Glacier d'Arolla measured by Hubbard et al. (1995).The hydraulic diffusivity D represents the velocity of a pressure pulse through the media and is given as: values, Eq. ( 30) and the physical parameters needed to evaluate S i (given in Table 2), the IDS transmitivity T i is estimated to range from 1.4 × 10 −4 to 8 × 10 −2 m 2 s −1 .The choice of the EPL transmitivity is more complicated.It cannot be directly measured as it represents the mean behaviour of a number of discrete channels.However, based on the previous measurements, the lowest value of the EPL transmitivity is set to 4 × 10 −3 m 2 s −1 , which corresponds to the value measured closest to the channel margin.The higher limit for the EPL transmitivity is then fixed at 8 × 10 −1 m 2 s −1 , ten times larger than the maximum value of the IDS transmitivity.This scaling is consistent with the results of Nienow et al. (1998) which describe the differences between the mean flow velocity of the distributed and of the channelized drainage systems.Finally, the leakage factor cannot be constrained by measurements and so, a large range of values, from 1 m up to 50 m, is adopted.Assuming these broad ranges for the three unknown hydrological parameters (see values in Table 3), two configurations of the hydrological system are then constructed and compared to measurements.The first configuration is characteristic of an end of winter system whereas the second reproduces the development of the draining system during summer.The comparison is done using two metrics constructed using large scale features of the hydrological system.
The first metric is defined as the maximum length of the active EPL, which represents the development of the efficient drainage system.The EPL length can be compared to the maximum channel length estimated by dye tracing experiments performed at different dates during the summer season (Nienow et al., 1998).The control observations are extracted from a number of dye tracing experiments that were undertaken during summers 1990, 1991 and 1995 (Nienow et al., 1998;Mair et al., 2002b).The earliest dye tracing measurements performed on Haut Glacier d'Arolla are dated around the 10 June and show channel lengths of slightly less than 700 m.At this time of the year, the discharge at the snout of the glacier is already ten times higher than the base winter discharge used for the forcing of our simulations and we would therefore expect to model a maximum length of the EPL substantially lower than the recorded 700 m.The Introduction

Conclusions References
Tables Figures

Back Close
Full other constraint given by these observations is the timing of the upglacier migration of the head of the channel system throughout the melt season.This dynamic aspect is compared to the evolution of the maximum EPL length during the transient summer simulation.
The second metric gives the mean water head of the IDS in the borehole array shown in Fig. 4. It gives a good overview of the general state of the coupled draining system and will be further referenced as to IDS water head.However, this last metric is more difficult to compare with actual data because such measurements show large spatial variability even between two boreholes spaced a few tens of meters apart (Willis et al., 2003).The IDS water head should then be regarded more as a metric used to compare different simulations than for a proper validation of the model against measurements.
Due to the difficulty of interpretation and comparison of these water head measurements, these data will not be presented in the current work.

End of winter configuration
The end of winter configuration is achieved by distributing the observed winter discharge at the snout of the glacier (5 × 10 −2 m 3 s −1 ) over the whole glacier surface, giving an input of ∼ 8 × 10 −9 m s −1 m −2 of water.This constant water flux is maintained until the water head of both the IDS and the EPL reach a steady state.Figure 5 presents the IDS water head on the whole glacier and the extent of the active EPL at the end of winter for three different values of the IDS transmitivity.The length of the EPL is depicted by the white thick line in Fig. 5. Comparison of steady state configurations indicates that an increasing IDS transmitivity leads to a decreasing water head and a less developed EPL.As expected, the large range of values for the IDS transmitivity leads to a large spread in the model results.For T i = 1.6×10 −4 m 2 s −1 the drainage system is overwhelmed by the EPL.Conversely, for T i = 1.6 × 10 −2 m 2 s −1 the drainage capacity of the IDS is such that it can drain all the input water and the development of the EPL is therefore not required.These two extreme cases indicate that the chosen range of IDS transmitivity covers all possible behaviour of the IDS.Nevertheless, Introduction

Conclusions References
Tables Figures

Back Close
Full excessively high values of the IDS transmitivity leads to unrealistic behaviour.Using the observed maximum channel extent at the beginning of spring (Nienow et al., 1998;Mair et al., 2002b), the range of the IDS transmitivity can be constrained.Thus, only transmitivity values that lead to an EPL maximum length lower than 700 m and larger than 200 m will be considered as admissible.Figure 6 shows the evolution of the EPL maximum length (upper panel) and of the IDS water head (lower panel) for IDS transmitivity ranging from 1.4 × 10 −4 to 8 × 10 −2 m 2 s −1 .As in Fig. 5, a decrease in IDS transmitivity leads to an increase in the EPL maximum length.The adopted restriction on the EPL length is represented by the grey zone on the upper panel.This constraint leads to a new range for the admissible values of the IDS transmitivity from 3 × 10 −4 to 3 × 10 −3 m 2 s −1 .The thickness of the lines in Fig. 6 represents the scattering of the metrics in response to a modification of the EPL transmitivity.The relatively small line thickness indicates that, in the case of the steady state configuration characterizing the end of winter, both draining systems are quite insensitive to the EPL transmitivity.
The IDS water head for its part increases with decreasing IDS transmitivity until the EPL extent is such that it can drain the borehole array (around 1200 m from the snout, as depicted in Fig. 4).For the configurations where the EPL reaches the borehole array, the water from the IDS can then be easily drained explaining the sharp decrease of the IDS water head.After reaching a minimum for T i = 1.5 × 10 −4 m 2 s −1 , the IDS water head increases again in response to the decrease in IDS transmitivity.The amplitude of this saw-tooth behaviour is a function of the drainage efficiency of the EPL but is not sensitive to the EPL transmitivity.It should then be driven by the leakage factor as discussed below.
Figure 7 shows the sensitivity experiments to the leakage factor ϕ. As explained in Sect.2.4, a large leakage factor implies low efficiency in the EPL and so a larger extent of the EPL and a higher water head in the IDS.The IDS water head is more sensitive to the leakage factor than is the EPL length metric.As shown in Fig. 7, the drop of IDS water head is amplified for smaller leakage factor.This amplitude variation is explained by the ability of the EPL to drain water from the IDS.The smaller the leakage factors, the Introduction

Conclusions References
Tables Figures

Back Close
Full easier the transfer of water from the IDS system to the EPL system.For unrealistically large leakage factors (i.e.ϕ ≥ 20 m), even if the EPL is activated the water head in the IDS returns quickly to the value it has before the opening of the EPL.The range for the leakage factor is therefore restricted from 1 to 20 m.In summary, the end of winter configuration allows the range of both the IDS transmitivity and the leakage factor (values given in the second column of Table 3) to be decrease, but not the EPL transmitivity.To continue, the second configuration corresponding to the development of the hydrological system during the summer is then used.

Transient summer configuration
The previous steady-state configuration corresponding to an end of winter is subsequently used as the initial state of the transient summer simulations.The transient response of the model is obtained by imposing a time dependent water flux at prescribed moulin positions.To this aim, we use the moulin positions recorded during the 1993 melt season and the associated influx modelled by Arnold et al. (1998) for the 1993 summer season.Each moulin is assumed to be perfectly vertical and is represented by a single node of the mesh.Unfortunately, dye tracing experiments were not performed during the 1993 season and the evolution of the channel drainage system is therefore compared against the 1990 measurements.Nevertheless, the comparison of the 1990, 1991 and 1995 melt seasons indicates that the changing extent of the hydrological system during the summer seasons develops at a similar rate and follows similar structures despite some variations in the spatial and temporal distribution of the water sources.Moreover, Fig. 10 shows that the moulin positions in 1993 are very similar to the 1990 positions.Comparing the model results to the 1990 melt season drainage system evolution is therefore a reasonable assumption, especially given the uncertainties in modelled moulin influxes.
Starting from the poorly developed EPL observed at the end of winter, a fully developed EPL is obtained by the end of the melt season.As for the end of winter Introduction

Conclusions References
Tables Figures

Back Close
Full configuration, the model results are strongly dependent on the IDS transmitivity value but rather insensitive to variations in EPL transmitivity.Because of this lack of sensitivity, only the simulation performed for T i = 7.9 × 10 −2 m 2 s −1 will be presented in the following.
Figure 8 shows for various IDS transmitivities the evolution during the summer of the EPL maximum length and IDS water head.In the range of applied transmitivity, an increase in IDS transmitivity leads to a delay in the spreading of the EPL during the melt season.This delay is induced by the postponing of the EPL opening induced by lower water head at the end of winter for large transmitivity, as can be shown in the lower panel of Fig. 8. From this sensitivity study, an IDS transmitivity of 1.6 × 10 3 m 2 s −1 seems to best reproduce the observed development of the drainage system during summer.The model's sensitivity to the leakage factor is presented in Fig. 9.The lower panel indicates that a higher leakage factor leads to a higher IDS water head because of the less efficient transfer from the IDS system to the EPL system.As shown on the upper panel, the EPL length metric show little sensitivity to varying the leakage factor except for very low leakage factor values (ϕ < 10 m).As for the end of winter configuration, the IDS water head is more sensitive to the leakage factor than is the EPL length.Gordon et al. (1998) have reported a decrease of the water head by 70 m at the opening of the EPL, in good agreement with the modelled drop for leakage factor ϕ ≥ 10 m.
Figure 10 compares the modelled EPL extend for ϕ = 10 and ϕ = 20 m to the reconstructed channel system at the end of the 1990 melt season (Sharp et al., 1993).This comparison shows a good agreement between the extent of the modelled active EPL and the channel system observed at the end of summer for both values of ϕ.Therefore, even if the two modelled EPL extents show some differences, they are too similar to identify an optimum leakage factor value, bearing in mind that the observed channel system is itself a reconstruction from dye tracing measurments.In the following, the value ϕ = 10 m is therefore adopted.Introduction

Conclusions References
Tables Figures

Back Close
Full  and 9 for the adopted values of the three hydrological parameters.The simulation is compared to the estimated length of the channelized drainage system (black dots) but also the up glacier extent of channelised drainage as interpreted by Nienow et al. (1998) (black line).Our results seem to indicate that the EPL extent is less smooth than the one proposed by Nienow et al. (1998) and it evolves by steps driven by the moulin positions and the bedrock topography which is consistent with the interpretation of Mair et al. (2002a) In summary, comparison between model results and observations for two distinct configurations of the hydrological system has allowed selection of the most reasonable values of the three hydrological parameters constrained by observations and independent interpretation.The adopted values are given in Table 3.This set of parameters is now used in the following section to model the coupling between the ice flow and the hydrological system throughout the melt season.

Modelling of spring speed-up events
Adopting the newly defined set of parameters, the ice flow and hydrological models are coupled with the aim of modelling the spring speed-up observed at Haut Glacier d'Arolla.Speed-up events were recorded during four melt seasons on Arolla glacier, in 1994 (Mair et al., 2001), 1995(Mair et al., 2002a), 1998and 1999(Mair et al., 2003).Again, unfortunately, no velocity measurements are available from the 1993 melt season for which we have water inputs (Arnold et al., 1998).However, from these various speed-up observations, we can characterise them as being short lived (three to four days) periods during which the surface velocities show a two to four-fold increase with respect to their average values.
The hydrological and ice flow models are coupled through the friction law (Eq.28) which depends on 4 parameters.As stated in Sect.3.2, the post-peak decrease exponent m and the friction law exponent n are, respectively fixed to 1 and 3.The latter is fixed at the usually accepted value of the Glen's exponent whereas the for-Introduction

Conclusions References
Tables Figures

Back Close
Full mer is chosen to achieve a better numerical stability.The sliding parameter in the absence of cavities is assumed to be uniform at the base of the glacier.The value A s = 1.6 × 10 −23 m Pa −3 s −1 is adopted to reproduce the observed winter velocity when the water pressure is low.In addition, the effective pressure is not computed for bed elevations above 3000 m and is fixed to N = 1.2 MPa.This limitation was necessary to stabilise the ice flow model on these parts of the glacier where bedrock slopes are high.Comparison between limited and not limited simulations shows that this upstream limitation does not impact the velocities on the lower part of the glacier on which we now focus.The fourth parameter in the friction law represents the maximum value reached by τ bi /N and should be smaller than the maximum value of the slope of the local obstacles m max .In Pimentel and Flowers (2010), the value C = 0.84m max for a sinusoidal roughness distribution was adopted, with m max = 0.5.Here we will test the model sensitivity to values of C from 0.5 to 0.9.
The results of the simulation performed with this set of parameters are presented in Fig. 12 for the reference point shown in Fig. 4. The effective pressure computed by the hydrological model controls the variability in the surface velocity throughout the melt season.A spring speed-up occurs between days 185 and 190 before the activation of the efficient system.The speed-up continues until the efficient system is activated and the water pressure drops.In the 12 days prior to this major event, two minor speed-ups are modelled.After day 205, the effective pressure decreases again in response to a heavy water input which triggers a new increase in the glacier speed with a peak around day 235.In comparison to the first one, this second speed-up period is characterized by higher daily variability.By this stage of the season, the subglacial hydrological system has reached its maximum capacity and any melt occurring during the day induces a quasi-instantaneous increase in water pressure, explaining this higher daily variability.After this second speed-up period, the glacier enters a quieter regime due to the relatively low water input during fall.Introduction

Conclusions References
Tables Figures

Back Close
Full Winter observations cannot be used to constrain the value of the parameter C which has a very negligible impact on the velocity when the water pressure is low.However, during spring, the model is quite sensitive to C due to higher water pressure.Figure 13 presents the evolution of the surface velocities obtained for different values of C during the spring speed-up event.The larger the value of C, the less marked the acceleration during the spring speed-up event.With C = 0.9, the speed-up is hardly distinguishable from the background velocity, whereas C = 0.5 gives the highest speed-up event with velocities that tend to stay higher in between the speed-up events.
The surface longitudinal velocity pattern for C = 0.5 during the first speed-up is presented in Fig. 14 for the lower part of the glacier.The evolution of the pattern during the spring speed-up matches the one that is observed on the glacier with a two-fold increase in the velocities during the speed up.For comparison purposes, the velocities measured during the first 1998 spring speed-up event (Mair et al., 2003) are superimposed on the modelled velocities.The comparison with this specific event shows that the modelled speed-up is less pronounced than the observed one and that the maximum speed is shifted downstream by approximately 400 m.Considering the various hypotheses in the model and the non synchronous data-sets used for velocity and water input, a complete agreement with the observations is not to be expected.

Conclusions
We have presented a new hydrological model especially designed to be coupled to an ice flow model.This hydrological model is based on a double continuum approach and solves the same set of equations using different parameters for both the inefficient and efficient drainage systems.The two systems are coupled so that the total amount of water is conserved and an ad-hoc scheme is proposed to activate the efficient drainage system where the water pressure exceeds the overburden ice pressure.In our approach, the channels are not represented individually but in a continuous manner, presenting the advantage of not requiring a very fine description of the basal topography.Introduction

Conclusions References
Tables Figures

Back Close
Full In this paper, the hydrological model and its coupling to an ice flow model are validated by performing a series of three applications of increasing complexity using the dataset of Haut Glacier d'Arolla.A first application aiming for a steady state configuration, corresponding to the winter state, is used to decrease the range of possible values for the three most poorly constrained parameters of the model, i.e. the transmitivity of both drainage systems and the leakage factor.In the second experiment, the evolution of the drainage system during the spring season is studied, starting from the previously obtained winter steady-state.Again, the model sensitivity to the three most poorly constrained parameters is tested.The third and last application couples the hydrological and ice flow models and results are compared to observed glacier speeds.Despite the use of non-synchronous data-sets and the number of simplifying assumptions in the model, good agreement is obtained both in terms of the temporal and spatial drainage  1. Definition of the different variables, constants and parameters in the model.As stated in Sect.2.1, the subscript j is used to refer to the porous media in general, and j = i for the sediment layer IDS and j = e for the EPL.

Conclusions References
Tables Figures

Back Close
Full Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | pressure p w : N = −σ nn − p w .
Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Ice flow is governed by the Stokes equations, that consist of the conservation of mass for incompressible fluids tr ε = divu = 0, (25) and the conservation of linear momentum divσ + ρ ice g = 0. (26Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | in the same area gives a range of values for the sediment hydraulic conductivity in the vicinity.Some other studies, involving dye tracing experiments (e.g., Mair et al., 2002a; Nienow et al., Discussion Paper | Discussion Paper | Discussion Paper |

)
Hubbard et al. (1995) measured hydraulic diffusivities ranging from 4 × 10 1 m 2 s −1 near an efficient drainage zone up to 7 × 10 −2 m 2 s −1 70 m away from this zone.From these Introduction Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Figure11is an enlargement of the upper panel of Figs.8 and 9for the adopted values Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | system evolution during the spring season and of the magnitude and duration of the observed speed-ups.The largest uncertainties in the model's velocity response are in the values of the leakage factor and the parameter C in the friction law.The other parameters in the model have proven to be relatively well constrained by observations or empirical interpretations (as for the IDS transmitivity or A s ), or the model has proven to be quite insensitive to them (as for the EPL transmitivity)Discussion Paper | Discussion Paper | Discussion Paper | Teutsch, G. and Sauter, M.: Groundwater modeling in karst terranes: scale effects, data acquisition and field validation, in: Proc.Third Conf.Hydrogeology, Ecology, Monitoring, and Management of Ground Water in Karst Terranes, Nashville, TN, 17-35, 1991Discussion Paper | Discussion Paper | Discussion Paper | Table

a
mass balance [m s −1 ] A fluidity parameter [Pa −n s −1 ] A s sliding parameter without cavity [m Pa −n s −1 ] C friction law maximum value e j thickness of the layer [m] g gravitational acceleration [m s −2 ] h j water head of the porous media [m] k j intrinsic permeability of the porous media [m 2 ] K j permeability of the porous media [m s −1 ] m friction law exponent n flow law exponent N effective pressure [Pa] p j water pressure in the media [Pa] q j volumic sink/source term [s −1 ] Q j water flux by unit of surface [m s −1 ] Q t water transfer between the two layers [m s water with respect to the porous media [m s −1 ] U j filtration velocity of porous media [m s −1 ] U w filtration velocity of water [m s −1 ] z vertical coordinate [m] α compressibility of the solid [Pa −1 ] β s compressibility of the sediment [Pa −1 ] β w compressibility of the water [Pa −1 ] ε strain-rate tensor [s −1 ] εe strain-rate invariant [s −1 ] η effective viscosity of ice [Pa s] µ w viscosity of water [Pa s] ρ w density of water [kg m −3 ] ρ ice density of ice [kg m −3 ] ρ j density of the porous media [kg m −3 ] σ Cauchy stress tensor [Pa] τ deviatoric stress tensor [Pa] τ bi mean basal drag [Pa] ϕ leakage factor [m] ω j porosity of the media Introduction Discussion Paper | Discussion Paper | Discussion Paper | Table 2. Values of the parameters used in the hydrological model along with their sources.Well known parameters are referenced as wkp and poorly known parameters are labeled pkp.Discussion Paper | Discussion Paper | Discussion Paper |

Fig. 1 .
Fig. 1.Description of the coupling between the two layers of the hydrological model.The top panels represent the water load in the IDS (solid line), in the EPL (dotted line) and the flotation limit (dashed line).The lower panels show the routing of water, (a) when the EPL is not active, (b) when the EPL is active in a transitional state (grey and white chess) and (c) when the EPL is active and efficient.

Fig. 2 .Fig. 3 .Fig. 4 .Fig. 4 .
Fig. 2. Description of the evolution of the boundaries of the EPL.Panel (a) shows the transitional phase where Ω is the IDS domain with Γ 1 and Γ 2 the Neumann and Dirichlet boundary conditions, respectively, ω is the domain in which the EPL equation is solved with a zero flux boundary condition γ.In panel (b), the EPL is efficient, in this case the boundary is such that γ = γ 1 + γ 2 where γ 2 = γ Γ 2 is a Dirichlet boundary condition and γ 1 is a zero flux boundary condition.

Fig. 5 .Fig. 5 .Fig. 6 .
Fig. 5. Maps of Arolla Glacier showing the water head of the IDS (color scale in meters) and the development of the EPL (hatched zone)at the end of the winter season for four different IDS transmitivity values (Ti).For the highest IDS transmitivity, all the produced water is drained by the sediment layer explaining the very low water head and the non-development of the EPL.The white thick line indicates how the maximum length of the EPL is determined.

Fig. 7 .
Fig. 7. Maximum active length of the EPL (upper panel) and IDS water head in the borehole array (lower panel) as a function of the IDS transmitivity T i for four different values of the leakage factor ϕ and an EPL transmitivity of T e = 7.9 × 10 −2 m 2 s −1 .

Fig. 8 .
Fig. 8. Evolution with time of the maximum length of the modelled EPL (lines) compared to the computed position of the head of the channelized drainage system (black dots) in the upper panel and of the IDS water head in the lower panel for different values of the IDS transmitivity.The dashed line in the lower panel represents the flotation limit.Simulations are performed with a constant EPL transmitivity T e = 7.9 × 10 −2 m 2 s −1 and a constant leakage factor ϕ = 10 m.

Fig. 9 .
Fig. 9. Same as Fig. 8 but with varying leakage factors.The grey line of the lower panel is the IDS water head for a simulation without the EPL.Simulations are performed with a constant IDS transmitivity T i = 1.6 × 10 −3 m 2 s −1 and a constant EPL transmitivity T e = 7.9 × 10 −2 m 2 s −1 .

Fig. 10 .
Fig. 10.Comparison between the modelled active EPL (black dashed zone) and the observed channel system (blue line) at the end of the melt season for two different values of the leakage factor ϕ. The observed channel system and the corresponding moulin positions (red circles) are reproduced from Sharp et al. (1993).The moulins used for the simulations are marked by yellow squares.The color scale represents the water head of the IDS in meters.Parameters of the simulation are T i = 1.6 × 10 −3 m 2 s −1 and T e = 7.9 × 10 −2 m 2 s −1 .

Fig. 11 .
Fig. 11.Comparison of the evolution of the maximum length of the modelled EPL (red line), the computed position of the head of the channelized drainage system (black dots) and the interpretation of the channel spreading by Nienow et al. (1998) (black line).The modeling is performed using the water input of the 1993 melt season (Arnold et al., 1998) with T i = 1.6 × 10 −3 m 2 s −1 , T e = 7.9 × 10 −2 m 2 s −1 and ϕ = 10 m.The position of the head of the channelized system is computed following Nienow et al. (1998) using dye tracing data from the 1990 melt season.

Fig. 12 .
Fig. 12. Evolution of the effective pressure (upper panel, left axis, black curve), water input (upper panel, right axis, blue curve) and longitudinal velocity of the reference point, depicted in Fig. 4, (lower panel) throughout the melt season.The simulation is performed with A s = 1.6 × 10 −23 m Pa −3 s −1 and C = 0.5.

Fig. 13 .
Fig. 13.Evolution of the effective pressure (upper panel) and surface longitudinal velocity of the reference point, depicted in Fig. 4, (lower panel) during the spring speed-up event for different values of the parameter C.

Fig. 14 .
Fig. 14.Spatial pattern of modelled surface velocities on the tongue of Haut Glacier d'Arolla.Left panel is before the spring event (days 182 to 185), center panel is during the speed-up (days 186 to 188) and right panel is after the event (days 188 to 190).The surface velocities are given by the color scale in m d −1 and contoured every 0.005 m d −1 .The dots on the three panels present the measured velocities of the first 1998 speed-up event as documented in Mair et al. (2003) with the same colorscale as the one of the model results.

Table 3 .
Values of the hydrological parameters for the different steps of the selection procedure.