Sensitivity of centennial mass loss projections of the Amundsen basin to the friction law

Abstract. Reliable projections of ice sheets' future contributions to sea-level rise require models that are able to accurately simulate grounding-line
dynamics, starting from initial states consistent with observations. Here, we
simulate the centennial evolution of the Amundsen Sea Embayment in response
to a prescribed perturbation in order to assess the sensitivity of mass loss
projections to the chosen friction law, depending on the initialisation
strategy. To this end, three different model states are constructed by
inferring both the initial basal shear stress and viscosity fields with
various relative weights. Then, starting from each of these model states,
prognostic simulations are carried out using a Weertman, a Schoof and a Budd
friction law, with different parameter values. Although the sensitivity of
projections to the chosen friction law tends to decrease when more weight is
put on viscosity during initialisation, it remains significant for the most
physically acceptable of the constructed model states. Independently of the
considered model state, the Weertman law systematically predicts the lowest
mass losses. In addition, because of its particular dependence on effective
pressure, the Budd friction law induces significantly different grounding-line retreat patterns than the other laws and predicts significantly higher
mass losses.



Introduction
The West Antarctic Ice Sheet mean annual contribution to global sea-level rise (SLR) has tripled over the last 25 years as a consequence of a growing imbalance between the mass it receives as snowfall and that which is discharged to the ocean by ice streams (Shepherd et al., 2018). The most ac-tive basin of this region is the Amundsen Sea Embayment (ASE), where marine-terminating outlet glaciers draining ice to the oceans have shown sustained acceleration and thinning over the last decades, with their grounding lines (GLs), i.e. the limit between the grounded ice sheet and the floating ice shelf, retreating at rates higher than 1 km a −1 since 1992 Rignot et al., 2014). These observations raise concerns regarding the near future of the ASE as they suggest that this sector of West Antarctica is undergoing a marine ice-sheet instability (MISI), which would imply a significant additional global SLR in the coming decades (Joughin et al., 2014;Favier et al., 2014;Cornford et al., 2015). Producing reliable estimates of this future contribution requires accurate modelling of GL dynamics on subcentennial timescales as well as the ability to produce model states that are as close as possible to observations. Using a synthetic flow line geometry, Brondex et al. (2017) have shown that the GL dynamics depends critically on the choice of the friction law, i.e. the mathematical relationship linking basal shear stress to other parameters including sliding velocity. The ice-bed interface being usually out of reach, the formulation of a friction law has been a long-standing problem in glaciology. Various laws intended to describe different physical processes at the roots of basal motion have been developed over the years (Weertman, 1957;Budd et al., 1979;Schoof, 2005;Tsai et al., 2015). Most ice flow modelling studies published so far use the Weertman friction law, which aims to describe the basal motion of ice over and around obstacles of a rigid bedrock by a combination of viscous creep and regelation (Weertman, 1957). However, many rapid ice streams of Antarctica are known to be lying on soft beds made of water-laden till, the deformation of which explains most of the motion observed at the surface Alley et al., 1986;Kavanaugh and Clarke, 2006). Numerous laboratory studies on till samples and in situ measurements have shown that, at large strain, the till rheology is plastic with a critical strength τ * depending on effective pressure N , i.e. the difference between ice overburden pressure and water pressure (Boulton and Jones, 1979;Blankenship et al., 1986;Alley et al., 1986). To account for both the cases of rigid and soft beds, Tsai et al. (2015) proposed a law inducing a Coulomb friction regime at low N, which instantaneously switches to a Weertman friction regime at higher N . By construction, this law induces an upper-bound function of N of the basal shear stress. Although it was originally intended to describe the ice flowing over a rigid bed with the opening of water-filled cavities, the law proposed by Schoof (2005) behaves similarly, except that the transition between the Coulomb and Weertman regimes is continuous, which makes it easier to handle numerically (Brondex et al., 2017).
Although a growing amount of information about the current state of the ice-sheet (e.g. surface velocities, surface elevation, surface elevation rates of change) is made available by the rapid development of satellite observations, several model parameters remain poorly constrained (e.g. bedrock elevation, ice viscosity, friction law coefficients). Gradientbased optimisation methods are routinely used to estimate uncertain model parameter fields and boundary conditions so that the initial ice-sheet geometry and surface velocity field are as close as possible to the observations (e.g. Morlighem et al., 2010Morlighem et al., , 2013Gillet-Chaulet et al., 2012;Cornford et al., 2015). Although such methods enable us to deduce the current basal shear stress field from observed surface velocities, the form of the friction law cannot be discriminated with a unique set of observations (Joughin et al., 2004;Gillet-Chaulet et al., 2016). Furthermore, when several uncertain fields are simultaneously inferred, several different initial states consistent with observations can potentially be constructed. Adhalgeirsdóttir et al. (2014) have shown that SLR projections on decadal timescales are sensitive to the initial state of the model, which can account for an important source of uncertainty in the model response.
In the present study, we aim to assess the relative sensitivity of centennial mass loss projections of the ASE to the chosen friction law and initialisation strategy. Our work being based on a prescribed perturbation scenario, the results presented here should not be considered as actual projections of the future contribution of the ASE to SLR. The first step consists of building three different model states of the ASE by simultaneously inferring the basal shear stress and the ice viscosity, with various relative weights attributed to each one of these two fields during the inversion. Then, for each of these states, we follow the same procedure as in Brondex et al. (2017) to identify the distributions of the friction coefficients of three commonly used friction laws that lead to the same model initial states. Finally, we apply a synthetic per- turbation to the basal melting rate under floating ice, to the different initial states and compare the dynamical responses obtained with the various friction laws. In Sect. 2, we give a precise description of the model used to conduct this study and describe the experimental setup. The results obtained at each step of the experiments are presented in Sect. 3 and discussed in the last section.

