Drivers of Pine Island Glacier retreat from 1996 to 2016

Pine Island Glacier in West Antarctica is among the fastest changing glaciers worldwide. Over the last two 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 commonly attributed to ocean-induced thinning of its floating ice shelf and the associated reduction in buttressing forces. However, other drivers of change such as large-scale calving, 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 5 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 two decades, that (2) simulated changes in ice flow over a viscously deforming bed are only compatible with observations if large and widespread changes 10 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.


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, 2008;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 20 poorly constrained variations in physical ice properties and basal sliding. Evidence from indirect observations have 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 two decades, thinning of the grounded ice (Shepherd et al., 2001;Pritchard et al., 2009;Bamber and Dawson, 2020), intermittent retreat of the grounding line , changes in calving front position (Arndt et al., 2018) and the partial loss of ice shelf integrity (Alley et al., 2019) have all been 25 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 Payne et al., 2004;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), model studies have primarily focused on the important 30 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 changes in ice dynamics. Increased damage in the shear margins of the ice shelf, for example, has 35 been reported by Alley et al. (2019) and is known to reduce the buttressing capacity of an ice shelf (Sun et al., 2017). Moreover, a series of recent calving events has caused 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 emphasised in numerical studies by e.g. Schmeltz et al. (2002); Payne et al. (2004) Brondex et al. (2019). In all cases, some combination of thickness changes, ice softening, a reduction in ice shelf buttressing or variations in basal shear stress were 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 45 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 between years 1996 and 2016, this study brings together the latest observations and modelling techniques. We consider how 50 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 tributaries. Initial observations indicated that the speed-up of PIG was primarily confined to its fast-flowing central trunk Rignot, 2008), though more complex, spatiotemporal 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 55 to changes in the force balance, and it is imperative that numerical ice flow models are capable of reproducing this complex behavior in response to the correct forcing. In general, the driving stress (τ d ), which depends on the ice thickness distribution, is balanced by resistive stresses that include the basal drag (τ b ), side drag through horizontal shear (τ W ), longitudinal resistive forces (τ L ) and back forces by the ice shelf (τ IS ): It is conceivable that each of the terms in Eq. 1 has changed considerably in recent decades, whilst maintaining a balance at all times. This interplay between different changing forces, in combination with the appropriate boundary conditions, underlie the observed dynamical changes of PIG, and form the backbone of any numerical model simulation. In response to changes in the stress balance, modeled changes in ice velocity between time t 0 and t 1 can be expanded as follows: where terms on the right hand side indicate different contributions to the changes in ice flow, caused by variations in calving front position (∆U Calv ), changes in ice thickness (∆U Thin ), changes in ice properties commonly parameterized by a rate factor A (∆U A ), and changes in basal slipperiness C (∆U C ). Note that these contributions are not generally independent due to feedbacks within the system, and that only the total sum, ∆U , can be observed directly. The velocity component related to changes in ice thickness, ∆U Thin , generally consists of two contributions: an instantaneous response due to ice shelf thinning, 70 as investigated by e.g. Gudmundsson et al. (2019), and the induced dynamic loss and redistribution of mass upstream of the grounding line. We will separately assess the impact of both components. Present-day observations of ∆U are generally assumed to be dominated by ∆U Thin , whereas other contributions remain unquantified and are not generally included in model simulations of future ice flow at decadal to centennial timescales. In particular, temporal changes in ice viscosity and basal sliding are ignored such that ∆U A + ∆U C = 0, whereas only a minority of ice flow models include (simple) parameterizations 75 of calving. 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, e.g., τ IS in an attempt to reproduce observed values of ∆U .
In this study we used a regional configuration of the shallow ice stream flow model, Úa (Gudmundsson, 2020), for PIG to estimate the individual components of ∆U in Eq. 2 between t 0 = 1996 and t 1 = 2016. Results enabled us to quantify the relative importance of each driver of change for the contemporary evolution of PIG, and validate the ability of current-80 generation ice flow models to reproduce the complex response of PIG to a range of realistic forcings. It is important to note that results were not derived from transient simulations of glacier flow based on (uncertain) estimates of the initial model state and external forcings. Instead, the diagnostic model response to a series of prescribed changes in ice geometry was analysed, based on the latest observations of ice thinning rates and calving between 1996 and 2016. For each perturbation, changes in the stress balance (Eq. 1) and the associated ice flow response (∆U Calv and ∆U Thin ) were computed. Assuming closure of 85 the velocity budget in Eq. 2 and observed values for ∆U , an estimate for ∆U A + ∆U C was obtained. Knowledge about the magnitude and spatial distribution of each contribution in Eq. 2 allowed us to verify whether common model assumptions such as ∆U Calv = 0 and ∆U A + ∆U C = 0 are indeed justified, at least for contemporary flow conditions.
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, the partitioning between different components about the form of the basal sliding law are likely to precondition the partitioning of ∆U . Indeed, previous model studies have shown that different forms of the sliding law can produce a distinctly different simulated response of PIG to changes in geometry (Joughin et al., 2010;Gillet-Chaulet et al., 2016;Joughin et al., 2019;Brondex et al., 2019). Based on the assumption that ∆U ≈ ∆U Thin , Joughin et al. (2019) showed that a regularized Coulomb law or the plastic limit of a Weertman power-law 95 provide a better fit between modeled and observed changes in surface velocity along the central flowline of PIG, compared to a commonly used viscous Weertman law. 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.

