Drivers of Pine Island Glacier speed-up between 1996 and 2016

Pine Island Glacier in West Antarctica is among the fastest changing glaciers worldwide. Over the last 2 decades, the glacier has lost in excess of a trillion tons of ice, or the equivalent of 3 mm of sea level rise. The ongoing changes are thought to have been triggered by ocean-induced thinning of its floating ice shelf, grounding line retreat, and the associated reduction in buttressing forces. However, other drivers of change, such as large-scale calving and changes in ice rheology and basal slipperiness, could play a vital, yet unquantified, role in controlling the ongoing and future evolution of the glacier. In addition, recent studies have shown that mechanical properties of the bed are key to explaining the observed speed-up. Here we used a combination of the latest remote sensing datasets between 1996 and 2016, data assimilation tools, and numerical perturbation experiments to quantify the relative importance of all processes in driving the recent changes in Pine Island Glacier dynamics. We show that (1) calving and ice shelf thinning have caused a comparable reduction in ice shelf buttressing over the past 2 decades; that (2) simulated changes in ice flow over a viscously deforming bed are only compatible with observations if large and widespread changes in ice viscosity and/or basal slipperiness are taken into account; and that (3) a spatially varying, predominantly plastic bed rheology can closely reproduce observed changes in flow without marked variations in ice-internal and basal properties. Our results demonstrate that, in addition to its evolving ice thickness, calving processes and a heterogeneous bed rheology play a key role in the contemporary evolution of Pine Island Glacier.

1 Introduction and motivation

Since the 1990s, satellite measurements have comprehensively documented the sustained acceleration in ice discharge across the grounding line of Pine Island Glacier (PIG, Fig. 1) in West Antarctica (Rignot et al., 2002; Rignot, 2008; Rignot et al., 2011; Mouginot et al., 2014; Gardner et al., 2018; Mouginot et al., 2019b). The changes in flow speed are an observable manifestation of the glacier’s dynamic response to both measurable perturbations, such as calving and ice shelf thinning, and poorly constrained variations in physical ice properties and basal sliding. Evidence from indirect observations has indicated that changes in ice shelf thickness have occurred since at least some decades before the 1970s (Jenkins et al., 2010; Smith et al., 2017; Shepherd et al., 2004; Pritchard et al., 2012). Within the last 2 decades, thinning of the grounded ice (Shepherd et al., 2001; Pritchard et al., 2009; Bamber and Dawson, 2020), intermittent retreat of the grounding line (Rignot et al., 2014), changes in calving front position (Arndt et al., 2018), and the partial loss of ice shelf integrity (Alley et al., 2019) have all been reported in considerable detail. At the same time, numerical simulations of ice flow have confirmed the strong link between ice shelf thinning, which reduces the buttressing forces, and the increased discharge across the grounding line (Schmeltz et al., 2002; Payne et al., 2004; Joughin et al., 2010; Seroussi et al., 2014; Favier et al., 2014; Arthern and Williams, 2017; Reese et al., 2018; Gudmundsson et al., 2019). Due to the dynamic connection between ocean-driven ice shelf melt rates and tropical climate variability (Steig et al., 2012; Dutrieux et al., 2014; Jenkins et al., 2016; Paolo et al., 2018), several model studies have focused on the important problem of simulating the response of PIG to a potential anthropogenic intensification of melt. Such external perturbations, in combination with ice-internal feedbacks including the marine ice sheet instability, can force PIG along an unstable and potentially irreversible trajectory of mass loss (Favier et al., 2014; Rosier et al., 2020). Whereas significant progress has been made in simulating the melt-driven retreat of PIG, less attention has been given to other processes that could affect the force balance and thereby inhibit or foster changes in ice dynamics. Increased damage in the shear margins of the ice shelf, for example, has been reported by Alley et al. (2019) and Lhermitte et al. (2020), and is known to reduce the buttressing capacity of an ice shelf (Sun et al., 2017). Moreover, a series of recent calving events has led to a sizeable reduction in the extent of the ice shelf and caused a potential loss of contact with pinning points along the eastern shear margin (Arndt et al., 2018).

The relative impact of changes in ice geometry, basal shear stress, and/or ice rheology on the dynamics of PIG has previously been emphasized in numerical studies by, for example, Schmeltz et al. (2002), Payne et al. (2004), Gillet-Chaulet et al. (2016), Joughin et al. (2019), and Brondex et al. (2019). In all cases, some combination of thickness changes, ice softening, a reduction in ice shelf buttressing, and variations in basal shear stress was required to attain an increase in flow speed comparable to observations. Similar conclusions were reached for other Antarctic glaciers. Based on a comprehensive series of model perturbation experiments, Vieli et al. (2007) suggested that the acceleration of the Larsen B Ice Shelf prior to its collapse in 2002 could not solely be explained by the retreat of the ice shelf front or ice shelf thinning, but required a further significant weakening of the shear margins. Complementary conclusions were reached by Khazendar et al. (2007), who demonstrated the important interdependence of the calving front geometry, a variable ice rheology, and flow acceleration based on data assimilation and model experiments for the Larsen B Ice Shelf.

In order to comprehensively diagnose the importance of all processes that have contributed to the acceleration of PIG over the period 1996 to 2016, this study brings together the latest observations and modelling techniques. We consider how calving, ice shelf thinning, the induced dynamic thinning upstream of the grounding line, and potential changes in ice-internal and basal properties have caused a different dynamic response across the ice shelf, the glacier's main trunk, the margins, and the tributaries. Initial observations indicated that the speed-up of PIG was primarily confined to its fast-flowing central trunk (Rignot et al., 2002; Rignot, 2008), though more complex spatio-temporal patterns of change have emerged more recently (Bamber and Dawson, 2020). The rapid and spatially diverse acceleration of the flow is an expression of the glacier's dynamic response to changes in the force balance, and it is imperative that numerical ice flow models are capable of reproducing this complex behaviour in response to the correct forcing. In general, the driving stress (τ_d), which depends on the ice thickness distribution, is balanced by resistive stresses, which include the basal drag (τ_b), side drag through horizontal shear (τ_W), longitudinal resistive forces (τ_L), and back forces by the ice shelf (τ_IS):

$$\begin{array}{}\text{(1)}& {\mathit{\tau }}_{\mathrm{d}}={\mathit{\tau }}_{\mathrm{b}}+{\mathit{\tau }}_{\mathrm{W}}+{\mathit{\tau }}_{\mathrm{L}}+{\mathit{\tau }}_{\mathrm{IS}}.\end{array}$$

It is conceivable that each of the terms in Eq. (1) has changed considerably in recent decades in response to changes in calving front position, ice thickness, ice properties, and/or basal slipperiness. The interplay between different changing forces, in combination with the appropriate boundary conditions, underlies the observed speed-up of PIG (ΔU). Present-day observations of ΔU are generally assumed to be dominated by ice shelf thinning and induced dynamic loss and redistribution of mass upstream of the grounding line, which includes grounding line retreat and the associated loss of basal traction. Other possible contributions to ΔU, such as ice front retreat or changes in ice viscosity (including damage) or basal slipperiness, remain unquantified and are not generally included in model simulations of future ice flow at decadal to centennial timescales. These missing processes, if important, could lead to a systematic bias in model projections of future ice loss, or could prompt the use of unrealistically large perturbations in, for example, τ_IS in an attempt to reproduce observed values of ΔU.

In this study we used a regional configuration of the shallow ice stream (SSA) flow model Úa (Gudmundsson, 2020) for PIG to diagnose how individual processes (calving, ice thinning and associated grounding line movement, and changes in ice viscosity and basal slipperiness) have contributed to ΔU over the period 1996 to 2016. The diagnostic model response to prescribed changes in ice geometry was analysed, based on the latest observations of calving and ice thinning rates between 1996 and 2016. For each perturbation, changes in the stress balance (Eq. 1) and the associated ice flow response were computed. Any further discrepancies between modelled and observed changes in velocity were attributed to variations in ice properties, commonly parameterized by a rate factor A, and changes in basal slipperiness C. Results enabled us to validate the ability of current-generation ice flow models to reproduce the complex response of PIG to a range of realistic forcings, and to verify whether common model assumptions such as a static calving front and fixed ice viscosity and basal slipperiness are indeed justified.

Although the aforementioned method provides insights into the individual contribution of geometrical perturbations and changes in ice viscosity and basal slipperiness to overall changes in ice flow, results will likely depend on a number of structural assumptions within the ice flow model. Previous studies have shown that different forms of the sliding law, for example, can produce a distinctly different simulated response of PIG to changes in ice thickness (Joughin et al., 2010; Gillet-Chaulet et al., 2016; Joughin et al., 2019; Brondex et al., 2019). Joughin et al. (2010, 2019) showed that a regularized Coulomb law or the plastic limit of a non-linear viscous power law provides a better fit between modelled and observed changes in surface velocity along the central flow line of PIG compared to a commonly used Weertman law with a cubic dependency of the sliding velocity on the basal shear stress. Motivated by the above considerations, we explore new ways to derive spatially variable constraints on the form of the sliding law and thereby provide the first comprehensive, spatially distributed map of basal rheology beneath PIG.