Model description
The modelled domain is represented in Fig. 1. For the stress balance, we solve the two-dimensional Shelfy-Stream Approximation (SSA) equations (MacAyeal, 1989) The Cryosphere, 13, 177-195, 2019 www.the-cryosphere.net/13/177/2019/ where ρ i is the ice density, g is the gravity norm and H = z s − z b is the ice thickness, with z s and z b the top and bottom surface elevations. The vertically averaged effective viscosity η reads as follows: where D e is the second invariant of the strain-rate tensor, n is the Glen's law exponent and η 0 is given by In Eq. (3), A is the rate factor. It is related to the temperature relative to the pressure-melting point T ( • C) through an Arrhenius law: where A 0 is the pre-exponential factor, Q is an activation energy and R is the gas constant (Cuffey and Paterson, 2010). For the basal shear stress τ b in Eq. (), we consider three different friction laws: where C W , C B and C S are friction parameters. Equation (5) corresponds to the widely used Weertman law (Weertman, 1957), where m is a positive exponent, often related to the creep exponent n of the Glen's law as m = 1/n. Equations (6) and (7) correspond to the Budd and Schoof friction laws, respectively (Budd et al., 1979;Schoof, 2005). Note that, on the contrary to the two other laws, the Schoof law induces an upper bound of the ratio |τ b |/N known as Iken's bound that is equal to the value of the parameter C max (Iken, 1981;Schoof, 2005;Gagliardini et al., 2007). Although it is not its original purpose, the use of a Schoof law -which can be seen as the equivalent of the Tsai law (Brondex et al., 2017) -is more and more justified by its ability to represent the deformation of till, as it induces a Coulomb friction regime at low N. In that case, the parameter C max ought to be seen as a Coulomb friction parameter related to the rheological properties of the till. Laboratory measurements have shown that the values of C max should then range between 0.17 and 0.84 (Cuffey and Paterson, 2010). Both the Budd and Schoof laws make explicit use of the effective pressure N. Using an effective pressure-dependent friction law when running transient simulations proves necessary in order to account for the enhancement of basal sliding due to the presence of subglacial drainage systems connected to the ocean, which have been reported in several studies (Gray et al., 2005;Fricker et al., 2007;Le Brocq et al., 2013). However, unlike the Schoof law, for which the dependence of |τ b | on N is limited to the close vicinity of the GLs, the Budd law induces a dependence of |τ b | on N over the whole grounded part of the domain, even far into its interior where such a drainage system may not exist. Ideally, N should be computed by a subglacial hydrology model but the few available models of that kind usually rely on a multitude of poorly constrained parameters (Schoof, 2010;Hewitt et al., 2012;Werder et al., 2013;de Fleurian et al., 2014) and the attempts to couple these models to ice flow models are relatively recent (Hewitt, 2013;Bougamont et al., 2014;Bueler and van Pelt, 2015;Gagliardini and Werder, 2018). Therefore, in the present study, N is calculated assuming a perfect hydrological connection between the subglacial drainage system and the ocean, so that where ρ w is the water density. A systematic comparison of the dependence of τ b on N and u b implied by these various friction laws is available in Brondex et al. (2017). Note that, since floating ice does not experience any friction, τ b is set to zero wherever ice is afloat, independently of the chosen friction law. The temporal evolution of ice thickness is governed by a two-dimensional mass transport equation: where a s is the surface mass balance and a b the oceanic melt rate applied to the bottom surface of the ice shelves only, basal melt at the ice-bed interface being neglected. Ice is assumed to be in hydrostatic equilibrium and the bottom surface elevation can be deduced from the bedrock topography b(x, y) by applying the no-penetration condition and the floating condition. Here we assume a constant sea level, z sl = 0, so that The GL being the limit beyond which grounded ice starts floating, its position (x G , y G ) can directly be deduced from Eq. (10) by solving the following: Therefore, the GL can be located at any point of the domain and its position is free to evolve over time as a result of the evolving geometry. The positions of the boundaries of the model domain are set based on the observations made by the satellite ICESat for the IMBIE2 effort (Zwally et al., 2012). In the following, n = (n x , n y ) is the normal unit vector to the considered boundary.
www.the-cryosphere.net/13/177/2019/ The Cryosphere, 13, 177-195, 2019 At the ice divides, there is no ice flux entering the domain, so the following Dirichlet condition applies: and is completed by a free slip condition in the tangential direction. At the ice shelves and glaciers fronts, which remain fixed in time, the following Neumann condition applies: where H sub | f is the submerged height at the considered front, which is null in the very few places where glacier fronts are grounded and relates to H | f through the floating condition at floating fronts. Similarly to Pollard and DeConto (2012), the oceanic melt rate a b (m a −1 ), prescribed on the bottom surface of ice shelves, is parameterised as follows: where T 0 and T f are the ocean water temperature and freezing point at the considered depth, K T is a transfer coefficient for sub-ice oceanic melting, C W the specific heat of ocean water and L f the latent heat of fusion of ice. Following Pollard and DeConto (2012), the temperature difference T 0 − T f is given by T 0 − T f = 0.5 • C for 0 > z > −170 m, 3.5 • C for z < −680 m and linearly interpolated for z between −680 and −170 m. Values of parameters prescribed in this study are presented in Table 1. For the surface mass balance a s , we use outputs of simulations performed with the MAR model and averaged over the 1979-2015 period (Cécile Agosta, personal communication, 2016). The bedrock topography is taken from Bedmap2 (Fretwell et al., 2013), except that we include two pinning points in contact with the bottom surface of Thwaites ice shelf using the bathymetry of Millan et al. (2017). The mechanical role of these two pinning points has indeed been shown to be critical because of the buttressing effect that they exert on the upstream ice stream (Fürst et al., 2016).
All the equations presented above are solved using the open-source finite element code Elmer/Ice (Gagliardini et al., 2013). The mesh is generated using the anisotropic meshadaptation technique described in Gillet-Chaulet et al. (2012) so that the distance between two nodes at the GL is of the order of ∼ 50 m. Finally, we end up with a mesh made of 274 678 nodes covering the whole model domain, the total surface of which is about 4.5 × 10 5 km 2 . The same mesh is

Parameters
Values Units used for all the successive steps of all the experiments presented in the following. In addition, we take advantage of the fact that the GL position is determined at subelement precision through Eq. (11) to make use of the subelement parameterisation SEP3 as proposed by Seroussi et al. (2014). In our case, the number of quadrature points is raised to 20 in the elements containing the GL.

Experimental setup
In the present study, we aim to compare the results in terms of GL dynamics and volume losses produced by the following friction laws: (1) a linear Weertman law given by Eq. (5) with m = 1, (2) a non-linear Weertman law given by Eq. (5) with m = 1/3, (3) a non-linear Schoof law given by Eq. (7) with m = 1/3 and C max = 0.4, (4) a non-linear Schoof law given by Eq. (7) with m = 1/3 and C max = 0.6, (5) a linear Budd law given by Eq. (6) with m = 1, and (6) a non-linear Budd law given by Eq. (6) with m = 1/3. To reach this goal, we start with an initialisation step in which three model states -hereinafter referred to as inferred states -with different initial basal shear stress and viscosity fields are constructed from available observations, before undergoing a 15-year relaxation period. Then, following the same procedure as the one described in Brondex et al. (2017), we identify the distributions of the friction coefficients of the aforementioned friction laws so that the basal shear stress field computed with these laws is the same as the one obtained at the end of the 15-year relaxation period. For the two Weertman laws and the two Schoof laws, this procedure is repeated for each of the three inferred states. In contrast, because they imply a dependence of |τ b | on N over the whole grounded part of the domain, we consider the two Budd laws as being less physically acceptable so that, for these laws, the identification procedure is performed for one of the inferred states only due to limited computational resources. Therefore, after the identification step, we are left with 14 new model states, hereinafter referred to as initial states, corresponding to four The Cryosphere, 13, 177-195, 2019 www.the-cryosphere.net/13/177/2019/ friction laws -the linear and non-linear Weertman laws and the non-linear Schoof laws with C max = 0.4 and C max = 0.6 -applied to the three inferred states as well as the linear and non-linear Budd laws applied to one of the inferred states only. These 14 initial states constitute the starting points of a set of two 100-plus-5 years of prognostic simulations, i.e. an unperturbed control run and a run, for which the oceanic melt rate is perturbed. Given the size of the domain and the mesh refinement required to capture essential flow features, each prognostic simulation has a numerical cost of around 11 000 CPU h. Figure 2 summarises all the consecutive steps of the experimental setup, which are described in detail below.

