1. INTRODUCTION
The Amundsen Sea Embayment (ASE) in West Antarctica is the dominant region of ice mass loss from the Antarctic ice sheet (Shepherd and Wingham, Reference Shepherd and Wingham2007). Ice drains through a number of ice streams, including Pine Island, Thwaites and Smith glaciers (Fig. 1). These ice streams rest on bedrock that is below sea level and generally slopes downwards towards the interior of the ice sheet (which is referred to as retrograde; prograde refers to bedrock sloping upwards towards the interior) (Holt and others, Reference Holt2006; Vaughan and others, Reference Vaughan2006). This configuration could be especially susceptible to grounding line retreat, through marine ice-sheet instability (Weertman, Reference Weertman1974).
The Amundsen Sea Embayment ice streams, with grounding line positions at t = 50 a overlain on the initial (t = 0 a) velocity field of the default r bm simulation. The thick black curve indicates the grounding line position at the start of the simulations. Each coloured curve is an ensemble member from each of the four groups: Lr′ (grey), Lr bm (orange), Pr′ (blue) and Pr bm (yellow). The polygons represent the mask areas used in the regression analysis. The inset shows location of the model domain (dashed red) and this plot area (solid red).
Studies using satellite radar (Flament and Rémy, Reference Flament and Rémy2012) and laser (Pritchard and others, Reference Pritchard, Arthern, Vaughan and Edwards2009) altimetry have shown that the ASE ice streams (predominantly Pine Island Glacier) are thinning at accelerated rates, and that the thinning is propagating up-glacier (Wingham and others, Reference Wingham, Wallis and Shepherd2009). Thinning occurred simultaneously with ice acceleration (Mouginot and others, Reference Mouginot, Rignot and Scheuchl2014). Rignot (Reference Rignot2008) found that Pine Island Glacier accelerated by 42% between 1996 and 2007, and Smith Glacier by 83%. For this same period, the fast flowing portion of Thwaites Glacier did not undergo significant velocity changes but did experience widening, encroaching on the slower eastern portion (Rignot, Reference Rignot2008). Since 2006, the fast flowing portion of Thwaites has sped up by 33% (Mouginot and others, Reference Mouginot, Rignot and Scheuchl2014). The acceleration of the ice streams in the ASE, like the thinning, is not confined to the areas close to the grounding line, but can be detected far inland (up to 250 km inland from the grounding line in the case of Pine Island) (Mouginot and others, Reference Mouginot, Rignot and Scheuchl2014). Observations have shown that the Pine Island Glacier grounding line retreated at an average rate of 0.95 ± 0.09 km a−1 between 1992 and 2011 (Park and others, Reference Park2013).
The observed thinning, acceleration and grounding line retreat indicate that the ASE region is experiencing dynamic changes in ice flow, which have been inferred to be caused by external forcing (Wingham and others, Reference Wingham, Wallis and Shepherd2009). Observational and modelling studies suggest that warmer waters in the sub-ice shelf cavity are the cause of these dynamic changes (Payne, Reference Payne2004; Shepherd, Reference Shepherd2004; Payne and others, Reference Payne2007; Jacobs and others, Reference Jacobs, Jenkins, Giulivi and Dutrieux2011). Enhanced sub-ice shelf melting results in thinning and a reduction in the buttressing effect that the ice shelves have on the ice streams draining into them (Dupont and Alley, Reference Dupont and Alley2005; Pritchard and others, Reference Pritchard2012; Paolo and others, Reference Paolo, Fricker and Padman2015). This can lead to sustained acceleration, and in turn thinning and grounding line retreat.
The contrasting changes in the ASE suggest that the geometric features of the individual ice streams, such as bed topography and ice-shelf configuration, might be important in determining how they respond to external forcings. The Crosson (which Smith Glacier flows into) and Pine Island ice shelves are situated in embayments, so lateral stresses are likely to play an important role in the stability of the upstream grounded ice. If the buttressing effect of the ice shelves is reduced, due to increased melt, these ice streams could accelerate and retreat (Dupont and Alley, Reference Dupont and Alley2005). However, Thwaites Glacier has a less constrained ice shelf, which may offer less buttressing to the grounded ice stream (Parizek and others, Reference Parizek2013) compared with the Pine Island and Crosson ice shelves. Despite this, the results of Joughin and others (Reference Joughin, Smith and Medley2014) indicate that ocean driven melt has the potential to sustain widespread retreat of Thwaites once initiated.
Much of the modelling research on ice stream stability in the ASE has focused on Pine Island Glacier (e.g. Joughin and others, Reference Joughin, Smith and Holland2010; Gladstone and others, Reference Gladstone2012; Favier and others, Reference Favier2014). More recently there has been effort to model Thwaites Glacier (Parizek and others, Reference Parizek2013; Docquier and others, Reference Docquier, Pollard and Pattyn2014; Joughin and others, Reference Joughin, Smith and Medley2014) as it has a large drainage basin and if unstable retreat occurs, it has the potential to cause considerable sea-level rise (Joughin and others, Reference Joughin, Smith and Medley2014). Cornford and others (Reference Cornford2015) found that the onset of grounding line retreat in Thwaites results in the greatest cause of uncertainty in sea-level rise from the ASE in the coming centuries.
In this paper we investigate the sensitivity of Pine Island and Thwaites, along with the group including Smith Glacier that drain in to the Crosson Ice Shelf, to changes in model parameters, in order to understand how and why the ice streams differ in their behaviour. To do this we run an ensemble of ice-sheet model simulations, in which we perturb parameters related to basal traction, ice viscosity stiffening factor and ice shelf melt rate. We also explore the effect of basal topography and sliding laws on the behaviour of the ASE.
2. METHODS
We constructed an ensemble of ice-sheet simulations, whose parameters vary around two optimal parameter sets. In the following section we describe the ice-sheet model, BISICLES. Next we describe the two optimal parameter sets: the first is optimal in the sense that the basal traction coefficient and viscosity stiffening factor are chosen to minimise the mismatch between model and observed velocity; the second additionally modifies the bed geometry to reduce the large flux discrepancies, which have been attributed to incorrect thickness (Morlighem and others, Reference Morlighem2011). Finally, we describe our experimental design used to create and run the collection of unique model simulations that formed our ensemble.
2.1. Ice-sheet model
We used BISICLES, a vertically-integrated ice flow model, based on the ‘L1L2’ model devised by Schoof and Hindmarsh (Reference Schoof and Hindmarsh2010) and described by Cornford and others (Reference Cornford2013). In the model, we assume ice is in hydrostatic equilibrium, so that the upper surface elevation s is related to the ice thickness h and bedrock topography r through
$$s = {\rm max}\left[ {h + r,\;\left( {1 - \displaystyle{{\rho _{\rm i}} \over {\rho _{\rm w}}}} \right)h} \right],$$
where ρ
i and ρ
w are ice and water density, respectively. Ice thickness and horizontal velocity
$\vec u$
satisfy a conservation of mass equation
$$\displaystyle{{\partial h} \over {\partial t}} + \nabla. (\vec uh) = M_{\rm s} - M_{\rm b}, $$
where M s is the surface mass balance and M b is the sub-shelf melt rate; and a stress balance equation
$$\nabla. \left[ {\varphi h\bar \mu (2 \dot{\bf \epsilon} + 2tr(\dot {\bf \epsilon} ){\bf I})} \right] + \vec \tau _{\rm b} = \rho _{\rm i} gh\nabla s,$$
alongside appropriate boundary conditions.
$\dot{\bf \epsilon} $
is the horizontal rate-of-strain tensor
$$\dot{\bf \epsilon} = \displaystyle{1 \over 2}\left[ {\nabla \vec u + (\nabla \vec u)^T} \right]$$
and
I
is the identity tensor. The vertically-integrated effective viscosity
$\varphi h\bar \mu $
is computed from a stiffening factor φ and a vertically varying effective viscosity μ derived from Glen's flow law. μ includes a contribution from vertical shear strains and, given that the flow rate exponent n = 3, satisfies
$$2\mu A(T)(4\mu ^2 \dot{\bf \epsilon} ^2 + \vert \rho _i g(s - z)\nabla s \vert ^2 ) = 1,$$
with a temperature dependent rate factor A(T) computed from the formula of Cuffey and Paterson (Reference Cuffey and Paterson2010). φ accounts for uncertainty in temperature T and the rate factor A(T), as well as macroscopic damage and variation in ice fabric. It is found by solving an inverse problem described by Cornford and others (Reference Cornford2015). Finally, depending on the model configuration, the basal traction
$\vec \tau _{\rm b} $
is assumed to satisfy either a linear viscous relation, where m = 1, or a non-linear power law, where m = 1/3,
$$\vec \tau _{\rm b} = \left\{ {\matrix{ { - C \vert \mu \vert ^{m - 1} \vec u,} \hfill & {h\displaystyle{{\rho _{\rm i}} \over {\rho _{\rm w}}} \gt - r} \hfill \cr {0,} \hfill & {{\rm otherwise}} \hfill \cr}} \right..$$
Like φ, the traction coefficient C is found by solving an inverse problem (see below).
BISICLES uses block-structured adaptive mesh refinement (AMR) to provide fine resolution in the region of the grounding line and shear margins, and coarse resolution elsewhere. Meshes are built up from square cells with spacing Δx = 4, 2, 1, 0.5 and 0.25 km and evolve as the simulations progress. A previous study of Pine Island Glacier showed these meshes are sufficient to track the movement of the grounding line and shear margins (Favier and others, Reference Favier2014).
2.2. Initialising the model
We required numerous data to solve the model equations described above, some of which we took directly from published sources, while others had to be determined indirectly. We took maps of surface elevation (s) and surface mass balance (M s) directly from BedMap2 (Fretwell and others, Reference Fretwell2013), and 3-D temperature from the results of a higher order model (Pattyn, Reference Pattyn2010), covering our domain Ω (Fig. 1 inset). We found the remaining fields (C, φ, h, r and M b) by solving an inverse problem, and applying perturbations to the result.
Values for the basal traction coefficient C and the stiffening factor φ (or equivalently, an enhancement factor) can be found by solving an optimisation problem. We minimised an objective function
$$J = \int_\Omega ( \vert \vec u \vert - \vert \vec u_{{\rm obs}} \vert )^2 {\rm d}\Omega + \alpha _C \int_\Omega (\nabla C)^2 {\rm d}\Omega + \alpha _\varphi \int_\Omega (\nabla \varphi )^2 {\rm d}\Omega,$$
with respect to C(x, y) and φ(x, y), where the first term measures the mismatch between ice-surface modelled speed
$\vert\vec u\vert$
and the observed speed
$\vert\vec u_{{\rm obs}} \vert$
from Rignot and others (Reference Rignot, Mouginot and Scheuchl2011), and the remaining two terms serve to regularise the solution. Methods to solve this inverse problem are well established (MacAyeal, Reference MacAyeal1992; Joughin and others, Reference Joughin2009; Morlighem and others, Reference Morlighem2010).
In principle, it is possible to take thickness data h
bm and bedrock data r
bm from BedMap2, solve the inverse problem described above, and take the results to be the initial state of the ice sheet. In practice, the observed speed and the thickness data are incompatible in the sense that
$\nabla. (\vec uh)$
(and hence ∂h/∂t) is dominated by large amplitude, short wavelength components (e.g. Morlighem and others, Reference Morlighem2011). In particular, a thickening tendency of order 100 m a−1 develops in the region of Pine Island Glacier's grounding line. We could have compensated for this by imposing large amplitude, short wavelength surface mass balance and/or melt rates, but instead we preferred to modify the thickness and bedrock elevation along the lines of Morlighem and others (Reference Morlighem, Rignot, Mouginot, Seroussi and Larour2014) to find a smoother
$\nabla. (\vec uh)$
.
Over grounded ice, we assumed the surface elevation and speed are well-known, and the ice thickness and bedrock elevation are less well-known. For floating ice, we could not adjust the thickness without also altering the surface elevation, but on the other hand, the melt-rate could have large amplitude. Lacking any other knowledge, we determined the melt rate from
$\nabla. (\vec uh)$
but parameterised it to make it spatially smooth and able to evolve to concentrate melt rates close to the grounding line. We set
$$M_{\rm b} (x,y,t) = \left\{ {\matrix{ {M_{\rm G} (x,y)p(x,y,t)} \hfill\cr{\quad+ M_{\rm A} (x,y)(1 - p(x,y,t)),} \hfill & {{\rm floating}} \hfill \cr {0,} \hfill & {{\rm grounded}} \hfill \cr}} \right.$$
where p(x, y, t) = 1 at the grounding line and decays exponentially with distance from it, by solving,
$$p - \lambda ^2 \nabla ^2 p = \left\{ {\matrix{ {1,} \hfill & {{\rm grounded}} \hfill \cr {0,} \hfill & {{\rm elsewhere}} \hfill \cr}} \right.$$
with boundary condition,
$$\nabla p.n = 0.$$
We computed the grounding line localised and ambient coefficients M
G and M
A by taking values of
$\nabla. (\vec uh)$
from regions where p > (1/100) and <(1/100) respectively, smoothing to remove short-wavelength features, and extrapolating over the whole domain. p is recomputed at every time step. A grid cell is only considered to be floating if its centre is floating. This method was employed by Cornford and others (Reference Cornford2015).
We employed an iterative procedure to find compatible C, φ, h, r, s, M G and M A:
-
1. Set h ← h bm, r ← r bm.
-
2. Find C, φ by minimising J.
-
3. Compute M b and M A.
-
4. Evolve h by integrating an equation that represents a compromise between mass conservation and consistency with h bm,
$\displaystyle{{\partial h} \over {\partial t}} + \nabla. (\vec uh) = \alpha _{\rm h} (h_{{\rm bm}} - h) - M_{\rm b} $
, forward in time for 1 a. -
5. Set r ← r − (h(t) − h(t − 1)).
-
6. Repeat 2–5 until
$\nabla. (\vec uh)$
is smooth.
Note that steps 1–5 were not repeated enough to bring the ice sheet into steady state with any given surface mass balance. This iterative procedure primarily modified the bedrock in the region of very large ∂h/∂t (10–199 m a−1) and left the basin as a whole losing mass at a rate of 0.25 mm a−1 sea level equivalent, in line with current observations (McMillan and others, Reference McMillan2014; Medley and others, Reference Medley2014). We will refer to the resulting thickness and bedrock data as h′ and r′.
We found that the bed geometries required some additional adjustments. Firstly, BedMap2 does not provide an observed bathymetry beneath the Crosson Ice Shelf. Given this limitation, grounding line advance experiments could not be performed, so the bedrock was therefore set to −900 m below sea level below the shelf. Secondly, we added a topographic rise that acts as a pinning point to the Thwaites Ice Shelf. This rise is evident in the velocity data (Rignot and others, Reference Rignot, Mouginot, Morlighem, Seroussi and Scheuchl2014), and coincides with the eastern peak in Tinto and Bell (Reference Tinto and Bell2011).
2.3. Construction of the ensemble
We used Latin hypercube sampling to create 64 parameter vectors, each one with a distinct basal traction coefficient, ice viscosity stiffening factor and sub-ice shelf melt rate, which all affect ice flux across the grounding line (Schoof, Reference Schoof2007). While ice flux may also be dependent on surface mass balance on millennial time scales, transport times mean that it is unlikely to play a role on the decadal time scales considered here (Seroussi and others, Reference Seroussi2014). We used a Latin hypercube sampling method in order to capture any non-linear interactions between the parameters, while maximising coverage of the sub-spaces perpendicular to any null space.
The parameters were scaled between a halving or a doubling, i.e. they were adjusted as:
$$F(x,y) = 2^{2(P_F - (1/2))} F_{\rm o} (x,y)$$
where F is the parameter coefficient (traction, C; viscosity, φ; melt rate, M b); P F is a scalar between 0 and 1 (obtained from the Latin hypercube sampling method); and F o is the parameter coefficient found from the inverse problem. We also ran end members, which involved halving (P F = 0) and doubling (P F = 1) each parameter one at a time, while keeping all other parameters at their default values (P F = 0.5). On top of that we created a default member, with P F = 0.5 for all the three parameters, to see how the model responds given no extra perturbation to the system. We made the parameter ranges large to attempt to capture the realistic response somewhere within the parameter space, and did not intend for the extremes of the ranges to necessarily represent plausible scenarios, caused by well-defined physical processes. Note that P F is uniform in space, but that does not necessarily result in a spatially uniform response. For example, decreasing C has the greatest effect on velocity in the fast flowing regions (i.e. the channels, where increased lubrication is more likely to occur). φ affects stress near the grounding line and shear margins, where strain rates are high, and therefore, a uniform factor of change in φ disproportionately affects the velocity in these regions. Applying a uniform factor to M b effectively only results in a change near the grounding line, because M b is only large near the grounding line (Eqns (8)–(10)).
We duplicated the set of viscosity, traction and melt rate parameters across the two sets of geometry data: Bedmap2 (h bm, r bm) and the modified geometry (h′, r′). Our intention here was to investigate the effect of bedrock topography on the ice dynamics, and to that end we could have constructed (h, r) in a similar fashion to the other parameters. However, in preliminary experiments, where each (h, r) was a linear combination of (h bm, r bm) and (h′, r′), we observed a categorical response, with, for example, an additional stable position for the Pine Island Glacier grounding line when (h, r) was closer to (h′, r′). That being the case, response to variations in the bedrock would not be amenable to the same kind of empirical modelling as the response to the variations in the other parameters, so we chose to treat the bedrock as a categorical parameter from the outset.
We also treat the sliding law as a categorical parameter. For every combination of viscosity, traction, melt rate and geometry parameters, we carry out simulations with a linear sliding law (m = 1) and a shear-thinning power sliding law (m = 1/3). Again, it would be possible to construct an ensemble that varied m continuously, but the existing literature tends to treat it as categorical (Joughin and others, Reference Joughin, Smith and Holland2010), not least because the different laws arise from different physical pictures of the dynamics at the base of the ice.
We will refer to the four groups of experiments as Lr′ and Lr bm for the linear sliding law experiments, and Pr′ and Pr bm for the non-linear power sliding law experiments. For each of these groups, we ran the inversion (described in the previous section) once for each geometry and sliding law, to give the initial condition of the ‘default’ simulation. Then we used each F(x, y) defined by the Latin hypercube to perturb the model. We did not rerun the inversion using each F(x, y). In addition, we ran the end members of the parameter perturbations with both max[r bm, r′] and min[r bm, r′] (using the linear sliding law only). In total, this resulted in an ensemble of 298 individual model configurations (Table 1), which were each used to simulate the ASE for 50 a. We ran these simulations using the University of Bristol's BlueCrystal Phase 3 high-performance computer facility. Each run took between 24 and 60 h to run on 64 processors.
Summary of model simulations in the ensemble
The first three columns refer to the naming structure and the corresponding sliding law and bed geometry used, and the subsequent columns give details about the number of simulation in each group. The default simulations refer to those run forward in time from the initial conditions. The perturbed parameter column refers to those simulations that perturbed C, φ and M b using Latin hypercube sampling. The end members are the simulations with each parameter altered separately to the extreme bounds of the parameter ranges, while the other parameters are held at the default value.
3. RESULTS
In this section, we report the sea-level contributions produced by the ensembles and how the grounding line retreat varied depending on the ensemble parameters. We also explore the influence of bed topography and sliding law. In some cases, the non-linear power sliding law resulted in very high and unrealistic simulated velocities, which therefore caused slow running speeds due to the short time steps required. We checked the mean square error of each velocity field
$(\vert\vec u\vert - \vert\vec u_{{\rm obs}} \vert)^2 $
, and terminated the ensemble members that had a velocity field with a mean square error >1 std dev. away from the mean of the ensemble. This amounted to a total of 26 of the 298 simulations, which tended to have a combination of low traction (C) and viscosity (φ) coefficients.
3.1. Sea-level contribution
We calculated sea level equivalent mass imbalance from the change in the ice volume above flotation (VAF). For the linear sliding law simulations, the mean rate of sea-level rise over 50 a from the ASE across the ensemble is 0.37 mm a−1 for the BedMap2 (Lr bm) simulations and 0.30 mm a−1 for the modified geometry (Lr′) simulations. For the non-linear sliding law, the mean rate is 0.51 mm a−1 for the Pr bm simulations and 0.38 mm a−1 for the Pr′ simulations. While we do not intend this to be a predictive study, we note that many simulations are close to current trends: recent satellite altimetry observations (Cryosat-2, 2010–2013) indicate that the current rate is 0.33 ± 0.05 mm a−1 (McMillan and others, Reference McMillan2014). Medley and others (Reference Medley2014) used an input-output method to determine the ice flux from the ASE and find the rate of sea-level contribution has increased over the last two decades, from 0.09 ± 0.04 mm a−1 in 1994–96 to 0.27 ± 0.04 mm a−1 in 2010. The default simulations (where P F = 0.5) result in an average rate over the 50 a model period of 0.27 mm a−1 for Lr′ and Lr bm, and 0.26 mm a−1 for Pr′ and Pr bm.
The rate of sea-level rise tends to remains constant over the 50 a period for each group of simulations. The Pr′ and Pr
bm simulations have wider distributions than the Lr′ and Lr
bm simulations (Fig. 2). The distribution of rates for the Lr′ and Pr′ simulations become narrower between the start of the 50 a simulations and the end. However, the results of the linear and non-linear sliding law are not directly comparable, as C was perturbed in the same way between the two laws, resulting in different
$\vec \tau _{\rm b} $
depending on the value of m used. The ensemble member with the most extreme ice loss results in 1.62 mm a−1 sea-level rise, averaged over the 50 a simulation, and belongs to the Pr
bm ensemble group. Approximately 4% of the linear sliding law simulations produced a slight increase in VAF, with the most extreme case resulting in a fall in sea level of 0.03 mm a−1 averaged over 50 a. For the non-linear sliding law simulations, this proportion is larger (17%) and the most extreme case results in 0.19 mm a−1 of sea level fall.
Histograms showing the distribution of the rate of sea-level rise produced by the four geometry and sliding law experimental groups during the first 10 a and last 10 a of the 50 a model period. The dashed lines represent observed rates. The colours represent the different ensemble groups: Lr bm (orange), Lr′ (grey), Pr bm (yellow) and Pr′ (blue).
Overall the range of ice loss rates produced by this ensemble spans a far larger range than the current observed rates. The simulations with the highest sea-level contributions are characterised by low P C and P φ , i.e. simulations with a slippery bed and low viscosity ice. In contrast, the simulations that have the lowest sea-level contributions tend to be stiff ice on a sticky bed.
The experiments with the max[r bm, r′] and min[r bm, r′] beds show that the minimum bed tends to result in a greater contribution to sea-level rise than the maximum bed, because thicker ice at the grounding line (e.g. for the minimum bed, given the same ice surface elevation) results in a higher flux across the grounding line (Schoof, Reference Schoof2007). However, the minimum bed does not tend to result in a higher sea-level contribution than Lr bm and Lr′, despite the greater ice thickness. The maximum and minimum beds are rougher than r bm and r′ due to the lack of spatial correlation between neighbouring cells. This roughness may provide stabilising features, which act to reduce grounding line retreat and ice loss.
3.2. Grounding line retreat
For Pine Island and Thwaites glaciers, grounding line retreat occurs within 50 a in most simulations (Fig. 1). Across the ASE domain, ~70% of the Lr′ and Lr bm simulations and 73% of the Pr′ and Pr bm simulations result in retreat, with the remainder resulting in a slight advance. As shown in Figure 1, the r bm experiments tend to retreat more than the r′ experiments. This difference was confirmed by using a Wilcoxon signed rank test, which shows that the change in grounded area for each ice stream is significantly different depending on which of the two geometries was used. Thwaites experiences the most dramatic retreat of almost 100 km in the most extreme Lr bm simulation. For the Lr′ simulations, the maximum retreat, averaged across the width of Thwaites is ~40 km. For Pine Island, the difference between Lr bm and Lr′ is smaller, with maximum retreat of ~50 km and 40 km, respectively. Grounding line advance occurs in some cases, but due to the model set up whereby C = 0 in the sub-ice shelf areas, this is limited to a width-averaged advance of ~12 km for Pine Island and 8 km for Thwaites. The grounding line positions shown in Figure 1 are not necessarily in their steady-state positions, rather they are just the locations reached after 50 a. For Pine Island and Thwaites, the grounding lines from several ensemble members tend to cluster together at certain locations, indicating regions of stability. However, grounding line stability in these regions may only be transient, where topographic features cause a temporary slowdown of grounding line retreat. For Smith Glacier, due to the modification of the sub-ice shelf bathymetry, grounding line advance does not occur. Patterns in the retreat of Smith are less clear than for Pine Island and Thwaites. The glaciers, including Smith, that make up the group that flow into the Crosson Ice Shelf all experience retreat, up to a maximum of ~60 km.
The relationship between change in grounded area and VAF varies between the three ice stream areas (Fig. 3). For the same retreat in grounded area, Pine Island Glacier loses more VAF compared with the other ice streams. The VAF loss from the ice streams can be partitioned between the loss due to grounding line retreat itself and loss due to dynamic thinning upstream. To distinguish between these two sources of VAF loss, we calculated the VAF in the area between the initial grounding line and the simulated grounding line at 50 a, and divided the result by the total change in VAF, to give the proportion of VAF change that is due to the change in grounded area alone. The proportion tends to be low (<0.1), showing that the majority of ice loss from the ice streams is due to widespread dynamic thinning, rather than localised thinning at the grounding line. Those linear sliding law simulations (Lr bm and Lr′) that exhibit grounding line advance, show no corresponding increase in VAF (within the ice stream areas analysed). This indicates that despite grounding line advance, increased velocities result in the flux of ice across the grounding line being greater than the supply of ice from upstream. Pr bm and Pr′ do include some simulations that result in an increase in VAF within the ice stream areas.
Change in grounded area plotted against change in volume above flotation (VAF), for the three ice streams (according to mask areas delineated in Fig. 1). Each colour represents the four ensemble groups: Lr′ (grey), Lr bm (orange), Pr′ (blue) and Pr bm (yellow). The dashed box delineates the boundaries of the inset.
3.3. Regression analysis
We performed regression analysis in order to quantify the extent that the three parameters sampled using the Latin hypercube influence the simulated retreat (Table 2). We analysed each of the three ice streams separately, with the aim to determine the variation in importance of the various parameters between the ice streams. For each simulation, we calculated the change in grounded area for three areas of interest, corresponding to the main trunk of Pine Island, Thwaites and the collection of ice streams including Smith Glacier, that flow into the Crosson Ice Shelf (mask areas in Fig. 1).
Multiple regression results, showing the intercept (ΔA 0), and the coefficients of the regression equation (Eqn (12)) normalised by the maximum change in grounded area (max ΔA) for each of the three ice streams across all the four groups of the ensemble
The coefficients represent the relationship between the three parameters and the change in grounded area after the 50 a model period. The exponent of the p-value for each coefficient is indicated in parentheses. This gives a measure of the significance of the various coefficients, with more negative exponents having a higher significance, i.e. the more likely the parameter had a true influence on the grounded area response.
Initial exploration of the grounded area response to the parameters revealed that grounded area does not respond linearly to P φ . Simulations with a high stiffening factor (P φ > 0.5) are more likely to advance, due to stiffer ice causing thickening. However, in our model set up, grounding line advance was limited by the lack of knowledge of the basal traction coefficient for newly grounded areas (which we therefore set to 0), and so we regard these members with some scepticism. For the purposes of this study, in which we are more interested in retreat because that is the current observed direction of change, we explored the low viscosity half of the ensemble further, and excluded the high viscosity half from the following regression analysis. For each ice stream area, we created a linear empirical model of grounding line retreat,
$$\eqalign{\Delta A &= \Delta A_0 + \alpha _C \left[ {2^{2(P_C - (1/2))} - 1} \right] + \alpha _\varphi \left[ {2^{2(P_\varphi - (1/2))} - 1} \right] \cr & \quad+ \alpha _{M_{\rm b}} \left[ {2^{2\left( {P_{M_{\rm b}} - (1/2)} \right)}} - 1 \right],}$$
where the change in grounded area (ΔA) is presented as a function of the three parameters (P F , where F = C, φ or M b) and an intercept (ΔA 0). The intercept represents the change in grounded area that occurs when the parameters are unchanged (i.e. if P F = 0.5). The magnitude of the coefficient (α) of each parameter, given in Table 2, indicates the relative importance of that parameter for determining the change in grounded area. Also in Table 2, the R 2 value indicates the proportion of variation in ΔA that can be explained by the regression model. Note that α is normalised by the maximum change in grounded area for each ice stream area in Table 2, in order to allow comparison of α between the ice streams.
For each ice stream, the viscosity parameter (P
φ
) has a significant effect on change in grounded area, with lower viscosity resulting in greater retreat. Out of the three parameters, halving and doubling φ exerts the greatest control on grounded area change for all three ice streams. The sub-ice shelf melt rate parameter (
$P_{M_{\rm b}} $
) has a negative relationship with the change in grounded area, meaning the higher the melt rate the greater the grounding line retreat. However, this relationship has a higher significance for Pine Island and Smith glaciers, than Thwaites (Table 2).
$\alpha _{M_{\rm b}} $
, normalised to account for the size of the ice stream, shows that for each simulation group, the degree to which
$P_{M_{\rm b}} $
influences grounding line retreat is less pronounced in Thwaites Glacier than Pine Island and Smith. The traction parameter has a negative relationship with change in grounded area. However, as shown in Figure 3, most of the simulations that experienced grounding line advance still experience VAF loss (within the polygon areas). The differences in grounding line retreat between the two geometries, which are evident in Figure 1, are also shown in Table 2. For example, for Thwaites the coefficients α tend to be greater for r
bm than r′, indicating that, given the same perturbations, r
bm is more sensitive to changes in the parameters.
4. DISCUSSION
The ensemble experiments show that all the three ice streams are particularly sensitive to changes in the coefficient associated with ice viscosity (φ). Changes in φ are responsible for the greatest change in grounded area for all the three ice streams. The simulations with low viscosity also lose the most VAF across the ASE domain, caused by thinning across the main trunks of the ice streams. In most cases, the basal traction coefficient also has a significant effect on the grounded area, with lower traction resulting in grounding line advance, but also VAF loss. Slippery beds, due to low C, cause increased velocities across the domain. This results in more ice being delivered to the grounding line, which can cause grounding line advance in some cases. However, due to the greater velocities, the flux of ice across the grounding line also increases, and therefore, more VAF is lost. We would expect this imbalance to eventually lead to grounding line retreat due to dynamic thinning (Pattyn and others, Reference Pattyn2013), however this does not occur during the 50 a time scale used in this study. In contrast, lowering φ causes increased velocities predominantly in the downstream areas of the ice streams, near the grounding line. This causes more localised thinning and grounding line retreat. The sub-ice shelf melt rate influences the grounding line retreat of Pine Island and Smith glaciers, however it does not have such a significant effect on the grounded area of Thwaites.
We should note, this is a model sensitivity study and the importance of the three parameters cannot be compared directly within the same ice stream because the range in which we varied the parameters (between a halving and a doubling) represents something different for each parameter. For example, doubling M
b is plausible, as melt rates are driven by an oceanographic forcing, and changes to this forcing have been observed in the Pine Island Glacier ice shelf cavity (Jacobs and others, Reference Jacobs, Jenkins, Giulivi and Dutrieux2011). In contrast, C and φ are unlikely to vary in the same way. Indeed, the P
φ
= 0 and P
C
= 0 end members result in maximum velocities of >15 km a−1 for the non-linear sliding law, compared with velocities of <10 km a−1 for the
$P_{M_{\rm b}} = 1$
end member. However, our results do show that the three ice streams have varying responses to the same parameters. These differences in the parameter sensitivities show that the ice streams' dynamics are controlled by different factors, most likely because of variation in their bed and ice-shelf geometries. We will now examine each of the ice streams in turn.
4.1. Pine Island Glacier
Pine Island Glacier is situated in a narrow bedrock channel with a retrograde slope. The ice shelf is laterally constrained in a relatively narrow embayment (~50 km wide). Its flow is therefore slower than that of an equivalent unconstrained ice shelf, due to the effect of lateral drag (Dupont and Alley, Reference Dupont and Alley2005), although this is limited by weak shear margins. Our regression results show that the sub-ice shelf melt rate is a significant control on the change in grounded area, due to ice shelf buttressing. High melt rates cause the ice shelf to thin and therefore reduce the buttressing of the grounded ice, which in turn causes the ice streams to speed up. Ice shelf thinning in the ASE has been observed using satellite altimetry (Pritchard and others, Reference Pritchard2012). This coincides with the observed acceleration and grounding line retreat. Seroussi and others (Reference Seroussi2014) test the sensitivity of Pine Island to various climate forcings and they find that the sub-ice shelf melt rate has the strongest influence on the ice dynamics.
The pattern of grounding line retreat along the main trunk of Pine Island Glacier indicates there are regions of relative stability (Fig. 1), provided by topographic rises where the bed is locally prograde (sloping downwards in the direction of flow), within the overall retrograde bed. These regions of slower retreat vary between the two bed geometries. For the modified geometry r′, Figure 4 shows there is a cluster of grounding line positions ~15 km up-glacier from the initial grounding line position. The retreating grounding lines experience a slow down of retreat in this area, resulting in a cluster of grounding lines at 50 a. This coincides with a topographic high, relative to r bm, which has a flat bed at this location. The simulations with grounding lines in this cluster are not in steady state, rather the retreat only slows temporarily over topographic rises, so that after 50 a the grounding lines of each simulation are more likely to be found on these rises. If the other parameters are sufficiently favourable to retreat, the grounding line is able to retreat beyond these features. The results of Joughin and others (Reference Joughin, Smith and Holland2010) also indicate that Pine Island's grounding line may jump from one transient position to another, as it passes over regions of retrograde slope to the next bedrock high.
Pine Island Glacier grounding line retreat over time. In the top two plots, each curve represents the movement in grounding line position over the 50 a period, relative to the initial grounding line position (0 km). Retreat from the initial grounding line is represented by positions <0 km and advance by positions >0 km. The dashed lines represent linear retreat rates. The top plot shows the linear sliding law results and the middle plot shows the non-linear sliding law results. The bottom plot shows the geometry of these cross sections with the black vertical line showing the initial grounding line, and the coloured vertical lines showing the grounding line positions after 50 a for the default simulations of each ensemble group: Lr′ (grey), Lr bm (orange), Pr′ (blue) and Pr bm (yellow).
In contrast, r bm has a topographic rise near the initial grounding line, relative to r′, and so the grounding lines tend to stay in this region for a longer period than the r′ simulations. Sun and others (Reference Sun, Cornford, Liu and Moore2014) find that errors on the order of 10 m in the height of this ridge can alter the timing of the onset of retreat. However, the r bm simulations that are more favourable for retreat, tend to rapidly retreat ~20 km in <10 a, as a result of the smooth retrograde slope down from the initial grounding line, before a locally prograde slope ~25 km upstream.
The narrow channel of Pine Island Glacier is the principal gate from which ice is lost from this catchment. Upstream acceleration results in more ice being delivered to the grounding line, which in turn limits the grounding line retreat, in a negative feedback. This can be seen in Figure 3, where, despite high VAF loss in some simulations, loss of grounded area is more limited, compared with Thwaites and Smith glaciers.
4.2. Thwaites Glacier
Thwaites Glacier is situated in a large, wide bedrock basin. Its outlet is divided into the slow-flowing eastern ice shelf and the fast-flowing western ice tongue, and a similar division is evident upstream from the grounding line (Fig. 5). In comparison with Pine Island Glacier, the Thwaites ice shelves are unconstrained by slower flowing margins or fjord walls, and therefore it is believed that they provide less buttressing to the grounded ice (MacGregor and others, Reference MacGregor, Catania, Markowski and Andrews2012). Our results (Table 2) support this as they show that grounding line retreat of Thwaites is less sensitive to changes in the sub-ice shelf melt rate compared with Pine Island and Smith glaciers. Parizek and others (Reference Parizek2013) also find that the Thwaites Glacier tongue provides limited buttressing to the grounded ice.
The bed topography under Thwaites Glacier for (a) BedMap2 (r bm) and (b) the modified bed (r′). The grey dashed contours represent initial velocity of the default simulations and the solid black contour show the initial grounding line. For presentational purposes we have masked out the non-ice shelf area.
Tinto and Bell (Reference Tinto and Bell2011) created a new bathymetric model, based on the NASA IceBridge data, which shows an offshore ridge, with two peaks (eastern and western). These peaks are found in both r bm and r′ (Fig. 5). However, in r′, the eastern peak reaches a higher elevation (just 60 m below the surface), than in r bm, which is ~300 m below the surface. The eastern ice shelf of Thwaites is currently pinned by the eastern peak, while the Thwaites Glacier tongue is lightly grounded, perhaps only ephemerally, on the western peak (Schmeltz and others, Reference Schmeltz, Rignot and MacAyeal2001; Tinto and Bell, Reference Tinto and Bell2011). So while the ice shelves of Thwaites are unconstrained laterally, these pinning points may contribute to the buttressing effect acting on the grounded ice. In particular this appears to be the case for the slow flowing eastern ice shelf. The eastern pinning point does not completely unpin in any of the Lr′ simulations, whereas over 80% of the Lr bm simulations have unpinned by the end of the 50 a.
The unpinning of the eastern peak correlates with the simulations with greatest VAF loss (Fig. 6). However, it is unclear whether the VAF loss is the result of the reduction in basal traction from unpinning, or whether both unpinning and high VAF loss are due to changes in the grounded ice caused by the ensemble parameters. To explore this further, we performed a series of additional experiments, in which we removed the eastern peak and ran four simulations: one in the default configuration and three additional members of the linear sliding law ensemble that did not experience high VAF loss. We ran the simulations with both bed geometries for 50 a. The experiments, represented by the crosses in Figure 6, show that removing the pinning point has little effect on the VAF loss from Thwaites Glacier. This suggests that the ensemble members where unpinning occurs have a high VAF loss due to a reduction in
$\vec \tau _{\rm b} $
over the whole area (i.e. from modifying C), rather than it being a direct consequence of the reduction in
$\vec \tau _{\rm b} $
integrated over the pinning point alone.
The time averaged basal traction,
$\vec \tau _{\rm b} $
integrated over the area of Thwaites Glacier's eastern peak, plotted against the total change in VAF over the 50 a period for each ensemble member (circles). The crosses represent the experiments where the peak was removed, with the dashed lines connecting them to the original ensemble members. Grey symbols denote the Lr′ simulations and orange, the Lr
bm simulations.
Bed geometry has a significant influence on the change in grounded area. The r bm simulations tend to have greater loss of grounded area than the r′. In r′ there is a deep trough upstream of the grounding line in the centre of the fast flowing western portion of the ice stream (Fig. 5). A similar feature was identified by Tinto and Bell (Reference Tinto and Bell2011). As well as this trough, there are a series of rises and smaller troughs to the west along the grounding line. These smaller scale features, which BedMap2 would not resolve if present, result in greater grounding line stability compared with the smoothed features of r bm. This supports the work of Durand and others (Reference Durand, Gagliardini, Favier, Zwinger and le Meur2011), who investigate the relationship between modelled grounding line retreat and roughness of the bed topography. They find that grounding line retreat is slower when the bed topography is rougher (i.e. contains more small scale features). This is because small-scale topographical features contain stabilising prograde portions, on which the grounding line retreat slows down. Sun and others (Reference Sun, Cornford, Liu and Moore2014) demonstrate that the pattern of grounding line retreat is more sensitive to low frequency noise applied to the bed geometry, compared with high frequency noise. Our results also suggest that the amount of retreat is more variable over a smoother bed (i.e. r bm) compared with one with higher frequency features (i.e. r′). Durand and others (Reference Durand, Gagliardini, Favier, Zwinger and le Meur2011) recommend that bed elevation models have a resolution of the order of 1 km. BedMap2 is derived from radar surveys with a between-track distance of ~5 km in the downstream areas of the ASE ice streams and ~15 km further inland (Fretwell and others, Reference Fretwell2013), so bed roughness is likely to be under represented; whereas the procedure used to derive r′ is able to generate features at scales <5 km. These results suggest that subglacial topographic features, such as troughs and rises, could control grounding line position and the dynamics of Thwaites Glacier if they exist, rather than the ice shelf pinning points.
The distribution of sea-level contribution from each ensemble group (Fig. 2) provides further evidence that the stabilising features in r′ play a role in the dynamics of the ASE ice streams. Over the 50 a model period, the distribution of sea level contribution rates becomes narrower for the r′ ensembles, compared with r bm, suggesting that the simulations become more consistent with one another. This could be because the grounding lines in these simulations tend to cluster on topographic rises, which are not present in r bm.
4.3. Smith Glacier
Smith Glacier is located in a narrow, deep bedrock channel. Along with Pope and Kohler glaciers, Smith Glacier is buttressed by the Crosson Ice Shelf, which is confined by the Bear Peninsula to the northwest and Mount Murphy to the south east (MacGregor and others, Reference MacGregor, Catania, Markowski and Andrews2012). Similar to Pine Island Glacier, the grounding line retreat was limited to the narrow channels (Fig. 6), due to these topographic constraints. The basal melt rate experienced by the Crosson Ice Shelf is highly influential on the change in grounded area (Table 2). This indicates that, similarly to Pine Island, the Crosson Ice Shelf buttresses the upstream ice and ice shelf thinning is therefore linked to grounding line retreat.
5. CONCLUSIONS
We employed an ice-sheet model to investigate the sensitivity of the ASE ice streams to changes in key physical parameters. We used Latin hypercube sampling to generate an ensemble of model simulations with varying basal traction, viscosity and sub-ice shelf melt rate. We tested the influence of bed geometry, by replicating the ensemble over two different bedrock topographies. We also replicated the ensemble, using a linear and non-linear (m = 1/3) viscous sliding law. Given the linear sliding law, the average rate of contribution to sea level from the ASE over the 50 a period is 0.37 mm a−1 for the simulations using BedMap2 geometry (Lr bm) and 0.30 mm a−1 using a geometry modified to better conserve mass (Lr′). For the non-linear sliding law, the mean rate is 0.51 mm a−1 for the Pr bm simulations and 0.38 mm a−1 for the Pr′ simulations. The recently observed rate is 0.33 ± 0.05 mm a−1 (McMillan and others, Reference McMillan2014).
The simulated grounding line retreat of the ice streams in the ASE is highly dependent on the viscosity and traction. Ice with a low viscosity results in increased velocity and thinning, especially near the grounding line, resulting both in grounding line retreat and VAF loss. On the other hand, reducing the basal traction has the counter-intuitive effect of both increasing the grounded area and decreasing the VAF. The position of the grounding line advances with lower traction, because more ice is delivered to the grounding line. However, the flux across the grounding line also increases for the more slippery simulations, resulting in increased VAF loss. This imbalance suggests that on time scales >50 a the grounding line would eventually retreat due to the continued mass loss and related dynamic thinning.
Our work confirms that there is a requirement for bed geometry to include km-scale resolution features (Durand and others, Reference Durand, Gagliardini, Favier, Zwinger and le Meur2011). The grounding line dynamics are sensitive to these features, as shown by the Pine Island and Thwaites glacier results, because even small topographic rises can result in a localised slow down in grounding line retreat. Sub-ice shelf topography is also important for modelling grounding line dynamics, unless retreat can be guaranteed. This is especially the case under the Crosson and Thwaites ice shelves, which have not benefited from the more detailed observations made for the Pine Island Glacier ice shelf cavity (Jenkins and others, Reference Jenkins2010).
The ice shelves of Pine Island and Smith glaciers are laterally constrained in bedrock channels, buttressing the grounded ice. Reducing this buttressing through increased sub-ice shelf melt rate, and hence ice shelf thinning, results in increased grounding line retreat and ice loss. However, our results suggest that the pinning points on the Thwaites Glacier ice shelf do not provide substantial buttressing to the grounded ice because imposed unpinning of the ice shelf does not affect the amount of VAF loss from Thwaites.
In summary, we have seen that while our models of Pine Island, Smith and Thwaites glaciers respond similarly to changes in basal traction and englacial viscosity, Thwaites Glacier appears to be somewhat less sensitive to likely near-future ice shelf thinning than its neighbours, although its larger size means that its contribution to future sea-level rise may well dominate. Ice dynamics in both Pine Island and Thwaites glaciers are also strongly affected by perturbations to the bedrock, which if not directly observed, are suggested by the observed velocity. Finally, although our parameter ranges are large – halving to doubling the initial basal traction, effective viscosity and melt rate – the maximum rate of sea-level rise is no >~1.62 mm a−1 (about five times the current rate). In future work, we will use observations of present day thinning rates to further constrain the results.
ACKNOWLEDGEMENTS
This work was supported by funding from the UK Natural Environment Research Council (NERC) through the iSTAR-C programme (grant number NE/J005738/1). BISICLES simulations were carried out on the University of Bristol's Blue Crystal Phase 3 supercomputer. BISICLES development is led by D. F. Martin at Lawrence Berkeley National Laboratory, California, USA, and S. L. Cornford at the University of Bristol, UK, with financial support provided by the US Department of Energy and the UK Natural Environment Research Council, through the NERC Centre for Polar Observation and Modelling (CPOM). We thank the reviewers and editor for their helpful and constructive comments.
Open access
