The predictive power of ice sheet models and the regional sensitivity of ice loss to basal sliding parameterisations: a case study of Pine Island and Thwaites glaciers, West Antarctica

Ice sheet models use a wide range of sliding laws to define a relationship between ice velocity and basal drag, generally comprising some combination of a Weertman-style power law and Coulomb friction. The exact nature of basal sliding is not known from observational data, making assessment of the suitability of different sliding laws difficult. The question of how much this choice could affect predictions of future ice sheet evolution is an important one. Here we conduct a model study of a large sector of the West Antarctic Ice Sheet (WAIS), a particularly critical component of the cryosphere, using a range of sliding parameterisations, and we provide an assessment of the sensitivity of ice loss to the choice of sliding law. We show that, after initialisation, various sliding laws result in broadly similar ranges of sea level contribution over 100 years, with the range primarily dependent on exact parameter values used in each sliding law. Comparing mass loss from Thwaites and Pine Island glaciers and the neighbouring regions reveals significant qualitative geographical differences in the relationship between sliding parameters and the modelled response to changes in forcing. We show that the responses do not necessarily follow universal systematic patterns, and, in particular, higher values of the sliding exponent m do not necessarily imply larger rates of mass loss. Despite differences in the magnitudes of ice loss and rates of change in the system, all of our experiments display broad similarities in behaviour which serve to reinforce the decade-to-century-scale predictive power of ice sheet models, regardless of the choice of basal sliding.

1 Introduction

For fast-flowing areas of the Antarctic Ice Sheet, a significant portion of the forward motion is due to sliding along the bedrock beneath them. Ice sheet models must mathematically represent this sliding in order to simulate the physical system, but the exact nature of the relationship between a glacier and its bedrock cannot be determined from observational data and is not fully understood. This process enters models in the form of a sliding law describing the form of basal drag. In general, these laws express basal drag in terms of the basal velocity and/or effective pressure, the difference between the hydrostatic overburden of ice and the basal water pressure. Sliding laws contain some parameters whose values are unknown and are prescribed as part of setting up a model. As a result, they can vary between models and between simulations, giving different results in future predictions. While there are some common choices for these parameters, it is unknown which are the most representative of the physics of the ice sheet system.

An early sliding law, and one which is still widely used in modelling, is that of Weertman (1957). This is a power law of the form

$${\mathit{v}}_{\mathrm{b}}=C{\mathit{\tau }}_{\mathrm{b}}^m,$$

where τ_b is the basal drag, C is a sliding parameter, v_b is the basal velocity and m is a sliding exponent. We refer to this as the Weertman sliding law.

Throughout its history, the Weertman law has been criticised for neglecting the effects of cavitation, in which water under the ice can affect the basal drag. Lliboutry (1968) proposed the idea that a sliding law should depend on effective pressure. This can be done using Coulomb's law of solid friction, giving the relationship

$${\mathit{\tau }}_{\mathrm{b}}={\mathit{\mu }}_{k}N,$$

where μ_k is a coefficient of friction and N is the effective pressure.

The idea of including effective pressure in the sliding law has been explored further. Budd et al. (1979) proposed a law of the form

$${\mathit{v}}_{\mathrm{b}}=N^{-q}{C}_B{\mathit{\tau }}_{\mathrm{b}}^m,$$

where C_B is a sliding parameter and q is a positive constant. Fowler (1987) produced a sliding law with broadly similar behaviour to that of Budd but which approached Coulomb friction at high velocities. Schoof (2005) asserts that laws taking the form of the Budd law are unphysical as they allow arbitrarily large stresses at the bed. Motivated by the work of Iken (1981), Schoof (2005) derived a regularised Coulomb sliding law, stated to be essentially the same as that of Fowler (1987). Variations of this have been used in several models since (e.g. Gagliardini et al., 2007; Leguy et al., 2014; Cornford et al., 2020).

Another approach, suggested by Tsai et al. (2015), is to compute the basal drag using both a Weertman sliding law and a Coulomb friction law and to take the minimum value of the two at each point. We refer to this as the Tsai sliding law. The Coulomb friction is applicable to softer beds and is more representative of sliding near the grounding line. Taking the minimum value of this and the Weertman law effectively applies Coulomb friction near the grounding line and Weertman sliding elsewhere.

With several sliding laws available and parameter options within them, it is useful for modellers to know more about the possible implications of their choice. Previous studies have compared sliding laws in different ways. Brondex et al. (2017) compare the Weertman, Budd, regularised Coulomb and Tsai sliding laws in an idealised setting, finding different levels of transient response in each case. This was followed up in Brondex et al. (2019) by simulations on a West Antarctic domain using linear and non-linear versions of the Weertman and Budd laws and two friction coefficients in the regularised Coulomb law. In both studies, the Weertman law was found to produce the smallest changes in volume above flotation (VAF), while the Budd law caused the largest changes. Simulations using the non-linear versions of the Budd and Weertman sliding laws lost ice faster than those using the linear versions.

