Earth and Planetary Science Letters The evolving instability of the remnant Larsen B Ice Shelf and its tributary glaciers

Following the 2002 disintegration of the northern and central parts of the Larsen B Ice Shelf, the tributary glaciers of the southern surviving part initially appeared relatively unchanged and hence assumed to be buttressed suﬃciently by the remnant ice shelf. Here, we modify this perception with observations from IceBridge altimetry and InSAR-inferred ice ﬂow speeds. Our analyses show that the surfaces of Leppard and Flask glaciers directly upstream from their grounding lines lowered by 15 to 20 m in the period 2002–2011. The thinning appears to be dynamic as the ﬂow of both glaciers and the remnant ice shelf accelerated in the same period. Flask Glacier started accelerating even before the 2002 disintegration, increasing its ﬂow speed by ∼ 55% between 1997 and 2012. Starbuck Glacier meanwhile did not change much. We hypothesize that the different evolutions of the three glaciers are related to their dissimilar bed topographies and degrees of grounding. We apply numerical modeling and data assimilation that show these changes to be accompanied by a reduction in the buttressing afforded by the remnant ice shelf, a weakening of the shear zones between its ﬂow units and an increase in its fracture. The fast ﬂowing northwestern part of the remnant ice shelf exhibits increasing fragmentation, while the stagnant southeastern part seems to be prone to the formation of large rifts, some of which we show have delimited successive calving events. A large rift only 12 km downstream from the grounding line is currently traversing the stagnant part of the ice shelf, deﬁning the likely front of the next large calving event. We propose that the ﬂow acceleration, ice front retreat and enhanced fracture of the remnant Larsen B Ice Shelf presage its approaching demise.


