Investigating changes in basal conditions of Variegated Glacier prior to and during its 1982-1983 surge

The Variegated Glacier (Alaska) is known to surge periodically after a su ﬃ cient amount of cumulative mass balance is reached, but this observation is di ﬃ cult to link with changes in the basal conditions. Here, using a 10-year dataset, consisting in surface topography and surface velocity observations along a ﬂow line for 25 dates, we 5 have reconstructed the evolution of the basal conditions prior and during the 1982– 1983 surge. The model solves the full-Stokes problem along the central ﬂow line using the ﬁnite element method. For the 25 dates of the dataset, the basal friction parameter distribution is inferred using the inverse method proposed by Arthern and Gudmunds-son (2010). This method is here slightly modiﬁed by incorporating a regularisation term 10 in the cost function to avoid short wave length changes in the friction parameter. Our results indicate that dramatic changes in the basal conditions occurred between 1973 to 1983. Prior to the surge, periodical changes can be observed between winter and summer, with a regular increase of the sliding from 1973 to 1982. During the surge, the basal friction decreased dramatically and an area of very low friction moved from 15 the upper part of the glacier to its terminus. Using a more complex friction law, these changes in basal sliding are then interpreted in terms of basal water pressure. It con-ﬁrms that dramatic changes took place in the subglacial drainage system of Variegated Glacier, moving from a relatively e ﬃ cient drainage system prior to the surge to an in-e ﬃ cient one during the surge. By reconstructing the water pressure evolution at the 20 base of

