Journal cover
Journal topic
**The Cryosphere**
An interactive open-access journal of the European Geosciences Union

Journal topic

- About
- Editorial board
- Articles
- Special issues
- Highlight articles
- Manuscript tracking
- Subscribe to alerts
- Peer review
- For authors
- For reviewers
- EGU publications
- Imprint
- Data protection

- About
- Editorial board
- Articles
- Special issues
- Highlight articles
- Manuscript tracking
- Subscribe to alerts
- Peer review
- For authors
- For reviewers
- EGU publications
- Imprint
- Data protection

**Research article**
22 Jan 2019

**Research article** | 22 Jan 2019

Sensitivity of Amundsen basin centennial evolution to friction law

- Univ. Grenoble Alpes, CNRS, IRD, IGE, 38000 Grenoble, France

Abstract

Back to toptop
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.

How to cite

Back to top
top
How to cite.

Brondex, J., Gillet-Chaulet, F., and Gagliardini, O.: Sensitivity of centennial mass loss projections of the Amundsen basin to the friction law, The Cryosphere, 13, 177-195, https://doi.org/10.5194/tc-13-177-2019, 2019.

1 Introduction

Back to toptop
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 active
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 (Mouginot et al., 2014; 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
(Cornford et al., 2015; Favier et al., 2014; Joughin et al., 2014). Producing reliable
estimates of this future contribution requires accurate modelling of GL
dynamics on sub-centennial 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
(Budd et al., 1979; Schoof, 2005; Tsai et al., 2015; Weertman, 1957). 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; Blankenship 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 (Alley et al., 1986; Blankenship et al., 1986; Boulton and Jones, 1979). 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). Gradient-based 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 (Cornford et al., 2015; Gillet-Chaulet et al., 2012; Morlighem et al., 2010, 2013). 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 (Gillet-Chaulet et al., 2016; Joughin et al., 2004). 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 perturbation 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.

2 Methods

Back to toptop
The modelled domain is represented in Fig. 1. For the stress balance, we solve the two-dimensional Shelfy-Stream Approximation (SSA) equations (MacAyeal, 1989), for which the horizontal velocity field $\mathbf{u}=(u,v)$ is a solution of

$$\begin{array}{ll}{\displaystyle}\mathrm{4}{\displaystyle \frac{\partial}{\partial x}}\left(H\stackrel{\mathrm{\u203e}}{\mathit{\eta}}{\displaystyle \frac{\partial u}{\partial x}}\right)+\mathrm{2}{\displaystyle \frac{\partial}{\partial x}}\left(H\stackrel{\mathrm{\u203e}}{\mathit{\eta}}{\displaystyle \frac{\partial v}{\partial y}}\right)& {\displaystyle}+{\displaystyle \frac{\partial}{\partial y}}\left(H\stackrel{\mathrm{\u203e}}{\mathit{\eta}}\left({\displaystyle \frac{\partial u}{\partial y}}+{\displaystyle \frac{\partial v}{\partial x}}\right)\right)\\ {\displaystyle}& {\displaystyle}+{\mathit{\tau}}_{\text{b},x}={\mathit{\rho}}_{\text{i}}gH{\displaystyle \frac{\partial {z}_{\text{s}}}{\partial x}}\end{array}$$

$$\begin{array}{ll}{\displaystyle}\mathrm{4}{\displaystyle \frac{\partial}{\partial y}}\left(H\stackrel{\mathrm{\u203e}}{\mathit{\eta}}{\displaystyle \frac{\partial v}{\partial y}}\right)+\mathrm{2}{\displaystyle \frac{\partial}{\partial y}}\left(H\stackrel{\mathrm{\u203e}}{\mathit{\eta}}{\displaystyle \frac{\partial u}{\partial x}}\right)& {\displaystyle}+{\displaystyle \frac{\partial}{\partial x}}\left(H\stackrel{\mathrm{\u203e}}{\mathit{\eta}}\left({\displaystyle \frac{\partial u}{\partial y}}+{\displaystyle \frac{\partial v}{\partial x}}\right)\right)\\ \text{(1)}& {\displaystyle}& {\displaystyle}+{\mathit{\tau}}_{\text{b},y}={\mathit{\rho}}_{\text{i}}gH{\displaystyle \frac{\partial {z}_{\text{s}}}{\partial y}},\end{array}$$

where *ρ*_{i} is the ice density, *g* is the gravity norm and $H={z}_{\text{s}}-{z}_{\text{b}}$
is the ice thickness, with *z*_{s} and *z*_{b} the top and bottom surface
elevations. The vertically averaged effective viscosity
$\stackrel{\mathrm{\u203e}}{\mathit{\eta}}$ reads as follows:

$$\begin{array}{}\text{(2)}& \stackrel{\mathrm{\u203e}}{\mathit{\eta}}={\stackrel{\mathrm{\u203e}}{\mathit{\eta}}}_{\mathrm{0}}{D}_{e}^{(\mathrm{1}-n)/n},\end{array}$$

where *D*_{e} is the second invariant of the strain-rate tensor, *n* is the
Glen's law exponent and ${\stackrel{\mathrm{\u203e}}{\mathit{\eta}}}_{\mathrm{0}}$ is given by

