Towards a better ice sheet model initialisation and basal knowledge using data assimilation

Ice flow models are now routinely used to forecast the ice-sheets contribution to 21st century sea-level rise. For such short term simulations, the model response is greatly affected by the initial conditions. Data assimilation algorithms have been developed to invert for the friction of the ice on its bedrock using observed surface velocities. A drawback of these methods is that remaining uncertainties, especially in the bedrock elevation, lead to non-physical ice flux divergence anoma5 lies resulting in undesirable transient effects. In this study, we compare two different assimilation algorithms based on adjoints and nudging to constrain both bedrock friction and elevation. Using synthetic twin experiments with realistic observation errors, we show that the two algorithms lead to similar performances in reconstructing both variables and allow the flux divergence anomalies to be significantly reduced. 10


Introduction
Robustly reproducing the responsible mechanisms and forecasting the ice-sheets contribution to 21st century sea-level rise is one of the major challenges in ice-sheet and ice flow modelling as highlighted by community-organized efforts such as SeaRISE (Sea-level Response to Ice Sheet Evolution) (Bindschadler et al., 2013;Nowicki et al., 2013a, b) or ice2sea (e.g., Gillet-Chaulet et al., 2012; 15 Shannon et al., 2013;Edwards et al., 2014).Such projections on decadal timescales are sensitive to the model initial state which can account for an important source of uncertainty in the model response (Aðalgeirsdóttir et al., 2014).Improving the reliability of the model projections requires the model initial state to be better constrained from observations.The problem is that observations are often uncertain, sparse in time and space et al., 2012).This allows admissible flux divergence rates to be reached but at the expense of the accuracy of the modelled surface elevation and surface velocities which can then depart significantly from observations after the relaxation step.
Among the remaining uncertainties, one of the most important is the uncertainty related to the bedrock elevation.The basal topography is derived from ice thickness measurements, mostly ob-35 tained from airborne ice penetrating radars.These measurements can have large uncertainties and are usually at a lower resolution than required model grids (Durand et al., 2011).Standard bedrock elevation maps for Antarctica and Greenland are then produced by interpolation or Kriging and report standard errors ranging from few tens of meters to several hundreds of meters depending on the distance to observations and local topographic variability (Fretwell et al., 2013;Bamber et al., 40 2013).For comparison the uncertainty on the surface elevation is usually one order of magnitude lower (Fretwell et al., 2013).
Because of theses large uncertainties, several methods have been proposed to consider the bedrock elevation as an optimisation variable.For example, Morlighem et al. (2011) derive the adjoint of the continuity equation for the ice thickness.The depth-averaged velocities and surface mass balance are 45 the optimised to minimise the mismatch between modelled and measured ice thicknesses.Surface velocity measurements are used as initial guess for depth-averaged velocities, and by construction the flux divergence produced by this approach is in equilibrium with the prescribed surface mass balance.However, there is no constrain that the optimised velocities are solution of the stress equilibrium equations, so that there is no guarantee that a mechanical flow model can match the optimised 50 divergence field.
van Pelt et al. (2013) develop an iterative algorithm where the discrepancy between the surface elevation predicted by the model and the observations is used to correct the bedrock elevation.So, the method does not rely on the accurate computation of the derivative of a cost function as in a control method and is then more similar to nudging methods that have been widely studied in 55 the past decades in meteorology (e.g., Hoke and Anthes, 1976) and later in oceanography (Verron, 1992;Blayo et al., 1994).However, the method does not use observed surface velocities to control the model parameters.
(2003) developed a least-squares inversion using analytical solutions for the transmission of small scale basal perturbations to the ice surface.This method has been extended in a non-linear Bayesian framework by Raymond and Gudmundsson (2009) and applied to an Antarctic ice stream by Pralong and Gudmundsson (2011).Bonan et al. (2014) have tested the performances of an ensemble Kalman Filter on twin experiments using a shallow ice flow line model.The adjoint method has been tested 65 by Goldberg and Heimbach (2013) and Perego et al. (2014) with models of different complexity.
All these methods usually show good performances in reconstructing both basal friction and basal topography when using observations of both surface elevation and surface velocities, so that mixing between the two variables does not seem to be too problematic for realistic applications (Gudmundsson and Raymond, 2008).In addition, Pralong and Gudmundsson (2011) and Perego et al. (2014)