Initialisation
A number of observations are used in order to initialise the model. The initial ice thickness is taken from Bedmap2 (Fretwell et al., 2013). A two-dimensional reference field of η 0 , denoted by η 0,ref , is computed using a three-dimensional ice temperature field of Antarctica (Brice Van Liefferinge, personal communication, 2016) obtained using a thermomechanical model forced with a mean geothermal heat flux (Van Liefferinge and Pattyn, 2013). In addition, we use the observed surface velocities of  to invert for basal shear stress and ice rheology. These velocities are given on each node of a regular grid with a 450 m resolution, except in some regions, which are sometimes very wide, over which the information is missing. The estimated error in the norm of observed velocities ranges, depending on regions, between 1 and 17 m a −1 . These observations have been collected over periods spanning several years, which could potentially induce inconsistencies in regions where the flow features are evolving rapidly. The flow solution u of the SSA Eq. () depends on both the basal shear stress field and the effective viscosity field. For simplicity, the initial basal shear stress field is computed with a linear Weertman law of which we infer the friction coefficient distributionĈ W . Although η 0,ref enables a first approximation of the effective viscosity, there are significant uncertainties regarding the temperature field used to derive this field as well as the parameters linking the rate factor of ice to its temperature. In addition, ice viscosity does not depend only on temperature but also on several other parameters including damage, anisotropy, water content, density, grain size and impurity content (Cuffey and Paterson, 2010). Therefore, the question is whether the local discrepancies between observed and modelled velocities at any given point of the domain should be attributed to an inappropriate basal shear stress field or to errors on the reference field η 0,ref . Indeed, several model states consistent with observations can be constructed, for which the fit between modelled and observed velocities is obtained by adjusting the basal shear stress and the viscosity. In order to assess the sensitivity of the results to the initialisation strategy, we construct three inferred states -denoted by I R γ ,∞ , I R γ ,100 and I R γ ,1 -by means of the control method, building on the approach described in Fürst et al. (2015).
The total cost function J tot used for minimisation includes two cost functions and two Tikhonov regularisation terms: The misfit between modelled (u) and observed (u obs ) velocities is comprised in the first cost term J v . To avoid errors due to interpolation of observed velocities at the model mesh nodes, J v is a discrete sum evaluated directly at the observation points. It reads as follows: where N obs is the total number of available observations, u i,obs is the velocity observed at point i, and u i is the modelled velocity, which is interpolated at the observation point i. The second cost function J div is intended to penalise the large gaps between ice flux divergence and mass balance, leading to inferred states closer to steady states. It reads as follows: where is the model domain. During the inversion, we actually optimise the variables α and γ , which are related to the linear Weertman law coefficient and η 0 , respectively, as follows: and where R γ is a constant. The variable changes (Eq. 18) and (Eq. 19) prevent non-physical negative values of basal shear stress and viscosity, respectively. In addition, the variable change (Eq. 19) enables us to tune the relative weight that will be put on basal shear stress and viscosity during inversion. Indeed, the minimisation of the cost functions relies on the calculation of their gradients with respect to the variables to optimise (MacAyeal, 1993). As a consequence, the higher the value attributed to R γ , the lower the gradients of J v and J div calculated with respect to γ relative to the ones calculated with respect to α. In such a case, the distribution of α -and thus of the basal shear stress -will be more affected (on the grounded part only as τ b is forced to zero wherever ice is floating) than for lower values of R γ . The two regularisation functions J reg,α and J reg,γ in turn penalise first spatial derivatives of α and γ , respectively. They are meant to avoid www.the-cryosphere.net/13/177/2019/ The Cryosphere, 13, 177-195, 2019 Figure 2. Flow chart summarising the consecutive steps of the present study for a given inferred state. The simulations EXP_CONTROL and EXP_ABMB are run for each of the five friction laws. The same procedure is repeated for each of the three inferred states, I R γ ,∞ , I R γ ,100 and I R γ ,1 , except that the non-linear and linear Budd friction laws are applied only to the inferred state I R γ ,100 .
overfitting the velocity observations and thus, improve the conditioning of the problem. Among the three inferred states considered in this study, the two states I R γ ,100 and I R γ ,1 are constructed by optimising both α and γ with R γ = 100 and R γ = 1, respectively, in Eq. (19). Hence, more weight is put on basal shear stress to the detriment of viscosity for I R γ ,100 than for I R γ ,1 . In contrast, the inferred state I R γ ,∞ is obtained by optimising α only, while η 0 = η 0,ref . The gradients of J tot are derived following the adjoint method (MacAyeal, 1993), while the minimisation itself is done using the quasi-Newton routine M1QN3 (Gilbert and Lemaréchal, 1989). The minimisation method being an iterative method, we need to provide initial guesses for the distributions of the variables to optimise. Following Morlighem et al. (2013), the initial guess for α is found by considering that, as a first approximation, the driving stress is exactly compensated by the basal shear stress at any given point of the ice-bed interface. Regarding the initial field of η 0 , it is simply set to η 0,ref for both I R γ ,100 and I R γ ,1 .
Although strong discrepancies between mass balance and ice flux are penalised through J div in Eq. (15), the three inferred states are not exactly steady states. Seroussi et al. (2011) and Gillet-Chaulet et al. (2012) have reported the occurrence of non-physical ice thickness rates of change in the first years following the inversion. These ice flux anomalies are due to uncertainties in the model initial conditions, in particular regarding the prescribed topography and the values attributed to the various parameters, and are usually dissipated within a few years. For this reason, each of the three inferred states undergoes a relaxation period, the duration of which is arbitrarily fixed to 15 years. These relaxations are performed with the direct model as described in Sect. 2.1, with the basal shear stress field being computed through the linear Weertman law (see Fig. 2).