From the well-studied 1982-1983 surge, it seems that Variegated Glacier is characterised by a two-phase surge, each phase with a reasonably distinct termination separated by one year (Eisen et al., 2005).One other characteristic is the seasonal timing of Variegated surges, with an onset in late autumn or winter and termination in late spring or early summer.
As shown by Eisen et al. (2001), the duration of the quiescent phase in between two surges is very well correlated with the total cumulative mass balance at a point located at the altitude of 1500 m in the accumulation area.The Variegated Glacier is found to surge each time the ice-equivalent cumulative balance at this particular point reaches the threshold value of 43.5 ± 1.2 m.This relation is not fulfilled for the 2003-2004 surge, for which the cumulative mass balance was only half of that required for previous surges (Harrison et al., 2008).As anticipated by Eisen et al. (2005), this loss of correlation might be explained by the early termination of the one-phase 1995 surge and its unusual post-surge surface topography corresponding to a relatively small mass transfer from the upper part to the lower part of the glacier.Because the 2003-2004 was a normal two-phase surge, Harrison et al. (2008) have predicted that the mass balance correlation will hold for the next surge.Nevertheless, the causality of this mass balance surface observation cannot easly be linked to the basal processes controlling the surge.Surges are initiated by a change in the basal hydrological system, which moves from a discrete efficient system with low water pressure and high water discharge to a distributed inefficient system with high water pressure (e.g., Raymond, 1987).A discrete efficient system is usually formed by a few large channels and its influence on the ice flow is relatively low, whereas an inefficient system consists in small linked cavities strongly influencing the basal velocity (Kamb, 1987) (2005), Variegated surges are certainly governed by a change in the glacier's geometry which affects the internal drainage system.When a threshold amount of geometry change is reached, the discrete system closes at the end of the melting season when the amount of water is insufficient to keep it open.Then, subsequent rain or meltwater from the surface, even in small volume, will progressively contribute to increase the basal water pressure, finally leading to the glacier surge.The following spring, when the amount of water is again sufficient, the discrete efficient system opens again and the surge stops (Harrison and Post, 2003;Lingle and Fatland, 2003).Note that this interpretation is consistent with the observed timing of Variegated surges, which start during the winter and end during the summer.
During the surge, short-term variations (hours to days) of ice velocity, water pressure and outflow stream at the glacier terminus have been observed.These observations indicate the predominant contribution of basal sliding during the surge phase.Measurements of the internal deformation in a borehole during the surge show that 95 % of the surface velocity is due to sliding (Kamb et al., 1985).Velocities as high as 50 m day −1 were measured during the second phase of the 1982-1983 surge.Simultaneous records of water pressure from borehole measurements indicate the strong correlation between water pressure and velocity.Pulses in surge movement do indeed correspond to peaks in pressure.Conversely, the increase of the outflow stream at the terminus is closely correlated with a rapid slowdown of the glacier (Kamb et al., 1985).This last observation indicates that a large amount of water is stored in subglacial cavities, inducing an increase in water pressure and a consequent increase in ice sliding velocities.But when a threshold pressure is reached, the subglacial water storage purges, leading to flooding at the terminus outflow and to a slowdown of the ice sliding.
Extensive measurements of the surface topography and surface velocities were carried out during the 1973-1983 decade (Bindschadler et al., 1977;Kamb et al., 1985;Raymond and Harrison, 1988).This measurement period covers the quiescent phase which follows the [1964][1965]  period, surface elevation and horizontal surface velocity were measured at 25 different dates, twice a year prior to the surge and 8 times during the 2 years of the surge.At each date, the dataset is composed of the horizontal surface velocity and the surface elevation every 250 m along the 20 km of the central flow line.Most of the datasets are incomplete, mainly in the upper and lower parts, but also where the glacier was too crevassed to be accessible.Few attempts have been made to reconstruct the basal condition history below Variegated Glacier from these datasets.Raymond and Harrison (1988), using a very simple flow line model, determined that basal sliding increased from 1973 to 1981 and concluded that by 1981, basal sliding might be 50 % or more of the total surface velocity in the upper part of the glacier.Again with a relatively simple hydrological model, Eisen et al. (2005) showed that the flux required to keep the efficient drainage system open is very sensitive to the basal shear stress.Combining the model and the observations, they determined a critical basal shear stress along the flow line which initiates a surge.
In this paper, we propose to use the very well documented period from 1973 to 1983 to reconstruct the history of the basal conditions below the Variegated Glacier using a full-Stokes model.Determining the optimal basal conditions from the glacier topography and the surface velocities is an inverse problem.Recently, three methods have been proposed to solve this particular inverse problem using a full-Stokes direct model.The first one, is a Bayesian method developed by Gudmundsson and Raymond (2008) and further applied to the Rutford ice stream (West Antarctica, Raymond Pralong and Gudmundsson, 2011).Note that for this application to real data, both the basal friction and the bedrock topography were inverted.The two others, which belong in the class of the the variational methods, are a control method using the adjoint model of the linear Stokes equations (Morlighem et al., 2010) and a Robin inverse method (Arthern and Gudmundsson, 2010).These two variational methods rely on the minimisation of a cost function that measures the mismatch between the model and the observations.In each case, the gradient of the cost function with respect to the basal drag coefficient is obtained analytically Introduction assuming a linear flow law and a linear sliding law.Theoretically, these results could be extended to non-linear laws but this would require further analytical and numerical developments.In their applications, Morlighem et al. (2010) and Arthern and Gudmundsson (2010) show that even by using the gradient derived in the linear case, it is possible to minimise the cost function with non-linear laws, but this could fail for some applications (Goldberg and Sergienko, 2011).These two methods should lead to very similar solutions for the basal drag coefficient and both have advantages and drawbacks.The control method needs the derivation of the adjoint model but it is easy to modify the cost function to take into account the error on the observed velocities.
The Robin inverse method can be easily implemented using the direct model only, but does not integrate the observation errors in the cost function.The two methods were implemented in the finite element code Elmer/Ice and compared based on some of the datasets.They showed very good agreements, and in this paper we present only the results obtained with the Robin inverse method (Arthern and Gudmundsson, 2010) extended with a regularisation term.To our knowledge this is the first application of this method to real data in glaciology.
The direct full-Stokes flow line model is presented in the first Section.The associated inverse model (Arthern and Gudmundsson, 2010) and its extension is presented in the second Section.In the third Section, the inverse model is used to infer the basal friction distribution along the flow line at each measurement date.In the fourth Section, following the idea proposed by Flowers et al. (2011), changes in the basal friction parameter are interpreted in terms of changes in basal water pressure through the use of the water pressure dependent friction law proposed by Schoof (2005) and Gagliardini et al. (2007).Finally, using the basal friction parameter distributions inferred from the inverse method, a transient simulation is run over the 10-year data period to compare Therefore, the modelling is limited to a two-dimensional flow line geometry, delimited by the bedrock b(x) and the upper surface z s (x).We further assume a Cartesian coordinate system such that x is the horizontal direction and z the up-oriented vertical one.For a given geometry, the ice flow is governed by the Stokes equations, i.e. the mass and momentum conservation equations in which the acceleration terms are neglected.The Stokes equations write: divσ Here u = (u x ,0,u z ) is the velocity vector, σ = τ −pI is the Cauchy stress tensor and p the isotropic pressure, ρ the ice density and g = (0,0,−g) the gravity vector.The body force f is added to account for the friction arising on the lateral side of a real glacier in the flow line model.To this end, the concept of shape factor (Nye, 1965) is here extended to the full-Stokes formulation by defining the body force f l as where the shape factor f = f (x) is a scalar function of the transversal shape of the glacier and t is the unit vector tangent to the upper surface.As shown by this equation, the concept of shape factor adds a resistive body force parallel to the tangent to the upper surface.When f = 1, the limit case of an infinitely large glacier is obtained, whereas small f stands for narrow and/or deep transverse sections.
Here, we evaluate f (x) by assuming that the transverse shape of the bedrock is a parabola of the form b(x,y) = b(x) + a( x in 1973 (Raymond and Harrison, 1988).This approach accounts for variations with time of the shape factor induced by changes in ice thickness.Following the approach of Nye (1965), the relation between the shape factor and the ice thickness in the central flow line is inferred from three-dimensional full-Stokes simulations of an infinitely long glacier flowing over a parabola-type bedrock using different values of the friction paramter.All these three-dimensional simulations (not shown here), are well reproduced with a two-dimensional flow line model using the following empirical estimate of the shape factor: where is the ice thickness.Figure 1 shows the evolution of the shape factor f (x) along the flow line for the 1973 geometry.
The ice rheology is described through a power-type flow law, known as Glen's law in glaciology, linking the strain-rate tensor ε to the deviatoric stress tensor τ such as: where τ 2 e = τ i j τ i j /2 is the square of the second invariant of the deviatoric stress and A a rheological parameter, which depends on the ice temperature via an Arhenius law.Since the Variegated Glacier is temperate, the constant value A = 100 MPa −3 a −1 is adopted (close to the ones proposed in Cuffey and Paterson, 2010).

Boundary conditions
The upper surface Γ s , i.e. z = z s , is a stress-free surface and the following Neumanntype boundary condition applies: At the bedrock interface Γ b , i.e. z = b, zero basal melting is assumed (u • n = 0) as well as a linear friction law (Robin type boundary condition).This linear friction law relates the basal drag τ nt to the sliding velocity u t such as: where n and t are the normal and tangent unit vectors to the bedrock surface, and β ≥ 0 is the basal friction parameter.All these equations are solved using the finite element method with the code Elmer/Ice.More details on the numerics can be found in Gagliardini et al. (2007) and Gagliardini and Zwinger (2008).

Robin method
The inverse problem is, for each dataset (surface geometry and velocities), to determine the basal friction parameter β that gives the smallest mismatch between observed and modelled surface velocities.
We use the inverse Robin method adapted to glaciology by Arthern and Gudmunds- The cost function that expresses the mismatch between the solution of the two models on the upper surface Γ s is given by where superscripts N and D refer to the Neumann and Dirichlet problem solutions, respectively.
The G âteaux derivative of the cost function J o with respect to the friction parameter β for a perturbation β ′ is given by (Arthern and Gudmundsson, 2010): where the symbol |.| defines the norm of the velocity vector.
In this paper, to avoid unphysical negative values of the friction parameter, β is expressed as The optimisation is now done with respect to α and the G âteaux derivative of J o with respect to α is obtained as follows: In the presence of noise in the observed velocities, the method can lead to spurious small wavelength oscillations of the inferred friction parameter.Arthern and Gudmundsson (2010) suggest to stop the minimisation when the cost function starts to stagnate at a certain level.Furthermore, the authors show that this is in agreement with a heuristic stopping criterion based on the observation errors.One drawback of this approach is that on a glacier, the magnitude of the velocities and the observation Introduction errors could vary strongly from one place to another, but also from one dataset to another, so that the stopping criterion should be different for each area and each dataset.
Here, an additional Tikhonov regularisation term that penalises the small wavelength oscillations of the friction parameter β, taken as is added to the cost function J o .The total cost function now writes where ūobs is the mean value of the observed surface velocities and λ is a weighting parameter used to adjust the influence of the added regularisation with respect to the initial cost function.The term ūobs takes into account the large changes in velocity observed along the 10-year dataset and allows to use a unique value of the λ parameter for all the datasets.Regularisation is classical in data assimilation: the minimisation of J o alone is an ill-posed problem, and the addition of a regularisation term ensures existence of a global minimum.If Λ is large enough, the problem becomes well-posed, with a unique minimum, and therefore the minimisation algorithm shows improved convergence properties.The form of the additional term ensures that the optimal α is smooth.The effect of this regularisation term and the sensitivity of the β distribution to the parameter λ is discussed in the results Section.The minimisation of the cost function J tot with respect to β is done using the limited memory quasi-Newton routine M1QN3 (Gilbert and Lemar échal, 1989) implemented in Elmer/Ice in reverse communication.This method uses an approximation of the second derivatives of the cost function and is thus more efficient than a fixed-step gradient descent.

Technical aspects
The Stokes equations and the Robin problem are solved using the finite element code Elmer/Ice.For each date, a regular mesh is constructed using 80 horizontal times 20 vertical layers of quadrangle elements, between the bedrock and upper surface.For non-icy areas, a minimal thickness of 3 m is imposed to avoid zero volume elements.
Each of the 25 datasets is composed of the surface elevation and the horizontal surface velocity at 81 points regularly spaced every 250 m along the 20 km of the glacier length.Topographic measurements are representative of a given date whereas velocity measurements refer to the period in-between two measurement dates (Raymond and Harrison, 1988).For the quiescent period, the same surface topography is used for the summer and for the following winter.For the surge, because of the fast changing topography, for a given velocity measurement, the surface topography is taken as the median of the two surface topographies corresponding to the surface velocity measurement dates.To construct the 25 geometries corresponding to the 25 datasets, the surface elevation must be defined along the whole glacier.Where surface topography measurements are missing, the elevation is estimated from the other datasets using a linear adjustment to fulfil the current surface elevation continuity.For the velocity, the mesh is constructed so that point measurements and mesh nodes coincide.When solving the Dirichlet problem, measured velocities are imposed only where measurements are available and no interpolation is used to complete missing data.We verified that a finer mesh does not change significantly the results of the inversion of the friction parameter.surface and sliding velocities obtained for different values of the regularisation parameter λ are shown in Fig. 2. For the case with no regularisation, the iterations were stopped after 54 iterations when the cost function started to stagnate.The case where λ = 0 allows the best match with the observed velocities with a relative mean error equal to 4.9 %, but shows large variations of the friction coefficient (Fig. 2a) which result in large variations of the sliding velocity (Fig. 2c).The wavelength of the oscillations is larger than the mesh spacing and thus does not correspond to the short wavelength oscillations associated with the error noise presented by Arthern and Gudmundsson (2010).

TCD
With a non-zero regularisation term, at the scale of the glacier, the mean value of the inferred friction parameter is very similar for all the values of λ but the oscillations are smoothed and the smoothing is greater as λ increases.The difference between modelled and observed surface velocities remains small, but the short wavelength oscillations of the observed velocities are less well resolved when the regularisation term increases.As a result, the corresponding sliding velocities, illustrated in Fig. 2c, are also smoothed by the regularisation and some points show very large differences compared with no regularisation (λ = 0).The comparison between surface and basal velocities shows that all these small wavelength oscillations arising at the base have almost no visible influence at the glacier surface.The relative mean error increases from 5.0 % to 9.1 % when λ is increased from 10 4 to 10 6 .If in reality high and low friction areas may alternate, other approximations or uncertainties (such as the topography or 3-D effects) may explain the oscillations obtained with no regularisation.A smoothed distribution is thus preferred.Nevertheless, the optimal regularisation parameter λ can be estimated using a L-curve analysis (Hansen, 2001).The L-curve method uses the plot of the norm of the regularised solution J reg given by Eq. ( 13) versus the norm of the initial cost function J o (9) to choose the optimal regularisation parameter.As shown in Figure 3, for this application, a relatively large range of values of the regularisation parameter can be deduced from the curvature analysis of the L-curve, varying from λ ≈ 10 4 to λ ≈ 10 6 .Also, by one order on magnitude when λ is varied from 0 to 10 9 , which is relatively low in comparison to other applications.Note that the L-curve analysis should be conducted for all dataset and might conduct to different values of the regularisation parameter despite the ponderation by ūobs of the regularisation term in ( 14).Because this analysis necessitate a large number of simulations, it was not reasonably feasible for the 25 datasets, and in what follow, the value λ = 10 5 is then adopted for all the datasets.

Inferred basal friction parameter distributions
The distribution of the friction parameter was inferred using the same method for the 25 datasets available during the quiescent phase and the surge.Results are shown in Fig. 4. Seasonal changes between summer and winter can be observed.For a given year, the winter always presents higher friction than the previous summer.But, during the eight years of the quiescent phase, the friction parameter β regularly decreases, so that the last winter values are smaller than the first summer ones.Another remarkable feature is that the friction decrease is more pronounced in the upper part of the glacier than in the lower part during the quiescent phase.Contradictory to the assumption made by Bindschadler (1982), our results indicate that sliding already contributes for a large part of the total motion of the glacier during the quiescent phase.As shown in to September 1982 and ends with a punctual increase of the basal friction.The second phase starts in January 1983 and spreads down the glacier until July 1983 with a dramatic decrease of the basal friction.During the second part of the surge, we observe the propagation of a low basal friction area from the middle part of the glacier down to its terminus.At the end of the surge, the basal sliding accounts for more than 90 % of the observed surface velocities everywhere on the glacier, whereas it is only the case in the upper part during the first phase of the surge.Fig. 5 shows that in some places along the flow line, basal velocities are greater than the surface ones.This is possible because of the stress transmission when solving the Stokes system with no simplification.
The inversion procedure allows gives a good representation of the observed velocities of each date of the 10-year dataset.The observed changes in surface velocities during the quiescent and surge phases can be explained by changes in the basal sliding velocity and thus are clearly visible in the inferred distributions of the friction parameter β.Here, the simplest linear friction law ( 7) is assumed, in which the friction parameter β encompasses all the complexity of basal friction.In the next section, following the idea proposed by Flowers et al. (2011), the inferred friction parameter distribution β(x i ,t j ) (i = 1,81 points and j = 1,25 dates) is interpreted in terms of changes of the effective pressure at the base of the Variegated Glacier from 1973 to 1983, using a more complex friction law.Introduction from mathematical developments, and Gagliardini et al. (2007) from finite element simulations have both proposed a similar friction law for the flow of clean ice over a rigid bedrock in presence of cavitation.In its simplest form, where the basal drag tends asymptotically to its maximum value (post-peak exponent q = 1 in Gagliardini et al. (2007)), this friction law is of the form: In the above equation, A s is the sliding parameter in the absence of cavitation and n Glen's law exponent, resulting in a non-linear relation between the basal drag τ nt and the basal sliding velocity u t .Note that in the limit case where n = 1 and N ≫ 0, the sliding parameter A s and the friction parameter β are inversely proportional.As shown by Schoof (2005), the coefficient C is lower than the maximum local positive slope of the bedrock topography at a decimetre to meter scale, so that the ratio τ nt /N ≤ C fulfils Iken's bound (Iken, 1981).The friction law ( 15) is strongly related to the water pressure p w through the effective pressure N = −σ nn − p w .When p w = 0, the effective pressure is equal to the normal compressive Cauchy stress, and increasing water pressure leads to a decrease of N toward zero.The two parameters A s = A s (x) and C = C(x) are only a function of space whereas the time-dependent changes are due to changes in the effective pressure N = N(x,t).Note that effective pressure changes reflect changes in water pressure and/or basal normal stress, as discussed below.
Assuming A s and C are known, the effective pressure, and thus the water pressure, can be evaluated from the previous inversion results for all points x i and all dates t j , as follows: Physically, the effective pressure is bounded, i.e. 0 < τ nt /C ≤ N ≤ −σ nn , but the evaluation of N(x i ,t j ) from ( 16) can be out of these bounds due to the assumed A s and C distributions.The effective pressure being inversely proportional to the coefficient C > 0, we will hereafter assume a uniform C (C = 0.5) and discuss the bound issue only in terms of A s .The upper bound N ≤ −σ nn (or p w ≥ 0) is violated where the sliding parameter A s is too high, so that even with zero water pressure the sliding velocity is too great.The bound N ≥ τ nt /C (or p w ≤ −σ nn − τ nt /C, reminding σ nn < 0 and τ nt > 0) is never violated since it corresponds to infinitely great sliding velocity.
By considering the upper bound it is possible to estimate the sliding parameter in absence of cavitation A s so that N ≤ −σ nn is always fulfilled.Assuming zero water pressure, the distribution of the minimal sliding parameter A min s can be evaluated from Eqs. (7 and 15) in which p w = 0, such as: where β, σ nn and u t are defined at each date t j and each point x i where data are available.It is found that the minimal value of A s inferred from ( 17) is everywhere obtained for the winter 1973 dataset, except for the upper and lower points for which no measurements were performed in 1973.This indicates that basal sliding during winter 1973 represents the lowest values in the following 10-year period.

Modelled change in basal water pressure
Using the inferred sliding parameter distribution A min s (x), the water pressure for the 25 dates is then obtained from Eq. ( 16). Figure 6 shows the evolution with time (vertical axis) and space (horizontal axis) of the ratio of the water pressure over the normal stress −p w /σ nn .This ratio is plotted instead of the effective pressure N because it is visually easy to interpret.When this ratio tends toward 1, the effective pressure tends toward zero and the sliding is increased.The first noticeable result is that the large changes in the friction parameter β shown in Fig. 4 are induced by relatively small changes in terms of water pressure.In fact, as shown in Fig. 7, the basal normal stress is almost constant during the 10-year period, despite strong variations along the flow line.Therefore, changes with time of the ratio −p w /σ nn are mostly due to changes in water pressure.As shown in Fig. 6, the ratio −p w /σ nn only evolves between 0.7 to 1, if we except the winter 1973 for which it is zero due to the definition of the A s parameter.This can be explained by the shape of the friction law used and the fact that it is bounded (see Fig. 8 in Gagliardini et al., 2007).For low sliding velocity corresponding to great effective pressure, a great increase in water pressure is needed to increase the velocity.For great sliding velocity, it is the opposite, due to the asymptotical behaviour of the friction law, i.e. a small increase in water pressure leads to great increase in velocity.

TCD
Note, that even if the normal stress σ nn is almost constant, a slight and progressive increase of σ nn in the upper part of the glacier is visible in the years before the surge, induced by an increase of the observed ice thickness.During the surge, this increase propagates down the glacier attesting displacement of the ice mass.
Again, the seasonality of the sliding is visible in Fig. 6, as well as the two phases of the surge.Also, as was already inferred from the β inversion, the geatest changes are observed in the upper part of the glacier during the quiescent phase and the first surge phase.For the second phase of the surge, we observe the propagation of a high water pressure area from the upper part to the lower part of the glacier, while the pressure in the upper part still remains significant.During this last stage, the increase of water pressure, even though relatively small (2-6 %), leads to very large increase of the ice flow.

Effect of basal topography on basal water pressure
The topography clearly affects the basal water runoff below the Variegated Glacier.First, each bump in the bedrock induces a higher normal stress σ nn on its upstream face (Fig. 7).Surprisingly, it is also the places where the ratio −p that these are the areas where the increase of the water pressure is the highest.These bedrock bumps, and more particularly the one located at x = 10 km, seem to restrain the water in an upstream catchment.This interpretation is consistent with the surge hypothesis formulated by Lingle and Fatland (2003).Indeed, the progressive evolution of the surface topography of Variegated Glacier leads to an increasingly constricted water catchment upstream x = 10 km, which increases slightly the water pressure up to a threshold value for which the surge occurs.During the second phase of the surge, the propagation of the high water pressure area downstream x = 10 km (Fig. 6) is probably a consequence of the destruction of the water catchment during the first surge stage, leading to the opening of the initially constricted water catchment.

Direct prognostic simulations
From the previous analysis, the friction parameter β is reconstructed for the 25 datasets from summers 1973 to 1983.To test the sensitivity of the model to the surface geometry, from the inferred history of the friction parameter, we run a prognostic timedependent simulation assuming a constant mass balance over these 10 years.
Starting from the summer 1973 surface geometry, the upper free surface is evolved using the following kinematic boundary condition: where a s is the accumulation/ablation function.This function is estimated from the average mass balance measured on Variegated Glacier between 1972to 1976(Bindschadler, 1982).The mass balance is supposed to be linearly dependent on the surface altitude z s and the equilibrium line located at 1050 m.a.s.l., such as: The Stokes Eqs.(1 and 2) using the basal boundary condition ( 7), and the free surface evolution Eq. ( 18) are coupled and solved iteratively using a time step of 0.1 a.At each time step, the basal friction parameter β is interpolated linearly from the two closest datasets in the time-series.
The modelled surface is compared with the observed surface at four different dates just before, during and after the surge in on the modelled topography can be easily explained by the lack of data (topography, velocities and mass balance) and three-dimensional effects.

Conclusions
We have presented the first application to a real case of the inverse method proposed by Arthern and Gudmundsson (2010).It demonstrates the strong relevance of this inverse method, which allowed the reconstruction with a high accuracy of the basal conditions below Variegated Glacier along a 10-year period.From this reconstruction of the friction parameter, the water pressure was inferred using a water pressure dependent friction law.As an important result, we showed that very large changes in the basal friction parameter are induced by relatively small changes in basal water pressure.This is mainly due to the asymptotical behaviour of the friction law for great sliding velocity.Our results confirm the presence of a subglacial water storage in the mid-upper part of the glacier as proposed by Lingle and Fatland (2003).
To take this study further, our reconstruction of the basal water pressure over the 10year period should be used to constrain a hydrological model coupled with an ice flow model, as in Pimentel et al. (2010) andde Fleurian et al. (2011).This would, however, need to be conducted using a three-dimensional modelling approach to overcome the limitations of a flowline hydrological model.This would be the last step to fully relate the surge behaviour of Variegated Glacier to its surface mass balance over the decade [1973][1974][1975][1976][1977][1978][1979][1980][1981][1982][1983] is a temperate glacier located in the coastal St Elias Mountains in Alaska (USA).It is approximately 20 km long and 1 km wide, with ice flowing from the altitude of 2000 m.a.s.l down to the sea.Due to its surging behaviour, the Variegated Introduction Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | et al.
surge and includes the 1982-1983 surge.During this Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Available data for the Variegated Glacier are limited to the central flow line of the glacier.
) • y 2 , where the parabola coefficient a(x) is constant in time and estimated from the thickness and width measurements performed Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | son (2010).The method consists in solving alternately the Neumann-type problem defined by Eqs.(1, 2) and the surface boundary conditions (6), and the associated Dirichlet-type problem defined by the same equations excepted that the Neumann upper-surface condition (6) is replaced by a Dirichlet condition, such as: u(z s ) = u obs , (8) where u(z s ) and u obs are the model and observed surface horizontal velocities, respectively.This condition is enforced for each location where a surface velocity was measured.The natural Neumann condition is imposed in the vertical direction and where no observations are available.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 |