The effect of the value of m in the Weertman sliding law has been investigated in a few cases. In these studies, the sliding parameter C is set for each value of m, usually by using inversion methods, in order to match an initial velocity field. Gillet-Chaulet et al. (2016) test different values and find that higher values of m generally produce higher velocities in forward simulations and a closer match to an observed surface velocity field of Pine Island Glacier (PIG), with m=20 suggested as the best value in this case. Joughin et al. (2019) compare experiments using the Weertman law with different values of m to a regularised Coulomb law, focusing on a profile along PIG. The regularised Coulomb law was found to reproduce observations with greater success, and higher values of m in the Weertman law approached a better result. De Rydt et al. (2021) explore the idea that m in the Weertman sliding law does not necessarily need to take a uniform value, presenting a map of optimal m values over PIG which come closest to reproducing observed velocity changes. In general, this involves higher values of m on the fastest flowing areas. Sun et al. (2020) compare the response of several models to ice shelf collapse. It is shown that there is a general tendency for models using higher values of m to respond faster to the changes in forcing. This is in agreement with results from within single models in Brondex et al. (2019), Joughin et al. (2019) and De Rydt et al. (2021).

Some examples can be found of the opposite being true, indicating a more complex relationship between the value of m and the rates of change to a system. For Pine Island Glacier in particular, Nias et al. (2018) show that higher values of m can produce higher or lower changes in VAF depending on the bathymetry product being used. This sensitivity to topography is further demonstrated by Wernecke et al. (2022).

The focus of the work presented in this paper is on the effects of parameter choices in sliding laws on the forward evolution of an ice sheet system and the response to changes in ice shelf forcing. The values of m, μ_k, and q are varied, and we seek to quantify the difference these choices can make in the modelled response when different changes in the buttressing of an ice sheet system are applied, as well as whether there are universal systematic patterns when varying sliding parameters. We also set out to assess the degree of confidence with which predictions can be made by an ice sheet model, despite the lack of empirical knowledge of basal conditions to inform a sliding law.

The chosen domain for the experiments presented in this paper is the Amundsen Sea Embayment (ASE) in West Antarctica (Fig. 1), the region of the continent which is experiencing the fastest ice mass loss (Shepherd et al., 2018; Rignot et al., 2019), at a rate which is accelerating (Sutterley et al., 2014). Understanding how the West Antarctic Ice Sheet (WAIS) is evolving is seen as one of the most important questions in Antarctic research (National Academies of Sciences, Engineering, and Medicine, 2015). A lot of past work has been centred around PIG, including the previously mentioned studies. Recently, there has also been a focus on Thwaites Glacier due to the rapid changes observed in recent years (e.g. Mouginot et al., 2014; Milillo et al., 2019), and it is the subject of a large multidisciplinary investigation, the International Thwaites Glacier Collaboration (Scambos et al., 2017). Many models have been used to make predictions of the future evolution of Thwaites Glacier and the ASE region, and they generally agree that the trend of rapid retreat in the ASE will continue into the future (e.g. Favier et al., 2014; Joughin et al., 2014; DeConto and Pollard, 2016; Yu et al., 2018). Models tend to differ on the rates of change they calculate for a variety of reasons (e.g. Cornford et al., 2020), and as such efforts are ongoing to improve the understanding and functionality of all aspects of ice sheet models and to reduce uncertainty in their predictions.

Figure 1 The catchment basins of the Amundsen Sea Embayment (ASE) using MEaSUREs Antarctic boundaries (Rignot et al., 2013) and, inset, the location of the ASE within Antarctica. The different colours indicate the four regions we focus on. The thicker black outline is the boundary of our model domain, and the initial grounding line is indicated in blue.

2 Experimental design

2.1 Model setup and forcing

We use the ice sheet model Úa (Gudmundsson, 2020), which implements the vertically integrated shallow shelf approximation (SSA) (MacAyeal, 1989), and Glen's flow law (Glen, 1958) with exponent n=3. Our domain covers the ASE region including PIG, Thwaites Glacier, and the glaciers flowing into the Dotson and Crosson ice shelves. The inland boundary was created based on the MEaSUREs Antarctic boundaries (Rignot et al., 2013), with smoothing applied. We use the geometry from BedMachine Antarctica v2 (Morlighem et al., 2020), from which we also derive the calving front boundary for our domain. A Dirichlet condition is used to set all velocities along the inland boundary to zero, since the domain boundary generally follows the edges of drainage basins. The ice-front stress boundary condition is applied along the seaward boundaries, which are fixed at the calving front. The mesh resolution is 1.5 km at the grounding line and becomes coarser further upstream. The resolution was chosen for computational efficiency, and test cases were compared to finer meshes to ensure this resolution is sufficient to demonstrate the relative behaviours of sliding laws required for this study.