Identification of friction laws coefficients
For each of the three inferred states, the 15-year relaxation period leads to a new state characterised by a new velocity field, which we denote byû b , as well as a new basal shear stress field, denoted byτ b . These two fields relate to each other through the linear Weertman law, i.e. Eq. (5) with the inferred friction coefficientĈ W (which depends on the considered inferred state) and m = 1. As thoroughly described in Brondex et al. (2017), these two fields can be used to analytically identify the distribution of the friction coefficient of other laws, which would lead to the same reference basal shear stress fieldτ b . Thus, the distributions of the friction coefficient of the non-linear Weertman law, the two Budd laws and the two Schoof laws can be identified by respectively The Cryosphere, 13, [177][178][179][180][181][182][183][184][185][186][187][188][189][190][191][192][193][194][195]2019 www.the-cryosphere.net/13/177/2019/ solving the following: and, For the Budd and Schoof laws, the effective pressure field N used for the identification is calculated through Eq. (8) from the geometry obtained at the end of the relaxation period. Note that the fieldsĈ W ,û b ,τ b andN being dependent on the inferred state, these identifications are repeated for each one of the three inferred states, except for the linear and non-linear Budd laws for which the identifications have been done only for the case I R γ ,100 (see Fig. 2). Equations (20) and (21) enable us to identify the friction coefficients of the non-linear Weertman law and the two Budd laws at every grounded node covered with ice. At floating nodes or in the very few places of the domain which are ice-free when identification is performed, these coefficients are arbitrarily fixed to C W,nl = 10 −6 MPa m −1/3 a 1/3 and C B = 10 −6 m −1/3 a 1/3 , respectively. This choice does not appear to be critical as all the following prognostic experiments lead exclusively to a GL retreat. In contrast, because the Schoof law was constructed so that Iken's bound is satisfied (Schoof, 2005), Eq. (22) has a solution only in places where |τ b |/N < C max . Figure 3 shows, for each of the three inferred states, the percentage of grounded nodes where the value of C S cannot be identified directly from Eq. (22), depending on the value attributed to the parameter C max . Although we do not expect the parameter C max to be uniform all over the model domain, there is currently no way to constrain its spatial distribution. As a consequence, two different uniform values are considered: C max = 0.4 and C max = 0.6. Note that the identification of C S then has to be done for each of the two values of C max . Because the sensitivity of basal shear stress to C S is very small in places close to flotation, the identification of C S from Eq. (22) is only done at nodes where |τ b |/N ≤ 0.8C max . The values of C S at nodes where |τ b |/N > 0.8C max are linearly interpolated from the closest neighbouring nodes. In addition, we arbitrarily set C S = 10 −3 MPa m −1/3 a 1/3 in places which are ice-free or where ice is afloat when the identification is performed. Finally, we end up with six different distributions of C S corresponding to the two values of C max applied to the three inferred states I R γ ,∞ , I Rγ ,100 and I Rγ ,1 . For these six distributions, the percentage of grounded nodes at which C S has to be interpolated is comprised between ∼ 4 % (for the initial state I R γ ,∞ associated with the parameter C max = 0.6) and ∼ 8 % (for the initial state I Rγ ,1 associated with the parameter C max = 0.4). As expected, the nodes at which an interpolation is required are the closest to flotation and, therefore, are mostly located directly upstream of the GL.
As stated previously, the identification steps lead to 14 different initial states, which are used as the starting points of the prognostic simulations (blue dots in Fig. 2).

Prognostic simulations
The 14 initial states are then run forward in time for 105 years under two different scenarios: (1) a control run (EXP_CONTROL) and (2) a basal melt anomaly run (EXP_ABMB). The EXP_CONTROL run is an unforced forward experiment aiming to characterise model drift, which depends both on the initial state and the friction law. Therefore, all model parameters and forcing in EXP_CONTROL are the same as those used for initialisation and presented in Sect. 2.1.
The EXP_ABMB run consists of applying a synthetic anomaly of basal melting rate under floating ice, all other model parameters and forcing being kept the same as in EXP_CONTROL: after a first 5-year period free of any perturbation, a uniform basal melting rate anomaly "abmb" of 13.2 m a −1 is progressively added to the initial basal melting rate a b , given by Eq. (14), over the following 40 years of simulation. Thus, the basal melting rate a b,ABMB for this run is given by www.the-cryosphere.net/13/177/2019/ The Cryosphere, 13, 177-195, 2019 The basal melting rate anomaly abmb corresponds to the one prescribed for the ASE within the InitMIP-Antarctica framework (Nowicki et al., 2016).

Initialisation
The three inferred states I R γ ,∞ , I R γ ,100 and I R γ ,1 are compared in Fig. 4. The absolute difference between modelled u and observed u obs velocities turns out to shrink when both basal shear stress and viscosity are inferred (Fig. 4a-c). This is particularly true for the ice shelves, which do not feel any basal shear stress: although the basal shear stress directly upstream of the GL do influence velocities within the downstream ice shelf, the most efficient way to obtain a better match between modelled and observed velocities in floating areas is through a local adjustment of viscosity. Thus, the inferred state I R γ ,∞ shows important errors on Pine Island, Thwaites and, to a lesser extent, Crosson ice shelves (Fig. 4a). The root mean squared (rms) errors obtained for the inferred states I R γ ,∞ , I R γ ,100 and I R γ ,1 are 85.5, 52.4 and 37.5 m a −1 . The ratio between the norm of the basal shear stress and the norm of the driving stress is shown in Fig. 4d-f. In addition, plots of the norm of the basal shear stress obtained for the three inferred states are included in the Supplement (Fig. S1) to allow for a comparison with previous work in which inversions of basal shear stress in the regions of PIG and Thwaites were performed (Joughin et al., 2009;Morlighem et al., 2010). At the end of the inversions, there are several regions where the local driving stress is not entirely compensated by the local basal shear stress (blue regions). As expected, these regions correspond to regions of rapid flow; in particular, Pine Island Ice Stream can be easily distinguished. The regions of low τ b become narrower as more weight is put on viscosity during inversion. Indeed, when the viscosity is inferred, another way for the inversion algorithm to increase the modelled velocities in areas where they would be too low otherwise is to soften the ice locally: this is the case, for example, in the higher part of PIG, in particular for the inferred state I R γ ,1 for which the inversion algorithm induces a local reduction in viscosity rather than of the basal stress to increase the modelled velocities (panels f and i of Fig. 4). It is also this same mechanism which is behind the low viscosity bands, easily distinguishable in Fig. 4h (note that they are also present in Fig. 4i but less visible due to the wider colour scale). These soft bands correspond to shear margins, where high shear stresses induce a local reduction of ice viscosity through different physical processes, including damage, strain heating and the development of crystalline fabric (e.g. Borstad et al., 2013;Bondzio et al., 2017;Minchew et al., 2018). On the contrary, the inversion algorithm can render the ice locally stiffer in order to decrease the modelled velocities in areas where they are too high compared to the observations, which is typically the case for ice shelves.
Because the model states obtained after the initialisation step are not steady states, they drift toward new states during the following 15-year relaxation period. Overall, the three inferred states show similar evolution patterns with relatively moderate ice thickness changes, except for the Pine Island and Thwaites ice shelves, where ice thickness increases by about 100 m over the 15 years. In contrast, ice streams become slightly thinner, losing a few tens of metres of ice in some places. At the end of the relaxation period, the thickness rates of change have decreased to physically acceptable values for all the inferred states, with, at most, 30 m a −1 of absolute thickness rate of change, concentrated within very localised areas; elsewhere, ice thickness is nearly at equilibrium with, at most, a few m a −1 of absolute thickness rate of change. Although the modelled velocity structure remains similar to observations, the rms errors between observed velocities and velocities computed at the end of the relaxation period have risen to 124.8, 119.4 and 96.4 m a −1 for I R γ ,∞ , I R γ ,100 and I R γ ,1 , respectively.