Observed changes of Pine Island Glacier between 1996 and 2016
Our study area and model domain encompasses the 135,000 km 2 grounded catchment  and seaward floating extension of PIG in West Antarctica, as depicted in Fig. 1a. To investigate the physical processes that forced the contemporary  The surface velocity measurements used in this study were taken from the MEaSUREs database (Mouginot et al., 2019a, b Fig. 1b, and we refer to e.g. Rignot et al. (2014) and Gardner et al. (2018) for a more comprehensive description of these observations.
Recent estimates of ice thickness were obtained from the BedMachine Antarctica dataset (Morlighem et al., 2020), which provides both high-resolution surface topography based on the REMA mosaic (Howat et al., 2019)  combination of existing CPOM measurements of thickness change rates (Shepherd et al., 2016) for areas upstream of the 2016 grounding line, and newly analyzed data for the floating ice shelf. Detailed information about the latter can be found in App. A.
The resulting values for ∆H, linearly interpolated across the grounding line and in data sparse areas, are shown in Fig. 1c, and provide the most comprehensive observation-based ice shelf and grounded ice thickness changes for PIG to date.

Experimental design
We discuss the numerical experiments required to obtain an optimal model configuration for

Inverse experiments
We explicitly solved the stress balance in year 1996 (an analogous routine was applied for 2016) by assimilating the known ice 155 thickness (H 96 ), calving front position and surface velocity (U 96 ) in the shallow ice stream (SSA) model Úa (Gudmundsson et al., 2012;Gudmundsson, 2020). This 'data assimilation' or 'inverse' step is commonly adopted in glaciology (see MacAyeal  Glen's law, Eq. 3, relates the strain rates˙ to the deviatoric stress tensor τ . A creep exponent n = 3 was used throughout this study. Equation 4 is known as a Weertman sliding law (Weertman, 1957), and describes a linear, viscous or close-to plastic 165 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, which caused small variations in τ b between cases, but produced an optimal fit between modeled and observed surface velocities in each case. This method differs from other studies, e.g. Joughin et al. (2019), who performed a single inversion for m = 1, and obtained C for different values of m by solving Eq. 4 under the assumption that τ b remains constant. We consider our approach to be more appropriate for this study, 170 as our focus is primarily on an accurate model representation of the surface flow (e.g. Eq. 2). Results for A and C for m = 3 are provided in Appendix B. The outcome of the inverse step is a best estimate for each term in Eq. 1, based on observations of geometry and surface velocity of PIG in year 1996. Analogous results were obtained for 2016.

Geometric perturbation experiments
In the second step we carried out a series of numerical perturbation experiments, starting from the 1996 model configuration,

175
to simulate the impact of observed changes in geometry on the flow of PIG. The rate factor and basal slipperiness were kept fixed to their 1996 values. For each perturbation, the modified force balance (Eq. 1) and corresponding surface velocities, U * , were diagnosed within Úa. Experiments are referred to as E m * with a variable subscript to indicate the type of perturbation and a superscript to specify the value of the sliding exponent m. Experiments were carried out for a range of exponents so we leave m unspecified for now.
Calv . Changes in the calving front location were prescribed to reflect the loss of ice shelf between 1996 and 2016 (see -E m CalvThin . Combined changes in calving front position (as in E m Calv ), and thinning (as in E m Thin ) were prescribed.

190
A schematic overview of the experiments is provided in Fig. 2. While E m Calv allows us to assess the time-integrated impact of calving between 1996 and 2016 (∆U Calv ), and E m ISThin simulates the instantaneous response to total changes in ice thickness between 1996 and 2016 (∆U ISThin ), both experiments ignore the time-dependent, dynamic response of the upstream grounded ice. These separate perturbations make it possible to disentangle the changes in ice shelf buttressing caused by each process, and hence their relative importance for driving the transient evolution of the flow. Dynamic thinning of grounded ice, as well 195 as migration of the grounding line, is included in the experiments E m Thin , which allows us to determine the full response to changes in ice thickness (∆U Thin ). Finally, the experiment E m CalvThin accounts for all geometric perturbation, and provides a spatial distribution of ∆U Thin + U Calv .

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 200 the PIG catchment, i.e. ∆U = ∆U CalvThin . 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, i.e. ∆U A + ∆U C = 0. 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 inverse simulations can provide valuable insights into their evolution. We used the inverse method as described in 205 Sect. 2.2.1 and App. 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.
-E m A . The aim of this experiment is to determine possible changes in the rate factor between 1996 (A 96 ) and 2016 (A 16 ).
A 96 was previously obtained in part 1 (inverse step) of the experimental design. To estimate A 16 , an inverse optimization 210 problem was solved for the 2016 PIG geometry (H 16 ) and velocities (U 16 ), but using a cost function that was minimized with respect to A only. The slipperiness C was kept fixed to its 1996 solution.
-E m C . This experiment is analogous to E m A , 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 215 3.1 Ice dynamics 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 section 2.2.2, perturbations are split between four separate cases: calving (E 3 Calv ), thinning of the ice shelf (E 3 ISThin ), thinning of the ice shelf and grounded ice (E 3 Thin ), and the combined impact of all observed geometrical changes (E 3 CalvThin ). We did not previously specify the value of the sliding exponent, however, here we set m = 3, which is 220 a commonly adopted value in ice flow modeling and describes a viscous, rate-strengthening bed rheology. We will explore results for different values of m in Sect. 3.3. Results for the relative change in surface speed for each of the above perturbations are presented in Fig. 3a-d.
Calving as simulated in E 3 Calv causes changes in flow speed that are predominantly restricted to the 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 225 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 central, fast-flowing 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. The only area with negative relative changes is the western shear margin of the ice shelf, where modeled and observed changes in flow speed have the opposite sign. Extensive damage has caused this margin to migrate and significant interannual variations in flow speed have 230 been reported by Christianson et al. (2016), a process that is not captured by this experiment.
Flux gates provide an alternative, aggregated way to convey the above results. We present flux calculations for two gates perpendicular to the flow within the central part of PIG, as displayed in  Figure 3e shows that calving accounts for 2% and 13% of the observed flux changes through Gate 1 and 2 respectively. This supports the earlier conclusion that the retreat of the PIG calving front between 1996 and 2016 has caused only minor instantaneous changes to the flow upstream of the grounding line.
Thinning of the ice shelf as simulated in experiment E 3 Thin 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 similar perturbation in the buttressing forces. The 240 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 explain the magnitude of 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, E 3 CalvThin , the combined effect of calving and changes in ice thickness were simulated. Modeled versus observed changes in surface speed are shown in Fig. 3d. The spatial pattern is consistent with previous ex- 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 270 experimental design. Indeed, results show that for a viscous bed rheology described by a Weertman sliding law with constant sliding coefficient m = 3, changes in ice geometry alone cannot explain the complex and spatially variable pattern of speed-up over the observational period, i.e. ∆U = ∆U CalvThin . In the remainder of this study, two possible hypotheses are analyzed that enable to close the gap between geometry-induced changes in ice flow and the observed speed-up of PIG. The first hypotheses, which is considered in section 3.2, assumes that bed deformation can indeed be described by a viscous power 275 law with m = 3, but further temporal variations in ice viscosity and/or basal slipperiness are required in addition to changes in geometry: ∆U = ∆U CalvThin + ∆U A + ∆U C . The second, alternative hypotheses, discussed in section 3.3, assumes that internal properties of the ice and bed have not significantly changed between years 1996 and 2016, i.e. ∆U A + ∆U C ≈ 0, but a different physical description of the basal rheology is required instead.

Changes in the rate factor and basal slipperiness between years 1996 and 2016 280
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, exert a second-order control on changes in ice flow. As such, temporal variability in A and C are 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 E 3 A and E 3 C , as outlined in section 2.2.3, is to establish whether this is a 285 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 E 3 CalvThin . Experiment E 3 A assumes that, in addition to changes in geometry, temporal variations in A alone are able to explain the significant increases in flux that were unaccounted for in previous experiments. Alternatively, E 3 C assumes that, in addition to changes in geometry, temporal variations in C alone are able to explain the discrepancy in section 3.1 between the modeled and observed speed- increasingly complex and damaged morphology, as is apparent from satellite images (Alley et al., 2019). Weakening of the ice in these areas accounts for the remaining 50% of observed changes in ice-shelf speed-up that could not previously be explained 300 by calving and ice shelf thinning alone (experiment E 3 CalThin ). Projected changes in A along the flanks of the upstream glacier, on the other hand, are more ambiguous. Values in excess of 10 −7 yr −1 kPa −3 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 https://doi.org/10.5194/tc-2020-160 Preprint. Discussion started: 7 July 2020 c Author(s) 2020. CC BY 4.0 License.
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 wide-spread changes in the thermal and mechanical properties have 305 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 explain the discrepancies between modeled and observed changes in surface speed between 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 310 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 two decades. Further information, such as a timeseries 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. 315 We note that in the E 3 C 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 explain the changes in ice dynamics along the slow-moving flanks of PIG (Fig. 4a), whereas only small changes in C less than 10 −3 yr −1 kPa −3 m were required to explain this behaviour. It is therefore conceivable that, in addition to PIG's evolving geometry, an intricate combination of changes 320 in both the rate factor and basal slipperiness are required to explain the glacier's complex and spatially-diverse patterns of speed-up over the last two decades. It is however not straightforward to disentangle these processes in the current modeling framework.

Evidence for a heterogeneous bed rheology
The relationship between changes in geometry and the dynamic response of a glacier crucially depends on the mechanical 325 properties of the underlying bed and subglacial hydrology. So far, we have assumed that basal sliding can be represented by a viscous power law with spatially uniform stress exponent m = 3 (see Eq. 4). A viscous rheology is particularly suitable for the description of hard-bedded sliding without cavitation, 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 330 et al., 2016) or laboratory experiments (Zoet and Iverson, 2020)  14 https://doi.org/10.5194/tc-2020-160 Preprint. Discussion started: 7 July 2020 c Author(s) 2020. CC BY 4.0 License.
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 Weertman 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 : The transfer amplitude |T U S | contains complicated positive functions f 1 and f 2 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 App. C. Despite the simplifying assumptions that underlie the analytical expression of |T U S | obtained by Gudmundsson (2008), results from our simulations E mi CalvThin , m i ∈ {1, 3, · · · , 21}, indicate that Eq. 5 is also applicable to the more complex 355 setting of PIG. Indeed, as explained in detail in App. C, we found that across a large portion of the PIG catchment, the transfer amplitude |T U S | provides a suitable model to describe the dependency of the relative velocity changes ∆U CalvThin /∆U on m. The parameters f 1 and f 2 were treated as spatially variable fields, and best estimates for f 1 (x) and f 2 (x) were obtained by minimizing the misfit between f1(x)m m+f2(x) and Given the non-linear dependency of ∆U CalvThin /∆U on m with known fields f 1 (x) and f 2 (x), one can derive an 'optimal' 360 spatial distribution of the sliding exponent, m optimal (x), such that ∆U CalvThin /∆U = 100% everywhere, namely By construction, the variable sliding exponent m optimal (x) enables to reproduce 100% of the observed speed-up of PIG in response to calving and ice thickness changes. The results, depicted in Fig. 6a Two interesting properties of the regression model in Eq. 5 are worth noting. Firstly, for m → ∞, the function |T U S | approaches a horizontal asymptote with limit equal to f 1 . As a consequence, the associated solution for m optimal diverges to ∞

370
for locations x where f 1 (x) = 100, and becomes negative where f 1 (x) < 100. In these areas, indicated by black dots in Fig. 6a, no non-negative, finite value of m exists such that ∆U CalvThin (x)/∆U (x) = 100%, 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 either contains significant measurement 375 uncertainties, or approaches the limit ∆U − → 0. In these areas, the non-linear regression was generally found to be poor, with R 2 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 used regression method crucially relies on non-trivial measurements of changes in surface velocity (∆U = 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 1 (x) and f 2 (x) 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 1 (x 1 ), f 1 (x 2 ) = 0 and f 2 (x 1 ), f 2 (x 2 ) = 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 App. D.
In order to demonstrate the improved model response to thinning and calving for a spatially variable sliding exponent 385 m optimal (x), we performed a new inversion with m optimal (x), and subsequently repeated the geometric perturbation experiments E 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 modeled 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 390 Gate 1 and 2 respectively, compared to 28% and 64% for a uniform 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 years 1996 and 2016. The remaining mismatch between the observed and modeled 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 395 regression method in Eq. 5 did not provide a reliable or finite estimate for m optimal , and where Weertman theory of sliding could break down all together (Iverson et al., 1998;Schoof, 2006).

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 analyzed the relative importance of ice shelf 400 thinning, calving and grounding line retreat for the speed-up of Pine Island Glacier between years 1996 and 2016. The detailed comparison between simulated and observed changes in flow speed has provided unprecedented 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 modeled changes in flow were found, and were addressed by either allowing changes in ice viscosity and basal slipperiness, or by varying the mechanical properties of the ice-bed interface. For viscous sliding at the 405 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 two 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 410 tributaries.
Code and data availability. The open-source ice flow model Úa is available from Gudmundsson (2020). All model configurations files specific to this study, as well as model output and plotting routines for each figure are available from DOI TBC. Ice shelf thinning rates are available upon request from FP (fernando.serrano.paolo@jpl.nasa.gov).
Appendix A: Observations of Pine Island Ice Shelf thickness changes between 1996 and 2016 415 We derived a new ice-shelf height time series from measurements acquired by four overlapping ESA satellite radar altimetry (RA) missions: ERS-1 (1991ERS-1 ( -1996, ERS-2 (1995-2003), Envisat (2002), and CryoSat-2 (2010. For this study, we constructed a record of ice-shelf height spanning 20 years , with a temporal sampling of 3 months. We integrated all measurements along the satellite ground tracks and gridded the solution on a 3 by 3 km grid. Our adopted processing steps for RA data are a modification/improvement from Paolo et al. (2016) and Nilsson et al. 420 (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); we 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, and calculating mean values along the satellite reference tracks; we removed height estimates less than 2 m above the Eigen-425 6C4 geoid (Chuter and Bamber, 2015) to account for ice-shelf mask imperfections near the calving front; applied all of the standard corrections to altimeter data over ice shelves (for example, removed gross outliers, and residual heights with respect to mean topography > 15 m; ran an iterative three-sigma filter; minimized the effect of variations in backscatter (Paolo et al., 2016); 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 3 km × 3 km × 3 month cube, for each mission independently. We 430 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 std. 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 435 mitigate floating-grounded mask imperfections, and the limitation of geophysical corrections within the ice-shelf flexural zone.

Appendix B: Model configuration and inverse methodology
The open source ice flow model Úa (Gudmundsson, 2020) uses finite element methods to solve the shallow ice stream equations, commonly referred to as SSA or SSTREAM, 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 440 1.2 km and local mesh refinement down to 500 m in areas with above-average horizontal shear, 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).
The inverse capabilities of Úa follow commonly applied techniques in ice flow modeling to optimize uncertain model parameters, p i , based on prior information,p i , and a range of observations with associated measurement errors. Úa uses an 445 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 ∈ {A, C} were obtained as a solution to the minimization problem d p J with the cost function J defined as the sum of the misfit term I and Tikhonov regularization R: J = I + R, with Eq. 3) and non-linear viscous Weertman sliding (arbitrary m in Eq. 4) 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: T U S = τ d mγ (1 + ψ) + ηH j 2 ψ + k 2 + 4l 2 Hmγ 2 + γηH 2 [l 2 (4 + m) + k 2 (1 + 4m)] + 4H 3 j 4 η 2 , where H is the local ice thickness, α is the local bed slope, ρ is the ice viscosity, τ d = ρgHsinα is the driving stress, η is the 475 effective viscosity and γ = τ 1−m d mC , ψ = ikHcotα and j 2 = k 2 +l 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 T U S 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. Following Gudmundsson (2008), physical quantities can be rescaled to obtain the non-dimensional form of the transfer 480 function. After substitution of the scalings H → 1, η → 1/2, τ d → 1 into Eq. C1 and some reordering, one obtains The resulting transfer amplitude takes the form |T U S | = f1m m+f2 as in Eq. 5, where functions f 1 and f 2 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 485 complicated perturbation composed of a range of wavelengths and amplitudes, and Eq. C2 does not automatically hold. Based Figure C1. (a) Goodness of fit between f 1 m m+f 2 and model simulations ∆U m i CalvThin /∆U, mi ∈ {1, 3, · · · , 21}. Red areas correspond to R 2 ≥ 0.9 and fitting parameter f1 ≥ 100. An example of the fit at location 1 and resulting m optimal (Eq. 6) are shown in panel b. Black areas in (a) correspond to R 2 ≥ 0.9 and fitting parameter f1 < 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 are shown in panel c. Grey areas in (a) correspond to R 2 < 0.9, indicating a poor fit between f 1 m m+f 2 and ∆U m i CalvThin /∆U, mi ∈ {1, 3, · · · , 21}. An example at location 3 is shown in panel d.
on experiments E mi CalvThin , m i ∈ {1, 3, · · · , 21} 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 f1m m+f2 , but more complicated expressions for f 1 and f 2 are required that do not exist in analytical form. A best-estimate for the spatially varying fields f 1 and f 2 was obtained by minimizing the misfit between ∆U mi CalvThin /∆U, m i ∈ {1, 3, · · · , 21} and f1m m+f2 . The resulting misfit, quantified by R 2 values, 490 is summarized in Fig. C1a. Red and black areas indicate a good fit with R 2 ≥ 0.9, though an important distinction was made between solutions with f 1 ≥ 100 (red) and f 1 < 100 (black). The difference between both 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 2 < 0.9) and the dependency of ∆U m CalvThin /∆U on m cannot be adequately described by the function f1m m+f2 . Possible reasons for this discrepancy are discussed in Sect. 3.3. 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 multiplier γ s in Eq. B2. One method to choose an 'optimal' value for γ s is the L-curve approach presented in App. B. However, this is an ad-hoc method and it remains to be shown that results are robust for a range of γ s values.
In case of the perturbation experiments E 3 * , which were designed to simulate the velocity response to a series of prescribed 505 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 = 25000 m, identical perturbation experiments were carried out for γ s = 10000 m and γ s = 50000 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.
Experiments E 3 A and E 3 C were also repeated for γ s = 10000 m and γ s = 50000 m. Maps of A and C (not shown) were 510 compared to the default results for γ s = 25000 m shown in Fig. 4, and no significant qualitative differences were found.
Perturbation experiments E m * for a range of sliding law exponents 1 ≤ m ≤ 21 were repeated for γ s = 10000 m and γ s = 50000 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 exist for m optimal (shaded in black) is reduced in size and 515 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 . Author contributions. JDR and RR designed and initiated the project and prepared the manuscript; FP processed the ice shelf thickness data; JDR performed the model simulations, carried out the analysis and produced the figures; FP and GHG reviewed and edited the paper.
Competing interests. JDR serves as topical editor for The Cryosphere. and (c) γs = 50000 m. White dots indicate areas where results for the non-linear regression method were poor, with a R 2 -value smaller than 0.9. Black dots indicate areas where the value of f1 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 = 25000 m was used throughout the main part of this study. Table D1. Sensitivity of the relative flux changes in the E 3 * experiments (see Fig. 3) with respect to the choice of regularization multiplier γs. The optimal value, γs = 25000 m, used throughout this study was based on the L-curve presented in Fig. B1