$$\begin{array}{}\text{(3)}& {\stackrel{\mathrm{\u203e}}{\mathit{\eta}}}_{\mathrm{0}}={\displaystyle \frac{\mathrm{1}}{H}}\underset{{z}_{\text{b}}}{\overset{{z}_{\text{s}}}{\int}}{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{A}^{-\mathrm{1}/n}\text{d}z.\end{array}$$

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:

$$\begin{array}{}\text{(4)}& A={A}_{\mathrm{0}}{e}^{\left(-Q/\left[R\right(\mathrm{273.15}+{T}^{\prime}]\right)},\end{array}$$

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:

$$\begin{array}{}\text{(5)}& {\displaystyle}& {\displaystyle}{\mathit{\tau}}_{\text{b}}+{C}_{\text{W}}|{\mathbf{u}}_{\text{b}}{|}^{m-\mathrm{1}}{\mathbf{u}}_{\text{b}}=\mathrm{0},\text{(6)}& {\displaystyle}& {\displaystyle}{\mathit{\tau}}_{\text{b}}+{C}_{\text{B}}N|{\mathbf{u}}_{\text{b}}{|}^{m-\mathrm{1}}{\mathbf{u}}_{\text{b}}=\mathrm{0},\text{(7)}& {\displaystyle}& {\displaystyle}{\mathit{\tau}}_{\text{b}}+{\displaystyle \frac{{C}_{\text{S}}|{\mathbf{u}}_{\text{b}}{|}^{m-\mathrm{1}}{\mathbf{u}}_{\text{b}}}{{\left(\mathrm{1}+{\left(\frac{{C}_{\text{S}}}{{C}_{\text{max}}N}\right)}^{\mathrm{1}/m}\left|{\mathbf{u}}_{\text{b}}\right|\right)}^{m}}}=\mathrm{0},\end{array}$$

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=\mathrm{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 $\left|{\mathit{\tau}}_{\text{b}}\right|/N$ known as Iken's bound that is equal to the value
of the parameter *C*_{max} (Gagliardini et al., 2007; Iken, 1981; Schoof, 2005).
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
(Fricker et al., 2007; Gray et al., 2005; Le Brocq et al., 2013). However, unlike the Schoof law,
for which the dependence of $\left|{\mathit{\tau}}_{\text{b}}\right|$ on *N* is limited to the
close vicinity of the GLs, the Budd law induces a dependence of
$\left|{\mathit{\tau}}_{\text{b}}\right|$ 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 (Hewitt et al., 2012; Schoof, 2010; Werder et al., 2013; de Fleurian et al., 2014)
and the attempts to couple these models to ice flow models are relatively
recent (Bougamont et al., 2014; Bueler and van Pelt, 2015; Gagliardini and Werder, 2018; Hewitt, 2013).
Therefore, in the present study, *N* is calculated assuming a perfect
hydrological connection between the subglacial drainage system and the ocean,
so that

$$\begin{array}{}\text{(8)}& N=\left\{\begin{array}{ll}{\mathit{\rho}}_{\text{i}}gH+{\mathit{\rho}}_{\text{w}}g{z}_{\text{b}}& \text{if}\phantom{\rule{0.25em}{0ex}}{z}_{\text{b}}<\mathrm{0},\\ {\mathit{\rho}}_{\text{i}}gH& \text{if}\phantom{\rule{0.25em}{0ex}}{z}_{\text{b}}\ge \mathrm{0},\end{array}\right.\end{array}$$

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:

$$\begin{array}{}\text{(9)}& {\displaystyle \frac{\partial H}{\partial t}}+{\displaystyle \frac{\partial \left(uH\right)}{\partial x}}+{\displaystyle \frac{\partial \left(vH\right)}{\partial y}}={a}_{\text{s}}-{a}_{\text{b}},\end{array}$$

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

$$\begin{array}{}\text{(10)}& \left\{\begin{array}{ll}{z}_{\text{b}}(x,y,t)=b(x,y)& \text{for grounded ice,}\\ {z}_{\text{b}}(x,y,t)=-H{\displaystyle \frac{{\mathit{\rho}}_{\text{i}}}{{\mathit{\rho}}_{\text{w}}}}>b(x,y)& \text{for floating ice}.\end{array}\right.\end{array}$$

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:

$$\begin{array}{}\text{(11)}& H({x}_{\text{G}},{y}_{\text{G}})+b({x}_{\text{G}},{y}_{\text{G}}){\displaystyle \frac{{\mathit{\rho}}_{\text{w}}}{{\mathit{\rho}}_{\text{i}}}}=\mathrm{0}.\end{array}$$

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, $\mathbf{n}=({n}_{x},{n}_{y})$ is the normal unit vector to the considered boundary. At the ice divides, there is no ice flux entering the domain, so the following Dirichlet condition applies:

$$\begin{array}{}\text{(12)}& (\mathbf{u}.\mathbf{n}){|}_{\text{id}}=\mathrm{0},\end{array}$$

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:

$$\begin{array}{ll}{\displaystyle}\mathrm{4}H\stackrel{\mathrm{\u203e}}{\mathit{\eta}}{\left.{\displaystyle \frac{\partial u}{\partial x}}\right|}_{\text{f}}{n}_{x}& {\displaystyle}+\mathrm{2}H\stackrel{\mathrm{\u203e}}{\mathit{\eta}}{\left.{\displaystyle \frac{\partial v}{\partial y}}\right|}_{\text{f}}{n}_{x}+H\stackrel{\mathrm{\u203e}}{\mathit{\eta}}{\left.\left({\displaystyle \frac{\partial u}{\partial y}}+{\displaystyle \frac{\partial v}{\partial x}}\right)\right|}_{\text{f}}{n}_{y}\\ {\displaystyle}& {\displaystyle}={\displaystyle \frac{g}{\mathrm{2}}}({\mathit{\rho}}_{\text{i}}H{|}_{\text{f}}^{\mathrm{2}}-{\mathit{\rho}}_{\text{w}}{H}_{\text{sub}}{|}_{\text{f}}^{\mathrm{2}}){n}_{x}\\ {\displaystyle}\mathrm{2}H\stackrel{\mathrm{\u203e}}{\mathit{\eta}}{\left.{\displaystyle \frac{\partial u}{\partial x}}\right|}_{\text{f}}{n}_{y}& {\displaystyle}+\mathrm{4}H\stackrel{\mathrm{\u203e}}{\mathit{\eta}}{\left.{\displaystyle \frac{\partial v}{\partial y}}\right|}_{\text{f}}{n}_{y}+H\stackrel{\mathrm{\u203e}}{\mathit{\eta}}{\left.\left({\displaystyle \frac{\partial u}{\partial y}}+{\displaystyle \frac{\partial v}{\partial x}}\right)\right|}_{\text{f}}{n}_{x}\\ \text{(13)}& {\displaystyle}& {\displaystyle}={\displaystyle \frac{g}{\mathrm{2}}}({\mathit{\rho}}_{\text{i}}H{|}_{\text{f}}^{\mathrm{2}}-{\mathit{\rho}}_{\text{w}}{H}_{\text{sub}}{|}_{\text{f}}^{\mathrm{2}}){n}_{y},\end{array}$$

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:

$$\begin{array}{}\text{(14)}& {a}_{\text{b}}={\displaystyle \frac{\mathrm{8}{K}_{\text{T}}{\mathit{\rho}}_{\text{w}}{C}_{\text{W}}}{{\mathit{\rho}}_{\text{i}}{L}_{\text{f}}}}|{T}_{\mathrm{0}}-{T}_{\text{f}}|({T}_{\mathrm{0}}-{T}_{\text{f}}),\end{array}$$

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}_{\mathrm{0}}-{T}_{\text{f}}=\mathrm{0.5}$ ^{∘}C
for $\mathrm{0}>z>-\mathrm{170}$ m, 3.5 ^{∘}C for $z<-\mathrm{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 mesh-adaptation 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 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.

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=\mathrm{1}/\mathrm{3}$,
(3) a non-linear Schoof law given by Eq. (7) with $m=\mathrm{1}/\mathrm{3}$
and *C*_{max}=0.4, (4) a non-linear Schoof law given by Eq. (7)
with $m=\mathrm{1}/\mathrm{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=\mathrm{1}/\mathrm{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
$\left|{\mathit{\tau}}_{\text{b}}\right|$ 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 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.

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 ${\stackrel{\mathrm{\u203e}}{\mathit{\eta}}}_{\mathrm{0}}$, denoted by
${\stackrel{\mathrm{\u203e}}{\mathit{\eta}}}_{\mathrm{0},\text{ref}}$, is computed using a three-dimensional ice temperature
field of Antarctica (Brice Van Liefferinge, personal communication, 2016) obtained using
a thermo-mechanical model forced with a mean geothermal heat flux
(Van Liefferinge and Pattyn, 2013). In addition, we use the observed surface
velocities of Rignot et al. (2011) 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
${\widehat{C}}_{\text{W}}$. Although ${\stackrel{\mathrm{\u203e}}{\mathit{\eta}}}_{\mathrm{0},\text{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 ${\stackrel{\mathrm{\u203e}}{\mathit{\eta}}}_{\mathrm{0},\text{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}_{\mathit{\gamma},\mathrm{\infty}}}$,
${I}_{{R}_{\mathit{\gamma},\mathrm{100}}}$ and ${I}_{{R}_{\mathit{\gamma},\mathrm{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:

$$\begin{array}{}\text{(15)}& {J}_{\text{tot}}={J}_{v}+{\mathit{\lambda}}_{\text{div}}{J}_{\text{div}}+{\mathit{\lambda}}_{\text{reg},\mathit{\alpha}}{J}_{\text{reg},\mathit{\alpha}}+{\mathit{\lambda}}_{\text{reg},\mathit{\gamma}}{J}_{\text{reg},\mathit{\gamma}}.\end{array}$$

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:

$$\begin{array}{}\text{(16)}& {J}_{\text{v}}={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\sum _{i=\mathrm{1}}^{{N}_{\text{obs}}}\left(\right|{\mathbf{u}}_{i,\text{obs}}-{\mathbf{u}}_{i}|{)}^{\mathrm{2}},\end{array}$$

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:

$$\begin{array}{}\text{(17)}& {J}_{\text{div}}={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\underset{\mathrm{\Gamma}}{\int}{\left({\displaystyle \frac{\partial \left(uH\right)}{\partial x}}+{\displaystyle \frac{\partial \left(vH\right)}{\partial y}}-({a}_{\text{s}}-{a}_{\text{b}})\right)}^{\mathrm{2}}\mathrm{d}\mathrm{\Gamma},\end{array}$$

where Γ is the model domain. During the inversion, we actually
optimise the variables *α* and *γ*, which are related to the linear
Weertman law coefficient and ${\stackrel{\mathrm{\u203e}}{\mathit{\eta}}}_{\mathrm{0}}$, respectively, as follows:

$$\begin{array}{}\text{(18)}& {\widehat{C}}_{\text{W}}={\mathrm{10}}^{\mathit{\alpha}},\end{array}$$

and

$$\begin{array}{}\text{(19)}& {\stackrel{\mathrm{\u203e}}{\mathit{\eta}}}_{\mathrm{0}}={\displaystyle \frac{{\mathit{\gamma}}^{\mathrm{2}}}{{R}_{\mathit{\gamma}}}}{\stackrel{\mathrm{\u203e}}{\mathit{\eta}}}_{\mathrm{0},\text{ref}},\end{array}$$

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 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}_{\mathit{\gamma},\mathrm{100}}}$ and ${I}_{{R}_{\mathit{\gamma},\mathrm{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}_{\mathit{\gamma},\mathrm{100}}}$ than for
${I}_{{R}_{\mathit{\gamma},\mathrm{1}}}$. In contrast, the inferred state ${I}_{{R}_{\mathit{\gamma},\mathrm{\infty}}}$
is obtained by optimising *α* only, while
${\stackrel{\mathrm{\u203e}}{\mathit{\eta}}}_{\mathrm{0}}={\stackrel{\mathrm{\u203e}}{\mathit{\eta}}}_{\mathrm{0},\text{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 ${\stackrel{\mathrm{\u203e}}{\mathit{\eta}}}_{\mathrm{0}}$, it is simply set to ${\stackrel{\mathrm{\u203e}}{\mathit{\eta}}}_{\mathrm{0},\text{ref}}$
for both ${I}_{{R}_{\mathit{\gamma},\mathrm{100}}}$ and ${I}_{{R}_{\mathit{\gamma},\mathrm{1}}}$.

For each of the three inferred states, we follow a L-curve approach as
described in Fürst et al. (2015) to retrieve the optimal values of *λ*_{div}, *λ*_{reg,α} and *λ*_{reg,γ} (only for
${I}_{{R}_{\mathit{\gamma},\mathrm{100}}}$ and ${I}_{{R}_{\mathit{\gamma},\mathrm{1}}}$), leading to inferred states
showing good compromise between smooth distributions of *α* and
*γ*, low differences between ice fluxes and mass balance and a good fit
between observed and modelled velocities. We find $(\mathrm{5.1}\times {\mathrm{10}}^{-\mathrm{5}}$, 7.1×10^{5}) for ${I}_{{R}_{\mathit{\gamma},\mathrm{\infty}}}$, $(\mathrm{5.9}\times {\mathrm{10}}^{-\mathrm{5}}$,
2.5×10^{6}, 1.4×10^{6)} for ${I}_{{R}_{\mathit{\gamma},\mathrm{100}}}$, and $(\mathrm{1.8}\times {\mathrm{10}}^{-\mathrm{5}}$, 5.4×10^{6}, 4.7×10^{5}) for
${I}_{{R}_{\mathit{\gamma},\mathrm{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).

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
${\widehat{\mathbf{u}}}_{\text{b}}$, as well as a new basal shear stress field, denoted
by ${\widehat{\mathit{\tau}}}_{\text{b}}$. These two fields relate to each other through
the linear Weertman law, i.e. Eq. (5) with the inferred
friction coefficient ${\widehat{C}}_{\text{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 ${\widehat{\mathit{\tau}}}_{\text{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 solving the following:

$$\begin{array}{}\text{(20)}& {C}_{\text{W,nl}}={\widehat{C}}_{\text{W}}|{\widehat{\mathbf{u}}}_{\text{b}}{|}^{\mathrm{1}-m},\end{array}$$

$$\begin{array}{}\text{(21)}& {C}_{\text{B}}={\widehat{C}}_{\text{W}}{\displaystyle \frac{|{\widehat{\mathbf{u}}}_{\text{b}}{|}^{\mathrm{1}-m}}{\widehat{N}}},\end{array}$$

and,

$$\begin{array}{}\text{(22)}& {C}_{\text{S}}={\displaystyle \frac{\left|{\widehat{\mathit{\tau}}}_{\text{b}}\right|}{|{\widehat{\mathbf{u}}}_{\text{b}}{|}^{m}{\left(\mathrm{1}-{\left(\frac{\left|{\widehat{\mathit{\tau}}}_{\text{b}}\right|}{{C}_{\text{max}}\widehat{N}}\right)}^{\mathrm{1}/m}\right)}^{m}}}.\end{array}$$

For the Budd and Schoof laws, the effective pressure field $\widehat{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 ${\widehat{C}}_{\text{W}}$, ${\widehat{\mathbf{u}}}_{\text{b}}$, ${\widehat{\mathit{\tau}}}_{\text{b}}$ and
$\widehat{N}$ 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}_{\mathit{\gamma},\mathrm{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}_{\text{W,nl}}={\mathrm{10}}^{-\mathrm{6}}$ MPa m${}^{-\mathrm{1}/\mathrm{3}}$ a^{1∕3} and
${C}_{\text{B}}={\mathrm{10}}^{-\mathrm{6}}$ m${}^{-\mathrm{1}/\mathrm{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 $\left|{\widehat{\mathit{\tau}}}_{\text{b}}\right|/\widehat{N}<{C}_{\text{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 $\left|{\widehat{\mathit{\tau}}}_{\text{b}}\right|/\widehat{N}\le \mathrm{0.8}{C}_{\text{max}}$. The
values of *C*_{S} at nodes where $\left|{\widehat{\mathit{\tau}}}_{\text{b}}\right|/\widehat{N}>\mathrm{0.8}{C}_{\text{max}}$ are linearly interpolated from the closest neighbouring nodes. In
addition, we arbitrarily set ${C}_{\text{S}}={\mathrm{10}}^{-\mathrm{3}}$ MPa m${}^{-\mathrm{1}/\mathrm{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}_{\mathit{\gamma},\mathrm{\infty}}}$, *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}_{\mathit{\gamma},\mathrm{\infty}}}$ 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).

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

$$\begin{array}{}\text{(23)}& \left\{\begin{array}{ll}{a}_{\text{b,ABMB}}={a}_{\text{b}}& \text{if}\phantom{\rule{0.125em}{0ex}}\phantom{\rule{0.125em}{0ex}}\mathrm{0}\le t<\mathrm{5}\phantom{\rule{0.125em}{0ex}}\text{a}\\ {a}_{\text{b,ABMB}}={a}_{\text{b}}+\text{abmb}{\displaystyle \frac{\text{floor}\left(t\right)-\mathrm{5}}{\mathrm{40}}}& \text{if}\phantom{\rule{0.125em}{0ex}}\phantom{\rule{0.125em}{0ex}}\mathrm{5}\le t<\mathrm{45}\phantom{\rule{0.125em}{0ex}}\text{a}\\ {a}_{\text{b,ABMB}}={a}_{\text{b}}+\text{abmb}& \text{if}\phantom{\rule{0.125em}{0ex}}\phantom{\rule{0.125em}{0ex}}t\ge \mathrm{45}\phantom{\rule{0.125em}{0ex}}\text{a}.\end{array}\right.\end{array}$$

The basal melting rate anomaly abmb corresponds to the one prescribed for the ASE within the InitMIP-Antarctica framework (Nowicki et al., 2016).

3 Results

Back to toptop
The three inferred states ${I}_{{R}_{\mathit{\gamma},\mathrm{\infty}}}$, ${I}_{{R}_{\mathit{\gamma},\mathrm{100}}}$ and
${I}_{{R}_{\mathit{\gamma},\mathrm{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}_{\mathit{\gamma},\mathrm{\infty}}}$ 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}_{\mathit{\gamma},\mathrm{\infty}}}$, ${I}_{{R}_{\mathit{\gamma},\mathrm{100}}}$ and
${I}_{{R}_{\mathit{\gamma},\mathrm{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}_{\mathit{\gamma},\mathrm{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 (Bondzio et al., 2017; Borstad et al., 2013; 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}_{\mathit{\gamma},\mathrm{\infty}}}$,
${I}_{{R}_{\mathit{\gamma},\mathrm{100}}}$ and ${I}_{{R}_{\mathit{\gamma},\mathrm{1}}}$, respectively.

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 ${\widehat{\mathit{\tau}}}_{\text{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, ${\widehat{\mathit{\tau}}}_{\text{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
${\widehat{\mathit{\tau}}}_{\text{b}}$ are the highest. As expected, the highest
relative differences are obtained close to the ice shelves (shown in green in
Fig. 5a–c), 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
${\widehat{\mathit{\tau}}}_{\text{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
${\widehat{\mathit{\tau}}}_{\text{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 linear Weertman law and the Schoof law induce
significant differences on the recomputed velocity fields (Fig. 5d–f).
These differences propagate from the
regions showing the highest differences in terms of *τ*_{b} to
the floating areas located directly downstream (particularly distinguishable on
Thwaites Ice Shelf). Indeed, the decrease in the total basal shear stress due
to the interpolation of *C*_{S} close to the GL needs to be compensated
elsewhere so that the global stress balance keeps being satisfied. Because
the floating areas do not support any basal shear stress, the perturbation is
transmitted 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}_{\mathit{\gamma},\mathrm{\infty}}}$), 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 ${\widehat{\mathbf{u}}}_{\text{b}}$.

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}_{\mathit{\gamma},\mathrm{1}}}$ and 434 000 km^{2} for ${I}_{{R}_{\mathit{\gamma},\mathrm{\infty}}}$.

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}_{\mathit{\gamma},\mathrm{1}}}$ and 7700 km^{2} for ${I}_{{R}_{\mathit{\gamma},\mathrm{\infty}}}$. In
contrast, the reduction of the grounded ice area obtained with the non-linear
Weertman law ranges between 2400 km^{2} for ${I}_{{R}_{\mathit{\gamma},\mathrm{100}}}$ and
4300 km^{2} for ${I}_{{R}_{\mathit{\gamma},\mathrm{\infty}}}$. 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 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
${\widehat{\mathbf{u}}}_{\text{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}_{\mathit{\gamma},\mathrm{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}_{\mathit{\gamma},\mathrm{\infty}}}$,
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}_{\mathit{\gamma},\mathrm{100}}}$ and ${I}_{{R}_{\mathit{\gamma},\mathrm{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}_{\mathit{\gamma},\mathrm{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}_{\mathit{\gamma},\mathrm{\infty}}}$ and ${I}_{{R}_{\mathit{\gamma},\mathrm{100}}}$. On the contrary, for
the inferred state ${I}_{{R}_{\mathit{\gamma},\mathrm{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}_{\mathit{\gamma},\mathrm{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 final positions obtained with the
linear Weertman law (cyan solid lines) are often very close to the ones
obtained with the non-linear Weertman law (magenta solid lines), except for
the inferred state ${I}_{{R}_{\mathit{\gamma},\mathrm{\infty}}}$, 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}_{\mathit{\gamma},\mathrm{\infty}}}$, 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}_{\mathit{\gamma},\mathrm{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}_{\text{f}}=(\mathrm{1}-{\mathit{\rho}}_{w}/{\mathit{\rho}}_{\text{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 non-linear
Budd laws is available in the Supplement (Fig. S2).

4 Discussion

Back to toptop
Among the three inferred states, ${I}_{{R}_{\mathit{\gamma},\mathrm{100}}}$ appears to be the most
physically acceptable one, despite an rms error that is slightly higher than the one
obtained for ${I}_{{R}_{\mathit{\gamma},\mathrm{1}}}$. Indeed, for ${I}_{{R}_{\mathit{\gamma},\mathrm{100}}}$, the
inferred ${\stackrel{\mathrm{\u203e}}{\mathit{\eta}}}_{\mathrm{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 ${\stackrel{\mathrm{\u203e}}{\mathit{\eta}}}_{\mathrm{0},\text{ref}}$. In contrast, for ${I}_{{R}_{\mathit{\gamma},\mathrm{1}}}$,
${\stackrel{\mathrm{\u203e}}{\mathit{\eta}}}_{\mathrm{0}}$ is at least 3 times higher than ${\stackrel{\mathrm{\u203e}}{\mathit{\eta}}}_{\mathrm{0},\text{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}_{\mathit{\gamma},\mathrm{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}_{\mathit{\gamma},\mathrm{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}_{\mathit{\gamma},\mathrm{\infty}}}$.

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}_{\mathit{\gamma},\mathrm{\infty}}}$.

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}_{\mathit{\gamma},\mathrm{\infty}}}$ than for ${I}_{{R}_{\mathit{\gamma},\mathrm{100}}}$ or ${I}_{{R}_{\mathit{\gamma},\mathrm{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 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}_{\mathit{\gamma},\mathrm{\infty}}}$, 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}_{\mathit{\gamma},\mathrm{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}_{\mathit{\gamma},\mathrm{1}}}$, there are several regions of
unexpectedly high viscosity located a few tens of 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}_{\mathit{\gamma},\mathrm{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. $\left|{\mathit{\tau}}_{\text{b}}\right|\sim {C}_{\text{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.
$\left|{\mathit{\tau}}_{\text{b}}\right|\sim {C}_{\text{S}}|{\mathbf{u}}_{\text{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}_{\mathit{\gamma},\mathrm{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
(Brisbourne et al., 2017; Smith et al., 2013) 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, especially 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.

5 Conclusions

Back to toptop
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}_{\mathit{\gamma},\mathrm{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}_{\mathit{\gamma},\mathrm{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 discrepancies 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}_{\mathit{\gamma},\mathrm{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}_{\mathit{\gamma},\mathrm{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

Back to toptop
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.

Supplement

Back to toptop
Supplement.

The supplement related to this article is available online at: https://doi.org/10.5194/tc-13-177-2019-supplement.

Author contributions

Back to toptop
Author contributions.

JB, FGC and OG designed the study. FGC collected the data and prepared the experimental setup. JB conducted the inversions, the identifications and the prognostic simulations. JB wrote the paper with contributions from all co-authors.

Competing interests

Back to toptop
Competing interests.

Olivier Gagliardini is a member of the editorial board of the journal. All other authors declare that they have no competing interests.

Acknowledgements

Back to toptop
Acknowledgements.

We would like to thank the editor, Carlos Martin, as well as the two referees,
Mathieu Morlighem and Thomas Kleiner, for their insightful and helpful comments. This
work was funded by the French National Research Agency (ANR) through the
TROIS-AS (ANR-15-CE01-0005-01) project. Julien Brondex, Fabien
Gillet-Chaulet and Olivier Gagliardini are part of
Labex OSUG@2020 (ANR10 LABX56). Most of the computational resources were
provided by CINES through the gge6066 project. Some of the computations were
performed using the Froggy platform of the CIMENT infrastructure
(https://ciment.ujf-grenoble.fr, last access: 17 January 2019), which is supported by the Rhône-Alpes
region (GRANT CPER07-13 CIRA), the OSUG@2020 labex (reference ANR10 LABX56),
and the Equip@Meso project (reference ANR-10-EQPX-29-01) of the programme
Investissements d'Avenir supervised by the ANR.

Edited by: Carlos Martin

Reviewed by: Mathieu Morlighem and Thomas Kleiner

References

Back to toptop
Adhalgeirsdóttir, G., Aschwanden, A., Khroulev, C., Boberg, F., Mottram, R., Lucas-Picher, P., and Christensen, J.: Role of model initialization for projections of 21st-century Greenland ice sheet mass loss, J. Glaciol., 60, 782–794, 2014. a

Alley, R. B., Blankenship, D. D., Bentley, C. R., and Rooney, S.: Deformation of till beneath ice stream B, West Antarctica, Nature, 322, 57–59, 1986. a, b

Blankenship, D. D., Bentley, C. R., Rooney, S., and Alley, R. B.: Seismic measurements reveal a saturated porous layer beneath an active Antarctic ice stream, Nature, 322, 54–57, 1986. a, b

Bondzio, J. H., Morlighem, M., Seroussi, H., Kleiner, T., Rückamp, M., Mouginot, J., Moon, T., Larour, E. Y., and Humbert, A.: The mechanisms behind Jakobshavn Isbræ's acceleration and mass loss: a 3D thermomechanical model study, Geophys. Res. Lett., 44, 6252–6260, https://doi.org/10.1002/2017GL073309, 2017. a

Borstad, C. P., Rignot, E., Mouginot, J., and Schodlok, M. P.: Creep deformation and buttressing capacity of damaged ice shelves: theory and application to Larsen C ice shelf, The Cryosphere, 7, 1931–1947, https://doi.org/10.5194/tc-7-1931-2013, 2013. a

Bougamont, M., Christoffersen, P., A L, H., Fitzpatrick, A., Doyle, S., and Carter, S.: Sensitive response of the Greenland Ice Sheet to surface melt drainage over a soft bed, Nat. Commun., 5, 5052, https://doi.org/10.1038/ncomms6052, 2014. a

Boulton, G. and Jones, A.: Stability of temperate ice caps and ice sheets resting on beds of deformable sediment, J. Glaciol., 24, 29–43, 1979. a

Brisbourne, A. M., Smith, A. M., Vaughan, D. G., King, E. C., Davies, D., Bingham, R., Smith, E., Nias, I., and Rosier, S. H.: Bed conditions of Pine Island Glacier, West Antarctica, J. Geophys. Res.-Earth Surf., 122, 419–433, 2017. a

Brondex, J., Gagliardini, O., Gillet-Chaulet, F., and Durand, G.: Sensitivity of grounding line dynamics to the choice of the friction law, J. Glaciol., 63, 854–866, 2017. a, b, c, d, e, f, g, h, i, j, k, l, m

Budd, W., Keage, P., and Blundy, N.: Empirical studies of ice sliding, J. Glaciol., 23, 157–170, 1979. a, b

Bueler, E. and van Pelt, W.: Mass-conserving subglacial hydrology in the Parallel Ice Sheet Model version 0.6, Geosci. Model Dev., 8, 1613–1635, https://doi.org/10.5194/gmd-8-1613-2015, 2015. a

Colleoni, F., De Santis, L., Siddoway, C. S., Bergamasco, A., Golledge, N. R., Lohmann, G., Passchier, S., and Siegert, M. J.: Spatio-temporal variability of processes across Antarctic ice-bed–ocean interfaces, Nat. Commun., 9, 2289, https://doi.org/10.1038/s41467-018-04583-0, 2018. a

Cornford, S. L., Martin, D. F., Payne, A. J., Ng, E. G., Le Brocq, A. M., Gladstone, R. M., Edwards, T. L., Shannon, S. R., Agosta, C., van den Broeke, M. R., Hellmer, H. H., Krinner, G., Ligtenberg, S. R. M., Timmermann, R., and Vaughan, D. G.: Century-scale simulations of the response of the West Antarctic Ice Sheet to a warming climate, The Cryosphere, 9, 1579–1600, https://doi.org/10.5194/tc-9-1579-2015, 2015. a, b

Cuffey, K. M. and Paterson, W. S. B.: The physics of glaciers, Academic Press, 4th Edn., Butterworth-Heinemann, Oxford, 2010. a, b, c

de Fleurian, B., Gagliardini, O., Zwinger, T., Durand, G., Le Meur, E., Mair, D., and Råback, P.: A double continuum hydrological model for glacier applications, The Cryosphere, 8, 137–153, https://doi.org/10.5194/tc-8-137-2014, 2014. a

Favier, L., Durand, G., Cornford, S., Gudmundsson, G., Gagliardini, O., Gillet-Chaulet, F., Zwinger, T., Payne, A., and Le Brocq, A.: Retreat of Pine Island Glacier controlled by marine ice-sheet instability, Nat. Clim. Change, 4, 117–121, 2014. a

Fretwell, P., Pritchard, H. D., Vaughan, D. G. et al.: Bedmap2: improved ice bed, surface and thickness datasets for Antarctica, The Cryosphere, 7, 375–393, https://doi.org/10.5194/tc-7-375-2013, 2013. a, b

Fricker, H. A., Scambos, T., Bindschadler, R., and Padman, L.: An active subglacial water system in West Antarctica mapped from space, Science, 315, 1544–1548, 2007. a

Fürst, J. J., Durand, G., Gillet-Chaulet, F., Merino, N., Tavard, L., Mouginot, J., Gourmelen, N., and Gagliardini, O.: Assimilation of Antarctic velocity observations provides evidence for uncharted pinning points, The Cryosphere, 9, 1427–1443, https://doi.org/10.5194/tc-9-1427-2015, 2015. a, b

Fürst, J. J., Durand, G., Gillet-Chaulet, F., Tavard, L., Rankl, M., Braun, M., and Gagliardini, O.: The safety band of Antarctic ice shelves, Nat. Clim. Change, 6, 479–482, 2016. a

Gagliardini, O. and Werder, M. A.: Influence of increasing surface melt over decadal timescales on land-terminating Greenland-type outlet glaciers, J. Glaciol., 64, 1–11, https://doi.org/10.1017/jog.2018.59, 2018. a

Gagliardini, O., Cohen, D., Råback, P., and Zwinger, T.: Finite-element modeling of subglacial cavities and related friction law, J. Geophys. Res.-Earth Surf., (2003–2012), 112, F02027, https://doi.org/10.1029/2006JF000576, 2007. a

Gagliardini, O., Zwinger, T., Gillet-Chaulet, F., Durand, G., Favier, L., de Fleurian, B., Greve, R., Malinen, M., Martín, C., Råback, P., Ruokolainen, J., Sacchettini, M., Schäfer, M., Seddik, H., and Thies, J.: Capabilities and performance of Elmer/Ice, a new-generation ice sheet model, Geosci. Model Dev., 6, 1299–1318, https://doi.org/10.5194/gmd-6-1299-2013, 2013. a

Gilbert, J. C. and Lemaréchal, C.: Some numerical experiments with variable-storage quasi-Newton algorithms, Math. Programm., 45, 407–435, 1989. a

Gillet-Chaulet, F., Gagliardini, O., Seddik, H., Nodet, M., Durand, G., Ritz, C., Zwinger, T., Greve, R., and Vaughan, D. G.: Greenland ice sheet contribution to sea-level rise from a new-generation ice-sheet model, The Cryosphere, 6, 1561–1576, https://doi.org/10.5194/tc-6-1561-2012, 2012. a, b, c

Gillet-Chaulet, F., Durand, G., Gagliardini, O., Mosbeux, C., Mouginot, J., Rémy, F., and Ritz, C.: Assimilation of surface velocities acquired between 1996 and 2010 to constrain the form of the basal friction law under Pine Island Glacier, Geophys. Res. Lett., 43, 10311–10321, doi10.1002/2016GL069937, 2016. a

Gray, L., Joughin, I., Tulaczyk, S., Spikes, V. B., Bindschadler, R., and Jezek, K.: Evidence for subglacial water transport in the West Antarctic Ice Sheet through three-dimensional satellite radar interferometry, Geophys. Res. Lett., 32, L03501, https://doi.org/10.1029/2004GL021387, 2005. a

Gudmundsson, G. H., Krug, J., Durand, G., Favier, L., and Gagliardini, O.: The stability of grounding lines on retrograde slopes, The Cryosphere, 6, 1497–1505, https://doi.org/10.5194/tc-6-1497-2012, 2012. a, b

Haseloff, M. and Sergienko, O. V.: The effect of buttressing on grounding line dynamics, J. Glaciol., 64, 1–15, https://doi.org/10.1017/jog.2018.30, 2018. a

Hewitt, I.: Seasonal changes in ice sheet motion due to melt water lubrication, Earth Planet. Sci. Lett., 371, 16–25, 2013. a

Hewitt, I., Schoof, C., and Werder, M.: Flotation and free surface flow in a model for subglacial drainage. Part 2. Channel flow, J. Fluid Mechan., 702, 157–187, 2012. a

Iken, A.: The effect of the subglacial water pressure on the sliding velocity of a glacier in an idealized numerical model, J. Glaciol., 27, 407–421, 1981. a

Joughin, I., MacAyeal, D. R., and Tulaczyk, S.: Basal shear stress of the Ross ice streams from control method inversions, J. Geophys. Res.-Solid Earth, 109, 3 B09405, https://doi.org/10.1029/2003JB002960, 2004. a

Joughin, I., Tulaczyk, S., Bamber, J. L., Blankenship, D., Holt, J. W., Scambos, T., and Vaughan, D. G.: Basal conditions for Pine Island and Thwaites Glaciers, West Antarctica, determined using satellite and airborne data, J. Glaciol., 55, 245–257, 2009. a

Joughin, I., Smith, B. E., and Medley, B.: Marine ice sheet collapse potentially under way for the Thwaites Glacier Basin, West Antarctica, Science, 344, 735–738, 2014. a

Kavanaugh, J. L. and Clarke, G. K.: Discrimination of the flow law for subglacial sediment using in situ measurements and an interpretation model, J. Geophys. Res.-Earth Surf., 111, F01002, https://doi.org/10.1029/2005JF000346, 2006. a

Le Brocq, A. M., Ross, N., Griggs, J. A. et al.: Evidence from ice shelves for channelized meltwater flow beneath the Antarctic Ice Sheet, Nat. Geosci., 6, 945–948, 2013. a

MacAyeal, D. R.: Large-scale ice flow over a viscous basal sediment: Theory and application to ice stream B, Antarctica, J. Geophys. Res.-Solid Earth (1978–2012), 94, 4071–4087, 1989. a

MacAyeal, D. R.: A tutorial on the use of control methods in ice-sheet modeling, J. Glaciol., 39, 91–98, 1993. a, b

Millan, R., Rignot, E., Bernier, V., Morlighem, M., and Dutrieux, P.: Bathymetry of the Amundsen Sea Embayment sector of West Antarctica from Operation IceBridge gravity and other data, Geophys. Res. Lett., 44, 1360–1368, 2017. a

Minchew, B. M., Meyer, C. R., Robel, A. A., Gudmundsson, G. H., and Simons, M.: Processes controlling the downstream evolution of ice rheology in glacier shear margins: case study on Rutford Ice Stream, West Antarctica, J. Glaciol., 64, 1–12, https://doi.org/10.1017/jog.2018.47, 2018. a

Morlighem, M., Rignot, E., Seroussi, H., Larour, E., Ben Dhia, H., and Aubry, D.: Spatial patterns of basal drag inferred using control methods from a full-Stokes and simpler models for Pine Island Glacier, West Antarctica, Geophys. Res. Lett., 37, L14502, https://doi.org/10.1029/2010GL043853, 2010. a, b

Morlighem, M., Seroussi, H., Larour, E., and Rignot, E.: Inversion of basal friction in Antarctica using exact and incomplete adjoints of a higher-order model, J. Geophys. Res.-Earth Surf., 118, 1746–1753, 2013. a, b

Mouginot, J., Rignot, E., and Scheuchl, B.: Sustained increase in ice discharge from the Amundsen Sea Embayment, West Antarctica, from 1973 to 2013, Geophys. Res. Lett., 41, 1576–1584, 2014. a

Nias, I., Cornford, S., and Payne, A.: New mass-conserving bedrock topography for Pine Island Glacier impacts simulated decadal rates of mass loss, Geophys. Res. Lett., 45, 3173–3181, 2018. a

Nowicki, S. M. J., Payne, A., Larour, E., Seroussi, H., Goelzer, H., Lipscomb, W., Gregory, J., Abe-Ouchi, A., and Shepherd, A.: Ice Sheet Model Intercomparison Project (ISMIP6) contribution to CMIP6, Geosci. Model Dev., 9, 4521–4545, https://doi.org/10.5194/gmd-9-4521-2016, 2016. a

Pollard, D. and DeConto, R. M.: Description of a hybrid ice sheet-shelf model, and application to Antarctica, Geosci. Model Dev., 5, 1273–1295, https://doi.org/10.5194/gmd-5-1273-2012, 2012. a, b

Rignot, E., Mouginot, J., and Scheuchl, B.: Ice flow of the Antarctic ice sheet, Science, 333, 1427–1430, 2011. a

Rignot, E., Mouginot, J., Morlighem, M., Seroussi, H., and Scheuchl, B.: Widespread, rapid grounding line retreat of Pine Island, Thwaites, Smith, and Kohler glaciers, West Antarctica, from 1992 to 2011, Geophys. Res. Lett., 41, 3502–3509, 2014. a

Ritz, C., Edwards, T. L., Durand, G., Payne, A. J., Peyaud, V., and Hindmarsh, R. C.: Potential sea-level rise from Antarctic ice-sheet instability constrained by observations, Nature, 528, 115–118, 2015. a

Schoof, C.: The effect of cavitation on glacier sliding, P. Roy. Soc. A, 461, 609–627, 2005. a, b, c, d, e

Schoof, C.: Ice-sheet acceleration driven by melt supply variability, Nature, 468, 803–806, 2010. a

Seroussi, H., Morlighem, M., Rignot, E., Larour, E., Aubry, D., Ben Dhia, H., and Kristensen, S. S.: Ice flux divergence anomalies on 79north Glacier, Greenland, Geophys. Res. Lett., 38, L09501, https://doi.org/10.1029/2011GL047338, 2011. a

Seroussi, H., Morlighem, M., Larour, E., Rignot, E., and Khazendar, A.: Hydrostatic grounding line parameterization in ice sheet models, The Cryosphere, 8, 2075–2087, https://doi.org/10.5194/tc-8-2075-2014, 2014. a

Seroussi, H., Nakayama, Y., Larour, E., Menemenlis, D., Morlighem, M., Rignot, E., and Khazendar, A.: Continued retreat of Thwaites Glacier, West Antarctica, controlled by bed topography and ocean circulation, Geophys. Res. Lett., 44, GL072910, https://doi.org/10.1002/2017GL072910, 2017. a

Shepherd, A., Ivins, E., Rignot, E. et al.: Mass balance of the Antarctic ice sheet from 1992 to 2017, Nature, 558, 219–222, https://doi.org/10.1038/s41586-018-0179-y, 2018. a

Smith, A. M., Jordan, T. A., Ferraccioli, F., and Bingham, R. G.: Influence of subglacial conditions on ice stream dynamics: seismic and potential field data from Pine Island Glacier, West Antarctica, J. Geophys. Res.-Solid Earth, 118, 1471–1482, 2013. a

Tsai, V. C., Stewart, A. L., and Thompson, A. F.: Marine ice-sheet profiles and stability under Coulomb basal conditions, J. Glaciol., 61, 205–215, 2015. a, b

Van Liefferinge, B. and Pattyn, F.: Using ice-flow models to evaluate potential sites of million year-old ice in Antarctica, Clim. Past, 9, 2335–2345, https://doi.org/10.5194/cp-9-2335-2013, 2013. a

Weertman, J.: On the sliding of glaciers, J. Glaciol., 3, 33–38, 1957. a, b, c

Werder, M. A., Hewitt, I. J., Schoof, C. G., and Flowers, G. E.: Modeling channelized and distributed subglacial drainage in two dimensions, J. Geophys. Res.-Earth Surf., 118, 2140–2158, 2013. a

Zwally, H. J., Giovinetto, M. B., Beckley, M. A., and Saba, J. L.: Antarctic and Greenland drainage systems, GSFC Cryospheric Sciences Laboratory, available at: http://icesat4.gsfc.nasa.gov/cryo_data/ant_grn_drainage_systems.php (last access: 17 January 2019), 2012. a