Identification of friction laws coefficients
For each of the three inferred states, the basal shear stress fields computed with the non-linear Weertman law and with the two Budd laws using the friction coefficient distributions deduced through Eqs. (20) and (21), respectively, are identical to the basal shear stress fieldsτ b used for the identifications. In contrast, because the value of C S cannot be exactly identified through Eq. (22) at every grounded node but needs to be interpolated in some places,τ b is not perfectly reproduced with the Schoof law. Figure 5a-c show an enlargement of regions within which the relative differences between the basal shear stress field calculated with the non-linear Schoof law (7) associated with C max = 0.4 and the C S distributions identified through Eq. (22) or interpolated, which is denoted by τ bS04 , and the reference basal shear stress fieldsτ b are the highest. As expected, the highest relative differences are obtained close to the ice shelves (shown in green in Fig. 5ac), where the nodes are close to flotation and C S needs to be interpolated. Elsewhere, the fields of τ b produced by the two friction laws are numerically identical. The relative difference between τ bS04 andτ b at nodes where C S is interpolated rarely exceeds 10 %; higher differences (up to 100 %) occur very locally at some of the last grounded nodes. The results obtained with the parameter C max = 0.6 (not shown) are similar, except that the nodes where the interpolation of C S is required are fewer and, therefore, the regions within which the relative difference between τ bS06 andτ b is not null are less extended. Note that in places where the relative difference is not null, the Schoof law systematically produces the weakest basal shear stress.
Although they are moderate and very localised, the differences between the basal shear stress fields computed with The Cryosphere, 13, 177-195, 2019 www.the-cryosphere.net/13/177/2019/ mitted to the front of ice shelves, unless it can be compensated by an increase in the buttressing effect through a contact on a pinning point or an increase in lateral stresses at shear margins in the case of confined ice shelves. This latter mechanism is well illustrated in Fig. 5d-f, which shows an important increase in velocities at the shear margins of PIG when the linear Weertman law is replaced by the Schoof law. Furthermore, it turns out that the relative differences in the velocity field are more pronounced when more weight is put on basal shear stress during the initialisations. Indeed, when the inversion focuses on the basal shear stress field only (inferred state I R γ ,∞ ), the ice flow is more sensitive to any perturbation of this field. Once again, the results obtained with the parameter C max = 0.6 (not shown) are very similar with slightly lower relative differences between the velocity fields calculated with the two laws. Finally, the identification (Eq. 22) being performed at mesh nodes, the Gaussian integration over elements induces small numerical errors on the recomputed velocity field. These numerical errors are concentrated in coarse-mesh areas and correspond to the discrepancies of, at most, a few tens of percents observed further inland, in places where the differences in terms of τ b are null (Fig. 5d-f). Similar numerical errors are observed within the same areas when the velocity fields produced with the non-linear Weertman law and the two Budd laws after the identification step are compared toû b .