70
show better performance when the rates of surface elevation change are also constrained from observations.
In this paper, we explore two different algorithms to infer both the basal friction and the basal topography and initialise the model state using simultaneous observations at a given time.The first algorithm is in line with Goldberg and Heimbach (2013) and Perego et al. (2014) since it uses the 75 adjoint solution of the force balance equation.We use the Shallow Shelf Approximation to facilitate the derivation of the adjoint.Indeed, in this case the ice thickness appears as a state variable, while it changes the geometry of the domain for higher order approximations (Perego et al., 2014).In its simplest formulation, the algorithm minimises the misfit between model and observed surface velocities, but an additional constraint where the flux divergence is close to a given surface mass 80 balance can be added.The second method is an algorithm combining inversion of basal friction using the adjoint method and nudging of the bedrock topography.The control from the surface velocity observations is imposed by the adjoint step while the nudging step allows to decrease the discrepancy between the flux divergence and the surface mass balance.The main motivation of this second algorithm is its ease of implementation as no inversion of the model with respect to 85 the ice thickness is required.To our knowledge, this second method has never been explored to reconstruct the bedrock topography in a systematic manner.Both algorithms are implemented in the finite element ice-sheet / ice flow model Elmer/Ice (Gagliardini et al., 2013).The methods and algorithms are described in details in Sec. 2. To test their performances, we design a twin experiment in Sec. 3. The results are discussed in Sec. 4. 90 2 Methods

Direct Model
For the force balance, we use the standard vertically integrated Shallow-Shelf Approximation (SSA) equations (MacAyeal, 1989).This approximation neglects the effects of vertical shearing and is, hence, more adapted to model the flow in areas where the friction is low, resulting in an ice motion where a is the surface mass balance and accumulation/ablation at the bedrock interface is neglected.

Inverse methods
The objective of the methods is to produce a model state from Eq. ( 1) that best fits the observations of surface velocities and the rates of change of ice thickness.To minimise the discrepancy between the model and the observations, the optimisation parameter space p contains both the basal friction coefficient β and the bedrock elevation z b . 115

Cost Functions
The misfit between the model and the corresponding observations is evaluated using cost functions.
The first cost function measures the difference between modelled (u) and observed (u obs ) surface velocities : where Γ is the model domain.
The second cost function measures the misfit between modelled and observed thickness rates of change: The modelled rate of change of ice thickness ∂H/∂t is evaluated from Eq. (3) as the difference 125 between the flux divergence solution of Eq. ( 1) and the prescribed surface mass balance.Observed rate of change of ice thickness (∂H/∂t) obs can be estimated from surface elevation trends extracted from radar altimetry measurements (Flament and Rémy, 2012).
The objective is then to find the parameter space p that minimises J v and J div .This can be achieved in different ways as illustrated in the following sections.

Adjoint method
The two cost functions have an implicit dependence on the parameter space p through the model surface velocities u which are solution of Eq. ( 1).The gradient of the cost functions with respect to p can be computed efficiently using the adjoint equations of Eq. (1).The derivation of the continuous adjoint equations and the gradient of J v with respect to the friction coefficient β can be found in 135 Macayeal (1993).This can be easily extended for the computation of the gradient with respect to H.
The implementation in Elmer/Ice is carried out in a way that stays as close as possible to the differentiation of the discrete implementation of the direct equations.This method should lead to a better accuracy on the gradient computation than the discretisation of the continuous equations.It is similar to Martin and Monnier (2013) except that the iterative algorithm used to solve the non-140 linearity of the direct equations, arising from the explicit dependence of the ice viscosity on the velocity field, is not reversed in our adjoint computation.However as our direct solver is equipped with a Newton linearisation of the ice viscosity it remains self-adjoint (Petra et al., 2012).The other parts of the discrete source codes have been differentiated manually.
Inverse problems are often ill-posed leading to instabilities.It is then necessary to add regular-145 isation terms to the cost function to avoid over-fitting of data.This can be done in the form of a Tikhonov regularisation.Here we define two different regularisations.The first one measures the norm of the first spatial derivative of p, thus allowing to give preference to smooth solutions: The second forces the optimisation variables to stay close to a certain prior or background infor-150 mation p b .This background can be based on observations or on empirical knowledge.This second regularisation term is written as where σ p is a spatial parameter allowing to give more or less weight to the prior information p b .
The computation of the gradients of these two functional with respect to p is trivial.How these 155 regularisation terms are weighted with respect to the model-data misfit functionals Eq. ( 4) and Eq. ( 5) is described in more details with the description of the algorithms in Sec.2.3.
Using the gradients, the minimisation of the functionals is achieved using the quasi-Newton routine M1QN3 (Gilbert and Lemaréchal, 1989) implemented in Elmer/Ice.This method uses an approximation of the second derivatives of the cost function and is therefore more efficient than a 160 fixed-step gradient descent.