Fig. 5 ,
the contribution of the basal velocity continuously increases during the quiescent phase, from a mean value of 10 % in winter 1973 up to 60 % for summer 1981 just before the surge started, confirming the values inferred by Raymond and Harrison (1988) using a simple model.At the onset of the surge in winter 1981-1982, dramatic changes in the basal friction occurre.The inferred friction parameter drops by about one order of magnitude from 10 −3 to 10 −4 MPa m −1 a, causing a high increase of the basal sliding, as shown in Figs. 4 and 5.After this onset phase, the friction continues to decrease regularly until the end of the surge in July 1983.Both phases of the surge are visible in the results.The first phase occures only in the upper part of the glacier from March 1982 Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | ) and then p w = −N − σ nn .Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | w /σ nn increases, indicating Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper |

Fig. 8 .
The modelled and observed surface elevations before and during the surge relative to the measured surface elevation of summer 1973 are shown in Figs. 9 and 10.Before the surge we observe a thickening of the upper part of the glacier and a thinning of the lower part.The timing and the magnitude of the changes are well reproduced by the model except on the highest part of the glacier where the model leads to a thinning of the glacier.Part of the discrepancies between the model and observations can be explained by errors in the mass balance and/or by three-dimensional effects (ice convergence along the central flow line, Raymond, 1987) as the observed volume increases whereas, in the model, the integrated mass balance is slightly negative leading to a decrease of the modelled ice volume.The modelled surge occurs in phase with the observations when the friction parameter dramatically decreases in March 1982.The surge is characterised by a thinning of the upper part of the glacier and a thickening of the lower part which results in the advance of the ice front.As already observed for the quiescent phase, the upper part of the glacier is too thin when compared with summer 1973, but the timing of the mass transfer from the upper part to the lower part of the glacier is well captured by the model, and particularly the advance of the ice front position.The ability of our model to reproduce the main characteristics of the surge validates a posteriori the use of the diagnostic model to infer the basal friction distribution for each dataset of the time-series.It demonstrates that the results are not very sensitive to the surface geometry and that the Robin inverse method (Arthern and Gudmundsson, 2010) with an appropriate regularisation allows to retrieve a good order of magnitude of the friction parameter β where surface velocity observations are available.The errors Discussion Paper | Discussion Paper | Discussion Paper |

Fig. 1 .
Fig. 1.Evolution of the shape factor f along the central flow line for the 1973 surface topography.

Fig. 2 .Fig. 3 .Fig. 6 .Fig. 7 .
Fig. 2. Sensitivity of the results to the regularisation penalisation λ for the 1978 summer data.(a) Distribution of the inferred basal friction parameter β along the central flow line and corresponding (b) modelled and measured (cross) surface velocities and (c) modelled basal velocities.