Prognostic simulations
The evolution of the grounded ice area during the control run for the three inferred states and the various friction laws is shown in Fig. 6. A decrease (resp. an increase) in the grounded ice area reflects a retreat (resp. an advance) of the GL. After relaxation, the grounded ice area at t = 0 a differs depending on the considered inferred state: it ranges between 433 400 km 2 for I R γ ,1 and 434 000 km 2 for I R γ ,∞ .
Irrespective of the chosen inferred state or friction law, the natural (i.e. unperturbed) evolution of the domain is toward a reduction of its grounded ice area: the Schoof law with C max = 0.4 systematically gives the more pronounced GL retreat with a reduction of grounded ice area over the 105a of the control run ranging between 5000 km 2 for I R γ ,1 and 7700 km 2 for I R γ ,∞ . In contrast, the reduction of the grounded ice area obtained with the non-linear Weertman law ranges between 2400 km 2 for I R γ ,100 and 4300 km 2 for I R γ ,∞ . Note that the evolutions of the grounded ice area obtained with the two Weertman laws are very close to one another over the entire 105 a of the control run, while the two Schoof laws produce noticeably different evolutions from The Cryosphere, 13,[177][178][179][180][181][182][183][184][185][186][187][188][189][190][191][192][193][194][195]2019 www.the-cryosphere.net/13/177/2019/ the beginning of the run, with the case C max = 0.6 showing less GL retreat than the case C max = 0.4. This could be partly due to the differences betweenû b and the velocity field recomputed with the Schoof law after initialisation, which are lower for C max = 0.6 than for C max = 0.4 as reported in Sect. 3.2. For the inferred state I R γ ,100 , the Budd laws produce intermediate results with a reduction of the grounded ice area of 3000 km 2 for the non-linear Budd law and of 4200 km 2 for the linear Budd law. The results of the perturbation experiments EXP_ABMB are shown relative to the control simulation EXP_CONTROL (Fig. 7). We consider both the relative reductions of grounded ice area as well as the relative losses of volume above flotation (VAF), given in millimetres of sea-level equivalent (SLE). For the inferred state I R γ ,∞ , the grounded ice areas produced by the two Weertman laws at the end of EXP_ABMB are about 2000 km 2 larger than the ones produced by the two Schoof laws. The difference between the linear and non-linear Weertman laws or between the two Schoof laws (C max = 0.4 and C max = 0.6) is about 10 times smaller (∼ 200 km 2 ). However, the differences in terms of VAF change are more significant with 6mm SLE of difference between the two Weertman laws and 3mm SLE of difference between the two Schoof laws. For comparison, the linear Weertman law and the Schoof law with C max = 0.4 show 24 mm SLE of difference in terms of VAF loss at the end of EXP_ABMB. When the viscosity is adjusted during inversion (cases I R γ ,100 and I R γ ,1 ), the GL shows less retreat and the relative VAF losses are less important (Fig. 7b-c and e-f). In addition, the relative evolutions of VAF obtained with the Weertman laws are closer to the ones obtained with the Schoof laws, except for the experiment associated with the inferred state I R γ ,1 and the linear Weertman law (Fig. 7e-f). Note that the Schoof laws produce more GL retreat and more SLR contribution than the Weertman laws for the inferred states I R γ ,∞ and I R γ ,100 . On the contrary, for the inferred state I R γ ,1 , the grounded ice areas seem to stabilise around t = 60 a with the Schoof laws, while they keep decreasing for a few more years with the Weertman laws, so that the grounded ice areas obtained with the two latter at the end of EXP_ABMB are slightly smaller than the ones obtained with the two former. Despite a very similar decrease in terms of grounded ice areas to the ones obtained with the two Schoof laws, the linear Weertman law predicts 8 mm SLE less VAF loss (Fig. 7f). For the inferred state I R γ ,100 , the two Budd friction laws show much stronger responses to the basal melting rate perturbation than the four other laws: at the end of EXP_ABMB, it results in a grounded ice area reduction comprised between 13 000 km 2 for the non-linear Budd law and 16 600 km 2 for the linear Budd law relative to EXP_CONTROL, while the grounded ice area reduction ranges between 6300 and 8900 km 2 for the four other experiments. Likewise, the relative VAF loss obtained at the end of EXP_ABMB is 25 mm SLE for the linear Budd law and 33 mm SLE for the non-linear Budd law, whereas it ranges between 6 and 15 mm SLE for the other experiments. This dramatic response of the GL dynamics to the perturbation when a Budd law is used is in line with the results reported in Brondex et al. (2017).
The GL positions obtained with the various initial states at the beginning and at the end of experiment EXP_ABMB are represented in Fig. 8. The GL produced with the two Schoof laws are almost systematically more retreated than the ones produced with the two Weertman laws. On the other hand, the GL final positions obtained with the Schoof law associated with C max = 0.6 (brown solid lines) are generally very close to the ones obtained with the Schoof law associated with C max = 0.4 (green solid lines). Likewise, the GL www.the-cryosphere.net/13/177/2019/ The Cryosphere, 13, 177-195, 2019 Figure 7. Evolution of grounded ice area in ×10 3 km 2 (first row) and of volume above flotation in millimetre sea-level equivalent (second row) as a function of time during EXP_ABMB (coloured solid lines) relative to EXP_CONTROL for the inferred states I R γ ,∞ -panels (a) and (d), I R γ ,100 -panels (b) and (e), and I R γ ,1 -panels (c) and (f) -with the friction laws, linear Weertman (cyan), non-linear Weertman (magenta), Schoof with C max = 0.4 (green), Schoof with C max = 0.6 (brown), linear Budd (orange, for the inferred state I R γ ,100 only) and non-linear Budd (blue, for the inferred state I R γ ,100 only). The results obtained with the various friction laws for the inferred state I R γ ,∞ are reported on, respectively, (b-c) and (e-f) for ease of comparison (black dotted lines). The vertical black dotted lines mark the introduction of the perturbation at t = 5 a. , except for the inferred state I R γ ,∞ , where the former is sometimes more advanced and sometimes more retreated than the latter. Note that this does not contradict the previously mentioned result regarding the similar evolution of grounded ice areas obtained with the two Weertman laws for I R γ ,∞ , as evolutions represented in Fig. 7 are relative to EXP_CONTROL, while the GL positions reported in Fig. 8 are absolute. The spatial distribution of the GL retreat is primarily controlled by the bed topography, as already reported by Seroussi et al. (2017). As a consequence, the Schoof and Weertman laws all give the same retreat patterns: the most pronounced retreats are observed in regions of gentle slope (particularly visible upstream of the Thwaites Ice Shelf), while the GL tends to wrap around prominences. However, for I R γ ,100 , the patterns of GL retreat obtained with the Budd laws are significantly different from the ones obtained with the four other laws. Indeed, the final GL position obtained with the non-linear Budd law is the most advanced in the Thwaites region, where the other laws (especially the Schoof laws and the linear Budd law) induce important retreats. On the contrary, both the linear and non-linear Budd laws induce retreats of several tens of kilometres in the regions of Cosgrove and Dotson ice shelves, where the four other laws give very similar final GL positions with limited retreat (especially in the Cosgrove Ice Shelf region). The initial and final ice-sheet profiles obtained with the Schoof and non-linear Budd laws are represented in Fig. 9, along four selected flow lines (reported in white in Fig. 8b). By definition, wherever ice is grounded, the ice thickness comprised between z s and the flotation altitude z f , given by z f = (1 − ρ w /ρ i )b, constitutes the thickness above flotation. Because the GL position is deduced from Eq. (11), it is located at the exact vertical of the point where z s = z f . Looking at the initial profiles, it turns out that the initial ice thickness above flotation increases rapidly upstream of the GL for the PIG and Thwaites flow lines. On the contrary, it increases very progressively as we go further upstream within the grounded part for the Dotson flow line and even more so for the Cosgrove flow line. For all the considered flow lines, the bed profiles have sections of retrograde slope susceptible to giving rise to a MISI, unless this can be prevented by sufficient lateral buttressing (Gudmundsson et al., 2012). A MISI seems to occur for the Cosgrove and Dotson flow lines when the non-linear Budd law is used and for the PIG and Thwaites flow lines when the Schoof law is used. Note that the ice shelves of Cosgrove and Dotson have completely disappeared at the end of EXP_ABMB with both the Schoof and Budd laws. On the contrary, whatever the chosen law, PIG and Thwaites ice shelves are conserved until the end of the experiment, although they become thinner. A similar figure comparing the profiles obtained with the linear and nonlinear Budd laws is available in the Supplement (Fig. S2).