The densities are given the values of 917 kg m⁻³ for ice and 1027 kg m⁻³ for ocean water, consistent with the BedMachine geometry. The effective pressure is calculated by assuming a perfect hydrological connection to the ocean. The surface mass balance is taken from RACMO2.3p2 (Van Wessem et al., 2018). Where the ice is afloat, we use a simple depth-based melt parameterisation of ocean-induced melt:

$$\begin{array}{}\text{(1)}& M=\left\{\begin{array}{ll}\mathrm{0}& \text{if }z\ge \mathrm{0}\\ -\frac{M_{\max }}{\mathrm{500}}z& \text{if }\mathrm{0}>z>-\mathrm{500}\\ M_{\max }& \text{if }z\le -\mathrm{500},\end{array}\right.\end{array}$$

where M is the basal melt rate in metres per annum (m a⁻¹), M_max is a prescribed maximum melt rate and z is the vertical coordinate in metres, positive upwards with zero at sea level.

Our experiments test four different sliding laws and vary the parameters within them. We use the Weertman, Budd, and Tsai sliding laws introduced in Sect. 1 and a regularised Coulomb law as used in Cornford et al. (2020), which is a power-weighted reciprocal sum of the Weertman and Coulomb sliding laws (see Eq. 6). We present these sliding laws below in a form which makes clear how they relate to each other.

$$\begin{array}{}\text{(2)}& {\displaystyle }\text{Coulomb friction:}{\mathit{\tau }}_C={\mathit{\mu }}_{k}N{\displaystyle \frac{{\mathit{v}}_{\mathrm{b}}}{\Vert {\mathit{v}}_{\mathrm{b}}\Vert }}\text{(3)}& {\displaystyle }\text{Weertman sliding:}{\mathit{\tau }}_W=C^{-\frac{\mathrm{1}}{m}}\Vert {\mathit{v}}_{\mathrm{b}}{\Vert }^{\frac{\mathrm{1}}{m}}{\displaystyle \frac{{\mathit{v}}_{\mathrm{b}}}{\Vert {\mathit{v}}_{\mathrm{b}}\Vert }}\text{(4)}& {\displaystyle }\text{Budd sliding:}{\mathit{\tau }}_B=N^{\frac{q}{m}}{C}_B^{-\frac{\mathrm{1}}{m}}\Vert {\mathit{v}}_{\mathrm{b}}{\Vert }^{\frac{\mathrm{1}}{m}}{\displaystyle \frac{{\mathit{v}}_{\mathrm{b}}}{\Vert {\mathit{v}}_{\mathrm{b}}\Vert }}\text{(5)}& {\displaystyle }\text{Tsai sliding:}{\mathit{\tau }}_T=\min ({\mathit{\tau }}_C,{\mathit{\tau }}_W)\text{(6)}& {\displaystyle }\text{Regularised Coulomb sliding:}{\displaystyle \frac{\mathrm{1}}{{\mathit{\tau }}_R^m}}={\displaystyle \frac{\mathrm{1}}{{\mathit{\tau }}_W^m}}+{\displaystyle \frac{\mathrm{1}}{{\mathit{\tau }}_C^m}}\end{array}$$

Note that Coulomb friction alone is not used in any of our model simulations, meaning that all of our sliding laws contain the sliding parameter C, which incorporates the effects of various thermal and mechanical properties of ice as well as basal roughness (e.g. Cuffey and Paterson, 2010). The parameters are given a range of values, exploring a range around common values used in modelling studies. We use m=1, 3, 5, 7 and 9. Commonly, m=3 is used by ice sheet models, but higher values of m in the Weertman law could produce behaviour closer to that of regularised Coulomb laws. We use q=1 and 2. The value q=2 was originally suggested by Budd et al. (1984), but q=1 is a more common choice in recent studies (e.g. Brondex et al., 2019). A range of friction coefficients have been used in modelling and mostly fall within the range of our chosen values, μ_k=0.25, 0.5 and 0.75.

2.2 Overview of experiments

We use an adjoint inversion method to calculate fields of values for A and C independently for each set of parameters in each sliding law. The method involves referring to observed velocities u_obs and v_obs with observational errors u_err and v_err, from the 2016–2017 year of the dataset presented by Mouginot et al. (2017), to minimise a cost function $\mathcal{J}=I+R$, consisting of a misfit function,

$$\begin{array}{}\text{(7)}& \begin{array}{rl}I& ={\displaystyle \frac{\mathrm{1}}{\mathrm{2}\mathcal{A}}}\int \left(\right(u-u_{\mathrm{obs}})/u_{\mathrm{err}}{)}^{\mathrm{2}}\mathrm{d}\mathcal{A}\\ & +{\displaystyle \frac{\mathrm{1}}{\mathrm{2}\mathcal{A}}}\int \left(\right(v-v_{\mathrm{obs}})/v_{\mathrm{err}}{)}^{\mathrm{2}}\mathrm{d}\mathcal{A},\end{array}\end{array}$$

and a regularisation term,

$$\begin{array}{}\text{(8)}& \begin{array}{rl}R& =\sum _{k=\mathrm{1},\mathrm{2}}{\displaystyle \frac{\mathrm{1}}{\mathrm{2}\mathcal{A}}}\int \left({\mathit{\gamma }}_{s,k}^{\mathrm{2}}\right(\mathrm{\nabla }(p_k-p_{k,\mathrm{prior}}){)}^{\mathrm{2}}\\ & +{\mathit{\gamma }}_{a,k}^{\mathrm{2}}(p_k-p_{k,\mathrm{prior}}{)}^{\mathrm{2}})\mathrm{d}\mathcal{A},\end{array}\end{array}$$

where $\mathcal{A}=\int \mathrm{d}\mathcal{A}$ is the total area of the domain; u and v are the modelled horizontal x and y velocity components, respectively; γ_s,k and γ_a,k are slope and amplitude regularisation parameters; and p₁=log ₁₀A, p₂=log ₁₀C and p_k,prior are prior values, or initial estimates, for the parameters p_k. For the inversions in this study, ${\mathit{\gamma }}_{s,\mathrm{1}}={\mathit{\gamma }}_{s,\mathrm{2}}=\mathrm{1}\times {\mathrm{10}}^{\mathrm{4}}$ m, ${\mathit{\gamma }}_{a,\mathrm{1}}=\mathrm{100}$ and ${\mathit{\gamma }}_{a,\mathrm{2}}=\mathrm{1}$. Further details of the inversion method can be found in Barnes et al. (2021). In most cases, the priors used in these experiments are the outputs of a previous inversion for the Weertman law with m=3, adjusted using an analytical conversion to account for the value of m. The original Weertman inversion used the priors named “Priors2” in Barnes et al. (2021), which are spatially varying based on observed velocities and a temperature field. A change had to be made to the value of p_2,prior in Eq. (8) in order for the inversion process to run successfully in the m=9 cases for the regularised Coulomb law with μ_k=0.5 and 0.75 and for the Tsai law with μ_k=0.75. In these cases, the inverted C field for the Tsai law with m=9 and μ_k=0.5 is used as the prior. The mean absolute speed misfit over the domain is within the range 27–52 m a⁻¹ for all inversions.

After optimising, the sliding parameters are fixed and each simulation is run for 20 years using M_max=100. This relaxation period allows ample time to move past any initial transient effects, which are fairly common in ice flow models (e.g. Pattyn et al., 2012) and in Úa only have a significant influence on the output during the first few years of model evolution. We do not introduce different forcing from the beginning of the run, as we wish to assess the effects of a change introduced to a system which is already in motion. Changes are introduced after 20 years of simulation, and the experiments are split into four branches from this point.

• Control branch. The simulation continues with the same forcing.

• Increased melt branch. There is an immediate step change increase in the melt rate, with the value of M_max raised to 200.

• Decreased melt branch. Melting stops immediately, with M_max lowered to 0.

• Ice shelf removed branch. All the floating ice is immediately removed and no more is allowed to form for the remainder of the run. This is achieved by setting the melt rate equal to the ice shelf thickness instead of using Eq. (1).

These four branches are run for a further 100 years. These are not realistic perturbations but are designed to test the most extreme changes to the ice shelf buttressing, as well as one more moderate case, in order to clearly see the differences caused by altering the sliding parameterisation.

Figure 2 VAF loss in the control branch for every set of parameter choices in the sliding laws.

Table 1 The VAF loss (mm sea level rise equivalent) at the end of the control branch simulations for all experiments. Values marked with an asterisk (*) indicate where a different prior had to be used in initialisation.

3 Results

The experiments in the control branch all display loss of ice from the system, with the loss of volume above flotation (VAF) over the 120 years of simulation being between 37 and 104 mm sea level equivalent, depending on the sliding law and parameter choices (Fig. 2 and Table 1). The mean ice loss is 59.5 mm sea level equivalent, and the median is 59.3 mm. There is a clear pattern of lower values of m causing smaller VAF changes in all sliding laws other than Budd. The regularised Coulomb and Tsai laws produce very similar curves, becoming more similar with higher values of m. When using the same values of m, both produce larger mass loss than the Weertman law even for the highest friction coefficient in our ensemble (μ_k=0.75), although the gap between the two becomes significantly smaller with higher values of m. There is also a systematic pattern of higher values of μ_k (i.e. higher friction) causing less change in volume above flotation. In a few cases, the curves do not fit the pattern, due to the use of a different prior during initialisation. These particular cases should perhaps be discounted when assessing patterns in the results but are presented here to demonstrate a dependence of results on the initialisation.

The Budd sliding law shows less systematic behaviour with respect to the value of m, but a higher value of q produces a larger change in VAF, for values of m up to 5. When m becomes too high, the Budd sliding law produces very high ice loss, as demonstrated when m=7 and q=1 in Fig. 2. Raising m further soon makes the sliding law unusable. For both values of q in our ensemble, simulations using m=9 ran into convergence issues early on and could not be completed.

Figure 3 Additional VAF loss compared to the control branch after the changes in forcing are introduced for the entire ensemble of sliding law and parameter choices.

Table 2 The percentage differences in VAF at the end of simulation compared to the control branch for all experiments. Values marked with an asterisk (*) indicate where a different prior had to be used in initialisation.

Turning our attention to the responses to changes in ice shelf forcing, the differences in VAF loss for the entire ensemble of experiments are displayed in Fig. 3 and Table 2. The results are presented with reference to the control branch in each case. Since the control branch is run for 20 years before any changes are made to the forcing, the system state at the point of perturbation will be slightly different in each case. It is therefore useful to look at VAF differences in terms of a percentage change compared to the control branch rather than as a volume in square kilometres (km²). The Weertman sliding law displays the smallest response to all changes, while at a glance the other three sliding laws appear to respond quite similarly, when using the same values of m.

For the decreased melt branch, it is notable in most panels of Fig. 3 that a larger difference is seen between linear sliding (m=1) and the non-linear versions (m>1) than between any two non-linear cases, which all respond very similarly. For the Weertman law, the VAF loss difference at the end of the simulation compared to the control branch is −23.4 % for m=1 and within 1 percentage point of 33.6 % for all other values of m. For the Tsai and regularised Coulomb laws, this gap between the m=1 response and other values of m is larger for higher values of μ_k, but the difference from the control branch is larger for lower values of μ_k. The difference in response with varying friction coefficient is not linear. As an example, when m=3, the VAF loss difference for the regularised Coulomb law is −52.3 % with μ_k=0.25, −43.4 % for μ_k=0.5 and −39.9 % for μ_k=0.75. Taking μ_k=0.5 as a midpoint, lowering the friction causes a significantly larger difference in the response than raising it by the same amount. For the Budd law, there is more variability and a slightly larger percentage difference with q=2 compared to q=1. The effect of m is also reversed compared to the other laws, with m=1 responding more than higher values.

The increased melt branch, as the less extreme of the three forcing changes, generally causes the lowest differences in the responses. The VAF loss differences using the Weertman law range from 3.1 % to 4.8 %. The obvious outliers are in the regularised Coulomb and Tsai laws with low friction (μ_k=0.25), where the simulations using m=1 start to display an increased response towards the end of the run, and the Budd law for which there is more variability in the responses for different values of m, particularly when q=1. Using the regularised Coulomb law with m=3 as an example again, the response changes with μ_k in a similar way to that seen in the decreased melt branch. For μ_k=0.25, the VAF loss is 13.3 % greater than in the control branch, and increasing the friction to μ_k=0.5 and μ_k=0.75 decreases the difference to 9.3 % and 7.7 %, respectively.

The largest variability is found in the ice shelf removed branch. The Budd sliding law shows similar behaviour here as in the other two branches, with lower values of m and higher values of q causing larger responses. The other sliding laws all display a clearly new behaviour. In all cases, the general trend seen in other branches holds initially, with higher values of m causing greater changes to the VAF. However, the response for lower values of m in this branch becomes much larger at some point during each simulation. This point is reached earlier with lower values of μ_k in the Tsai and regularised Coulomb laws and earliest when using the Tsai law, which produces just over twice the VAF loss of the control branch when m=1. Within these two laws, the difference is far more evident with m=1. Using the Weertman law, the difference in response is lower overall and very similar for m=1 and m=3, which have a VAF loss 23.7 % and 22.3 % higher than the control branch, respectively, while the other values of m are in the range 11.4 %–12.9 %. The effects of altering μ_k are similar to those in the other branches, with the difference from the control branch higher with μ_k=0.25 than when μ_k=0.5 and the difference when μ_k=0.75 a little lower.

Figure 4 Grounding lines at 120 years in the ice shelf removed branch for (a) all sliding laws using m=3, q=1, and μ_k=0.5 and (b) different values of m in the Weertman sliding law. The initial grounding line (t=0) is also indicated.

The VAF values are just one metric by which one can assess the system. Another option is to look at the positions of the grounding lines at the end of simulations. Figure 4 shows the grounding lines for all sliding laws using the values m=3, q=1, and μ_k=0.5 and for all values of m within the Weertman sliding law in the ice shelf removed branch. Since this shows the response to the most extreme forcing change, this is the case for which differences are most evident. When comparing the sliding laws using the same parameter values, it is clear that the Weertman law causes the least retreat and that the Tsai and regularised Coulomb laws cause the most, ending up in very similar positions. For the most part, the grounding line produced by the Budd law lies somewhere between the two. When looking at the grounding lines for the Weertman law using different parameters, it can be seen that the grounding line retreats less for higher values of m. This trend is much more noticeable around the Pine Island and Kohler glaciers, while the grounding lines are located in very similar positions around Thwaites Glacier and the Crosson Ice Shelf.

Figure 5 Differences in final speed (a–d) and ice thickness (e–h) in the ice shelf removed branch for the Weertman sliding law. All other values of m are compared to the case where m=1. Areas which are not grounded in either case are not shaded as part of each comparison.

We can also look directly at model outputs of velocities and ice thickness. Figure 5 displays the differences in speed and thickness at the end of the ice shelf removed branch for each non-linear version of the Weertman law, compared to that of the linear case (m=1). Across much of the domain, the ice is thinner and moves faster when higher values of m are used in the sliding law. However in some areas, specifically Pine Island and Kohler glaciers, the opposite is true. The differences between experiments using successive values of m become less as m increases, so the largest difference is between the linear case and m=3, while the difference between the m=7 and m=9 cases is much smaller. We note that there is evidence of numerical artefacts in the thickness field upstream where the resolution is coarse. When compared to a finer mesh, this choice of resolution does not change our findings regarding the patterns of behaviour due to sliding law choices, but the numerical artefacts do indicate that the experiments in this study should not be treated as predictive.

4 Discussion

4.1 The whole domain

As outlined in the previous section, the results display fairly systematic behaviour in many respects. Some of the patterns are simple, such as higher friction coefficients μ_k resulting in smaller changes to the system. Other trends are more complicated. In particular, there is an apparent inconsistency in the effects of the value of m. When the melt is decreased, it is clear that experiments using higher values of m respond more to the change, with the system gaining more ice when m>1. However, with an increased melt rate, and especially with the removal of the entire ice shelf, low values of m cause a larger response.

The Budd sliding law responds in ways which set it apart as an outlier among those tested. This is most likely because it relies on the value of effective pressure throughout the domain, which is calculated assuming a perfect hydrological connection with the ocean. This is a reasonable assumption close to the grounding line but becomes unlikely further upstream (e.g. Cuffey and Paterson, 2010), especially where the bedrock is above sea level and subglacial water pressure is set to zero. With more accurate information about the effective pressure at the base of the ice sheet, the Budd law would likely produce results more consistent with the other sliding laws. However, in our case the assumption causes erratic behaviour, and the Budd law is not considered in detail during our discussion. The Tsai and regularised Coulomb sliding laws also include N, but in these laws the Coulomb component of the equation is only relevant in close proximity to the grounding line, with a transition into Weertman-style sliding upstream.

Excluding the Budd sliding law, all of the control branch experiments (Fig. 2) show a similar evolution in the volume of ice contained within the domain. After 120 years of transient simulation, the loss of volume above flotation ranges from about 37 to 71 mm sea level rise equivalent, with a mean value of 57.05 mm. The standard deviation is 8.68 mm, giving a coefficient of variation of 0.15, which implies moderate variation in rates of change in the system. But no matter which sliding law or parameters are used, the direction of evolution is always the same, suggesting that the choice of sliding parameterisation does not introduce a level of uncertainty which would affect predictions of the broader patterns of future ice movement in a system forced by ice shelf melt rates in a realistic range. Indeed, it has been shown by Hill et al. (2021) that far larger uncertainties are generally introduced by oceanic forcing than by parameters of the ice flow dynamics.

The difference in rates of change when varying basal sliding conditions can be seen more clearly when perturbations are introduced. Again, all sliding laws react in broadly the same way, but with different rates and magnitudes of change. From this, we can say that the choice of sliding law does not impact predictions of how an ice sheet system will change in terms of the direction of evolution, but it does affect how fast the changes occur. It is also worth noting that the large differences in behaviour for low values of m in the ice shelf removed branch are due to an unrealistic extreme forcing accelerating the rate of change in the system and that the results in the more moderate increased melt branch do not differ so much.

In the decreased melt branch, the grounding line advances beyond its initial position, where C was not inverted for and the prior value could have an influence on the results. However, due to the 20 years of model evolution before the perturbation is introduced, the initial responses to the perturbation are unaffected by this. With the exception of the Budd law, none of our experiments show a change in behaviour midway through the simulation.

The significantly larger percentage rises in VAF resulting from using lower values of m cannot be accounted for by differences in the system states after the first 20 years of simulation. They demonstrate differences between the rates at which the forward simulations respond after the perturbations are introduced. It is likely that the overall trend is driven by the large speed reductions of Pine Island as the value of m increases. The grounding line positions (Fig. 4) and differences in ice thickness and speed (Fig. 5) clearly suggest that the modelled response to different m values, as well as the sensitivity to basal sliding representation in general, is affected by regional differences.

Figure 6 Magnitudes of additional VAF loss compared to the control branch for the Weertman sliding law with different values of m. The contributions of each region within the ASE domain are shown.

4.2 Regional variation

To investigate differences within the ASE, we divide the domain into four regions, as indicated by the colours on Fig. 1. These regions are the catchment basins of three glaciers – Pine Island, Thwaites and Kohler – and the combined catchments of the small glaciers which flow out to the Crosson Ice Shelf. Figure 6 shows the magnitude of differences in VAF for these four regions. Pine Island Glacier has the largest overall ice loss and thus contributes heavily to the behaviour seen for the whole domain in the previous section.

Figure 7 shows the percentage differences in VAF loss compared to the control branch for all other branches using the Weertman sliding law. Figure 8 displays the same information in the case of using the Tsai sliding law with μ_k=0.5. The results from the regularised Coulomb law are very similar to those from the Tsai law. We can see a variety of responses across the domain, showing that regional factors can play an important role in how sliding laws react to changes. In terms of percentage increases, Kohler Glacier responds the most by far to changes involving a loss of buttressing, with up to 150 % greater VAF using the Tsai law. The differences in VAF for Thwaites and Crosson at the end of the simulations is very small – less than 10 % for the Weertman law and 30 % for the Tsai law.

Figure 7 Additional VAF loss compared to the control branch for the Weertman sliding law with different values of m, showing the contributions of (a) the entire ASE domain and (b–e) each region within the domain. Note that the y axes for each region use different scales.

Figure 8 Additional VAF loss compared to the control branch for the Tsai sliding law with μ_k=0.5 and different values of m, showing the contributions of (a) the entire ASE domain and (b–e) each region within the domain. Note that the y axes for each region use different scales.

The large changes when the ice shelf is removed and m has a low value are not universal. They occur in the Pine Island and Kohler glaciers – the same locations in which the largest grounding line retreat is seen in Fig. 4. Thwaites Glacier and those flowing into the Crosson Ice Shelf do not display the same behaviour. Looking at flowlines through Pine Island Glacier, we discover that 50 years into the simulations the glacier is grounded on a ridge in the bathymetry. With low values of m, the ice unpins from this point earlier, producing the sudden rise in VAF differences. The region-specific effect at play here could be the bathymetry itself. Nias et al. (2018) found that for Pine Island Glacier, using different bathymetry datasets can affect whether more or less ice is lost with increasing values of m. They conclude that the buttressing force of a bump in one bathymetry dataset is amplified by larger values of m, slowing down the flow. Our bathymetry is a different product to either of those tested in that study, but it is likely that there is a similar cause for the pattern observed in our results and that such bathymetric features do not exist around Thwaites Glacier and Crosson Ice Shelf, so the same effect is not seen in those regions.

The relationship between the value of m and the rates of retreat and ice loss is an issue which could be explored further. Nias et al. (2018) identify local bathymetry as a likely source of differences in behaviour, but there could be other factors to consider, such as variation in the basal sliding parameter C or flow rate factor A. We note that even a small change required for successful inversion in three of our cases caused the responses in those cases to fall out of line with the patterns seen in the rest, suggesting a high sensitivity to changes in the inverted fields. Whatever the cause, we can confidently say that the relationship contains complexities and can be affected by regional differences, so there is no simple statement to be made about what the general effect of increasing m will be on a system. We believe the regional variance in the effect of m supports the idea that m should perhaps be treated not as a universal constant but as a distributed field of values, as suggested by De Rydt et al. (2021). At the very least, different glaciers may need to use different values in order to reproduce their observed rates of change.

In general, whatever differences are seen when increasing the value of m, they are smaller between higher values. This can be seen in the control experiments (Fig. 2), in which the m=7 and m=9 cases are closer together than those of lower sequential values in our ensemble. It is also seen clearly in most panels of Fig. 3 for the decreased melt branch that the m=1 case is more different from the other cases, as well as in the differences in final speed and thickness presented in Fig. 5. This demonstrates that when increasing the value of m in sliding laws, the response converges. With large values of m, the Weertman law is often said to approach the behaviour of a Coulomb friction law in this way.

An observation which appears strange at first is that when a change in forcing increases the melt rate or removes the buttressing ice shelf entirely, higher values of m can cause larger VAF losses but at the same time less retreat of the grounding line. We see this around Thwaites glacier when using the Weertman sliding law, for which the m=9 experiment produces the largest VAF loss in the control branch (Fig. 2a), the largest percentage response in the ice shelf removed branch (Fig. 7c), and the greatest ice speed and thinning (Fig. 5), yet the least grounding line retreat (Fig. 4b). We suggest that when more ice is lost from upstream areas, it crosses the grounding line at a rate faster than the ice shelf is melting away, thus causing a buildup of ice in this region. Because the ice around the grounding line is thicker, the glaciers stay grounded. In this way, it is possible for higher values of m to cause glaciers to lose a greater volume of ice while at the same time retaining a larger grounded area.

In most cases, we see that the responses in each region follow the same patterns using the Weertman and Tsai laws. There are some differences, most notably in the response of Thwaites, for which a significantly different pattern can be seen. However, this is another effect of the greater speed of response using the Tsai law. As can be seen in the grounding lines of Fig. 4, for the Tsai and regularised Coulomb laws, the ice retreats far enough for the Thwaites and Pine Island ice shelves to merge. The two glaciers become much more closely linked, and this has an effect on the evolution of Thwaites. Even before the ice sheets merge, the retreat of the small western section of Pine Island starts cutting into the Thwaites region, as we have defined the regions by present-day boundaries. The same behaviour would be seen in the Weertman case if the model was run for longer. This is a significant example of the importance of differing rates of change caused by the choice of basal sliding parameterisation.

5 Conclusions

Parameter choices in sliding laws have an important effect on the speed and magnitude of responses to changes in the forcing of the ice system, although the exact nature of the effect can be highly region-dependent. The Weertman, Tsai and regularised Coulomb laws all display similar behaviours in the effects of changing the values of m and μ_k. The Budd sliding law responds differently to the other laws, due to the simplistic hydrological condition used in our model, but in general raising the value of q increases the magnitude of the response.

Changing the value of the friction coefficient μ_k produces a very structured response in all cases. Simulations with lower values of μ_k respond more to changes in forcing. The relationship is not linear, with the differences in response compared to the control branch converging as μ_k increases. This is to be expected, as there will be a convergence to the limit at which friction is sufficiently high to allow no basal sliding at all.

Changing the value of the sliding exponent m produces a more varied range of responses and is affected by regional differences in bathymetry as well as the inverted fields A and C. In the case of lowering the ice shelf melt rate and allowing increased buttressing, m=1 is a notable outlier, producing less difference in VAF compared to the control branch, while higher values of m respond very similarly to the change. The response to increased melting or removal of the ice shelf also systematically shows higher values of m responding more in terms of VAF loss, up to a point at which the trend switches. This is due to regional changes at Pine Island and Kohler glaciers, where lower values of m cause more rapid retreat along the main glacier trunks. Thwaites Glacier and those flowing out to the Crosson Ice Shelf do not display this behaviour.

Despite the differences in our results when using different parameterisations of basal sliding, with the highest VAF loss at the end of the control branch being almost twice that of the lowest, there is also a degree of similarity between them. In broad terms, the system always reacts in the same way to a particular change in forcing, no matter which sliding law or parameter values are used. The magnitude and speed of these reactions vary, but the response is clearly bounded. Our results suggest that ice sheet models, if initialised by inverting for sliding law parameters using present-day geometry and surface velocities and driven by forcing in a realistic range, follow approximately the same trajectory for at least several decades and up to a century, irrespective of the choice of sliding law. In this sense, the predictions of ice sheet models about near-future ice loss are reasonably well constrained despite our incomplete knowledge of basal sliding processes. The implication is that initialised ice sheet models therefore have predictive power over shorter timescales of several decades, at least for fast-flowing glaciers such as those in West Antarctica. Further work is required to determine if this holds in other areas or when other poorly constrained processes such as calving and ice damage are included.