The remainder of this paper is organized as follows. In Sect. 2.1 we introduce the observational datasets used to constrain and validate the ice flow model. Additional details about our data processing methods are provided in Appendix A. Section 2.2 outlines the experimental design and provides a summary of the main model components. Further technical details about the model set-up and a discussion about the sensitivity of our results to numerical model details are provided in Appendix B and Appendix D respectively. Results and an accompanying discussion of all experiments are provided in Sect. 3.1–3.3. Final conclusions are formulated in Sect. 4.

2 Data and methods

The first aim of this study is to simulate the dynamic response of PIG to a series of well-defined geometric perturbations over the period 1996 to 2016 and compare model output to observed changes in surface speed over the same time period. As detailed in Sect. 1, geometric perturbations are considered to be observed changes in the calving front position and observed changes in ice thickness of the ice shelf and grounded ice. We are primarily interested in the relative contribution of each perturbation to the observed speed-up of PIG between 1996 and 2016. Each contribution can be characterized by a relative change in velocity, $\left(U_{\mathrm{pert}}-U_{\mathrm{96}}\right)/\left(U_{\mathrm{16}}-U_{\mathrm{96}}\right)$, where U₉₆ and U₁₆ were obtained from model optimization experiments, as described in Sect. 2.2.1, and velocities of the perturbed states (U_pert) were obtained from a series of diagnostic model calculations, as described in Sect. 2.2.2. First, the data sources required for each of these experiments are listed.

2.1 Observed changes of Pine Island Glacier between 1996 and 2016