Discussion
Among the three inferred states, I R γ ,100 appears to be the most physically acceptable one, despite an rms error that is slightly higher than the one obtained for I R γ ,1 . Indeed, for I R γ ,100 , the inferred η 0 has been reduced, at most, by 84 % (at the shear margins) and increased, at most, by 144 % (very locally in the higher part of the Kohler Glacier, which feeds the Dotson Ice Shelf) compared to the reference field η 0,ref .
In contrast, for I R γ ,1 , η 0 is at least 3 times higher than η 0,ref on large parts of the domain and more than 15 times higher in some very localised areas close to the Kohler Glacier (regions in red in Fig. 4i). Without considering any of the other parameters affecting the ice viscosity, such an increase could only be explained by a reduction of the ice temperature of several tens of degrees Celsius (up to 50 • C in some regions) compared to the reference temperature field. On the other hand, even if ice was assumed to be at the melting point, the viscosity reduction observed in some other parts of the domain (regions in blue in Fig. 4i) is non-physical. Note that for I R γ ,100 , an increase in the ice temperature of a few degrees only relative to the reference temperature field is sufficient to explain the local reduction of ice viscosity observed in some regions (much less extended than for I R γ ,1 ) after the inversion, except at the shear margins where the low viscosity bands can be attributed to the presence of damaged ice and/or to strain heating and/or to anisotropy. Finally, the difference between modelled and observed velocities, especially on the ice shelves, are too high for the inferred state I R γ ,∞ .
The basal melting rate increase prescribed for the perturbation experiments EXP_ABMB induces a thinning of the ice shelves, which reduces their buttressing effect and causes a retreat of the GL. The amplitude of this retreat increases as more weight is put on basal shear stress to the detriment of viscosity during initialisation. Indeed, adjusting viscosity during initialisation leads to low viscosity bands at shear margins, which hamper the transfer of lateral stresses toward the interior of ice shelves and, from there, toward the ice streams which feed them. It follows that the model states for which viscosity is inferred are less sensitive to a reduction of the buttressing effect as the contribution of this effect to the initial global stress balance is less important than for I R γ ,∞ .
The different friction laws show different sensitivities to velocity and effective pressure -from no explicit dependence on N for the Weertman laws to an explicit dependence over the whole domain for the Budd laws -which leads to different evolutions of basal shear stress as the geometry of the domain evolves in response to the perturbation. Therefore, it is not surprising that, overall, grounded ice area evolutions and VAF losses show more sensitivity to the chosen friction law when more weight is put on basal shear stress to the detriment of viscosity during initialisation, i.e. for I R γ ,∞ than for I R γ ,100 or I R γ ,1 . The choice of the friction law affects not only the grounded ice area evolution, but also the evolution of the ice thickness. Thus, the grounded ice area shrinkage is www.the-cryosphere.net/13/177/2019/ The Cryosphere, 13, 177-195, 2019 Figure 9. Ice-sheet profiles obtained for the inferred state I R γ ,100 at t = 0 a (black solid line), t = 55 a (coloured dotted line) and t = 105 a (coloured solid line) of EXP_ABMB with the Schoof law associated with C max = 0.4 (green) and the non-linear Budd law (blue), along the flow lines reported in Fig. 8. The solid light-brown line is the bed elevation. The red solid line is the flotation altitude z f , given by more pronounced when using the linear Budd law rather than the non-linear Budd law, whereas the latter induces more VAF loss than the former. Similarly, although the linear and non-linear Weertman laws lead to similar evolutions of the grounded ice area for I R γ ,∞ , they produce significantly different VAF losses. The Weertman laws systematically predict the lowest VAF losses. In addition, as for the Budd laws, the linear Weertman law always predicts less VAF loss than the non-linear Weertman law, which is in line with the results of Ritz et al. (2015), showing a higher contribution of the Antarctic ice sheet to SLR as the Weertman law tends toward a perfectly plastic law, i.e. m → 0 in Eq. (5). The two Schoof laws lead to significantly higher VAF losses than the Weertman laws, with one exception: the non-linear Weertman law for the inferred state I R γ ,1 , which predicts a VAF loss very close to the ones obtained with the two Schoof laws. This is likely due to the fact that, for I R γ ,1 , there are several regions of unexpectedly high viscosity located a few tens of The Cryosphere, 13, 177-195, 2019 www.the-cryosphere.net/13/177/2019/ kilometres upstream of the initial GL (Fig. 4i). Once the GL has reached these regions, ice is stiffer and hampers further retreat, which explains the stabilisation of the grounded ice area seen in Fig. 7c. As a consequence, the rate of VAF losses obtained with the two Schoof laws decreases, which also occurs with the non-linear Weertman law but later on, so that, at the end of EXP_ABMB, the relative VAF losses predicted by the three laws are almost identical. On the other hand, for the Schoof law, the lower the value of C max , the higher the predicted VAF losses. This is not surprising as the regions over which friction is governed by a Coulomb regime when using a Schoof law widen when the value of C max decreases (Brondex et al., 2017). For I R γ ,100 , the VAF losses obtained with the two Budd laws are dramatically higher than the one obtained with the four other laws, which is in line with the results obtained in Brondex et al. (2017). The sensitivity of ice thickness rates of change to the chosen friction law explains why the Budd laws produce a much larger GL retreat than the four other laws in the regions of Cosgrove and Dotson and, in the case of the non-linear Budd law, less retreat in the Thwaites region. As shown in Brondex et al. (2017), because the Schoof law induces a Coulomb friction regime (i.e. |τ b | ∼ C max N ) within a narrow region of a few kilometres directly upstream of the GL, the buttressing loss due to basal melting rate increase cannot be compensated within this region, despite an increase in ice velocities. Therefore, the perturbation is transmitted upstream, in areas governed by a Weertman friction regime (i.e. |τ b | ∼ C S |u b | m ), where a velocity increase enables a rapid compensation of the buttressing loss. As a consequence, although the Schoof law induces an important peak of thinning, it stays very localised in the immediate vicinity of the GL. In contrast, for the Budd laws, the basal shear stress depends on both u b and N over the whole domain. In addition, N can be low on large distances upstream of the GL, especially where the bed profile is retrograde and induces a rapid increase in the water column height as we go further upstream, while the ice thickness does not increase as fast. As a consequence, the peak of thinning obtained at the GL with a Budd law is less pronounced than the one obtained with the Schoof law, but, depending on the bed profile, it can propagate much further upstream. This is particularly well illustrated in Fig. 3 of Brondex et al. (2017). In addition, as for the Weertman law, perturbations are transmitted farther and faster inland when using the non-linear Budd law rather than the linear Budd law. As a consequence, the peak of thinning obtained with the linear Budd law in the immediate vicinity of the GL is slightly more pronounced than the one obtained with the non-linear Budd law, but it propagates over slightly shorter distances (Fig. S2). In the regions of Cosgrove and Dotson, the initial surface profiles are close to the flotation altitudes over large distances upstream of the GL (Fig. 9). In these two regions, both the linear and non-linear Budd laws induce a rapid purge of the VAF over these large distances, leading to an important retreat of the GL. This retreat might be enhanced by a potential MISI as the GL ends up on retrograde bed slopes. Note, however, that the Cosgrove and Dotson ice shelves being both well confined, knowledge of the bed slopes alone is not sufficient to predict the occurrence of a MISI (Gudmundsson et al., 2012;Haseloff and Sergienko, 2018).
In contrast, Thwaites ice shelf is mostly unconfined and, at the beginning of the experiment, the GL lies at the downstream bound of a section of retrograde bed slope (solid black line in the bottom-left panel of Fig. 9). The initial thinning produced by the Schoof law in the very close vicinity of the GL is sufficient to induce a retreat within the reverse slope area, where a MISI likely initiates. The same goes for the linear Budd law (Fig. S2). On the contrary, the thinning peak produced by the non-linear Budd law at the GL is not sufficient to induce a retreat of the latter over the 105 years of the experiment. Finally, the GL retreats obtained with the Budd and Schoof laws for PIG are more similar. This could be a consequence of the fact that the GL already lies on a retrograde slope at t = 0 a.
Although the differences in terms of GL dynamics and VAF losses obtained with the different friction laws tend to decrease as more weight is put on viscosity during inversion, they remain significant for the most physically acceptable inferred state, I R γ ,100 . In particular, the commonly used linear Weertman law predicts about 9 mm SLE less VAF losses at the end of the 105 years than the Schoof law.
In light of these results, using observations to constrain the form of the friction law which is best suited for a given application would constitute a huge leap towards the production of reliable projections of future SLR. Studies reporting the presence of soft sediments beneath PIG (e.g. Smith et al., 2013;Brisbourne et al., 2017) tend to support the use of a Schoof law in this region, as such a law induces a Coulomb friction regime in the vicinity of the GL. However, this does not constitute a validation of this law and neither does it provide any information regarding the spatial distributions of C S and C max . Because the basal stress must satisfy the global stress balance, constraining the form of the friction law would require observations at different times with significant differences in basal velocities, basal stresses and water pressure at the ice-bed interface. Unfortunately, these multiple sets of observations are not available or incomplete. In particular, the water pressure at the ice-bed interface is largely unknown and the assumption of perfect hydrological connectivity to the ocean is too gross for the purpose of constraining the form of the friction law.
More generally, the lack of a physically based subglacial hydrological model used to compute effective pressure is the main limitation of the work presented here. Indeed, the water pressure calculated based on the assumption of perfect hydrological connectivity to the ocean might well be underestimated in some places as geothermal heat flux and frictional heating are known to cause basal melting, which is susceptible to build up water pressure at the ice-bed interface, eswww.the-cryosphere.net/13/177/2019/ The Cryosphere, 13, 177-195, 2019 pecially in regions where the subglacial drainage system is inefficient (Colleoni et al., 2018). In addition, beside basal shear stress, ice viscosity is not the only poorly constrained field. In particular, there are important uncertainties regarding the bed elevation. Recently, Nias et al. (2018) derived new bedrock elevation and ice thickness maps of the PIG using a method based on the principle of mass conservation and compared these maps to the Bedmap2 data: they found substantial differences between the two, in particular directly upstream of the 1996 GL, where a topographic rise has been removed with their new method. In addition, they ran the two obtained geometries forward in time 50 years using four different friction laws, i.e. three Weertman laws with various values for the exponent m and a Tsai law with a Coulomb friction parameter set to 0.5 and showed that SLR projections are at least as sensitive to the accuracy of the bedrock topography and initial ice thickness as to the choice of the friction law.

