A double continuum hydrological model for glacier applications

The flow of glaciers and ice streams is strongly influenced by the presence of water at the interface between ice and bed. In this paper, a hydrological model evaluating the subglacial water pressure is developed with the final aim of estimating the sliding velocities of glaciers. The global model fully couples the subglacial hydrology and the ice dynamics through a water-dependent friction law. The hydrological part of the model follows a double continuum approach which relies on the use of porous layers to compute water heads in inefficient and efficient drainage systems. This method has the advantage of a relatively low computational cost that would allow its application to large ice bodies such as Greenland or Antarctica ice streams. The hydrological model has been implemented in the finite element code Elmer/Ice, which simultaneously computes the ice flow. Herein, we present an application to the Haut Glacier d’Arolla for which we have a large number of observations, making it well suited to the purpose of validating both the hydrology and ice flow model components. The selection of hydrological, under-determined parameters from a wide range of values is guided by comparison of the model results with available glacier observations. Once this selection has been performed, the coupling between subglacial hydrology and ice dynamics is undertaken throughout a melt season. Results indicate that this new modelling approach for subglacial hydrology is able to reproduce the broad temporal and spatial patterns of the observed subglacial hydrological system. Furthermore, the coupling with the ice dynamics shows good agreement with the observed spring speed-up.


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 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.Iken and Bindschadler, 1986;Mair et al., 2003;Harper et al., 2007;Fudge et al., 2009).
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.Pimentel et al. (2010) and Hewitt (2013) 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 drainage systems, they can be classified into two main groups: (i) inefficient drainage 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, 1986) and (ii) efficient drainage 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, 1986;Fountain and Walder, 1998;Boulton et al., 2007;Schoof, 2010a).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 drainage 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 water out of the inefficient one, which in turn induces a decrease in the basal velocities (Björnsson, 1974;Magnusson et al., 2010).
Recently published subglacial hydrological models take into account inefficient and efficient components for the drainage system (Pimentel et al., 2010;Schoof, 2012;Hewitt et al., 2012;Werder et al., 2013).Following the work initiated by Flowers and Clarke (2002a), 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 motivate 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, 1999;Gulley et al., 2012).
The model is tuned and validated by performing a set 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.To our knowledge, this is the first time that results from a three-dimensional coupled hydrology-ice flow model are compared to a data set including both hydrology and ice flow data for a winter to summer transition.
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 glaciers 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 pressure p w : N = −σ nn − p w . (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 vector normal 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 defines in a continuous manner the location where the efficient drainage system is most likely to develop.This strategy has the advantage of requiring a lower spatial resolution than the discrete approaches describing each channel individually (Schoof, 2010b;Hewitt et al., 2012;Werder et al., 2013).Furthermore, the use of a diffusion equation to compute the water head distribution in both systems allows for the implementation of an implicit time-stepping scheme which yields to a rather stable system, which in turn reduces the computational cost.

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 at each point of the porous medium (Bear, 1988).Mass Table 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.
filtration velocity of water with respect to the porous media [m s For the porous medium, mass conservation reads where U w refers to the filtration velocity of the fluid with respect to the fixed referential, while U j takes into account the filtration velocity of the solid with respect to the same referential.The filtration velocity corresponds to the velocity that the material would have should it use all of the available section.This definition relies on the first experiments made by Darcy in 1856, where the velocities 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 is an inflow (outflow if negative) by unit of volume which is due to 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 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, and p j is water pressure in the porous medium j at the altitude z j .In Eq. ( 5), U d is computed as the velocity of the water with respect to the porous medium in the fixed referential frame, which is obtained from the composition of U w and U j .
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 drainage 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 where β s and α are the compressibility of the solid phase and the porous media, respectively, while β w is the compressibility of water.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 is 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 transmissivity 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 drainage of water through a sediment layer (e.g.Boulton and Dobbie, 1993;Fischer et al., 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 water flux generating water heads larger than h max is transferred to the efficient layer 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.

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 system usually considered as inefficient in glaciology (e.g.Boulton and Hindmarsh, 1987;Hubbard and Nienow, 1997) 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 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 characteristics of this system, namely, a low storing capacity and high conductivity.The physical parameters of the EPL (i.e.T e , S e ) are adjusted to account for the high hydraulic transmissivity and the low storage coefficient which characterize such efficient drainage systems (Hubbard and Nienow, 1997).Moreover, the solving 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 has stayed below h max for all the preceding times (see Fig. 1a) or (ii) the EPL is active where h i has reach h max at least one time in is history (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 The Cryosphere, 8, 137-153,  layer.Once h i reaches h max at some places, the efficient drainage system is activated and starts to fill up from its initial head, which is given by the water head at the snout of the glacier.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 sink.In glaciology, an infinite sink can be a large subglacial lake, the snout of the glacier or the ocean.At these places, the 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 icebedrock interface on which the hydrological system is defined), the hydrological boundary condition is of Neumann type on the sides ( 1 ) and of Dirichlet type for any infinite sink ( 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 in a transitional state if γ 2 = 0 (Fig. 2a), ii. the EPL is in an effective 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 on the EPL.
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 neighbouring closed EPL element with the lowest IDS water head is activated.The solver is then iterated on the new domain to control that the new EPL water head is below the upper limit, if this is not the case, a new element is activated further down the hydropotential gradient until the upper limit condition is verified on all nodes of the domain where the EPL is in a transitional state.
Once the EPL becomes effective, 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, the EPL parameters are fixed throughout the simulation, which does not allow taking into account the modification in the draining capacity of the efficient drainage system observed on the field.The evolution of the draining capacity of the EPL, ultimately leading to its closure would be required to perform pluri-annual simulation.This feature is not yet included in the current version of the model, and therefore the applications are restricted to the opening phase of the efficient drainage system.Nevertheless, future work will focus on making the capacity of the EPL non-constant.

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 transmissivity T i ) and of the leakage length scale ϕ, such that The leakage length scale ϕ is a characteristic distance that the water has to cross to pass from one to the other 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 IDS S e is used instead.Once the EPL is in an effective state, the establishment of the transfer flux allows the water head in the IDS to be lowered due to the highest transmissivity of the EPL which yields lowest water head in this system.

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 .Discretization of this system finally leads to a formulation of the problem such that where H j is the solution vector, 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 discretize the time derivative, with being 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 is 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 X i matrix is the same as the X i matrix except for the pth 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 pth 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 treated as a source term in 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 drainage 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 summarized before 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.
Ice flow is governed by the Stokes equations that consist of the conservation of mass for incompressible fluids and the conservation of linear momentum 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.This coupling is achieved by computing the effective pressure N from the water pressure in the IDS.The pressure in the EPL is not taken into account as it represents a local pressure which is not likely to have a strong effect on glacier sliding (Hewitt and Fowler, 2008).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 fulfils the upper limit on the basal drag for a finite sliding velocity known as Iken's bound (Iken, 1981).
In the simplified case where the post-peak decrease exponent is equal to one, this Coulomb-type friction law reads with In these relations, τ bi = t i • (σ n)| b is the basal shear stress in the tangential direction i (i = 1, 2), and ) 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.With the assumption of a postpeak exponent equal to 1, τ bi /N monotonously increases to its upper bound C.
The coupling between the two systems through this friction law and the input of the normal stress in the hydrological model needs further iteration between the two to achieve stability.Regarding the different timescale of the involved processes, the hydrological and ice flow component of the model are solved with the time step required for their specific needs.A typical simulation would then have time steps of the order of an hour to solve the hydrological equations, and the ice dynamics will then be solved on a daily time step.
www.the-cryosphere.net/8/137/2014/The Cryosphere, 8, 137-153, 2014 For each time step when the hydrological and ice flow components are solved, iterations are performed between the two to achieve stability of the coupled system.

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 on unconsolidated sediments (Copland et al., 1997;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 labelled borehole array in Fig. 4. A study by Hubbard et al. (1995) in the same area gives a range of values for the sediment hydraulic conductivity in the vicinity.Some other studies, involving dyetracing experiments (e.g.Mair et al., 2002a;Nienow et al., 1998) give a good insight to the evolution of the subglacial drainage 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), their thicknesses e j (j = i, e), and the leakage length scale ϕ.We further assume that the transmissivity 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 transmissivity 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 prevents us from directly using parameter values that would characterize a true discrete channel.The EPL transmissivity, leakage length scale and layer thickness then have to be estimated by comparing the model results directly to observations of the hydrological systems.
We therefore adopted the following strategy.As a first approximation, a broad range of values for the unknown parameters are 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 hydro- logical system of Haut Glacier d'Arolla.The comparison to large-scale features allows the local variability of water head observed in neighbouring boreholes to be discarded.The strategy for evaluating the layers thickness is rather different.The presented simulations are performed for layer thickness of e i = 20 m for the IDS and e e = 1 m for the EPL; for each experiment, simulations are performed with values around these references to assess the sensitivity of the model.
The sediment transmissivity 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 T j S j .
(30) Hubbard et al. (1995) measured hydraulic diffusivities ranging from 4×10 1 m 2 s −1 near an efficient drainage zone down to 7×10 −2 m 2 s −1 70 m away from this zone.From these values, Eq. ( 30) and the physical parameters needed to evaluate S i (given in Table 2), the IDS transmissivity T i is estimated to range from 1.4 × 10 −4 to 8 × 10 −2 m 2 s −1 .The choice of The Cryosphere, 8, 137-153, 2014 the EPL transmissivity 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 transmissivity 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 transmissivity is then fixed at 8 × 10 −1 m 2 s −1 , ten times larger than the maximum value of the IDS transmissivity.This scaling is consistent with the results of et al.
(1998) which describe the differences between the mean flow velocity of the distributed and of the channelized drainage systems.Finally, the leakage length scale 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 drainage system during summer.The comparison is done using a large-scale feature of the hydrological system which will be comparable to the results given by the double continuum approach.
The 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 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 an EPL length substantially lower than the recorded 700 m.The other constraint given by these observations is the timing of the up-glacier migration of the head of the channel system throughout the melt season.This dynamic aspect is compared to the evolution of the EPL length during the transient summer simulation.
The mean water head of the IDS in the borehole array shown in Fig. 4 further referenced as to IDS water head is also presented to help the comparison between simulation results.

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 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 pressure on the whole glacier and the extent of the active EPL at the end of winter for three different values of the IDS transmissivity.The length of the EPL, measured along the effective EPL from the snout of the glacier to the farthest source of the EPL, is depicted by the white thick line in Fig. 5. Comparison of steady-state configurations indicates that an increasing IDS transmissivity leads to a decreasing water pressure and a shorter EPL.As expected, the large range of values for the IDS transmissivity leads to a large spread in the model results.For T i = 1.6×10 −4 m 2 s −1 the drainage system is dominated 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 transmissivity covers all possible behaviour of the IDS.Nevertheless, excessively high values of the IDS transmissivity lead to unrealistic behaviour.Using the position of the head of the channelized component at the end of spring (Nienow et al., 1998;Mair et al., 2002b), the range of the IDS transmissivity can be constrained.
Thus, only transmissivity values that lead to an EPL length lower than 700 m will be considered as admissible.Figure 6 shows the evolution of the EPL length (a) and of the IDS water head (b) for IDS transmissivity ranging from 1.4 × 10 −4 to 8 × 10 −2 m 2 s −1 .As in Fig. 5, a decrease in IDS transmissivity leads to an increase in the EPL length.The adopted restriction on the EPL length is represented by the grey zone in Fig. 6a.This constraint leads to a new range for the admissible values of the IDS transmissivity from 3 × 10 −4 to 8 × 10 −2 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 transmissivity.The relatively small line thickness indicates that, in the case of the steady-state configuration characterizing the end of winter, both drainage systems are quite insensitive to the EPL transmissivity.

Parameter
Starting range Range after end-of-winter selection Final value 1-50 1-20 10 The IDS water head for its part increases with decreasing IDS transmissivity 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 decrease of the IDS water head.After reaching a minimum around T i = 1.5 × 10 −4 m 2 s −1 , the IDS water head increases again in response to the decrease in IDS transmissivity.The value of the local minimum is a function of the drainage efficiency of the EPL but is not sensitive to the EPL transmissivity.It should then be driven by the leakage length scale as discussed below.
Figure 7 shows the sensitivity experiments to the leakage length scale ϕ.As explained in Sect.2.4, a large leakage length scale implies low efficiency of the water transfer between the two layers.This weak transfer triggers a larger extent of the EPL due to the higher water head in the IDS.The IDS water head is more sensitive to the leakage length scale than is the EPL length metric.As shown in Fig. 7, the drop of IDS water head is amplified for smaller leakage length scale.This amplitude variation is explained by the ability of the EPL to drain water from the IDS.The smaller the leakage length scales, the easier the transfer of water from the IDS system to the EPL system.For unrealistically large leakage length scales (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 length scale is therefore restricted from 1 to 20 m.In summary, the end-of-winter configuration allows the range of both the IDS transmissivity and the leakage length scale to be decreased, but not the EPL transmissivity (values given in the second column of Table 3).Modifying the thickness of both layers in this experiment while keeping the transmitivities at the same values (an increase in thickness then leads to a decrease in conductivity) does not lead to any change in the observed results.However, if the same changes in thickness are done with a constant conductivity (an increase in thickness then leads to an increase of the transmissivity), then the response of the model is on the line of the one that is observed for varying transmitivities with a negligible impact of the thickness change.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. 4 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 draining the major part of the glacier bed is obtained by the end of the melt season.As for the end-of-winter configuration, the model results are strongly dependent on the IDS transmissivity value but rather insensitive to variations in EPL transmissivity.Because of this lack of sensitivity, only the simulation performed for T e = 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 length and IDS water head.In the range of applied transmissivity, an increase in  IDS transmissivity 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 transmissivity, as can be shown in Fig. 8b.From this sensitivity study, an IDS transmissivity of 1.6 × 10 3 m 2 s −1 seems to best reproduce the observed development of the drainage system during summer.
The sensitivity of the model to a variation in the IDS layer thickness is similar to the response which is observed when varying the transmissivity.If the transmissivity is kept constant while the IDS thickness is increased, the increase of the IDS storing coefficient leads to a time lag in the opening of the EPL.The value of 20 m taken here is the one that fits best with the chosen IDS transmissivity value.The thickness of the EPL is limited by the fact that its storing coefficient should remain below that of the IDS.Modifying the thickness of the EPL while taking this limitation into account does not lead to significant differences in the model results.The EPL thickness is then kept at e e = 1 m.
The model's sensitivity to the leakage length scale is presented in Fig. 9; panel b indicates that a higher leakage length scale leads to a higher IDS water head because of the less efficient transfer from the IDS system to the EPL system.As shown in Fig. 9a, the EPL length metric shows little sensitivity to varying the leakage length scale except for very low leakage length scale values (ϕ < 10 m).As for the end-ofwinter configuration, the IDS water head is more sensitive to the leakage length scale than is the EPL length.Gordon et al. (1998)  70 m at the opening of the efficient drainage system, in good agreement with the modelled drop for leakage length scale ϕ ≥ 10 m. Figure 10 compares the modelled EPL extent 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 length scale value, bearing in mind that the observed channel system is itself a reconstruction from dye-tracing measurements.In the following, the value ϕ = 10 m is therefore adopted.
Figure 11 presents the evolution of the EPL length for the adopted values of the three hydrological parameters for different grid resolutions.Specific lengths of the element ranging from 25 to 100 m show similar results in the spreading velocity of the EPL.As in Werder et al. (2013), the results are also impacted by the position of the nodes, which could change the activation point of the EPL and then its general pattern.However, this sensitivity does not affect the global variables of the model.The estimated length of the channelized drainage system (black dots) is presented as a reference.Our results seem to indicate that the EPL extent is less smooth than the one proposed by Nienow et al. (1998) (black line) 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 this 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 characterize them as being short-lived (three to four days) periods during which the surface velocities show a twoto four-fold increase with respect to their average values.We will further use these general characteristics for the comparison with the modelled velocities.The first spring speed-up of 1998, as documented by Mair et al. (2003), will be used as a representative speed-up event.
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 and the friction law exponent n are fixed to 1 and 3, respectively.The latter is fixed at the usually accepted value of Glen's exponent, whereas the former 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 sliding for bed elevation above 3000 m is computed with an effective pressure of N = 1.2 MPa rather than the computed one.This limitation was necessary to stabilize the ice flow model on these parts of the glacier where the ice flow takes the form of steep ice falls (Hubbard et al., 2000).Contrary to what was done in previous studies (Hubbard et al., 1998), this steep area was kept in the modelled domain to avoid the use of fictitious boundary conditions, but, on the other hand, this requires constraining the friction law with a fixed effective pressure.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 et al. (2010), the value C = 0.84 m 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 which is highlighted by the sudden increase of the EPL head.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.
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 and summer, 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 speedup 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 speedup.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 assumptions 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.

Discussion
The approach presented here gives a new alternative for the modelling of subglacial hydrology.The aim of this work is to propose a hydrological model which evaluates the water pressure at the base of large ice sheet outlet glaciers and allow the computation of sliding at their bases.  of this model to a small valley glacier such as the Haut Glacier d'Arolla might not be the best suited to present the advantages of our approach, but the data set available on this glacier allows the building of convincing validation experiments for this new approach.Compared to existing largescale models (e.g.Le Brocq et al., 2009;Goeller et al., 2013), our model effectively computes the effective pressure.This feature makes it comparable to more physically based models (e.g.Schoof, 2010b;Hewitt et al., 2012;Werder et al., 2013) in which the efficient drainage system is modelled by channel, or the approach of Flowers and Clarke (2002a) in which a single porous layer with varying capacity accounts for both efficient and inefficient drainage systems.
In comparison to the model incorporating a discrete description of the channels (e.g.Schoof, 2010b;Hewitt et al., 2012;Werder et al., 2013), in the proposed approach the precise location of the efficient subglacial drainage system is not achieved.This prevents validation of the model against punctual water pressure measurements in borehole or precise modelling of dye tracer return.This first limitation is intrinsic to the model formulation but is not an issue as far as the model results can be compared to global variables as is done in this study.
A second limitation of the model, not intrinsic to its formulation, is the current lack of a closing mechanism for the EPL.This type of mechanism would require the implementation of an evolving transmissivity for the EPL.In the present study, the lack of this closing mechanism is of small concern as we focus on the opening of the EPL.However, the end of the summer simulation would probably take advantage of the implementation of an evolving draining capacity of the EPL, which would probably allow accommodation of the inputs that occur later in the season (days 220 to 240 of the cou-  pled simulation).Introducing this mechanism will be part of future development of the model.

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.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 data set 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 transmissivity of both drainage systems and the leakage length scale.In the second experiment, the evolution of the drainage system during the spring and summer seasons 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 in terms both of the temporal and spatial drainage system evolution during the spring and summer seasons 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 length scale 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 transmissivity or A s ), or the model has proven to be quite insensitive to them (as for the EPL transmissivity).

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 ] n flow law exponent N effective pressure [Pa] p j water pressure in the media [Pa] 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 length scale [m] ω j porosity of the media conservation for the water is given as

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 effective.

Fig. 2 .
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 effective; 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. 3 .
Fig. 3. Schematic description of the iteration scheme of the hydrological model.The outer box represents the entire model and each of the inside boxes is a component of the system which is solved with information from the other components.The convergence or not of each system is indicated with its iterative loop.The red arrows represent the start and end of a hydrological time step.

Fig. 4 .
Fig. 4. Map of Haut Glacier d'Arolla with the position of the borehole and velocity stakes arrays.The star is the position of the reference point used in Sect. 5.The glacier surface elevation is contoured every 100 m.Red circles are the position of the moulins used for the dye-tracing experiments in 1989 and 1990; the yellow squares are the position of the moulins recorded in 1993 and used for the modelling.

Fig. 5 .
Fig. 5. Maps of Arolla Glacier showing the water pressure of the IDS and the development of the EPL (hatched zone) at the end of the winter season for three different IDS transmissivity values (T i ).For the highest IDS transmissivity, all the produced water is drained by the sediment layer, explaining the very low water pressure and the non-development of the EPL.The white thick line indicates how the length of the EPL is determined.

Fig. 6 .
Fig. 6.Length of the EPL (a) and IDS water head in the borehole array (b) as a function of the IDS transmissivity T i .The grey zone in (a) indicates the admissible values for the EPL maximum length.The dashed line in (b) represents the flotation limit.The spread of the curves represents the scattering due to EPL transmissivity ranging from 4×10 −3 to 8×10 −1 m 2 s −1 , with the higher transmissivity values leading to the lowest EPL length and IDS water head.

Fig. 7 .
Fig. 7. Length of the EPL (a) and IDS water head in the borehole array (b) as a function of the IDS transmissivity T i for four different values of the leakage length scale ϕ and an EPL transmissivity of T e = 7.9 × 10 −2 m 2 s −1 .The dashed line in (b) represents the flotation limit.

Fig. 8 .
Fig. 8. Evolution with time for different values of the IDS transmissivity of (a) the maximum length of the modelled EPL (lines) and (b) the IDS water head.The position of the head of the channelized drainage system derived from observations (black dots) is presented in (a) for comparison.The dashed line in (b) represents the flotation limit.Simulations are performed with a constant EPL transmissivity T e = 7.9 × 10 −2 m 2 s −1 and a constant leakage length scale ϕ = 10 m.

Fig. 9 .
Fig. 9. Same as Fig. 8 but with varying leakage length scales.The grey line in (b) is the IDS water head for a simulation without the EPL.Simulations are performed with a constant IDS transmissivity T i = 1.6 × 10 −3 m 2 s −1 and a constant EPL transmissivity 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 length scale ϕ.The observed channel system and the corresponding moulin positions (red circles) for the summer of 1989 and 1990 are reproduced from Sharp et al. (1993).The moulins observed during the 1993 season which are used for the simulations are marked by yellow squares.The colour scale represents the water head of the IDS in metres.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 for different grid specific length (coloured lines).The position of the head of the channelized drainage system derived from observations (black dots) and the interpretation of the channel spreading byNienow et al. (1998) (black line) are shown for comparison.The modelling is performed using the modelled 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 followingNienow et al. (1998) using dye-tracing data from the 1990 melt season; the specific length (S L ) for each simulation is given in metres.

Fig. 12 .
Fig. 12. Evolution throughout the melt season at the reference point depicted in Fig. 4 of (a) the effective pressure (left axis, black curve) and the water input (right axis, blue curve), (b) the surface longitudinal velocity.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 during the spring speed-up event for different values of the parameter C of the effective pressure (a) and surface longitudinal velocity (b) at the reference point, depicted in Fig. 4.

Fig. 14 .
Fig. 14.Spatial pattern of modelled surface velocities on the tongue of Haut Glacier d'Arolla for the year 1993.Before the spring event (days 182 to 185), during the speed-up (b) (days 186 to 188) and after the event (c) (days 188 to 190).The surface velocities are given by the colour 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 colour scale as that of the model results.

www.the-cryosphere.net/8/137/2014/ The Cryosphere, 8, 137-153, 2014
The compressibility β s defines the compressibility of the solid phase of the media (i.e., it can be assessed by a

www.the-cryosphere.net/8/137/2014/ The Cryosphere, 8, 137-153, 2014 142 B. de Fleurian et al.: Glacio-hydrological model drainage
system.The introduction of the storage coefficient S j (j = i, e) is needed to convert the water head into volume due to the characteristic of the confined aquifer in which the storage is due to a compression of water and porous medium.

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 labelled pkp.

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