Figure 1 Pine Island Glacier (PIG) and its location in West Antarctica. (a) Surface speed of PIG in 1996 in m yr⁻¹, as reported by the MEaSUREs programme (Mouginot et al., 2019a). Solid black outlines delineate the extent of the PIG catchment (Rignot et al., 2011) and 1996 grounding line position (Rignot et al., 2014). The white line along the central flow line indicates the location of the transect in Fig. 2. The dashed rectangle corresponds to the extent of panels (b) and (c). (b) Observed increase in surface speed (Mouginot et al. (2019a), colours and contours in m yr⁻¹ and loss of ice shelf extent (grey shaded area) between 1996 and 2016. The blue line indicates the 2011 grounding line (Rignot et al., 2014). (c) Total change in ice thickness between 1996 and 2016 (ΔH in m), based on a combination of CPOM data (Shepherd et al., 2016) for the grounded ice and newly analysed data for the ice shelf (Appendix A). The zero contour is shown in black; other contours in grey are spaced at 20 m intervals.

Our study area and model domain encompasses the 135 000 km² grounded catchment (Rignot et al., 2011) and seaward-floating extension of PIG in West Antarctica, as depicted in Fig. 1a. To investigate the physical processes that forced the contemporary speed-up of the glacier and its increase in grounding line flux between the years 1996 and 2016, we needed detailed observations of the surface velocity, ice thickness, and calving front position for both years.

The surface velocity measurements used in this study were taken from the MEaSUREs database (Mouginot et al., 2019a, b). For 1996, synthetic-aperture radar data from the ERS-1 and ERS-2 missions were processed using interferometry techniques and combined into a mosaic with an effective timestamp of 1 January 1996. The MEaSUREs velocities for 2016 were based on feature tracking of Landsat 8 imagery with an effective timestamp of 1 January 2016. The change in surface speed between the two years is shown in Fig. 1b, and we refer to, for example, Rignot et al. (2014) and Gardner et al. (2018) for a more comprehensive description of these observations.

To obtain an accurate estimate of the ice thickness distribution of PIG in 1996, we compiled a time series of surface height changes from a comprehensive set of overlapping satellite altimeter data between 1996 and 2016. The integrated altimeter trend over the 20-year time interval, shown in Fig. 1c, was subtracted from the recent BedMachine Antarctica reference thickness (Morlighem et al., 2020). As satellite altimeters are very precise instruments capable of detecting small changes in surface elevation, our approach is more robust compared to thickness estimates based on a snapshot digital elevation model, such as the ERS-1-derived product (Bamber and Bindschadler, 1997), which has poor vertical accuracy. The BedMachine Antarctica reference thickness is based on the high-resolution Reference Elevation Model of Antarctica (REMA; Howat et al., 2019), tied to CryoSat-2 data, and an improved estimate of bedrock topography using mass conservation methods. The nominal date of the BedMachine geometry corresponds to the date stamp of the REMA elevation model, which is spatially variable but largely between 2014 and 2018 for PIG. We denote the 1996 and BedMachine Antarctica ice thickness estimates as H₉₆ and H₁₆ respectively, where $H_{\mathrm{96}}=H_{\mathrm{16}}-\mathrm{\Delta }H$, with ΔH the integrated altimeter trend. Estimates of ΔH were based on a combination of existing Centre for Polar Observation and Modelling (CPOM) measurements of surface elevation changes for areas upstream of the 2016 grounding line (Shepherd et al., 2016) and newly analysed data for the floating ice shelf. A detailed description of our methods and a map of the data coverage for ΔH can be found in Appendix A and Fig. A1 respectively.

The grounding line location for H₁₆ (blue line in Fig. 1b–c) corresponds to the DInSAR-derived grounding line in 2011 from Rignot et al. (2014), since this is included as a constraint in the generation of the BedMachine Antarctica bed topography (Morlighem et al., 2020). The grounding line for $H_{\mathrm{96}}=H_{\mathrm{16}}-\mathrm{\Delta }H$ approximately follows the 1992–1996 DInSAR estimates (Rignot et al., 2014), as shown in Fig. A1. To further improve the agreement between the model and DInSAR grounding line in 1996, some localized adjustments less than 150 m were made to the bed topography. The final grounding line location for H₉₆ is depicted in Fig. 1a–c.

Alongside the above-listed observed changes in flow dynamics and ice thickness, the calving front of PIG retreated by up to 30 km between 1996 and 2016 during a succession of large-scale calving events; see, for example, Arndt et al. (2018). We traced the calving front positions in 1996 and 2016 from cloud-free Landsat 5 and Landsat 8 panchromatic band images with timestamps of 18 February 1997 and 25 December 2016 respectively. Both outlines are included in Fig. 1b–c, and the ice shelf area that was lost between 1996 and 2016 is shaded in grey.

2.2 Experimental design

2.2.1 Optimization experiments

To obtain an optimal model configuration for the state of PIG in 1996, we explicitly solved the stress balance by assimilating the estimated ice thickness (H₉₆), measured calving front position, and measured surface velocities in the shallow ice stream model Úa (Gudmundsson et al., 2012; Gudmundsson, 2020). An analogous routine was applied for 2016. This “data assimilation” or “optimization” step is commonly adopted in glaciology (see MacAyeal, 1992, for one of the earliest examples) to minimize the misfit between modelled and observed surface velocities through the optimization of uncertain physical parameters. The optimization capabilities of Úa (further details are provided in Appendix B) were used to optimize the uncertain spatial distribution of the rate factor, A, and basal slipperiness, C. These physical parameters define the constitutive model and the relationship between basal shear stress τ_b and basal sliding velocity U_b respectively:

$$\begin{array}{}\text{(2)}& {\displaystyle }\dot{\mathit{ϵ}}& {\displaystyle }=A{\mathit{\tau }}_{\mathrm{E}}^{n-\mathrm{1}}\mathit{\tau },\text{(3)}& {\displaystyle }{\mathit{\tau }}_{\mathrm{b}}& {\displaystyle }=C^{-\mathrm{1}/m}{∥{\mathit{U}}_{\mathrm{b}}∥}^{\frac{\mathrm{1}}{m}-\mathrm{1}}{\mathit{U}}_{\mathrm{b}}.\end{array}$$

Glen's law, Eq. (2), relates the strain rates $\dot{\mathit{ϵ}}$ to the deviatoric stress tensor τ. A creep exponent n=3 was used throughout this study. Equation (3) is known as a Weertman sliding law (Weertman, 1957) and describes a linear viscous, non-linear viscous, or close-to-plastic bed rheology for m=1, m>1, and m≫1 respectively. Throughout this study, a range of values for m are considered, as specified below. For each m we performed a new inversion for A and C; example results for m=3 are provided in Appendix B. The outcome of the optimization step is an estimate for A and C that best fits the stress balance in Eq. (1) for given observations of geometry and surface velocity, associated measurement errors, and assumptions about the prior values of A and C. Solutions for A and C are not generally unique but rather depend on the choice of optimization scheme and several poorly constrained optimization parameters. Further details about the optimization scheme used and a discussion about the robustness of our results with respect to uncertain optimization parameters are provided in Appendix B and D.

2.2.2 Geometric perturbation experiments

The optimal model configuration in 1996 was subsequently used as the reference state for a series of numerical perturbation experiments, aimed at simulating the impact of observed changes in geometry on the flow of PIG. For each perturbation, the modified force balance (Eq. 1) and corresponding perturbed velocities were diagnosed within Úa. The rate factor and basal slipperiness were kept fixed to their 1996 values, although the basal traction was reduced to zero and slipperiness values became irrelevant in areas that ungrounded due to ice thinning. Experiments will be referred to as ${\mathcal{E}}_{*}^m$, with a variable subscript to indicate the type of perturbation and a superscript to specify the value of the sliding exponent m.

Figure 2 Overview of changes along the Pine Island Glacier centreline from (a) year 1996 to (b) year 2016. Increased ice flow is driven by a combination of calving, ice shelf thinning, and dynamic thinning with movement of the grounding line, as well as changes in basal sliding and ice rheology. Transects of the geometry are based on observations along the flow line indicated in Fig. 1; black dots indicate the respective grounding line positions in both years. Crevasses are introduced for illustration purposes only and do not strictly correspond to observed features. The importance of each “driver of change” was investigated in a series of numerical perturbation experiments, denoted by ${\mathcal{E}}_{*}^m$ in panel (b), with m indicating the sliding exponent and * the respective experiment described in Sect. 2.2.

• ${\mathcal{E}}_{\mathrm{Calv}}^m$. Changes in the calving front location were prescribed to reflect the loss of ice shelf between 1996 and 2016 (see Fig. 1b–c). All model grid elements downstream of the 2016 calving front (grey shaded area in Fig. 1b) were deactivated, whilst elements upstream of the 2016 calving front remained fixed to avoid numerical interpolation errors. All other model variables were kept fixed. The difference between the 1996 surface velocity and the perturbed velocity will be denoted by $\mathrm{\Delta }{U}_{\mathrm{Calv}}=U_{\mathrm{Calv}}-U_{\mathrm{96}}$.

• ${\mathcal{E}}_{\mathrm{ISThin}}^m$. Changes in ice shelf thickness were prescribed, corresponding to observed thinning of the ice shelf between 1996 and 2016 (Fig. 1c). Note that the calving front and grounding line location did not change in this experiment, which is similar to previous studies by, for example, Reese et al. (2018) and Gudmundsson et al. (2019). The instantaneous change in surface velocity due to ice shelf thinning will be denoted by $\mathrm{\Delta }{U}_{\mathrm{ISThin}}=U_{\mathrm{ISThin}}-U_{\mathrm{96}}$.

• ${\mathcal{E}}_{\mathrm{Thin}}^m$. Observed changes in both the floating and grounded parts of PIG were prescribed. This caused the grounding line to move from its 1996 position (black line in Fig. 1b–c) to the 2016 position (blue line in Fig. 1b–c). Velocity changes caused by thinning of the floating and grounded ice will be denoted by $\mathrm{\Delta }{U}_{\mathrm{Thin}}=U_{\mathrm{Thin}}-U_{\mathrm{96}}$.

• ${\mathcal{E}}_{\mathrm{CalvThin}}^m$. Combined changes in calving front position (as in ${\mathcal{E}}_{\mathrm{Calv}}^m$) and thinning (as in ${\mathcal{E}}_{\mathrm{Thin}}^m$) were prescribed. Corresponding velocity changes will be denoted by $\mathrm{\Delta }{U}_{\mathrm{CalvThin}}=U_{\mathrm{CalvThin}}-U_{\mathrm{96}}$.

A schematic overview of the experiments is provided in Fig. 2. While ${\mathcal{E}}_{\mathrm{Calv}}^m$ allows us to assess the instantaneous impact of calving between 1996 and 2016, the experiments ${\mathcal{E}}_{\mathrm{ISThin}}^m$ simulate the instantaneous response to total changes in ice shelf thickness between 1996 and 2016. The separate perturbations make it possible to disentangle changes in ice shelf buttressing caused by each process, and hence their relative importance for driving the transient evolution of the flow. However, both experiments ignore the time-dependent, dynamic response of the upstream grounded ice and the associated loss of basal traction due to grounding line movement. Dynamic thinning of grounded ice, as well as migration of the grounding line, is included in the experiments ${\mathcal{E}}_{\mathrm{Thin}}^m$, which allows us to determine the full dynamic response to changes in ice thickness. Finally, the experiments ${\mathcal{E}}_{\mathrm{CalvThin}}^m$ combine both calving and ice thinning, and thereby accounts for all geometric perturbations.

2.2.3 Estimates of changes in A and C

Later on we show that geometric perturbations alone are not able to fully reproduce the observed patterns of speed-up across the PIG catchment. It is conceivable that, along with the evolving geometry, variations in ice and basal properties have contributed to the changes in flow between 1996 and 2016. Indeed, feedback mechanisms are likely to cause an important interdependence between geometry-induced changes in ice flow, shear softening, and/or changes in basal shear stress. Reliable observations of changes in rheology and basal properties are not available, but numerical optimization simulations can provide valuable insights into their evolution. We used the inverse method as described in Sect. 2.2.1 and Appendix B to estimate necessary bounds on the magnitude and spatial distribution of changes in A and C that are required, beside the geometrical changes already applied, to produce the speed-up of PIG between 1996 and 2016. Changes in A and C are treated separately.

• ${\mathcal{E}}_A^m$. The aim of these experiments is to determine possible changes in the rate factor between 1996 (A₉₆) and 2016 (A₁₆). A₉₆ was previously obtained in part 1 (optimization step) of the experimental design. To estimate A₁₆, an inverse optimization problem was solved for the 2016 PIG geometry (H₁₆) and velocities (U₁₆), but using a cost function that was minimized with respect to A only. The slipperiness C was kept fixed to its 1996 solution.

• ${\mathcal{E}}_C^m$. These experiments are analogous to ${\mathcal{E}}_A^m$, but the cost function in the inverse problem was optimized with respect to C only, whereas the rate factor A was kept fixed to its 1996 solution.

3 Results and discussion

3.1 Ice dynamic response to changes in geometry between 1996 and 2016

We present results for the first set of perturbation experiments, which simulate the impact of observed changes in geometry on the flow of PIG. As detailed in Sect. 2.2.2, perturbations are split between four separate cases: (1) calving (${\mathcal{E}}_{\mathrm{Calv}}^{\mathrm{3}}$); (2) thinning of the ice shelf (${\mathcal{E}}_{\mathrm{ISThin}}^{\mathrm{3}}$); (3) thinning of the ice shelf and grounded ice (${\mathcal{E}}_{\mathrm{Thin}}^{\mathrm{3}}$), which includes associated movement of the grounding line and changes in basal traction; and (4) the combined impact of all the above (${\mathcal{E}}_{\mathrm{CalvThin}}^{\mathrm{3}}$). We did not previously specify the value of the sliding exponent; however, here we set m=3, which is a commonly adopted value in ice flow modelling and describes a non-linear viscous (or Weertman) bed rheology. Results for different values of m will be explored in Sect. 3.3.

Results for the relative change in surface speed, $\left(U_{\mathrm{pert}}-U_{\mathrm{96}}\right)/\left(U_{\mathrm{16}}-U_{\mathrm{96}}\right)$, for each of the above perturbations are presented in Fig. 3a–d. In addition to spatial maps of relative velocity changes, we present flux calculations for two gates perpendicular to the flow within the central part of PIG, as displayed in Fig. 3a. Gate 1 is situated about 50 km upstream of the 2016 grounding line and captures the inland propagation of changes in ice flow. Gate 2 approximately coincides with the 2016 grounding line position and captures changes in grounding line flux, which is a direct measure for PIG's increasing contribution to sea level rise and an important indicator of change.

Figure 3 Modelled changes in surface speed compared to 1996 for prescribed perturbations of the Pine Island Glacier geometry. (a) Retreat of the calving front. (b) Thinning of the ice shelf. (c) Thinning of the ice shelf and grounded ice, including grounding line retreat. (d) Calving and thinning combined. For each perturbation, the modelled change in speed (U_pert−U₉₆) is expressed as a percentage of the observed speed-up between 1996 and 2016 (U₁₆−U₉₆). Dashed black lines correspond to the 50 % contour. Panel (e) shows the percentage of the observed flux changes through Gate 1 and 2 that can be explained by the respective perturbations. The simulated impact of calving and thinning in experiment ${\mathcal{E}}_{\mathrm{CalvThin}}^{\mathrm{3}}$ only represents 28 and 64 % of the measured flux changes respectively. Possible explanations for the unaccounted-for increase in flow speed are provided in Sect. 3.2 and 3.3 .

Calving as simulated in ${\mathcal{E}}_{\mathrm{Calv}}^{\mathrm{3}}$ causes changes in flow speed that are predominantly restricted to the outer ice shelf, where it accounts for up to 50 % of the observed speed-up between 1996 and 2016 (Fig. 3a). A smaller dynamical impact is also felt upstream of the grounding line, caused by the calving-induced reduction in ice shelf buttressing and mechanical coupling between the floating and grounded ice. Along the fast-flowing central trunk of PIG, calving typically accounts for less than 10 % of the observed speed-up, with little or no effect on the dynamics of the upstream tributaries. Our results are consistent with earlier work by Schmeltz et al. (2002), in particular their calving scenario “part 2”. The only area with negative relative changes in our simulation is the western shear margin of the ice shelf, where modelled and observed changes in flow speed have the opposite sign. Extensive damage, a process that is not captured by this experiment, has caused this margin to migrate, and significant interannual variations in flow speed have been reported by Christianson et al. (2016). Figure 3e shows that calving accounts for 2 and 13 % of the observed flux changes through Gate 1 and 2 respectively, which confirms the minor instantaneous changes to the flow upstream of the grounding line.

Thinning of the ice shelf as simulated in experiment ${\mathcal{E}}_{\mathrm{ISThin}}^{\mathrm{3}}$ induces a flow response that is similar to calving, as shown in Fig. 3b, and indicates that calving and ice shelf thinning have caused a comparable perturbation in the buttressing forces. The largest percentage changes are found on the ice shelf and are typically less than 25 %, while the relative flux changes through Gate 1 and 2 are identical to the calving experiment (Fig. 3e). Ice shelf thinning is generally accepted to be the main driver of ongoing mass loss of PIG, and patterns of ice shelf thinning elsewhere in Antarctica are strongly correlated to observed changes in grounding line flux (Reese et al., 2018; Gudmundsson et al., 2019). However, the force perturbations that result from ice shelf thinning alone, in particular the instantaneous reduction in back forces τ_IS, are not sufficient to reproduce the magnitude of observed changes in upstream flow, consistent with previous studies (Seroussi et al., 2014; Joughin et al., 2010, 2019). Indeed, experiment ${\mathcal{E}}_{\mathrm{ISThin}}^{\mathrm{3}}$ demonstrates that the direct and instantaneous contribution of ice shelf thinning to observed changes in grounding line flux is less than 25 %. Instead, time-evolving changes in geometry and mass redistribution upstream of the grounding line, which may cause grounding line retreat and associated loss of basal traction, play a significant role in increasing the dynamic response of the glacier. These dynamic changes, caused indirectly by changes in the calving front position and ice shelf thinning, were not captured by the experiments ${\mathcal{E}}_{\mathrm{Calv}}^{\mathrm{3}}$ and ${\mathcal{E}}_{\mathrm{ISThin}}^{\mathrm{3}}$ but are considered in experiment ${\mathcal{E}}_{\mathrm{Thin}}^{\mathrm{3}}$.

In experiment ${\mathcal{E}}_{\mathrm{Thin}}^{\mathrm{3}}$ we prescribed the time-integrated change in ice thickness between 1996 and 2016 for both the floating ice shelf and upstream grounded ice. This perturbation incorporates the observed recession of the PIG grounding line between 1996 and 2016. The combined reduction in ice shelf buttressing, loss of basal friction due to grounding line retreat, and changes in driving stress caused a significant and far-reaching impact on the flow, as displayed in Fig. 3c. Modelled changes on the ice shelf are consistent with and similar in amplitude to ${\mathcal{E}}_{\mathrm{ISThin}}^{\mathrm{3}}$. Upstream of the grounding line, modelled changes relative to observations are between 25 and 50 % along the central trunk and up to 100 % along the tributaries. In addition, results demonstrate that glacier-wide changes in ice thickness account for 26 and 45 % of the observed changes in ice flux through Gate 1 and 2 respectively (Fig. 3e).

In the final perturbation experiment, ${\mathcal{E}}_{\mathrm{CalvThin}}^{\mathrm{3}}$, the combined effect of calving and changes in ice thickness was simulated. Modelled versus observed changes in surface speed are shown in Fig. 3d. The spatial pattern is consistent with previous experiments, and the amplitude of the response is approximately equal to the added response of experiments ${\mathcal{E}}_{\mathrm{Calv}}^{\mathrm{3}}$ and ${\mathcal{E}}_{\mathrm{Thin}}^{\mathrm{3}}$, i.e. $\mathrm{\Delta }{U}_{\mathrm{CalvThin}}\approx \mathrm{\Delta }{U}_{\mathrm{Calv}}+\mathrm{\Delta }{U}_{\mathrm{Thin}}$. The corresponding percentage changes in ice flux through Gate 1 and 2 are 28 and 64 % respectively, whereas modelled changes in flow across the actual grounding line account for about 75 % of the observed increase in flux between the years 1996 and 2016. Although this experiment prescribes all observed changes in PIG geometry over the observational period, model simulations are unable to capture a significant percentage of the observed speed-up. This is most noticeable along the fast-flowing central trunk upstream of the grounding line, whereas discrepancies decrease along the slow-flowing tributaries in the high catchment. We also note that, in one area between Gate 1 and 2, modelled and observed changes in surface speed have opposite signs.

Although it is not unexpected to find differences between diagnostic model output and observations, the consistently suppressed response of the model to realistic perturbations in ice geometry is indicative of a structural shortcoming within our experimental design. Indeed, results show that, for a non-linear viscous bed rheology described by a Weertman sliding law with constant sliding coefficient m=3, changes in ice geometry alone cannot account for the complex and spatially variable pattern of speed-up over the observational period, i.e. $U_{\mathrm{16}}-U_{\mathrm{96}}\ne \mathrm{\Delta }{U}_{\mathrm{CalvThin}}$. In the remainder of this study, two possible hypotheses are analysed that enable the gap to be closed between geometry-induced changes in ice flow and the observed speed-up of PIG. The first hypothesis, which is considered in Sect. 3.2, assumes that bed deformation can indeed be described by a non-linear viscous power law with m=3, but further temporal variations in ice viscosity and/or basal slipperiness are required in addition to changes in geometry: $U_{\mathrm{16}}-U_{\mathrm{96}}=\mathrm{\Delta }{U}_{\mathrm{CalvThin}}+\mathrm{\Delta }{U}_A+\mathrm{\Delta }{U}_C$. The second, alternative hypotheses, discussed in Sect. 3.3, assumes that internal properties of the ice and bed have not significantly changed between the years 1996 and 2016, i.e. ΔU_A≈0 and ΔU_C≈0, but a different physical description of the basal rheology is required instead.

3.2 Changes in the rate factor and basal slipperiness between 1996 and 2016

In transient model simulations of large ice masses such as Antarctica's glaciers and ice streams, it is common to assume that the advection of A with the ice, or changes due to temperature variations and fracture as well as changes in basal slipperiness C, exerts a second-order control on changes in ice flow. As such, temporal variability in A and C is often ignored, based on the assumption that these changes are sufficiently slow and do not significantly affect the flow on typical decadal to centennial timescales under consideration. The aim of experiments ${\mathcal{E}}_A^{\mathrm{3}}$ and ${\mathcal{E}}_C^{\mathrm{3}}$, as outlined in Sect. 2.2.3, is to establish whether this is a valid assumption, or whether previously ignored changes in A and/or C can provide a realistic explanation for the discrepancy between simulated and observed changes in the surface speed of PIG in the geometric experiment ${\mathcal{E}}_{\mathrm{CalvThin}}^{\mathrm{3}}$. Experiment ${\mathcal{E}}_A^{\mathrm{3}}$ assumes that, in addition to changes in geometry, temporal variations in A alone are able to account for the significant increase in flux that was unaccounted for in previous experiments. Alternatively, ${\mathcal{E}}_C^{\mathrm{3}}$ assumes that, in addition to changes in geometry, temporal variations in C alone are able to resolve the discrepancy in Sect. 3.1 between the modelled and observed speed-up. In line with previous experiments we assume a Weertman sliding law with m=3. The results for both experiments are summarized in Fig. 4.

Figure 4 (a) Results for the ${\mathcal{E}}_A^{\mathrm{3}}$ experiment: changes in the rate factor A required to fully reproduce observed changes in surface speed of the ice shelf and grounded ice between the years 1996 and 2016. The sliding exponent m=3 and basal slipperiness C were kept fixed for grounded areas. Magenta contours (in m yr⁻¹) correspond to the surface speed in 2016. (b) Results for the ${\mathcal{E}}_C^{\mathrm{3}}$ experiment: changes in the basal slipperiness C required to reproduce the observed increase in surface speed of the grounded ice between 1996 and 2016. The rate factor A is assumed constant between 1996 and 2016.

Changes in A (Fig. 4a) needed to fully reproduce the speed-up of PIG between the years 1996 and 2016 are spatially coherent and predominantly positive. This suggests a reduction in ice viscosity between 1996 and 2016, as a result of localized heating, enhanced damage within the ice column, or changes in anisotropy. The largest changes are found in distinct geographical areas: a localized increase within the shear margins of the ice shelf and a more widespread increase along the slower-moving flanks (magenta contours in Fig. 4a indicate surface speed in 2016) of the main glacier and westernmost tributary, about 20 km upstream of the 2016 grounding line. Changes within the ice shelf shear margins are consistent with their increasingly complex and damaged morphology, as is apparent from satellite images (Alley et al., 2019). Weakening of the ice in these areas is sufficient to account for the remaining 50 % of observed changes in ice shelf speed-up that could not previously be reproduced by calving and ice shelf thinning alone (experiment ${\mathcal{E}}_{\mathrm{CalThin}}^{\mathrm{3}}$). Projected changes in A along the flanks of the upstream glacier, on the other hand, are more ambiguous. Values in excess of 10⁻⁷ yr⁻¹ kPa⁻³ correspond to an equivalent increase in “ice” temperature by up to 40 ^∘C. This is nonphysical unless (part of) the change is attributed to damage or evolution of the ice fabric. Based on our analysis of Sentinel and Landsat satellite images, there is no obvious indication of recent changes in the surface morphology in these areas. Either significant and widespread changes in the thermal and mechanical properties have occurred beneath the surface or the observed speed-up and thinning in these areas, as previously reported by Bamber and Dawson (2020), cannot be convincingly attributed to changes in the rate factor.

Alternatively, temporal changes in C can be invoked to reproduce the discrepancies between modelled and observed changes in surface speed between the years 1996 and 2016. Results presented in Fig. 4b suggest that a complex and widespread pattern of changes in the slipperiness is required across an extensive portion of PIG's central basin and its upper catchment. Despite the complex and poorly understood relationship between C and quantifiable physical properties of the ice–bed interface, it is difficult to understand how any single process or combination of physical processes could be responsible for the large and widespread changes in C over a time period of 2 decades. Further information, such as a time series of maps similar to Fig. 4b, can potentially be used to test the robustness of this result and provide further insights into the physical processes that could control such changes. This is the subject of future research.

We note that, in the ${\mathcal{E}}_C^{\mathrm{3}}$ experiment, velocities on the floating ice shelf were largely unaffected by changes in C and remained significantly slower than observations (not shown). In contrast, changes in the rate factor were able to fully account for the speed-up of the ice shelf. On the other hand, large variations in A were needed to reproduce the changes in ice dynamics along the slow-moving flanks of PIG (Fig. 4a), whereas only small changes in C less than 10⁻³ yr⁻¹ kPa⁻³ m were required to explain this behaviour. It is therefore conceivable that, in addition to PIG's evolving geometry, an intricate combination of changes in both the rate factor and basal slipperiness are required to reproduce the glacier's complex and spatially diverse patterns of speed-up over the last 2 decades. It is however not straightforward to disentangle these processes in the current modelling framework.

3.3 Evidence for a heterogeneous bed rheology

The relationship between changes in geometry and the dynamic response of a glacier crucially depends on the mechanical properties of the underlying bed and subglacial hydrology. So far, we have assumed that basal sliding can be represented by a non-linear viscous power law with spatially uniform stress exponent m=3 (see Eq. 3). A power-law rheology is particularly suitable for the description of hard-bedded sliding without cavitation (Weertman, 1957), but missing processes such as variations in effective pressure or the deformation of a subglacial till layer with a maximum shear (yield) stress could be important limitations. Some evidence has been provided for plastic bed properties underneath ice streams either from observations (Tulaczyk et al., 2000; Minchew et al., 2016) or from laboratory experiments (Zoet and Iverson, 2020). Most recently, Gillet-Chaulet et al. (2016), Brondex et al. (2019), and Joughin et al. (2019) used numerical simulations to show that different sliding laws can cause a distinctly different dynamical response of PIG to changes in geometry, and observed changes in surface velocity were best reproduced for sliding exponents m≫1 or using a hybrid law that combines power law with Coulomb sliding. Although the results are compatible with a plastic bed underlying the central trunk of PIG, no constraints on the spatial variability in basal rheology were derived.

Figure 5 Dependency of simulated-versus-observed changes in surface speed on the sliding-law exponent: (a) m=1, (b) m=7, and (c) m=13. Dashed black lines correspond to the 50 % contour. Larger values of m cause an increased response of the modelled surface speed to geometrical changes (calving, thinning, and grounding line retreat). For m>3, the modelled response of slow-flowing ice in the upstream catchment exceeds observed changes by more than a factor of 2, whereas for m=13 modelled changes of the fast-flowing central trunk are still smaller than observed changes. (d) Changes in flux through Gate 1 and 2 as a percentage of observed changes for m=1, 7, and 13.

In order to quantify how different values of the sliding exponent affect the sensitivity of PIG to changes in geometry across the catchment, we repeated perturbation experiments ${\mathcal{E}}_{\mathrm{CalvThin}}^m$ for a range of sliding-law exponents, from m=1 to m=21 at increments of 2. Results for m=1, 7, and 13 are shown in Fig. 5. A linear rheology induces a simulated response to calving and thinning that accounts for less than 50 % of the observed changes everywhere. For m=7, relative changes in flow speed exceed 100 % along significant portions of the slower-flowing tributaries. For m=13, which effectively corresponds to a plastic rheology, the modelled response overshoots observations by more than 100 % in most areas, except along the main glacier, where the response approaches 100 %. Across the model domain, a significant positive correlation exists between m and relative velocity changes, indicating a stronger dynamic response to perturbations in geometry with increasing values of m. This finding is in agreement with Gillet-Chaulet et al. (2016) and Joughin et al. (2019); however our maps show that no single, spatially uniform value of the sliding exponent is able to produce a good match between model output and observations across the entire catchment.

The positive correlation between the flow response and m is an inherent property of the adopted physical description of glacier dynamics. For the shallow ice stream approximation with a non-linear viscous sliding law, the first-order response of the surface velocity, δU, to small perturbations in surface elevation, δS, was previously determined by Gudmundsson (2008) and depends on m in the following non-linear way:

$$\begin{array}{}\text{(4)}& \mathit{\delta }U\equiv \left|{\mathcal{T}}_{\mathrm{US}}\right(m\left)\right|\mathit{\delta }S={\displaystyle \frac{f_{\mathrm{1}}m}{m+f_{\mathrm{2}}}}\mathit{\delta }S.\end{array}$$

The transfer amplitude $\left|{\mathcal{T}}_{\mathrm{US}}\right|$ contains complicated positive functions f₁ and f₂ that generally depend on the wavelength of the surface perturbation, geometrical factors such as the local bed slope, and the basal slipperiness C. Further details are provided in Appendix C. Despite the simplifying assumptions that underlie the analytical expression of $\left|{\mathcal{T}}_{\mathrm{US}}\right|$ obtained by Gudmundsson (2008), results from our simulations ${\mathcal{E}}_{\mathrm{CalvThin}}^{m_i},m_i\in \mathit{\{}\mathrm{1},\mathrm{3},\mathrm{\cdots },\mathrm{21}\mathit{\}}$, indicate that Eq. (4) is also applicable to the more complex setting of PIG. Indeed, as explained in detail in Appendix C, we found that across a large portion of the PIG catchment the transfer amplitude $\left|{\mathcal{T}}_{\mathrm{US}}\right|$ provides a suitable model to describe the dependency of the relative velocity changes $\mathrm{\Delta }{U}_{\mathrm{CalvThin}}/\left(U_{\mathrm{16}}-U_{\mathrm{96}}\right)$ on m. The parameters f₁ and f₂ were treated as spatially variable fields, and best estimates $f_{\mathrm{1}}^{*}\left(\mathit{x}\right)$ and $f_{\mathrm{2}}^{*}\left(\mathit{x}\right)$ were obtained as a solution of the minimization problem

$$\begin{array}{}\text{(5)}& \begin{array}{rl}& \left(f_{\mathrm{1}}^{*}\left(\mathit{x}\right),f_{\mathrm{2}}^{*}\left(\mathit{x}\right)\right)=\underset{f_{\mathrm{1}},f_{\mathrm{2}}}{min}\left({\displaystyle \frac{f_{\mathrm{1}}\left(\mathit{x}\right)m}{m+f_{\mathrm{2}}\left(\mathit{x}\right)}}-{\displaystyle \frac{\mathrm{\Delta }{U}_{\mathrm{CalvThin}}^m\left(\mathit{x}\right)}{U_{\mathrm{16}}-U_{\mathrm{96}}}}\right),\\ & \text{ with}m\in \left\{\mathrm{1},\mathrm{3},\mathrm{\cdots },\mathrm{21}\right\}.\end{array}\end{array}$$

The non-linear dependency of $\mathrm{\Delta }{U}_{\mathrm{CalvThin}}^m/\left(U_{\mathrm{16}}-U_{\mathrm{96}}\right)$ on m can then be approximated by

$$\begin{array}{}\text{(6)}& {\displaystyle \frac{\mathrm{\Delta }{U}_{\mathrm{CalvThin}}^m}{U_{\mathrm{16}}-U_{\mathrm{96}}}}\approx {\displaystyle \frac{f_{\mathrm{1}}^{*}\left(\mathit{x}\right)m}{m+f_{\mathrm{2}}^{*}\left(\mathit{x}\right)}}.\end{array}$$

Using this dependency of the simulated velocity changes on m, one can derive an “optimal” spatial distribution of the sliding exponent, m_optimal(x), such that $\mathrm{\Delta }{U}_{\mathrm{CalvThin}}/\left(U_{\mathrm{16}}-U_{\mathrm{96}}\right)=\mathrm{100}\mathit{\%}$ everywhere, namely

$$\begin{array}{}\text{(7)}& m_{\mathrm{optimal}}\left(\mathit{x}\right)={\displaystyle \frac{f_{\mathrm{2}}^{*}\left(\mathit{x}\right)}{f_{\mathrm{1}}^{*}\left(\mathit{x}\right)-\mathrm{1}}}.\end{array}$$

By construction, the variable sliding exponent m_optimal(x) enables reproducing 100 % of the observed speed-up of PIG in response to calving and ice thickness changes. The results, depicted in Fig. 6a, indicate that plastic bed conditions (m≫1) prevail across most of the fast-flowing central valley and parts of the upstream tributaries. Values generally increase towards the grounding line, whilst linear or weakly non-linear bed conditions are consistently found in the slow-flowing inter-tributary areas. This finding is compatible with the presence of a weak, water-saturated till beneath fast-flowing areas of PIG and a hard bedrock or consolidated till between tributaries (Joughin et al., 2009). The transition to lower exponents in areas with slower flow (<600 m a⁻¹) is also consistent with results based on a Coulomb-limited sliding law, which produces Coulomb plastic behaviour at speeds >300 m a⁻¹ and weakly non-linear viscous sliding at slower speeds (Joughin et al., 2019).

Figure 6 (a) Optimal values of the sliding exponent, required to ensure close agreement between modelled and observed changes in flow velocity of Pine Island Glacier between the years 1996 and 2016. White and black dots mark areas where such an agreement cannot be achieved for different reasons: white dots indicate a poor fit between the transfer function $\left|{\mathcal{T}}_{\mathrm{US}}\right|$ and $\mathrm{\Delta }{U}_{\mathrm{CalvThin}}^{m_i}/\left(U_{\mathrm{16}}-U_{\mathrm{96}}\right),m_i\in \mathit{\{}\mathrm{1},\mathrm{3},\mathrm{\cdots },\mathrm{21}\mathit{\}}$, with R²<0.9; black dots indicate areas where a positive, finite solution for m_optimal in Eq. (7) does not exist and non-linear viscous sliding cannot reproduce observed changes in surface flow. (b) Same as Fig. 3d but for optimal values of the sliding-law exponent in panel (a). (c) Same as Fig. 3e but for optimal values of the sliding-law exponent in panel (a).

Two interesting properties of the regression model in Eq. (4) are worth noting. Firstly, for m→∞, the function $\left|{\mathcal{T}}_{\mathrm{US}}\right|$ approaches a horizontal asymptote with limit equal to f₁. As a consequence, the associated solution for m_optimal diverges to ∞ for locations x where $f_{\mathrm{1}}^{*}\left(\mathit{x}\right)=\mathrm{100}$ and becomes negative where $f_{\mathrm{1}}^{*}\left(\mathit{x}\right)<\mathrm{100}$. In these areas, indicated by black dots in Fig. 6a, no non-negative, finite value of m exists such that $\mathrm{\Delta }{U}_{\mathrm{CalvThin}}\left(\mathit{x}\right)/\left(U_{\mathrm{16}}-U_{\mathrm{96}}\right)=\mathrm{100}\mathit{\%}$, and conventional Weertman sliding is unable to fully reproduce the observed flow changes in response to thickness changes and calving. Either a different form of the sliding law is required or additional changes in the rate factor A and/or basal slipperiness C are needed. These findings are the subject of a forthcoming study. Our second observation concerns locations where U₁₆ or U₉₆ contains significant measurement uncertainties, or where no discernible changes in the surface velocity were measured, i.e. $U_{\mathrm{16}}-U_{\mathrm{96}}\approx \mathrm{0}$. In these areas, the non-linear regression was generally found to be poor, with R² values smaller than 0.9 as indicated by the white dots in Fig. 6a. As no reliable estimate for m_optimal could be obtained for areas shaded in white or black in Fig. 6a, values were instead based on a nearest-neighbour interpolation.

It is important to reiterate that the regression method used crucially relies on non-trivial measurements of changes in surface velocity ($U_{\mathrm{16}}-U_{\mathrm{96}}\ne \mathrm{0}$) and cannot be used to retrieve information about the basal rheology of ice bodies that are presently in steady state. It should also be noted that values of $f_{\mathrm{1}}^{*}\left(\mathit{x}\right)$ and $f_{\mathrm{2}}^{*}\left(\mathit{x}\right)$ were derived independently for each node of the computational mesh, whereas the continuum mechanical properties of glacier flow would suggest a non-zero spatial covariance $〈f_{\mathrm{1}}\left({\mathit{x}}_{\mathbf{1}}\right),f_{\mathrm{1}}\left({\mathit{x}}_{\mathbf{2}}\right)〉\ne \mathrm{0}$ and $〈f_{\mathrm{2}}\left({\mathit{x}}_{\mathbf{1}}\right),f_{\mathrm{2}}\left({\mathit{x}}_{\mathbf{2}}\right)〉\ne \mathrm{0}$. The optimal solution for m is therefore not automatically mesh independent or robust with respect to the amount of regularization in the inversion. This concern is discussed further in Appendix D.

In order to demonstrate the improved model response to thinning and calving for a spatially variable sliding exponent m_optimal(x), we performed a new inversion with m_optimal(x) and subsequently repeated the geometric perturbation experiments ${\mathcal{E}}_{*}^{\mathrm{optimal}}$. The results are presented in Fig. 6b and c. Compared to spatially uniform values of m (Fig. 3d and Fig. 5), a spatially variable basal rheology generally improves the fit between observed changes in flow and the modelled response across the entire basin. Based on the flux changes through Gate 1 and 2, we find that (1) calving and ice thickness changes in combination with a spatially variable, predominantly plastic bed rheology account for 67 and 105 % of flux changes through Gate 1 and 2 respectively, compared to 28 and 64 % for a uniform non-linear viscous sliding law with exponent m=3; that (2) calving and ice shelf thinning caused an almost identical response in ice dynamics upstream of the grounding line; and that (3) dynamic thinning and grounding line movement account for most of the flux changes between the years 1996 and 2016. The remaining mismatch between the observed and modelled response in Fig. 6b can, at least in part, be attributed to uncertainties in m_optimal(x). This is of particular relevance in the vicinity of the grounding line and for parts of the central trunk, where the non-linear regression method in Eq. (4) did not provide a reliable or finite estimate for m_optimal. Previous studies, for example by Gillet-Chaulet et al. (2016) and Joughin et al. (2019), have demonstrated a better agreement between modelled and observed speed-up using Coulomb-limited sliding laws, such as those proposed by Budd et al. (1984), Schoof (2006), and Tsai et al. (2015). Our results are consistent with these earlier studies and suggest that power-law sliding does not adequately capture the physical relationship between basal shear stress and sliding in the vicinity of the grounding line.

4 Conclusions

Based on the most comprehensive observations of ice shelf and grounded ice thickness changes to date, and a suite of diagnostic model experiments with the contemporary flow model Úa, we have analysed the relative importance of ice shelf thinning, calving, and grounding line retreat for the speed-up of Pine Island Glacier over the period 1996 to 2016. The detailed comparison between simulated and observed changes in flow speed has provided insights into the ability of a modern-day ice flow model to reproduce dynamic changes in response to prescribed geometric perturbations. Significant discrepancies between observed and modelled changes in flow were found and addressed either by allowing changes in ice viscosity and basal slipperiness or by varying the mechanical properties of the ice–bed interface. For non-linear viscous sliding at the bed, geometric perturbations could only account for 64 % of the observed flux increases close to the grounding line, whereas the remaining 36 % could be attributed to large and widespread changes in ice viscosity (including damage) and/or changes in basal slipperiness. Under the alternative assumption that ice viscosity and basal slipperiness did not change considerably over the last 2 decades, we found that the recent increase in flow speed of Pine Island Glacier is only compatible with observed patterns of thinning if a heterogeneous, predominantly plastic bed underlies large parts of the central glacier and its upstream tributaries, consistent with the earlier literature.

Appendix A: Observations of Pine Island Ice Shelf thickness changes between 1996 and 2016

We derived a new ice shelf height time series from measurements acquired by four overlapping ESA satellite radar altimetry (RA) missions: ERS-1 (1991–1996), ERS-2 (1995–2003), Envisat (2002–2012), and CryoSat-2 (2010–present). For this study, we constructed a record of ice shelf height spanning 20 years (1996–2016), with a temporal sampling of 3 months. We integrated all measurements along the satellite ground tracks and gridded the solution on a 3 km by 3 km grid.

Our adopted processing steps for RA data are a modification/improvement from Paolo et al. (2016) and Nilsson et al. (2016). Specifically for CryoSat-2, we retracked ESA’s SARIn L1B product over the Antarctic ice shelves using the approach by Nilsson et al. (2016), corrected for a 60 m range offset for data with surface types “land” or “closed sea”, and removed points with anomalous backscatter values (>30 dB). We estimated heights with a modified (from McMillan et al., 2014) surface-fit approach, with a variable rather than constant search radius to account for the RA heterogeneous spatial distribution, calculating mean values along the satellite reference tracks; we removed height estimates less than 2 m above the EIGEN-6C4 geoid (Chuter and Bamber, 2015) to account for ice shelf mask imperfections near the calving front; we applied all of the standard corrections to altimeter data over ice shelves (for example, removing gross outliers and residual heights with respect to mean topography >15 m); we ran an iterative 3σ filter; we minimized the effect of variations in backscatter (Paolo et al., 2016); and we corrected for ocean tides (Padman et al., 2002) and inverse barometer effects (Padman et al., 2004).

We then gridded the height data in space and time on a $\mathrm{3}\mathrm{km}\times \mathrm{3}\mathrm{km}\times \mathrm{3}$-month cube, for each mission independently. We merged the records (all four satellites) by only accepting time series that overlapped by at least three quarters of a year to ensure proper cross-calibration, and removed (and subsequently interpolated) anomalous data points that deviated from the trend by more than 5 SD. This removes data with, for example, satellite mispointing, anomalous backscatter fluctuations, grounded-ice contamination, high surface slopes, and geolocation errors. We fitted linear trends to the gridded product to obtain the ΔH field used in our model experiments (see Sect. 2.1). We also removed a 3 km buffer around the ice shelf boundaries to further mitigate floating–grounded mask imperfections and the limitation of geophysical corrections within the ice shelf flexural zone.

The thickness changes for the ice shelf were combined with existing data for thickness changes over the same time period on the grounded ice (Shepherd et al., 2016). The resulting dataset for ΔH, as used in the experiments described in Sect. 2.2, is shown in Fig. A1. The figure shows the data grids, including the 3 km buffer downstream of the 1996 grounding line, and other data-sparse areas along the central flow line. Here, thickness changes were obtained through linear interpolation from neighbouring data. The grounding line location associated with our 1996 thickness distribution was compared to independent measurements from DInSAR (Rignot et al., 2014), and both agree well (Fig.A1).

Figure A1 Ice thickness changes (ΔH) between 1996 and 2016, based on a comprehensive analysis of satellite altimeter data. The altimeter data coverage is represented by dots (ice shelf) and circles (grounded ice; Shepherd et al., 2016). The final 1996 ice thickness distribution was obtained by subtracting ΔH from the 2016 BedMachine ice thickness (Morlighem et al., 2020), as described in Sect.2.1. The associated 1996 grounding line location (blue line) compares well to independent DInSAR measurements (magenta line; Rignot et al., 2014).

Appendix B: Model configuration and optimization

The open source ice flow model Úa (Gudmundsson, 2020) uses the finite-element method to solve the shallow ice stream equations, commonly referred to as SSA or SSTREAM (Hutter, 1983; MacAyeal, 1989), on an irregular triangular mesh. The diagnostic velocity solver is based on an iterative Newton–Raphson method. A fixed mesh with 109 300 linear elements was used with a median nodal spacing of 1.2 km and local mesh refinement down to 500 m in areas with above-average horizontal shear, with strong gradients in ice thickness, and within a 10 km buffer around the grounding line. The mesh was generated using the open-source generator mesh2d (Engwirda, 2014).

Figure B1 (a) L curve used to determine the optimal value of the Tikhonov regularization multiplier γ_s, highlighted in red. (b) Misfit between modelled and observed surface speed in 1996 for γ_s=25 000 m. (c) Rate factor (A in Eq. 2) in 1996, obtained as a minimum of the cost function J in Eq. (B1) with γ_s=25 000 m. The equivalent depth-averaged ice temperature ranges from −35 ^∘C (grey) to 5 ^∘C (red). Colours are discretized at 5 ^∘C intervals, and the black lines indicate the 0 ^∘C contour. The white line corresponds to the 1996 grounding line position. (d) Optimal value of the basal slipperiness (C in Eq. 3) in 1996, estimated using the adjoint minimization approach.

The optimization capabilities of Úa follow commonly applied techniques in ice flow modelling to optimize uncertain model parameters, p_i, based on prior information, ${\widehat{p}}_i$, and a range of observations with associated measurement errors (MacAyeal, 1992). Úa uses an adjoint method to obtain a combined optimal estimate of the spatially varying rate factor A and basal slipperiness C across the full model domain, for given observations of surface velocity u_obs and measurement errors ε_u. Optimal values for $p_i\in \left\{A,C\right\}$ were obtained as a solution to the minimization problem d_pJ with the cost function J defined as the sum of the misfit term I and Tikhonov regularization R: $J=I+R$, with

$$\begin{array}{ll}\text{(B1)}& {\displaystyle }I& {\displaystyle }={\displaystyle \frac{\mathrm{1}}{\mathrm{2}\mathcal{A}}}\int \mathrm{d}\mathit{x}{\left(u_{\mathrm{model}}-u_{\mathrm{obs}}\right)}^{\mathrm{2}}/{\mathit{\epsilon }}_u^{\mathrm{2}},\text{(B2)}& {\displaystyle }R& {\displaystyle }={\displaystyle \frac{\mathrm{1}}{\mathrm{2}\mathcal{A}}}\int \mathrm{d}\mathit{x}\sum _i\left({\mathit{\gamma }}_{i,s}^{\mathrm{2}}{\left(\mathrm{\nabla }{\log }_{\mathrm{10}}(p_i/\widehat{p_i})\right)}^{\mathrm{2}}\right.{\displaystyle }& {\displaystyle }\left.+{\mathit{\gamma }}_{i,a}^{\mathrm{2}}{\left({\log }_{\mathrm{10}}\left(p_i/\widehat{p_i}\right)\right)}^{\mathrm{2}}\right),\end{array}$$

and $\mathcal{A}=\int \mathrm{d}\mathit{x}$ being the total area of the model domain. An iterative interior point optimization algorithm was used to calculate d_pJ and stopped after 10⁴ iterations, when fractional changes to the cost function were less than 10⁻⁵.

The gradient and amplitude contributions in the regularization term (Eq. B2) are multiplied by spatially constant Tikhonov regularization multipliers, γ_i,s and γ_i,a. Optimal values for γ_i,s and γ_i,a were determined using an L-curve approach. For γ_i,s results are shown in Fig. B1. The values ${\mathit{\gamma }}_{A,s}={\mathit{\gamma }}_{C,s}=\mathrm{25}\mathrm{000}$ m were used for all experiments in the main part of the text, as they produced the smallest misfit between observed and modelled surface velocities whilst limiting the risk of overfitting. The sensitivity of the main results with respect to the choice of γ_s is discussed in Appendix D. A similar L-curve approach was followed to determine optimal values for γ_i,a, with ${\mathit{\gamma }}_{A,a}={\mathit{\gamma }}_{C,a}=\mathrm{1}$ used throughout this study.

The pre-multipliers γ_i,s and γ_i,a are constant across the model domain. Some studies set ${\mathit{\gamma }}_{A,a}={\mathit{\gamma }}_{A,s}=\mathrm{0}$ for grounded areas and only optimize the rate factor on the ice shelves. This approach assumes perfect prior knowledge about the spatial distribution and magnitude of the rate factor upstream of the grounding line, often tied to (uncertain) estimates of ice temperature. Here we prefer to optimize A across the full domain to allow for spatial variations in ice temperature, damage, fabric, and other ice properties for both floating and grounded ice. The amplitude and gradient of A, relative to a spatially constant prior value, are controlled by γ_A,a and γ_A,s respectively, with optimal values given above.

All results presented here are based on optimization experiments with spatially constant a priori values for the rate factor and slipperiness: $\widehat{A}=\mathrm{5.04}\times {\mathrm{10}}^{-\mathrm{9}}$ kPa⁻³ yr⁻¹, which corresponds to a spatially uniform ice temperature of −15 ^∘C (Cuffey and Paterson, 2010), and $\widehat{C}=u_{\mathrm{b}}{\mathit{\tau }}^{-m}$, with u_b=750 m yr⁻¹ and τ=80 kPa and m being the sliding-law exponent. Different a priori values for the rate factor, equivalent to ice temperatures between −20 and −5 ^∘C, were tested but did not cause significant differences in the results. We did not consider optimization experiments with spatially variable $\widehat{A}$ based on independent estimates of ice temperature and/or damage, since such estimates contain significant uncertainties at the regional scales considered in this study.

Figures B1b–d summarize the results for an optimization with ${\mathit{\gamma }}_{i,s}=\mathrm{25}\mathrm{000}$ m, ${\mathit{\gamma }}_{i,a}=\mathrm{1}$, $\widehat{A}=\mathrm{5.04}\times {\mathrm{10}}^{-\mathrm{9}}$ kPa⁻³ yr⁻¹, and $\widehat{C}=\mathrm{1.46}\times {\mathrm{10}}^{-\mathrm{3}}$ m yr⁻¹ kPa⁻³. Modelled surface velocities are typically within 30 m per year or less of the observed values, with a mean misfit of $-\mathrm{1.68}\mathrm{m}{\mathrm{yr}}^{-\mathrm{1}}$ and standard deviation of 15.3 m yr⁻¹. The highest values of the rate factor are generally found within the shear margins, with positive equivalent ice temperatures suggesting the presence of a complex rheology or damage. The highest values of the slipperiness are consistently found in the fast-flowing central part of the glacier and along its upstream tributaries, with noticeably reduced values of C in an area between 5 and 40 km upstream of the 1996 grounding line. These results are broadly in agreement with previously published maps; see, for example, Arthern et al. (2015).

Appendix C: Non-linear dependency of the flow response on the sliding exponent

The transfer amplitude $\left|{\mathcal{T}}_{\mathrm{US}}\right|$, defined in Eq. (4), describes the linear response of the along-slope surface velocity to small harmonic perturbations in the surface elevation or, equivalently, ice thickness. Analytical solutions for the transfer function 𝒯_US (amplitude and phase) in the framework of the shallow ice stream approximation with a linear ice rheology (n=1 in Eq. 2) and a non-linear viscous sliding law (arbitrary m in Eq. 3) were previously obtained by Gudmundsson (2008). Note that the original expression (Eq. 29 in Gudmundsson, 2008) contained a printing error, so we repeat the correct form here:

$$\begin{array}{}\text{(C1)}& \begin{array}{rl}& {\mathcal{T}}_{\mathrm{US}}=\\ & {\displaystyle \frac{{\mathit{\tau }}_{\mathrm{d}}\left[m\mathit{\gamma }\left(\mathrm{1}+\mathit{\psi }\right)+\mathit{\eta }H\left(j^{\mathrm{2}}\mathit{\psi }+k^{\mathrm{2}}+\mathrm{4}{l}^{\mathrm{2}}\right)\right]}{Hm{\mathit{\gamma }}^{\mathrm{2}}+\mathit{\gamma }\mathit{\eta }{H}^{\mathrm{2}}\left[l^{\mathrm{2}}\left(\mathrm{4}+m\right)+k^{\mathrm{2}}\left(\mathrm{1}+\mathrm{4}m\right)\right]+\mathrm{4}{H}^{\mathrm{3}}{j}^{\mathrm{4}}{\mathit{\eta }}^{\mathrm{2}}}},\end{array}\end{array}$$

where H is the local ice thickness; α is the local bed slope; ρ is the ice viscosity; τ_d=ρgHsinα is the driving stress; η is the effective viscosity; and $\mathit{\gamma }=\frac{{\mathit{\tau }}_{\mathrm{d}}^{\mathrm{1}-m}}{mC}$, ψ=ikHcotα, and $j^{\mathrm{2}}=k^{\mathrm{2}}+l^{\mathrm{2}}$ are abbreviations, with k and l the along-slope and transverse wavelength respectively of the harmonic surface perturbation. Since we focus on the instantaneous response of the velocity to perturbations at the surface, the exponential decay of 𝒯_US with time has been omitted. An equivalent expression for the response of the transverse velocity component can be derived; we refer to Gudmundsson (2008) for more details.

Figure C1 (a) Goodness of fit between $\frac{f_{\mathrm{1}}m}{m+f_{\mathrm{2}}}$ and model simulations $\mathrm{\Delta }{U}_{\mathrm{CalvThin}}^{m_i}/\left(U_{\mathrm{16}}-U_{\mathrm{96}}\right),m_i\in \mathit{\{}\mathrm{1},\mathrm{3},\mathrm{\cdots },\mathrm{21}\mathit{\}}$. Red areas correspond to R²≥0.9 and fitting parameter f₁≥100. An example of the fit at location 1 and resulting m_optimal (Eq. 7) is shown in panel (b). Black areas in panel (a) correspond to R²≥0.9 and fitting parameter f₁<100. The horizontal asymptote with limit <100 indicates that a positive, finite solution m_optimal does not exist, and Weertman sliding cannot reproduce 100 % of the observed changes in surface velocity. An example of the fit and asymptote at location 2 is shown in panel (c). Grey areas in panel (a) correspond to R²<0.9, indicating a poor fit between $\frac{f_{\mathrm{1}}m}{m+f_{\mathrm{2}}}$ and $\mathrm{\Delta }{U}_{\mathrm{CalvThin}}^{m_i}/\left(U_{\mathrm{16}}-U_{\mathrm{96}}\right),m_i\in \mathit{\{}\mathrm{1},\mathrm{3},\mathrm{\cdots },\mathrm{21}\mathit{\}}$. An example at location 3 is shown in panel (d).

Following Gudmundsson (2008), physical quantities can be rescaled to obtain the non-dimensional form of the transfer function. After substitution of the scalings H→1, $\mathit{\eta }\to \mathrm{1}/\mathrm{2}$, and τ_d→1 into Eq. (C1) and some reordering, one obtains

$$\begin{array}{}\text{(C2)}& {\mathcal{T}}_{\mathrm{US}}={\displaystyle \frac{m\left[\frac{\mathrm{1}}{C}\left(\mathrm{1}+\mathit{\psi }\right)+\frac{\mathrm{1}}{\mathrm{2}}\left(j^{\mathrm{2}}\mathit{\psi }+k^{\mathrm{2}}+\mathrm{4}{l}^{\mathrm{2}}\right)\right]}{m\left[j^{\mathrm{4}}+\frac{\mathrm{1}}{\mathrm{2}C}\left(l^{\mathrm{2}}+\mathrm{4}{k}^{\mathrm{2}}\right)\right]+\frac{\mathrm{1}}{C^{\mathrm{2}}}+\frac{\mathrm{1}}{\mathrm{2}C}\left(\mathrm{4}{l}^{\mathrm{2}}+k^{\mathrm{2}}\right)}}.\end{array}$$

The resulting transfer amplitude takes the form $\left|{\mathcal{T}}_{\mathrm{US}}\right|=\frac{f_{\mathrm{1}}m}{m+f_{\mathrm{2}}}$ as in Eq. (4), where functions f₁ and f₂ depend on C, α, k, and l.

The analytical expression in Eq. (C2) describes the first-order response to small perturbations in ice thickness, δH≪1, for well-defined length scales characterized by k and l. However, in a realistic setting such as PIG, the system responds to a complicated perturbation composed of a range of wavelengths and amplitudes, and Eq. (C2) does not automatically hold. Based on experiments ${\mathcal{E}}_{\mathrm{CalvThin}}^{m_i},m_i\in \mathit{\{}\mathrm{1},\mathrm{3},\mathrm{\cdots },\mathrm{21}\mathit{\}}$, presented in Sect. 3.3, we found that the simulated surface response of PIG to observed geometrical perturbations retains it dependency on m of the form $\frac{f_{\mathrm{1}}m}{m+f_{\mathrm{2}}}$, but more complicated expressions for f₁ and f₂ are required that do not exist in analytical form. A best estimate for the spatially varying fields f₁ and f₂ was obtained by minimizing the misfit between $\mathrm{\Delta }{U}_{\mathrm{CalvThin}}^{m_i}/\left(U_{\mathrm{16}}-U_{\mathrm{96}}\right),m_i\in \mathit{\{}\mathrm{1},\mathrm{3},\mathrm{\cdots },\mathrm{21}\mathit{\}}$, and $\frac{f_{\mathrm{1}}m}{m+f_{\mathrm{2}}}$. The resulting misfit, quantified by R² values, is summarized in Fig. C1a. Red and black areas indicate a good fit with R²≥0.9, though an important distinction was made between solutions with f₁≥100 (red) and f₁<100 (black). The difference between the two cases is explained further in Sect. 3.3. Examples of the fit at locations 1 and 2 are shown in Fig. C1b and c respectively. Grey shading in Fig. C1a corresponds to a poor fit (R²<0.9), and the dependency of $\mathrm{\Delta }{U}_{\mathrm{CalvThin}}^m/\left(U_{\mathrm{16}}-U_{\mathrm{96}}\right)$ on m cannot be adequately described by the function $\frac{f_{\mathrm{1}}m}{m+f_{\mathrm{2}}}$. Possible reasons for this discrepancy are discussed in Sect. 3.3.

Appendix D: Dependency of the results on the regularization

The inverse problem of inferring information about the rate factor A and basal slipperiness C from uncertain observations of surface velocity is generally ill-posed. To remedy the ill-posedness of the problem, additional information in the form of a regularization term (Eq. B2) is commonly added to the cost function. The solution of the minimization problem generally depends on the choice of regularization. In the specific case of a Tikhonov regularization, which is used throughout this study, the solution for A and C will depend on the unknown multipliers γ_i,a and γ_i,s, and on the choice of prior information ${\widehat{p}}_i$ in Eq. (B2). One method of choosing an “optimal” value for the multipliers is the L-curve approach presented in Appendix B. However, this is an ad hoc method, and it remains to be shown that results are robust for a range of γ values. Below we discuss the robustness of our results for a range γ_s values. A similar analysis was carried out for a range of γ_a values and priors, but those results did not affect our conclusions and are not shown here.

In the case of the ${\mathcal{E}}_{*}^{\mathrm{3}}$ perturbation experiments, which were designed to simulate the velocity response to a series of prescribed changes in the PIG geometry, we are primarily interested in the γ_s dependency of the relative fluxes in Fig. 3e. In addition to the experiments with default value γ_s=25 000 m, identical perturbation experiments were carried out for γ_s=10 000 m and γ_s=50 000 m. The corresponding changes in flux, presented in Table D1, do not show any significant variability with γ_s, and results presented in Sect. 3.1 can be considered robust, at least across the range of tested γ_s values.

Figure D1 Optimal distribution of m, as in Fig. 6a, for different values of the regularization multiplier: (a) γ_s=10 000 m, (b) γ_s=25 000 m, and (c) γ_s=50 000 m. White dots indicate areas where results for the non-linear regression method were poor, with a R² value smaller than 0.9. Black dots indicate areas where the value of f₁ in the fit is less than 100, indicating that agreement between simulated and observed changes in surface velocity cannot be achieve for finite values of m. The value γ_s=25 000 m was used throughout the main part of this study.

Experiments ${\mathcal{E}}_A^{\mathrm{3}}$ and ${\mathcal{E}}_C^{\mathrm{3}}$ were also repeated for γ_s=10 000 m and γ_s=50 000 m. Maps of A and C (not shown) were compared to the default results for γ_s=25 000 m shown in Fig. 4, and no significant qualitative differences were found.

Perturbation experiments ${\mathcal{E}}_{*}^m$ for a range of sliding-law exponents $\mathrm{1}\le m\le \mathrm{21}$ were repeated for γ_s=10 000 m and γ_s=50 000 m. Following the approach outlined in Sect. 3.3, an optimal spatial distribution of the sliding exponent was computed for each γ_s. Results are presented in Fig. D1 and show a decreasing trend in m_optimal for increasing values of the regularization multiplier γ_s. In particular, the area where no positive, finite solution exists for m_optimal (shaded in black) is reduced in size and eventually disappears for increasing amounts of regularization. However, the spatial distribution of m_optimal is found to be in broad agreement across the considered range of γ_s.

Table D1 Sensitivity of the relative flux changes in the ${\mathcal{E}}_{*}^{\mathrm{3}}$ experiments (see Fig. 3) with respect to the choice of regularization multiplier γ_s. The optimal value, γ_s=25 000 m, used throughout this study was based on the L curve presented in Fig. B1.