Introduction
The disintegration of the northern and central parts of the Larsen B Ice shelf (LBIS) over six weeks in 2002 represented an invaluable large-scale natural experiment. It allowed a direct comparison between the responses of the glaciers terminating where the ice shelf no longer existed, with those glaciers still buttressed by the remnant part of the ice shelf that survived in the southern SCAR Inlet (Fig. 1a). The striking flow acceleration of the glaciers that lost buttressing, similar to the acceleration of the glaciers feeding the Larsen A Ice Shelf after its collapse in 1995, left no doubt about the importance of coupling between ice shelves and their tributaries, and the role of ice shelves in regulating the volume of ice that these glaciers discharge into the ocean (Rott et al., 2002;De Angelis and Skvarca, 2003;Rack and Rott, 2004;Rignot et al., 2004;Scambos et al., 2004). Several subsequent studies unsurprisingly focused on the fast changing tributary glaciers in the Larsen A and northern and central Larsen B embayments (e.g., Hulbe et al., 2008;Shuman et al., 2011;Rott et al., 2011;Berthier et al., 2012). In the meantime, earlier reports after the 2002 collapse emphasized the little change exhibited by the surviving southern LBIS remnant and its main tributary glaciers, Leppard, Flask and Starbuck. These reports detected no change in the surface elevation of Flask Glacier between 2003 and 2004 (Scambos et al., 2004), some acceleration in the flow of Flask between 1996 and 2003, and a 15% speeding up in the flow of Flask and Leppard by the end of 2003 (Rignot et al., 2004). Such early observations led to the suggestion that complete removal of an ice shelf might be necessary to initiate tributary glacier acceleration (Rignot et al., 2004;Scambos et al., 2004). More recent observations, however, are finding indications of change. Thus, surface lowering for both Leppard  . Panels 1b and 1c cover the same area as panel 1a. They also show the same latitude, longitude and grounding lines. (d) Location of the remnant LBIS in the northern Antarctic Peninsula. and Flask glaciers was reported for the period 2004-2006, but surface elevation gains or no change were found before and after that period (Shuman et al., 2011), and a negative mass balance for the three main tributary glaciers combined calculated for the period 2001-2010 . Also, progress was made in inferring the ice thickness and bed topography of Flask and Starbuck glaciers (Farinotti et al., 2013(Farinotti et al., , 2014. Here, we focus on the remnant LBIS and its main tributary glaciers, Leppard, Flask and Starbuck. To quantify changes in the surface elevations of the glaciers we analyze ATM laser altimetry from Operation IceBridge and pre-IceBridge campaigns. We find changes in ice flow speeds of these glaciers and the remnant ice shelf from InSAR data. We trace the front positions of the remnant LBIS using SAR mosaics. We finally apply numerical modeling to infer the rheology and backstress fields of the ice shelf to assess the changes in the prevalence of fracture and in buttressing. We present our findings in two parts. In Section 3 we describe the observations of the remnant LBIS and its tributary glaciers, and its modeled rheology and backstress fields. In Section 4 we explore the implications of these results to the question of ice shelf instability, noting that the concurrence of enhanced fracture, front retreat and ice flow acceleration being exhibited by the remnant LBIS is reminiscent of the events preceding the 2002 disintegration (Khazendar et al., 2007).

ATM laser altimetry
We find surface elevation changes from the measurements of the Airborne Topographic Mapper (ATM), which is a laser altimeter that has been flying since 2009 as part of Operation IceBridge but which had flown over Antarctica earlier during the years 2002, 2004 and 2008. Its measurements of surface elevations (Level-2 Icessn Elevation, Slope, and Roughness data) are condensed using the Icessn algorithm that fits a plane to blocks of points selected at regular intervals along track and several across track. The resulting data are presented with an along-track spacing of 50 m. Each nadir measurement is further complemented by 1 or 2 off-nadir measurements on either side, for a total of 3 or 5 points covering a swath of approximately a 100 m (Krabill et al., 2002). Elevations are provided with the values of south-to-north and west-to-east slopes. Our algorithm to find surface elevation changes between repeat ATM observations begins by finding for each point of the first set the spatially closest point in the second. Only points less than 25 m of each other are considered. Surface slopes associated with the identified point in the second set are then used to calculate the surface elevation corresponding to the location of the first point. The elevation difference is then found between the measured elevation of the first point and the elevation corresponding to its location in the second set.
We apply the following corrections to the measured surface elevations (Khazendar et al., 2013). Load tides are applied with the TPXO6.2 load medium-resolution, global inverse model of Egbert and Erofeeva (2002). Further corrections are needed for points on the ice shelf, as being afloat modulates their surface elevations by tidal movements and changes in atmospheric pressure. Hence, ocean tides are calculated with the CATS2008a_opt circum-Antarctic ocean regional model of Padman et al. (2002Padman et al. ( , 2008. The isostatic response of the ice shelf to changes in atmospheric pressure, known as the inverse barometer effect (Padman et al., 2003), is calculated using the value of 1.12 cm of surface elevation lowering for each 1 millibar of increase in atmospheric pressure based on Kwok et al. (2006).
The uncertainty of ATM measurements on grounded ice is assessed to be <9 cm (Krabill et al., 2002). Uncertainty associated with the model correction for load tide is estimated to be ∼0.5 cm (Phillips et al., 1999). On floating ice, the CATS2008a_opt tide model adds a further uncertainty of 7.8 cm (King et al., 2011) and the barometric correction another ∼2 cm, based on an error in surface atmospheric pressure values of ∼2 millibar (Bromwich et al., 2007). Assuming that these uncertainties are random and independent, the combined uncertainty of surface elevation measurements on grounded ice is hence ∼9 cm, and on ice shelves is ∼12 cm. When changes in surface elevation are calculated, propagation of error produces point-to-point uncertainties of ∼13 cm and ∼17 cm, respectively.
A larger source of uncertainty in observing elevation changes is firn compaction due to changes in accumulation, surface melting and temperature, among other factors (Ligtenberg et al., 2011). We constrain the contribution of these surface processes to observed elevation changes by seeking repeat ATM measurements of stagnant and grounded ice areas. One such area in the study region (Figs. 6 and 7) shows that for the years 2009, 2010 and 2011 surface elevation changes were mostly between 0 and 60 cm.
This range overlaps with firn compaction rates of −19 to 12 cm/yr modeled by Scambos et al. (2014) on grounded ice of the northern Antarctic Peninsula.

Velocity data
We derive ice flow velocities from spaceborne interferometric SAR (InSAR) data using speckle tracking (Michel and Rignot, 1999). For the study area, data are available from multiple sensors that cover the period 1997-2012: RADARSAT-1 (C-band;, 2000, ALOS PALSAR (L-band; 2006-2010) and TanDEM-X (X-band; 2012). We select interferometric data with the shortest time interval available and use a speckle tracking technique (Michel and Rignot, 1999;Rignot et al., 2011a) to calculate a two-dimensional velocity field from slant range and azimuth displacements. Tide and load correction are applied where necessary  using the CATS2008a_opt tide model (Padman et al., 2002(Padman et al., , 2008) and the TPXO6.2 load model (Egbert and Erofeeva, 2002). The absolute calibration of the velocity data is aided by the avail-ability of control points of known velocity (in most cases zero) such as mountain peaks, nunataks or domes. Measurement uncertainty largely depends on the sensor used as well as on the number of overlapping tracks and/or cycles available to form a map. Overlapping tracks and the use of multiple cycles, where available, reduce the measurement uncertainty. A more detailed description of our method is provided in Mouginot et al. (2012). To assess the uncertainty of the velocity data we identify a region of little to no motion in the scene and calculate the standard deviation for an area of about 500 km 2 . The one sigma uncertainties estimated for each sensor are: 10 m/yr for RADARSAT-1 in 2000 when a dedicated InSAR campaign was flown and orbits were better maintained, 33 m/yr for all other years, 6 m/yr for ALOS PALSAR and 10 m/yr for TanDEM-X. Some areas of noisier speeds still occur, in particular higher up the glaciers. This is because we compensate for sensor campaigns with lower data correlation by using larger correlation windows when estimating the offsets (about twice the size recommended by Mouginot et al., 2012). The drawback of this approach is that where glaciers are narrower the surrounding region will affect the speed estimate of the glacier and cause additional noise.

Glacier beds, ice-shelf front and grounding line
We use the Level-2 Multichannel Coherent Radar Depth Sounder (MCoRDS; Allen, 2010, 2014) airborne measurements to trace glacier profiles. Data are presented with an along-track resolution of ∼25 m and depth resolution of ∼18 m. Crossover analyses of thickness measurements by a system similar to MCoRDS flying at low altitude (Gogineni et al., 2001), and by MCoRDS itself at high latitude (Li et al., 2013), suggest an uncertainty of ±10 m, when the ice bottom is detected.
The ice fronts are hand digitized using SAR amplitude mosaics assembled for the years 2006 and 2010, viewed at a scale of approximately 1:100,000. At this scale, the uncertainty in ice front position is estimated to be ±1 pixel, or 150 m. In addition to the 2006 and 2010 fronts, the background MODIS image in Fig. 1a shows the front location in March 2003 .
We find grounding line locations by applying differential satellite synthetic-aperture radar interferometry to data from the Earth Remote Sensing Satellite 1 and 2 (ERS-1 and 2) from the years 1995 (October 31; November 1 and 7), 1996 (January 16,17,19,20,22,23;February 23,24,26,27) and 1999 (November 9). Two tandem interferograms, i.e., four acquisitions along the same satellite track, are used for each of four mappings of the grounding lines. The method is described in detail by Rignot et al. (2011b) who estimate the uncertainty to be ±100 m by comparing multiple mappings, multiple instruments and multiple epochs. This uncertainty refers to the spatial accuracy of the grounding line location, which is to be distinguished from cases of grounding line migration over several kilometers as we discuss in Section 4.3 below.

Modeling approach
The analytical and numerical models applied to infer the rheology and backstress of the ice shelf were described and tested in detail in previous work (Borstad et al., 2012(Borstad et al., , 2013. In summary, at a given point of the ice shelf backstress is defined as a scalar measure of resistance to flow relative to an unconfined ice shelf. Backstress is found from an analytical model of ice-shelf creep (Borstad et al., 2013) requiring ice rigidity as one of its inputs. The rheology field of the ice shelf is inferred with an inverse control method (MacAyeal, 1993;Rommelaere and MacAyeal, 1997) implemented in the Ice Sheet System Model (ISSM) (Morlighem et al., 2010;Larour et al., 2012). The method seeks a spatially variable rigidity field that minimizes the misfit between modeled and observed InSAR velocities obtained as described in Section 2.2 from the years 2000, 2006 and 2010 are used in the inverse method because of the near completeness of their spatial coverage. More data sets are used in inferring the rheology and backstress fields. Iceshelf thickness is from the Bedmap2 data set (Fretwell et al., 2013), which for ice shelves is based on the data of Griggs and Bamber (2011). A steady state temperature field for the ice shelf is calculated analytically (Borstad et al., 2012) as a function of surface and basal temperatures and basal melting rates, accounting for vertical advection and diffusion of heat within the ice column (Holland and Jenkins, 1999). Surface temperatures are specified from the Regional Atmospheric Climate Model (RACMO) mean temperature data for the period 1980(van den Broeke and van Lipzig, 2004. A constant basal temperature of −2 • C, the approximate pressure melting temperature, is assumed. Basal melting rates are estimated with the MIT general circulation model (MITgcm) with a three equation thermodynamic representation of the melting and freezing processes in the sub-ice-shelf cavity . The basal melting rates are the same as those used in a study of the Larsen B Ice Shelf in the year 2000 by Borstad et al. (2012). Moderate rates of melting, up to several meters per year, are calculated near the grounding line of the major tributaries of the ice shelf, with low (<1 m/yr) rates of marine ice accumulation calculated downstream of the major promontories of the shelf. Horizontal advection of colder ice from tributary glaciers into the ice shelf was not explicitly accounted for, nor was the possible influence of surface meltwater on the thermal regime of the ice shelf. The temperature field is only used in the calculation of backstress, and this calculation is not very sensitive to the specified ice temperature. For this reason it was considered that a more complicated thermodynamic treatment of the ice shelf thermal regime was beyond the scope of the present study.

Observed flow acceleration and front retreat of the remnant Larsen B Ice Shelf
The most prominent feature of the remnant LBIS is its persistent ice flow acceleration since the year 2002. Flow accelerated mostly in the sections of the ice shelf that are fed by the outflows of Leppard and Flask glaciers (Figs. 1b and 1c). To illustrate, ice flow speed at a point near the 2006 front of the ice shelf (Fig. 1a, red cross) increased from 475 m/yr in year 2003 to 750 m/yr in 2013. Between those two years, the change in ice flow speed was not monotonic, but the general tendency is that of flow acceleration. Meanwhile, the southeastern part of the ice shelf located between the southeastern edge of the outflow of the Leppard Glacier and the Jason Peninsula remained stagnant (Figs. 1b and 1c). The front of the ice shelf has been retreating since 2002 (Fig. 1a). The retreats are mostly through large calving events, the most prominent of which was the calving of iceberg A-54 in early 2006 (Shuman et al., 2011). Another large calving, that our examination of MODIS imagery suggests occurred in late 2007 or early 2008, pushed the southeastern part of the front further back (Fig. 1a).
Repeat ATM coverage of the remnant LBIS has been limited, and surface roughness of the most dynamic parts of the ice shelf made the ATM altimetry change signal too cluttered to draw conclusions about ice-shelf thickness changes there. The findings of Fricker and Padman (2012) from spaceborne radar altimetry, however, show an average surface lowering of 0.19 m/yr for the period 1992-2008. The surface of the stagnant southeastern part of the ice shelf is smoother, and our ATM analysis for the years 2009, 2010 and 2011 shows that there has been no significant change in its surface elevation during that period.

Observed flow acceleration and thinning of Leppard and Flask glaciers and the absence of change of Starbuck Glacier
Leppard and Flask glaciers, the two largest tributaries of the remnant LBIS, have been accelerating and thinning since at least 2004, the earliest year of change estimates being available from both InSAR-inferred velocities and ATM surface elevation data.
The along-flow transects (Figs. 2a and 3a) show that the increases in flow speeds of both glaciers were highest at the grounding line and diminished with distance upstream. Their evolutions, however, were not identical. In the case of Flask Glacier, change is observed already in the year 2000 when its flow speed directly upstream of the grounding zone (at latitude 297.56 in Fig. 3a) was ∼515 m/yr, up from ∼450 m/yr in 1997. The flow continued to accelerate afterwards until reaching ∼700 m/yr in 2012, 55% higher than its magnitude in 1997. The speed transects of Leppard Glacier are more noisy (Fig. 2a) (Figs. 2b and 3b). The two glaciers exhibited similar sur- face lowering spatial patterns and rates between 2002 and 2004, but they differed afterwards. The thinning rate of Leppard Glacier was higher between 2004 and 2008 after which it slowed, while Flask had a higher and more steady thinning rate between 2008 and 2011. By 2011, the cumulative outcome of these changes is a relatively more extensive thinning of Leppard. To illustrate, its surface between altitudes 500 and 600 m lowered by 8.5 to 5.0 m (Figs. 2b and 2c), and some lowering is still detectable at 800 m of altitude. By comparison, the surface of Flask lowered by 6.5 to 2.5 m between altitudes 500 and 600 m, and the lowering reached a maximum altitude of 700 m (Figs. 3b and 3c).
Compared with Leppard and Flask, Starbuck Glacier showed no clear flow speed changes (Fig. 4a) (Figs. 4b and 4c). The magnitude of these surface elevation changes diminished inland until becoming nil at altitude 350 m. The modest magnitude of these changes means that we cannot exclude the possibility that they are in part due to tidally induced flexure or surface processes that are more pronounced at lower altitudes. Above altitude 350 m, the signal is mostly of surface elevation gain of 0 to 1.5 m throughout the period.

Modeled changes in the rheology and backstress fields of the remnant Larsen B Ice Shelf
The ice rigidity fields of the remnant LBIS inferred with the inverse control method are presented in Figs. 5a, 5b and 5c. The spatial pattern of the modeled backstress for the year 2000 is shown in Fig. 5d, while Figs. 5e and 5f illustrate the changes in the backstress field calculated for the years 2006 and 2010 relative to the year 2000.
Values of inferred ice rigidity that are anomalously higher than would be expected for unconfined floating ice at the appropriate bulk temperature of the ice shelf indicate areas where backstress is present. Ice rigidity values that are anomalously lower are associated with fractured or damaged ice (Larour et al., 2005;Vieli et al., 2006Vieli et al., , 2007Khazendar et al., 2007Khazendar et al., , 2009Khazendar et al., , 2011. Those previous investigators found that the areas of anomalously low ice rigidity are the most robust features in the inferred rheology fields, and we focus on those features in the discussion below. The inverse calculations of ice rheology in the year 2000 show similar patterns to those in previous investigations (Vieli et al., 2007;Khazendar et al., 2007), with markedly weaker ice along the shear margins of the ice shelf, and a qualitatively similar pattern of ice rigidity elsewhere. This consistency lends confidence to the patterns of rigidity inferred in the years 2006 and 2010, especially the progressive weakening of the shear margins. On the other hand, in some areas near the grounding line the inferred ice rigidity and backstress values are more susceptible to model artefacts. Examples include the particularly high values of ice rigidity and backstress and their rapid change directly downstream from the Flask and Leppard grounding lines (Fig. 5). Such artefacts are probably the results of ice there not being in hydrostatic equilibrium as assumed in the model and/or ice thickness having abrupt jumps introduced by the interpolation techniques used in producing the thickness data set (Fretwell et al., 2013).

Ice shelf acceleration, front retreat and enhanced fracture
The evolving instability of the remnant LBIS is demonstrated by its flow acceleration, front retreat and increased fracture. Similar changes in LBIS were observed in the years preceding its 2002 disintegration (e.g., Skvarca et al., 2004;Khazendar et al., 2007;Vieli et al., 2007). The availability of more frequent observations in this study allows the changes in those properties to be better followed in time and for the links among them to be explored. Figs. 1b and 1c show that the acceleration of the remnant LBIS is not spatially uniform. The flow speed of the thinner, stagnant ice area between the Jason Peninsula and the outflow of Leppard Glacier did not change much. This is expected given the presence pre-2002 of an already well developed shear zone in the modeled rheology field (Fig. 5a) indicating a degree of mechanical decoupling from the rest of the ice shelf. The slight deceleration between 2006 and 2010 (Fig. 1c) including south of Rift 1, while the ice shelf on the other side of the shear zone continued to accelerate, could be the sign of further decoupling. The areas with the highest acceleration are those fed by the outflows of the Leppard and Flask glaciers. Of the two flow units, the one fed by Flask Glacier accelerated more as illustrated by Fig. 1c and Figs. 2a and 3a. The latter two figures show that the flow speed directly downstream of the grounding line of both glaciers in 1997 is between 320 and 390 m/yr, increasing by 2013 to ∼490-550 m/yr in the case of Leppard, and to ∼650-690 m/yr for Flask. The increased difference in the flow speeds of the two units manifests itself in the inferred rheology field as the extension and increased softening with time of the margin between them (Figs. 5a, 5b and 5c). That such a "suture zone" (Glasser and Scambos, 2008) can be weakened with time by enhanced shearing and fracture resulting from unequal accelerations along its flanks was previously hypothesized (Khazendar et al., 2007;Vieli et al., 2006Vieli et al., , 2007. Our analysis here of the temporal and spatial changes of inferred rheology provides strong evidence in support of the idea. Fig. 5c furthermore shows that within the Flask and Leppard flow units themselves several lines of less rigid ice have developed that were not present in the modeled rheology fields of the years 2000 and 2006 (Figs. 5a and 5b), indicating the susceptibility of this part of the ice shelf to future fragmentation.
The results also illustrate the connection between fracture and front retreat. The MOA image in Fig. 1a shows that the calving that produced the 2006 front occurred along Rift 3, and that the southeastern section of the 2010 front closely coincided with Rift 2. A precursor of the calving event that resulted in the 2010 front is already present in the 2006 modeled rheology field as a zone along Rift 2 of weaker ice (Fig. 5b) that did not exist in the rheology field of the year 2000 (Fig. 5a). Closer to the grounding line at only 12.5 km downstream from it, the rapid widening of Rift 1 is evident by comparing the MOA images of 2003-2004 (Fig. 1a) and [2008][2009] (Fig. 6). There are no repeat ATM measurements of the rift, but one trajectory crosses it nearly perpendicularly in 2009 (Fig. 6). The resulting ATM elevation profile (Fig. 8) reveals that the rift had a width of ∼1 km at the surface and a distinctive pattern of surface uplift on both flanks. The asymmetrical pattern, characterized by higher uplift on the seaward flank of the rift (Fig. 8), was proposed by Fricker et al. (2005) to be a sign of a recently opened rift. Figs. 5a, 5b and 5c unambiguously reveal a zone of weaker ice in the inferred rheology field along the location of Rift 1 progressively extending for more than half the distance from the southeastern margin of the Leppard flow unit towards the Jason Peninsula in less than 11 years. The 2003-2004 MOA image suggests that the rift opened sometime earlier and propagated into the stagnant ice zone before its tip arresting at a point not far from the crack tip location in the 2008-2009 MOA image. The extension of the weaker ice zone in the rheology field between 2006 and 2010 (Figs. 5b and 5c) was therefore mostly due to the flanks of the rift pulling apart rather than the tip of the crack propagating. We propose tentatively that, if the widening of Rift 1 in the stagnant part of the ice shelf is maintained at the rate indicated by Figs. 5a, 5b and 5c, Rift 1 is the likely location of the next large calving of the remnant LBIS, and that this calving would take place before the end of the decade. Furthermore, the same 2009 ATM profile (Fig. 8) also shows the presence of an opening that is ∼3 km wide directly downstream from the grounding line. The ATM data suggest that this feature is not a rift yet, but it too has been growing as can be seen in the MOA images ( Fig. 1 and Fig. 6). The development and enlargement of such an opening at the grounding line raises the intriguing possibility that the ice shelf is being pulled away from its margin at that location, and that would be consistent with the flanks of Rift 1 being pulled apart.
In summary, the stagnant southeastern part of the remnant LBIS appears more prone to the formation of large rifts that initiate either at the margin with Jason Peninsula or at the shear zone with the Leppard flow unit. On the other hand, the fast flowing northwestern part of the ice shelf exhibits extensive fragmentation reminiscent of that revealed by the 2002 collapse of LBIS (e.g., Rack and Rott, 2004).
Another factor possibly contributing to the instability of the remnant LBIS is the absence of marine ice in the observed rifts (e.g., Khazendar et al., 2007;Holland et al., 2009;Kulessa et al., 2014) as Fig. 8 shows Rift 1 and rifts near the front to be open nearly down to sea level. This absence could be due to the rifts having opened relatively recently (Khazendar and Jenkins, 2003), but it could also support speculation by Kulessa et al. (2014) that marine ice layers were relatively thin beneath LBIS or that they were weakened by prolonged basal melting. The values of ice rigidity inferred here in the shear margins of the shelf, however, are lower than could be accounted for even by ice at the melting temperature (Cuffey and Paterson, 2010), therefore mechanical weakening of the ice must be more important than thermal influences associated with the presence of warmer marine ice.

Early glacier acceleration and dynamic thinning
The enhanced fracture and weakening of the remnant LBIS contributes to the reduction of its backstress and the buttressing it affords its tributary glaciers (Figs. 5d, 5e and 5f). Diminished buttressing is consistent with the acceleration and thinning of Leppard and Flask glaciers. A dynamic origin of these changes is supported by the concurrence of the acceleration and thinning. It is also supported by the magnitude of the surface lowering which is too large to be explained by surface processes alone (Fig. 7), and its spatial pattern of decreasing with distance up-glacier (Figs. 2b and 3b).
The extent of ice flow acceleration and glacier surface lowering of Leppard and Flask glaciers is smaller than those of glaciers further north in the Larsen B embayment where the ice shelf disintegrated completely. Yet, the changes are significant, and they started as early as 2003, the year after the collapse of the ice shelf, or before. Moreover, the acceleration and dynamic thinning of Leppard and Flask glaciers occurred despite the continued presence of an ice shelf. Hence, the observations here contradict previous assumptions early after 2002 that the two glaciers continued to be largely unchanged due to the presence of the remnant ice shelf (Rignot et al., 2004;Scambos et al., 2004). Indeed, Fig. 3a shows that the flow of Flask Glacier has already been speeding up from the year 2000, suggesting that the glacier was reacting to the growing instability of LBIS even before its disintegration in 2002.

Differences in glacier evolutions
The earlier acceleration of Flask is one of the differences in the changes exhibited by the two glaciers. ATM surface elevation observations of this region started in 2002, so it is likely that Flask Glacier was already thinning before that year given its earlier acceleration. After 2002, while the two glaciers have similar overall patterns of surface lowering, that of Leppard seems to have reached up-glacier faster (Figs. 2b and 3b). On the other hand, the flow acceleration of Leppard slowed between the years 2006 and 2012 relative to earlier years in the study period (Fig. 2a), while that of Flask was more sustained until 2012 (Fig. 3a). One clue to explain these differences, at least in part, is the bed topography of the two glaciers shown in Figs. 2c and 3c. Both glaciers exhibit low surface slopes and surface elevations that are only 10-20 m higher than freeboard surface elevations would be in case of hydrostatic equilibrium, as can be seen in the case of Flask Glacier in Fig. 9. The two glaciers therefore probably have lightly grounded ice plains (Corr et al., 2001;Brunt et al., 2011) with Flask showing a characteristic break in slope at longitude 297.56 (Fig. 3c). Lightly grounded ice plains are those where the degree of grounding of a  glacier varies over the tidal cycle (Brunt et al., 2011), which seems to be the case of Leppard and Flask in view of the displacements of their inferred grounding lines by several kilometers (Figs. 1a,  2 and 3). On the other hand, the configurations of the multiple grounding lines inferred for the two glaciers are remarkably different. The grounding lines of Leppard are confined within its flow channel, while those of Flask project outwardly into the ice shelf reaching an unnamed island or area of grounded ice (Fig. 1a). This difference between the two glaciers suggests that Flask is more lightly grounded over a wider zone making it more susceptible to stress perturbations in the ice shelf. Furthermore, the bed of the Leppard grounding zone dips inland before rising again, while that of Flask shows no such topographic features. The acceleration and thinning of Leppard Glacier, therefore, could have proceeded faster when the grounding line was retreating over the inland deepening part of the bed, and slowed when reaching the rising, more firmly grounded, part. The Bed topographies of the two glaciers further upstream from the grounding line suggest that the higher susceptibility of Flask is likely to persist into the future if acceleration and thinning continue. The bedrock observations we analyzed here are unavailable upstream of longitude 297.52 (Fig. 3c), but the combination of modeling and observation applied by Farinotti et al. (2013) infers a Flask bedrock that slopes inland reaching as deep as 1200 m below sea level. On the other hand, the observed bedrock profile of Leppard presented here (Fig. 2c) continues to be mostly between 500 and 700 m below sea level beyond the Fig. 9. A cross section of Leppard and Flask glaciers along the 2011 ATM and MCoRDS flight path between the two red dots in Fig. 1a. The hydrostatic floatation line overlaps with the surface near the center of Flask Glacier in the cross section, which is grounded ice, indicating light grounding. Freeboard elevation was calculated at 0.12 of ice thickness (a value obtained from floating ice in Figs. 2c, 3c and 4c) and taking sea surface to be at 14.4 m above the WGS84 ellipsoid, as shown by Fig. 8. grounding zone. Glaciers with deeper bedrocks were found on average to exhibit higher mass losses following ice-shelf loss in the Antarctic Peninsula .
The different evolutions of Flask and Leppard could also be in part influenced by different spatial and temporal patterns of backstress change-a possibility that we cannot ascertain here. Explaining more rigorously the differences between the two glaciers in their reactions to changes in ice-shelf buttressing requires 3D numerical modeling of the coupled glacier-ice shelf system, which we do not attempt in this study but is the subject of future work.
The narrower, thinner and slower Starbuck Glacier exhibits much less change compared with its two larger neighbors. Three factors could have contributed to this lack of reaction to changes in the ice shelf. First, observations depict a glacier that is firmly grounded having grounding lines that overlap (Fig. 1a) and bedrock that slopes upwards upstream from the grounding line (Fig. 4c). The deepening parts of the glacier inferred by Farinotti et al. (2014) lie further inland. Second, the flow speed transects (Fig. 4a) show a high degree of decoupling between the ice flow a few kilometers downstream from the grounding line (at longitude 298.08 in Fig. 4) and the rest of the ice shelf. At the location of the decoupling in Figs. 5a, 5b and 5c, a distinct area of softer ice is shown to develop in the modeled rheology field, most likely as a result of the glacier flowing into a heavily fractured shear zone of the ice shelf. Finally, the relative narrowness of Starbuck Glacier indicates that lateral shear stresses possibly have a larger role in the stress regime than longitudinal stresses (Hulbe et al., 2008) hence making the glacier less sensitive to changes in the degree of buttressing offered by the remnant ice shelf.

Factors favoring the study of glaciological changes in the northern Antarctic Peninsula
Our investigation here focused on the one area in Antarctica, other than the Amundsen Sea sector, where there is a wealth of observational IceBridge and pre-IceBridge data collected over years. The glaciological changes in this area offer a large-scale natural experiment that has some advantages compared with the situation in the Amundsen Sea sector. These include the relatively large number of surveyed tributary glaciers that have been reacting to the same initial perturbation in different manners as we demonstrated here; the main perturbation itself is known (the removal of iceshelf buttressing or its weakening in the case of the remnant LBIS); available observations allow the calculation of pre-perturbation buttressing (e.g., Fig. 5d); and the elapsed time since perturbation occurred (abrupt ice-shelf disintegration) is also known. This is un-like the case of Pine Island or Thwaites where the changes initiated by warmer ocean water and enhanced ice-shelf basal melting have been occurring over an undetermined number of years. Indeed, Dupont and Alley (2005) point to the difficulty of diagnosing accurately the changes in the Amundsen Sea sector in the absence of sufficient information on pre-perturbation buttressing and the timing of buttressing loss.

Conclusions
The evolution of the remnant LBIS offers a rare opportunity to study an ice shelf as it weakens and retreats and the concurrent reactions of its tributary glaciers. Applying a combination of observation and modeling, we demonstrated in this work that the final phase of the demise of LBIS is most likely in progress. The weakening of the remnant ice shelf is manifested by its acceleration, front retreat, enhanced fracture including the rapid widening of a large rift close to the grounding line and possibly the detachment of the stagnant part of the ice shelf from neighboring grounded ice. Large rifts tend to open in the stagnant part of the ice shelf, some of which become the site of future calving events, while fragmentation is more prevalent in the fast flowing part of the shelf fed by Flask and Leppard glaciers. Shear zones between the different flow units that compose the ice shelf have weakened. Concurrently with the diminished buttressing provided by the remnant ice shelf, the main tributary glaciers of Flask and Leppard have accelerated and thinned, exhibiting change even before the 2002 disintegration in the case of Flask Glacier. Flask also appears to be the more susceptible to present and future mass loss due to its less firm grounding and its much deeper bedrock (Farinotti et al., 2013) upstream from the grounding zone compared with Leppard. Given that the two glaciers are already flowing out of balance, it is difficult without further modeling work to predict to what extent the eventual complete removal of the remnant ice shelf would affect their flow. Nonetheless, Flask Glacier alone is thicker than, and already flows at nearly half the speed of, Crane Glacier in 2003 (Rignot et al., 2004), which is among the glaciers with the larger mass losses to the ocean since the 2002 disintegration of the northern and central Larsen B Ice Shelf.