Uneven onset and pace of ice‐dynamical imbalance in the Amundsen Sea Embayment, West Antarctica

We combine measurements acquired by five satellite altimeter missions to obtain an uninterrupted record of ice sheet elevation change over the Amundsen Sea Embayment, West Antarctica, since 1992. Using these data, we examine the onset of surface lowering arising through ice‐dynamical imbalance, and the pace at which it has propagated inland, by tracking elevation changes along glacier flow lines. Surface lowering has spread slowest (<6 km/yr) along the Pope, Smith, and Kohler (PSK) Glaciers, due to their small extent. Pine Island Glacier (PIG) is characterized by a continuous inland spreading of surface lowering, notably fast at rates of 13 to 15 km/yr along tributaries draining the southeastern lobe, possibly due to basal conditions or tributary geometry. Surface lowering on Thwaites Glacier (THG) has been episodic and has spread inland fastest (10 to 12 km/yr) along its central flow lines. The current episodes of surface lowering started approximately 10 years before the first measurements on PSK, around 1990 on PIG, and around 2000 on THG. Ice‐dynamical imbalance across the sector has therefore been uneven during the satellite record.


Introduction
The glaciers of the Amundsen Sea Embayment (ASE) were identified as making high contributions to presentday and predicted future sea-level rise [Alley et al., 2015]. Pine Island Glacier (PIG), Thwaites Glacier (THG), and the smaller Pope, Smith, and Kohler glaciers (PSK) all flow over bedrock lying several hundred meters below sea level and deepening toward the interior of the West Antarctic Ice Sheet. All else being equal, marine glaciers flow faster when their grounding lines are located on deeper bedrock, so that grounding-line migration tends to accelerate over bedrock sloping up toward the ocean [Schoof , 2007]. While buttressing in the ice shelf can work against this tendency [Gudmundsson et al., 2012], it may also serve to support a stable state that will be prone to retreat if the ice shelf thins [Asay-Davis et al., 2016]. Contemporary ice-dynamical changes in the ASE have typically been attributed to the intrusion of warm ocean waters onto the continental shelf and into the ice-shelf cavities, melting the underside of the ice shelves Thoma et al., 2008;Jacobs et al., 2011;Alley et al., 2015]. The buttressing provided by ice shelves has weakened by their thinning, leading to speed-up [Joughin et al., 2003;Mouginot et al., 2014], and subsequent grounding-line retreat [Rignot, 1998;Park et al., 2013;Rignot et al., 2014] and upstream thinning [Shepherd et al., 2001[Shepherd et al., , 2002Flament and Rémy, 2012].
Grounding-line retreat also led to a loss of basal traction and, in turn, to an ongoing dynamical response [Joughin et al., 2010[Joughin et al., , 2014, which has resulted in a loss of ice above flotation growing to 120 Gt/yr in the period between 2010 and 2013 [McMillan et al., 2014]. Modeling studies have indicated that rates of ice discharge from this sector are expected to further grow in future, in part due to ongoing ice dynamics and in part due to projected ocean warming [Joughin et al., 2010[Joughin et al., , 2014Favier et al., 2014;Seroussi et al., 2014;Cornford et al., 2015;Goldberg et al., 2015].
Satellite radar and laser altimetry has been used extensively to assess the thinning of the ASE glaciers and their associated ice-dynamical imbalance [Shepherd et al., 2002;Zwally et al., 2005;Pritchard et al., 2009;Wingham et al., 2009;Helm et al., 2014]. Here we combine observations from the ERS-1 (1991ERS-1 ( -2000, ERS-2 (1995-2011), Envisat (2002, and CryoSat-2 (since 2010) satellite radar altimeter missions with measurements from the ICESat satellite laser altimeter (2003)(2004)(2005)(2006)(2007)(2008)(2009) over grounded ice to obtain a time series of elevation rates over the ASE covering the period 1992 to 2015. The main feature of these observations is the dynamic thinning and associated surface lowering, aspects of which have been reported elsewhere. However, the long altimetric record we assembled allows us to carry out more detailed analysis, because the dynamically induced surface lowering inland-which we refer to as drawdown-has been delayed relative to the onset of thinning near to the grounding line [e.g., Nye, 1963;Payne et al., 2004;Williams et al., 2012;van der Veen, 2013]. To study this, we track threshold values of surface lowering along a series of glacier flow lines and, by fitting a linear model to the propagation, quantify the onset and pace at which drawdown spread during the observational record. Finally, we discuss how differences in the pace of upstream migration of drawdown might have been influenced by the glaciers' individual geometric and glaciological settings.

Methods
We first constructed a series of maps of surface-lowering rate (positive values for a lowering surface) posted on a 10 km × 10 km grid covering the ASE and spanning the period 1992.5-2015 at 6-monthly intervals. For each cell with center x i , y j in polar stereographic coordinates, each time t k , and each satellite mission we selected data located at all x, y, t such that |x − x i | < 5 km, |y − x i | < 5 km, and |t − t i | < 2.5 years, then chose the coefficients of an empirical model which is linear in t and quadratic in x and y to give the least squares best fit to the data [McMillan et al., 2014]. The surface-lowering rate with respect to time at each x i , y j , t k was then simply the mean of the coefficients of t derived for each mission where data exist for at least half of the 5 year window surrounding t k . Note that the earliest data came from ERS-1 Phase-C in April 1992, hence the start of our series at 1992.5. The maximum mismatch of surface-lowering rate between two different satellite missions anywhere and at any time was 8 m/yr, but this value is an outlier as the median of mismatch is 0.09 m/yr, while the 95th percentile is 0.72 m/yr. Several refinements to the basic procedure outlined above were needed. A backscatter correction was applied to account for the effects of temporally correlated fluctuations in echo power and elevation [Wingham et al., 1998;McMillan et al., 2014], while biases in ascending and descending tracks were accounted for in the case of CryoSat-2 [Armitage et al., 2014]. Rates exceeding ±10 m/yr were considered unrealistic and discarded, given that McMillan et al. [2014] found maximum absolute values of 9 m/yr on Smith Glacier (see also the supporting information). Locally weighted scatterplot smoothing in time [Cleveland, 1979] of degree 1 with a 2.5 year window was applied, and the data were smoothed spatially using a Gaussian filter with = 3.5 km. Finally, since we were primarily concerned with grounded ice, we discarded all data seaward of the grounding line [Bindschadler et al., 2011].
We discuss the observed surface-elevation trends in the context of dynamic thickening or rather thinning. In principle, other mechanisms could be responsible: for example, anomalies in surface mass balance [Kuipers Munneke et al., 2015] or subglacial hydrology [Fricker et al., 2007] also lead to changes in surface elevation. Alternatively, changes in the snowpack properties, for example, due to strong melt events or surface mass balance anomalies, affect penetration depth of radar, which may lead to a bias in the apparent surface elevation [Arthern et al., 2001]. However, all of these mechanisms are likely to occur as episodic or periodic events, so that the length and continuity of the available record suggests an ice-dynamical cause.
In order to assess variations in the glaciers' dynamics both between different glaciers and between the tributaries of individual glaciers, we evaluated surface lowering along representative flow lines (Figures 1a and 3). The flow lines were defined on the basis of velocity observations  and run along the distinct tributaries of PIG (labels P1-P7), over the large lateral spread of the THG basin (labels T1-T6, HA for Haynes Glacier), and along the central flow lines of the PSK glaciers (labels PO for Pope, SM for Smith, and K1 and K2 for Kohler Glacier).
For each of the available time slices, we computed surface-lowering rates at 5 km intervals along these flow lines from the gridded data using nearest-neighbor interpolation. As a measure of drawdown migration, we then tracked two threshold values of surface-lowering rate, 0.5 m/yr and 1.0 m/yr, as they moved upstream. For each flow line and threshold value, we collected pairs of distance along the flow line, s, and time, t, for which the observed lowering rate was within 0.1 m/yr of the given threshold value and then fitted a straight line (label F1) to these pairs. The results were often affected by episodic events recorded in the data and thus did not necessarily reflect threshold migration. Therefore, in a second approach, we excluded those (s, t) pairs which were clearly related to episodic events, for example, close to the grounding line in the most recent past  [Bindschadler et al., 2011]; the dashed black lines are the outlines of the drainage basins (for PSK, Rignot et al. [2008], Zwally and Giovinetto [2011], else). Frame A also shows the flow lines P3, T3, and K1 in red, along which surface-lowering rates are shown in Figure 2. The situation of the ASE in Antarctica is shown in the map above frame A. on PIG, or any data before 2000 on most of the THG basin (so that our discussion is limited to the most recent phase of drawdown propagation). We also removed data outside one standard deviation, so that we found an alternative linear fit (label F2). The fitted parameters are referred to as "initiation" (intercept) and "spreading rate" (slope inverse) here. We consider the F2 fitted parameters as our best estimate and the difference between F1 and F2 as the uncertainty inherent in our approach. The quality of F1 and F2 fits is discussed in the supporting information. A straight-line approach does not necessarily represent the actual evolution at any point but rather the mean spreading rate along each flow line, so that differences in spreading rates of the 0.5 m/yr and 1.0 m/yr thresholds can be attributed to nonlinear spreading.
Both of our threshold values are relatively small compared to the maximum values reached at the grounding line and are hence suitable for discriminating between relatively steady surface conditions (before exceeding the threshold) and ice drawdown (after exceeding it). On PIG, where the pattern of inward migration of surface lowering was rather undisturbed by processes other than ice dynamics (see section 4), the fitted spreading rates did not depend on the considered threshold value. Therefore, we consider our choice of threshold values suitable. The grounding line retreated along all ASE glaciers . Our approach of delineating grounded ice with a steady grounding line over the entire 1992-2015 period does not take this into account. If the actual grounding-line position in the early 1990s was more seaward than inferred from data covering several years in the 2000s [Bindschadler et al., 2011], our approach would be biased toward later initiation times. However, the spreading rates would not be affected by this.

Results
PIG saw a steadily growing region of surface lowering throughout the period of observation (Figures 1 and 2a). In the earliest records, the surface close to the PIG grounding line lowered at moderate pace (approximately 1 m/yr over a 6 year period; Figure 1a), while the surface in the interior was either steady or slowly gained altitude. From the late 1990s to 2004, the region of lowering grew in extent but was largely confined to the main trunk where the distinct tributaries coalesce to flow along a deep bedrock trough [see also Shepherd et al., 2001] and only began to spread farther toward the interior afterward (Figure 2a Figure 1a. Distance is taken from the grounding line [Bindschadler et al., 2011]. The flow line on Kohler Glacier is only 150 km long (vertical black line). The data points obtained from upstream threshold tracking are filtered as explained in section 2 so that the subsets indicated by blue (0.5 m/yr threshold) and red (1.0 m/yr) markers remain as a basis for the best linear fit (F2). Initiation time and spreading rates for each threshold value and each flow line are obtained from this best fit as intercept and inverse slope, respectively; see Table 1.
Widespread surface lowering over THG occurred later than over PIG. The seaward portion of the basin exhibited an overall acceleration of surface lowering, from 1 to 2 m/yr close to the grounding line at the start of the satellite era toward an average value of 4 m/yr in the same location over the 2010-2015 period (Figure 1). Surface lowering did not propagate uniformly into the interior, but ceased around 2000 before recommencing around 2004 (Figure 2b).
The strongest signal of surface lowering in any part of our record occurred in the PSK basin (Figures 1 and  2c). Surface lowering rates exceeded 3 m/yr in the early 1990s and 7 m/yr in the most recent years. Surface lowering spread from the grounding line toward the interior in this basin over the 1992-2015 period, too; however, the increase in drawdown-affected area was far less pronounced than in the PIG and THG basins. Table 1 presents initiation dates and spreading rates of thinning-that is, the coefficients from the linear least squares fits-for each of the flow lines. According to that, the area around the PIG grounding line responded to oceanic forcing in the late 1980s to mid-1990s (Table 1), where the detection of initiation depends on the considered threshold value. Drawdown spread into the interior at rates ranging from ∼5 km/yr (flow line P7) to ≥12.9 km/yr (P4 and P5; Table 1 and Figure 3). Comparable changes in the THG basin occurred approximately 10 years later, albeit accompanied by larger uncertainties of the individual initiation dates which results from the occurrence of two separate episodes of surface lowering visible in the altimetry record (pre-and post-2000; Figure 2b and the supporting information). The spreading rates varied strongly between flow lines for the 0.5 m/yr threshold (6.5-17.6 km/yr) but less so for the 1.0 m/yr threshold (5.6-12.8 km/yr). In the latter case, high values above 10 km/yr were concentrated in the central part of the glacier (flow lines T3-T5). The initiation dates associated with the PSK glaciers predated the altimetry records, varied significantly, and were accompanied by large errors. Spreading of drawdown on these glaciers occurred at low rates below 6 km/yr, mostly around 3 km/yr. It should be noted that the applicability of our approach to representing the spatiotemporal evolution by linear fits varied between flow lines for THG and PSK and that the examples in Figures 2b and 2c are among the better representatives (supporting information).

Discussion
The steady evolution of surface lowering on PIG, visible in Figures 1 and 2a, suggests that our linear least squares description of upstream threshold propagation is an adequate representation with the exception of the recent reduced surface-lowering rates near the grounding line. These lower rates might be a consequence of inland diffusion of the dynamical response [Joughin et al., 2010] or could be due to further grounding-line retreat [Park et al., 2013] accompanied by smaller signals in newly floating areas, but our methods are unable to make the distinction. Apart from these recent lower rates at the grounding line, the evolution has been uniform: For a given time, surface-lowering rates decreased with distance from the grounding line; for a given site, they increased with time. This allowed for robust estimates of spreading rates, which fit well to those implicitly given by Payne et al. [2004] for the respective thresholds (7.5 km/yr to 12.5 km/yr) assuming diffusive propagation, which excludes propagation through membrane-stress gradients and which was found appropriate by Scott et al. [2009] and Williams et al. [2012]. Southwestern flow line P7 stood out with the lowest spreading rate for both threshold values (∼5 km/yr). The catchment of this tributary is small, which may explain the marginal drawdown spreading. The southeastern flow lines P4 and P5 exhibited the greatest rates (≥12.9 km/yr)-accompanied, however, by relatively large uncertainties mostly above 3 km/yr. A faster spreading of drawdown in this area was also evident in the most recent spatial pattern (Figure 1d; southward curvature of the area affected by surface lowering at 1.5 m/yr or faster). Earlier studies showed that inverting for subglacial conditions yielded higher basal shear stress relative to ice-flow speed for the southeastern tributaries represented by P4 and P5 compared to the eastern (P3) and northeastern tributaries (P1 and P2) [Joughin et al., 2009;Morlighem et al., 2010;Arthern et al., 2015]. This difference in effective drag, whether due to till properties or hydrology, could serve as an explanation for the faster spreading in the southeast if the drawdown was spread through membrane stress, as greater basal stress should lead to greater longitudinal strains and hence thinning rates in that case, but variation in lateral drag associated with differing tributary width might serve just as well. Note that although most of the flow lines run close to one another along the main trunk for ∼100 km, we do not automatically expect close agreement between the spreading rates in neighboring flow lines because the rates were computed over the entire length, including the clearly separated tributaries.
Surface lowering in the THG basin occurred in two episodes over the satellite altimetry era with a period of abatement around 2000 (Figure 2b). Notably, the ceasing surface lowering stands in apparent contradiction to a continuous retreat due to THG's inherent geometry-driven instability [Schoof , 2007]. Our calculated spreading rates and initiation dates reflect only the most recent episode from 2000 onward. Due to this episodic behavior, the area affected by drawdown in the most recent years (2010)(2011)(2012)(2013)(2014)(2015) was smaller than on PIG in the same period but rather similar to the area affected on PIG between the late 1990s and 2010 (Figure 1). Some localized patches of faster surface lowering might have been related to drainage of subglacial lakes as observed by Smith et al. [2016] at higher resolution. In contrast to PIG and the PSK glaciers, THG is not constrained by a relatively narrow bedrock trough, has a greater lateral extent, and the less confined floating glacier tongue and ice shelf provide less buttressing to the grounded ice, making THG less susceptible to an increase in subshelf melting [Parizek et al., 2013;Nias et al., 2016]. Consequently, if the episodic surface lowering was of ice-dynamical origin, it could have been due to temporally well-defined events of ungrounding and a relatively abrupt response of the fast-flowing sections. Due to this, THG is less well described by the linear empirical model than PIG, which is reflected in higher uncertainties for both initiation dates and spreading rates (Table 1) and the lower goodness of fit (supporting information). It also led to less consistency between the results using the two different thresholds: the eastern flow lines T1-T4 stood out with high spreading rates (≥15 km/yr) in the case of the 0.5 m/yr threshold, while T3-T5 exhibited greatest spreading rates (10.2-12.8 km/yr) for the 1.0 m/yr threshold. Lower spreading rates were evident toward the eastern and western margins of the basin (6 km/yr along the Haynes glacier flow line HA, ∼8 km/yr on the eastern flow lines T1 and T2). This pattern is in accordance with the deeper inland migration of drawdown along the center (Figure 1d), with higher velocities along the central section  and with the decay of longitudinal strain rates and associated elevation rates expected in the lateral direction outward from the center of an ice stream [Raymond, 1996], which should be particularly important for THG due to its wide lateral range. The migration of the higher 1.0 m/yr threshold value potentially better reflects the mean spreading as it may have been less affected by disturbances.
The PSK glaciers, while exhibiting the highest rates of surface lowering near the grounding line, have seen less inland spreading. The glaciers have small catchments and flow along short, narrow troughs in the bedrock topography, at least compared to PIG and THG. The slow spreading of drawdown in the PSK basin may have been a direct consequence of this. The observed grounding-line retreat in this area reported by Rignot et al. [2014], implying an open ice-shelf connection between Dotson and Crosson ice shelves, is obviously not reflected in the assumed steady grounding line, which led to an apparent shift of the maximum rates of surface lowering toward the interior (Figure 2c) in contrast to an expected peak at the grounding line.
A notable difference between glaciers is the variation in the initiation date. THG's current episode of drawdown began in the early 2000s, PIG's in the early 1990s, and the PSK glaciers likely before that. We surmise that the ice-dynamical response to warmer ocean waters has been more persistent on PIG than THG during the observed period. As discussed above, the geometry of THG and its ice tongue and ice shelf may have favored its rather episodic behavior, leading to the relatively small extent of drawdown until the most recent years. Apart from the PSK basin's smaller extent, the pattern at the start of the observational period was similar to PIG in the late 2000s, i.e., at least 15 years after the onset of thinning at the grounding line (Figure 1). Extrapolating the <6 km/yr spreading rates into the past indicates an onset of thinning between the 1960s and the mid-1980s, depending on where the grounding line was at that time. Episodic thinning was evident on THG, as two different episodes occur in the altimetry era. It is entirely possible that a sequence of such episodes might have occurred on PIG or the PSK glaciers before 1992, too. For example, Jenkins et al. [2010] reported evidence of grounding-line retreat of PIG between the 1970s and 1990s; at the same time the slow surface lowering at the start of the altimetry record suggests that the current thinning trend did not start long before 1990 so that any thinning related to earlier retreat episodes should have ceased by then.

Conclusions
We combined observations from five satellite altimetry missions over nearly 25 years to assess dynamical change in the Amundsen Sea Embayment. Pine Island Glacier (PIG), Thwaites Glacier (THG), and the Pope, Smith, and Kohler (PSK) glaciers that flow into the Dotson and Crosson ice shelves, all exhibited a widening extent of surface lowering over the observational period, but there are notable differences both between and within basins. The present episode of drawdown in PIG began at the grounding line around 1990, and its amplitude and extent have grown steadily since then, modulated by a recent reduction in thinning rates around the grounding line. In contrast, drawdown in THG did not begin to spread into the interior until 2000, although a prior episode of thinning near the grounding line ceased during the late 1990s. Surface lowering in the PSK basin must have begun before the altimetry record: simple extrapolation indicates onset before the mid-1980s, but the extent of drawdown has grown less quickly than in PIG or THG. Within basins, the extent of drawdown grew more quickly in the southeastern branches of PIG (∼13 km/yr), perhaps due to variable basal conditions or channel width. The picture in THG is not quite so simple, but the region of fast drawdown grew more quickly along the center of the ice stream than along the margins.
The nonuniform onset and spreading of drawdown during the altimetry era visible in our results can be used to calibrate and test numerical simulations of the ASE glaciers [Goldberg et al., 2015]. Under persistent future oceanic forcing, models indicate that thinning will expand to larger areas in the PIG and THG basins, which would contribute to accelerated ice-mass loss and associated sea-level rise [e.g., Joughin et al., 2010Joughin et al., , 2014Cornford et al., 2015]. However, the nonuniform pace of drawdown spreading gives rise to the possibility that the spreading may be decelerated [Joughin et al., 2010] or even cease in certain regions.