Conclusions
The present study constitutes an extension of the sensibility analysis of GL dynamics and VAF loss predictions to the choice of a friction law presented in Brondex et al. (2017). Whereas the latter was based on a synthetic 2-D flow line case, here we consider a real-world application, the Amundsen Sea Embayment, and therefore, have to address specific problems. First of all, the basal shear stress is not known and has to be deduced from observations of ice flow surface velocities through inverse methods. In addition, other model parameters are uncertain and can require adjustment based on observations. Thus, when ice viscosity is not inferred but deduced from available ice temperature fields, discrepancies between modelled and observed velocities are high, in particular within ice shelves. On the contrary, putting too much weight on viscosity to the detriment of basal shear stress during inversion leads to a non-physical viscosity field. The inferred state I R γ ,100 turns out to be the most physically acceptable of the constructed model states as it allows a good fit between observed and modelled velocities while showing a viscosity field closer to our expectations compared with I R γ ,1 . Once the basal shear stress has been inferred, the friction coefficient distributions of the Weertman and Budd laws producing the same basal shear stress field with these laws can easily be identified. In contrast, by construction the Schoof law induces an upper bound of the computed basal shear stress equal to the value of the parameter C max . This latter parameter accounts for till deformation and the values it can take are partly constrained by laboratory experiments. As a result, the value of C S cannot be identified directly at nodes too close to flotation. This difficulty is overcome by interpolating (or extrapolating) the value of C S at these nodes from the values at nodes where it can be directly identified. This procedure induces significant but very localised discrepan-cies between the recomputed velocity field and the reference velocity field used for the identification, in particular within ice shelves.
The perturbation experiments carried out following the identification step demonstrate a significant influence of the chosen friction law on GL dynamics and mass loss projections on a century timescale. In line with results obtained in our previous study, the VAF losses obtained with the commonly used Weertman law are systematically lower than the ones obtained with the Schoof law, which themselves are surpassed in the simulations carried out with the Budd laws (for I R γ ,100 ). In addition, the Budd laws being dependent on effective pressure over the whole grounded domain, they induce GL retreat in places where the other tested laws do not. As for the Weertman laws, the non-linear Budd law induces more VAF loss than the linear Budd law. Although the differences between the results produced with the various law tend to decrease as more weight is put on viscosity during inversion, they are still significant for the most physically acceptable model state constructed, i.e. I R γ ,100 .
In light of these results, we conclude that no reliable projection of future sea-level rise can be obtained without the use of a physically based friction law. Therefore, significant efforts still need to be made for a better understanding of the physical processes at play at the ice-bed interface in order to constrain the form of the friction law that needs to be used in models. In particular, the recognised importance of water pressure for basal sliding must be explicitly accounted for through the use of an effective pressure-dependent law. Because it implies a Coulomb friction regime at low effective pressure, which is best suited to model till deformation, we suggest using the Schoof law rather than the Budd law. However, the use of such a law adds the further difficulty of estimating the spatial distribution of water pressure at the ice-bed interface, which seems to be satisfactorily overcome only through a coupling of the ice flow model to a physically based subglacial hydrological model.
Code and data availability. Elmer/Ice code is publicly available through GitHub (https://github.com/ElmerCSC/elmerfem, last access: 17 January 2019; Gagliardini et al., 2013). All the simulations were performed with the version 8.2 (Rev: 997cb45) of Elmer/Ice. All scripts used for simulations and post-treatment as well as model output are available upon request from authors. The data used are listed in the references, except the surface mass balance, which is available upon request from the authors.