Nudging
By definition, the steady-state solution of Eq. ( 3) where a is replaced by the apparent mass balance a − (∂H/∂t) obs is the minimum for J div .Running the model forward in time with a constant forcing is then a simple way to minimise J div , equivalent to a relaxation step.Here we assume that 165 the surface elevation is known so that computed changes in ice thickness are used to correct the bedrock elevation z b .During this process, the ice thickness can substantially deviate from observations.Nudging methods, also called Newtonian relaxation, can remedy to this problem by constraining the thickness to fit observations through an additional call-back term in Eq. ( 3), which now writes: where the coefficient k defines the amplitude of the call-back in each node of the model.These methods imply a trade-off, adjustable through k, between model physics and observations.The call-back term can depend on many different criteria such as observation accuracy or distance to observation (Hoke and Anthes, 1976).Here, we take k as a Gaussian function of the distance to the closest ob-175 servation so that the call-back is maximum where an observation is available and decreases to zero far from all observations.

Algorithms
From the methods presented in the previous section we design two algorithms to infer simultaneously the friction coefficient β and the bedrock elevation z b .To ensure that the friction coefficient remains 180 positive during the inversion we use the following change of variable β = 10 α (9)

Adjoint method with two parameters (ATP )
This algorithm uses the gradients of the cost functions derived using the ajoint method to optimise both α and z b .For the regularisation, a constraint on the smoothness is imposed for α using Eq. ( 6 while a constraint on the background information is imposed for z b using Eq. ( 7).The total cost function then writes where γ is a constant fixed to give a similar weight to J v and J div while λ α and λ zb are two constants allowing to adjust the weight given to the regularisation terms.Following Fürst et al. (2015), several couples (λ α ,λ zb ) are tested using a L-curve approach, and optimal values are taken from the combinations that avoid two extremes: overfitting of the observations or excessive regularisation.

Adjoint-nudging coupling (ANC )
In this algoritm, the adjoint method is first used to optimise α only by minimising the following total cost function The bedrock elevation is then updated using the nudging method by solving Eq. ( 8) for a given time period T .T should be neither too short nor too long to allow to reduce J div without over-fitting observations.The sensitivity of the method to T is discussed in the results section.
These two steps are then repeated iteratively until changes in J v and J div between two iterations 200 are less than 1%.

Manufactured data sets
A twin experiment is design to investigate the ability of the two methods to reproduce simultaneously good estimates of the basal friction coefficient and the bedrock elevation.A flowline geometry is preferred to reduce the computational cost and easily test the method, however all the algorithms can 205 be applied to 2D plane view simulations.A reference experiment for which all the model parameters are prescribed is produced to generate synthetic observations.These observations are then used to test the performances of the two algorithms.

Reference experiment
A flowline of Jackobshavn Isbrae, Greenland, is used to test the two algorithms with realistic con-210 ditions.Jackobshavn Isbrae is one of Greenland's three largest outlet glaciers and has one of the largest drainage basin on the ice sheet's western margin (Bindschadler, 1984).It is also the fastest Greenland glacier with a terminus velocity greater than 13 km a −1 (Joughin et al., 2008(Joughin et al., , 2014)).The flowline is 550 km long and runs from the ice divide to the ice front.The surface and bedrock elevations are taken from available digital elevation models (Bamber et al., 2013).The basal friction The steady state solution is used as the reference of the twin experiment.

220
The geometry is discretised through a mesh of 500 linear elements, increasingly refined to the front of the glacier.The element size decreases from ∼ 2 km in the upper part of the glacier to ∼ 400 m down to the front.
Results will only be presented on the first 100 km upstream of the glacier front where velocities are above 100 m a −1 and where the SSA is more appropriate but the inversion is done all along the 225 flowline up to the ridge.

Synthetic observations
Synthetic observations are generated by sampling and/or adding noise to the reference simulation.
Details for each required field are given below.These synthetic observations and initial fields for the inverse methods are compared to the reference in Fig. 1.The surface mass balance, a, and thickness rate of change, (∂H/∂t) obs , in Eq. ( 5) are assumed to be perfectly observed.As the reference simulation corresponds to a steady-sate, (∂H/∂t) obs = 0.
However, the methods are also tested in Sec.4.3 for cases where (∂H/∂t) obs = 0 to show their ability to initialise the model when the flux divergence is not in equilibrium with the surface mass 240 balance.

Surface elevation
The surface elevation is assumed to be perfectly observed.

Bedrock elevation
For the bedrock elevation z b , we simulate observations representing airbone radar measurements observation points of the entire flowline.This error is similar to the errors given in practice on recent bedrock elevation maps (Fretwell et al., 2013;Bamber et al., 2013).For the mesh nodes between the observations, the bedrock is linearly interpolated as shown in Fig. 1a.This is used as the first guess 250 for the inverse methods and as the background information for the regularisation, Eq. ( 7).

Model parameters
The ice viscosity is assumed to be perfectly known and corresponds to the viscosity used in the reference experiment.
Assuming that no observation of the friction coefficient is available, an initial solution has to be 255 postulated.A good first guess for β is provided by using the driving stress to estimate the basal shear stress : where H(x), θ(x), and u(x) are respectively the ice thickness, the surface slope and the surface velocity at position x.The reference and initial values are shown in Fig. 1b.

260
The rms errors on the surface velocities and the rate of change of ice thickness between the initial state and the synthetic observations are, respectively, 761 m a −1 and 357 m a −1 , while the relative mean error on the basal shear stress is 394 %.The performances of the two algorithms in reducing these initial errors are presented in the following section and will be compared to these initial errors.

Adjoint method with two parameters (ATP )
A set of 255 couples (λ α ,λ zb ) is tested to adjust the weighting of the regularisation terms of Eq. ( 10).
The constant γ is fixed to 1 since J v and J div have the same order of magnitude.
The optimisation of both α and z b simultaneously allows to reach a rms misfit of 49.7 m a −1 on ( 2010), but with a greater risk of fitting noise since the relative observation error is higher in these regions.

285
In order to assess the influence of accounting for J div on the method, the optimisation is repeated without the J div term in the total cost function Eq. ( 10).The couple (λ α , λ zb ) is kept equal to the previous case since the optimum is hardly affected by the absence of J div in the total cost function.The result is given and compared to the previous one in Fig. 2. The friction coefficient β is again pretty well reconstructed, with a corresponding relative misfit of 31 % on basal shear stress to be compared 290 to the 25 % obtained with the optimisation of the J div term.However, z b shows non-consistent high frequencies involving a higher discrepancy with respect to the reference bedrock elevation that in the case of optimising J div .Therefore, the optimisation J div has a clear regularisation effect on the parameter z b , by reducing the non-consistent high frequency oscillations of the solution.
Introduction of a Gaussian noise on (∂H/∂t) obs has been experienced in order to assess its effect 295 on the optimisation of J div .Different levels of standard deviation σ have been tested.Experiences show that the optimisation is little affected by this noise even for standard deviations σ going up to the same order of magnitude than the surface accumulation a. Introduction of systematic bias on (∂H/∂t) obs in a physically acceptable range, i.e. of the same order of magnitude than surface accumulation a, also have little consequences on the optimisation.300

Adjoint-nudging coupling (ANC )
The steps consisting in the optimisation of α only are conducted with a value λ α = 5.10 9 and a value γ = 1.In order to find the optimal time period for the nudging steps, different T are tested looking at the resulting convergence on J v and J div after each cycle.The convergence is substantially similar for T from 1 to 4 years.Longer periods mainly involve a worse minimisation of J v since there is

Further sensitivity experiments
In order to evaluate the efficiency of both algorithms in transient states, we construct new reference cases where (∂H/∂t) obs = 0.This is achieved by multiplying β ref by a factor 2, 3, 4, 5 and 10.Therefore, these increases in the basal friction involve a disequilibrium of the glacier, an ice thickening and a decrease of ice flow velocities.

330
The time period for the glacier to come back to equilibrium, after this change of friction parameter, depends on the amplitude of the perturbation.Here, the perturbation is only applied during 5 years in order to keep the five cases in disequilibrium.Resulting thickness rates of change (∂H/∂t) obs = 0 are in the same order of magnitude than the tuned surface mass balance a.The five new reference cases are presented in Fig. 4.

335
The results of the optimisations for the five cases of perturbation are shown in Fig. 5 for ATP and Fig. 6 for ANC .The velocity misfit for ATP increases with the amplitude of perturbation with rms values between 47.8 m a −1 and 52.3 m a −1 while the rms misfit for thickness rate of change increase from 12.7 m a −1 to 21.8 m a −1 .ANC reaches rms misfit from 45.9 m a −1 to 47.2 m a −1 for velocities and 12.5 m a −1 to 21.5 m a −1 for thickness rate of change.The friction coefficient 340 β is well reconstructed for both methods.The corresponding average relative error (with respect to each reference) on basal basal shear stress varies from 22 to 30 % for ANC and 20 to 28 % for ATP , still according to the amplitude of the perturbation.Both algorithms also allow to improve the knowledge of the bedrock elevation z b with regard to the first guess.We notice a tendency to overestimate the amplitude of bumps and pits in some locations which generally corresponds to an 345 under estimation in the amplitude of variations of β.This last behaviour highlights the limits of the algorithms and the difficulty of distinguishing the effects of two basal parameters as closely linked as the friction and the bedrock topography.This behaviour had been already highlighted in Goldberg and Heimbach (2013) and Gudmundsson and Raymond (2008) where a higher ice thickness with respect to the reference is compensated by a higher basal friction and conversely.ANC and ATP initial states involve thickness rates of change much closer to zero than the optimi-365 sation of "β only".This also leads to a lower mismatch ε s on the surface elevation with respect to the reference after 10 years of simulation.Indeed, this mismatch is well below 20 m for both ANC and ATP, except on a few kilometres in the upstream region, whereas the optimisation of β only gives rise to a mismatch globally above 20 m with some regions exceeding 50 m.
In that way, the two algorithms implemented in this study show substantial improvements com-370 pared to the optimisation of "β only".We especially notice a better reproduction of low scale variations of the surface elevation due to the transfer of similar variations from the bedrock elevation z b (Fig. 8).These variations tend to disappear with the optimisation of the friction only, giving rise to a lower resolution of the surface.However, we should point out that this direct transfer of bedrock variations to the surface is a consequence of the SSA ice-flow model used and that a full-Stokes 375 model would produce a more diffusive transfer response .

Conclusions
The presented algorithms allow the reconstruction of two poorly known parameters: the bedrock topography z b and friction coefficient β at the same time.
215coefficient field is first adjusted so that the model velocities fit observed velocities(Joughin et al.,  Geosci.Model Dev.Discuss., doi:10.5194/gmd-2016-7,2016   Manuscript under review for journal Geosci.Model Dev.Published: 28 January 2016 c Author(s) 2016.CC-BY 3.0 License.2010).To have realistic thickness rates of change, the free surface is relaxed to steady state.The surface mass balance a in Eq. (3) has been calibrated so that the steady state is close to the initial geometry, and is meant to take into account the flow convergence or divergence along the flowline.
Surface velocities are assumed to be observed at the same resolution as the reference simulation but with a white Gaussian noise with a mean µ = 0 and a standard deviation σ = 50 m a −1 .This corresponds to a root mean squared error of 47.8 m a −1 for the entire flowline.The reference and noisy observed surface velocities are shown in Fig.1ctogether with their absolute difference.235 3.2.2Surface mass balance and thickness rate of change 245 crossing the flowline.Bedrock elevations are sampled every 10 km with a Gaussian noise centred on zero and with a standard deviation of σ = 50 m.This leads to a rms error of 62.4 m on the 55 Geosci.Model Dev.Discuss., doi:10.5194/gmd-2016-7,2016 Manuscript under review for journal Geosci.Model Dev.Published: 28 January 2016 c Author(s) 2016.CC-BY 3.0 License.
similar to the observation rms error, showing no overfitting of velocity data.The rate of change of ice thickness misfit is also largely decreased with a rms value of 19.2 m a −1 .The resulting β and z b as well as the misfit for the surface velocities are given in Fig. 2. The friction variability is accurately reproduced with a corresponding relative misfit of only 25 % for the entire flowline with respect to the reference basal shear stress, i.e. more than 10 times less than the misfit 275 obtained with β ini .We only notice local over-estimations of slipperiness in bedrock pits without significant impacts on the flow velocities.Indeed, under a defined value of β corresponding to a nearly perfectly sliding case, an additional reduction in friction has no impact on the flow.The same reasoning applies to a nearly perfectly sticky case, where an increased friction would not involve more decrease of the velocity.The bedrock elevation z b is well reconstructed in the first 280 Geosci.Model Dev.Discuss., doi:10.5194/gmd-2016-7,2016 Manuscript under review for journal Geosci.Model Dev.Published: 28 January 2016 c Author(s) 2016.CC-BY 3.0 License.50 km upstream of the glacier front.The discrepancy with respect to the reference bedrock is larger upstream where the cost function J v is less sensitive because of lower velocities.This could possibly be improved by using a cost function measuring the logarithm of the misfit as in Morlighem et al.

305
no control on velocities during nudging.Shorter relaxation times do not involve sufficient change of z b inducing a lower minimisation of J div for a given number of cycles.Therefore, a relaxation time T of 1 year is adopted, which seems sufficient to allow significant changes of z b without too much adaptation to the previous intermediate value of the friction coefficient.The algorithm is stopped after 10 cycles, corresponding to the stopping criterion of Sec.2.3.2.310 The model is in good agreement with observations with a rms misfit of 46.1 m a −1 in the range of observation noise, for velocities and 15.8 m a −1 for thickness rates of change.The friction coefficient β is close to the reference one despite exacerbated variations at some locations.The corresponding relative misfit on basal shear stress is 30 % for the entire flowline.The reconstructed bedrock elevation z b is also close to the reference on almost 100 km upstream of the front of the glacier.This 315 reflects, especially in fast flowing region, a real improvement of the basal knowledge with respect Geosci.Model Dev.Discuss., doi:10.5194/gmd-2016-7,2016 Manuscript under review for journal Geosci.Model Dev.Published: 28 January 2016 c Author(s) 2016.CC-BY 3.0 License.to the first guess.Moreover, the use of nudging, instead of adjoint method, does not show the same problem of non-sensitivity in region of slower flow velocities, as mentioned in Sec.4.1.Note however that z b significantly departures from the reference bedrock elevation in the 80 to 100 km to the front region, strongly linked to the poorer fit of β (see Fig.3).320As for ATP, introduction of a Gaussian noise in the observed thickness rate of change (∂H/∂t) obs has also been tested.Results show no significant impacts on the optimisation.Nevertheless, introduction of systematic bias in (∂H/∂t) obs has direct consequences on the nudging steps inducing an offset of z b of the range of the systematic bias cumulated on the nudging period T .ANC is therefore more sensitive to systematic bias than ATP. 325 350

4. 4
Flow divergence in transient modelIn this section, we assess the impact of our initialisation algorithms on the prognostic response of the model forward in time assuming the same constant forcing used to build the reference state.By doing so, if the initialisation was perfect, one would expect no change of the geometry and ice flow during this prognostic simulation.The experiment is performed from ATP and ANC initial states.355Athird initialisation state is constructed for which only the friction coefficient has been optimised, keeping z b equals to the a priori z bb .This third initialisation, called "β only" involves a rms misfit on velocities of 43.3 m a −1 and a relative mean error of 36 % on basal shear stress, similar to the ATP and ANC initial states.However, the rms misfit on the thickness rate of change is significantly higher, 147.8m a −1 .360Theprognostic simulations are conducted on a 10 years period in order to see how the initial thickness rate evolves during this time and how it impacts the final ice thickness and ice surface.The thickness rate of change after 1 and 10 years of simulation are shown in Fig.7aand b respectively, while the mismatch on the surface elevation after 10 years is shown in Fig 7c.

380
easier to measure: ice surface velocities and thickness rates of change.Some local measurements of bedrock elevation and associated errors are necessary in order to define a background z b .The two algorithms aim to infer the set of parameters which minimises the misfit between the model and the corresponding observations of ice surface velocities and thickness rates of change.If the optimisation of ice surface velocities is usually sufficient to infer β, the inference of a second parameter requires 385 12 Geosci.Model Dev.Discuss., doi:10.5194/gmd-2016-7,2016 Manuscript under review for journal Geosci.Model Dev.Published: 28 January 2016 c Author(s) 2016.CC-BY 3.0 License.more information to distinguish the effects of each parameters on the flow.Observations of rates of change of ice thickness are necessary to allow optimising z b as well.The two algorithms are based on the optimisation of the friction coefficient β with the adjoint method.The bedrock geometry z b is reconstructed in two different ways, again with adjoint method for the first algorithm (ATP) and with a nudging method based on mass conservation equation for 390 the second one (ANC).We have shown that ATP algorithm reproduces β well and the corresponding basal shear stresses while z b is particularly well reproduced in high velocities regions but does not depart from the background when the velocities become lower.The iterative algorithm coupling adjoint method and nudging (ANC) gives as good results.Moreover, ANC allows a better reconstruction of the bedrock 395 geometry z b in most regions.This is a very good sign for an adaptation of the method non depthintegrated flow models such as full-Stokes model.Indeed, there is no need to inverse a shape variable like bedrock topography which is a usual obstacle to adjoint-based methods.Furthermore, the transient simulations over 10 years from initial states reconstructed with the two algorithms developed give very encouraging results.The model divergence is clearly decreased with 400 respect to usual inversion methods of the friction coefficient only.The integration of observations like thickness rates variation through an optimisation of the divergence during inversion or nudging steps, allows to regularise the solution in a physical way and also clearly improves the results.Finally, the sensitivity experiments shows that the different algorithms can take into account the disequilibrium of mass balance, which is particularly interesting considering that a large amount of 405 outlet glaciers in both Greenland and Antarctica present this feature.13 Geosci.Model Dev.Discuss., doi:10.5194/gmd-2016-7,2016 Manuscript under review for journal Geosci.Model Dev.Published: 28 January 2016 c Author(s) 2016.CC-BY 3.0 License.Geosci.Model Dev.Discuss., doi:10.5194/gmd-2016-7,2016 Manuscript under review for journal Geosci.Model Dev.Published: 28 January 2016 c Author(s) 2016.CC-BY 3.0 License.

Figure 1 .
Figure 1.Reference (solid lines) and initial (dashed lines) state for (a) the bedrock elevations zb, (b) the friction coefficient β and (c) the surface velocities.In (a) Synthetic observations every 10 km are the plain black circles.In (c) the "observed" velocities are depicted by the circles and the shaded green curve is the absolute difference between observed and reference surface velocities (right axis).

Figure 2 .
Figure 2. Results of the ATP algorithm with (orange) and without (red) optimisation of Jdiv, i.e. γ = 1 or γ = 0, respectively, in Eq. (10): (a) absolute difference between observed and model velocities, (b) estimated friction parameter and (c) estimated bedrock elevation.The green shaded area is the difference between the noisy reference velocities and the true velocities.The green solid lines are the reference values and the black dashed line is the initial guess for the bedrock elevation.

Figure 3 .
Figure 3. Results of the ANC algorithm (purple): (a) absolute difference between observed and model velocities, (b) estimated friction parameter and (c) estimated bedrock elevation.The green shaded area is the difference between the noisy reference velocities and the true velocities.The green solid lines are the reference values and the black dashed line is the initial guess for the bedrock elevation.

Figure 5 .
Figure 5. Range of values for ATP algorithm for the 5 perturbations of the friction coefficient β.(a) and (b) minimum (dark orange shade) and maximum (light orange shade) of absolute difference between observed and model velocities and relative error for β respectively.(c) range of values for bedrock elevation zb (orange shade).The green solid line is the reference value and the black dashed line is the initial guess for the bedrock elevation.

Figure 6 .
Figure 6.Same as Fig. 5 but for the ANC algorithm.

Figure 7 .
Figure 7. Evolution of ∂H/∂t after 1 year (a) and 10 years (b) of prognostic simulation and the resulting missmatch after 10 years between surfaces obtained with three different initial states and reference surface (c).The orange and purple lines give the results for ATP and ANC .The red line gives the result for inversion of β only.

Figure 8 .
Figure 8. Ice surface elevation zs after 10 years of prognostic simulation for 3 different initial states : initialisation with ATP algorithm (orange line), with ANC algorithm (purple line) and with the inversion of β only (red line).The green line is the reference surface elevation.The figure focuses on the first 50 km next to the front of the glacier.
is the friction coefficient, ν the effective mean viscosity, ρ i the ice density, g the gravity, and H = z s − z b the thickness, with z s and z b the top and bottom surface elevations, respectively.