Ongoing grounding line retreat and fracturing initiated at the Petermann Glacier ice shelf, Greenland, after 2016

. The Petermann ice shelf is one of the largest in Greenland, buttressing 4 % of the total ice sheet discharge, and is considered dynamically stable. In this study, we use differential synthetic aperture radar interferometry to reconstruct the grounding line migration between 1992 and 2021. Over the last 30 years, we ﬁnd that the grounding line of Petermann retreated 4 km in the western and eastern sectors and 7 km in the central part. The majority of the retreat in the central sector took place between 2017 and 2021, where the glacier receded more than 5 km along a retrograde bed grounded 500 m below sea level. While the central sector sta-bilized on a sill, the eastern ﬂank is sitting on top of a down-sloping bed, which might enhance the glacier retreat in the coming years. This grounding line retreat followed a speedup of the glacier by 15 % in the period 2015–2018. Along with the glacier acceleration, two large fractures formed along ﬂow in 2015, splitting the ice shelf in three sections, with a partially decoupled ﬂow regime. While these series of events followed the warming of the ocean waters by 0.3 ◦ C in Nares Strait, the use of a simple grounding line model suggests that enhanced submarine melting may have been responsible for the recent grounding line migration of Petermann Glacier.


Introduction
The Greenland Ice Sheet (GrIS) is a main contributor to sea level rise at present (Reager et al., 2016;The IMBIE Team, 2020), and its contribution is expected to increase in the coming years (e.g., Choi et al., 2021). While the vast majority of Greenland marineterminating glaciers do not end in extensive floating extensions, ∼ 18 % of Greenland's total ice volume is buttressed by ice shelves in the north Morlighem et al., 2017). It is of prime importance to document the evolution of these ice shelves as their weakening may destabilize this sector and significantly increase the contribution of GrIS to sea level rise (Mouginot et al., 2015). Indeed, enhanced submarine melting (Wilson et al., 2017) and deep rheological damaging (Lhermitte et al., 2020) have the potential to reduce the buttressing effect of glaciers, triggering increased rates of mass losses through dynamic ice discharge.
Petermann Glacier is located on the northwestern part of the GrIS and ends in the second largest floating extension after Nioghalfjerdfjorden glacier, with a length of 70 km in the 1990s and a width of 20 km. Since 1995, this glacier seems to have increased its ice discharge from 10.6 up to 11.7 ± 1.2 Gt yr −1 in 2018 . Recent studies using Worldview digital elevation models (DEMs) re-vealed submarine melt rates near the glacier grounding line that exceed 50 m yr −1 (Wilson et al., 2017). The grounding line is one of the most sensitive regions of a glacier as its migration significantly modulates the ice shelf resistance to flow (Fürst et al., 2016;Reese et al., 2018;Thomas, 1979). Satellite radar interferometry data from the European Remote Sensing satellites (ERS-1 and ERS-2) were used to map the grounding line of Petermann Glacier in the 1990s and detect its migration between 1992 and 2011 (Rignot, 1996Hogg et al., 2016. From these observations, it was concluded that the long-term changes in grounding line locations on Petermann Glacier were indistinguishable from tidal fluctuations, despite the reported thinning of the ice shelf (Wilson et al., 2017;Rückamp et al., 2019;Münchow et al., 2014). Petermann Glacier also lost about half of its floating extension via two large calving events in 2010 and 2012. The glacier showed a limited velocity response to these calving events (Hill et al., 2021;Rückamp et al., 2019), suggesting it was still dynamically stable. Hill et al. (2021) proposed, however, that if future calving occurs within 12 km of the grounding line, the changes in the ice shelf buttressing would be sufficient for the glacier to accelerate substantially.
Since 2011, no extensive grounding line mapping of Petermann Glacier has been conducted (Hogg et al., 2016). Nevertheless, with the launches in 2014 and 2016 of a constellation of two Sentinel-1 (S1) C-band synthetic aperture radar satellites by the European Space Agency (ESA) for the European Copernicus initiative, the opportunity to document the dynamics and the grounding line position of Petermann exists (Milillo et al., 2019;Mohajerani et al., 2021). Following recommendations by the Polar Space Task Group under the umbrella of the World Meteorological Organization, Sentinel-1 has continuously acquired since June 2015 a set of six tracks that cover the entire coast of Greenland, hence providing a systematic and comprehensive coverage at a revisit time of 6 d. This continuous observation record has made it possible to map the grounding lines of ice shelves in a more comprehensive and systematic way. Indeed, S1 allows us to not only map the grounding line again since ERS-1, but on multiple occasions, several times a year, which provides information about the zone of short-term migration or grounding zone (Mohajerani et al., 2021) and its long-term evolution. In addition to S1, the Agenzia Spaziale Italiana's COSMO-SkyMed ® satellite constellation (X-band synthetic aperture radar) has acquired 1 d repeat pass data over Petermann Glacier in 2013 and 2020, hence providing a complementary and high-resolution coverage of Petermann's grounding line.
Here, we present the grounding line history of Petermann Glacier, spanning between 1992 and 2021, or 30 years. We combine satellite radar interferometry and optical imagery to track both the grounding line and the change in ice velocity of the glacier. We also examine concurrent changes in the ocean waters surrounding the glacier. Finally, we conclude on the evolution of the glacier grounding line and its consequences for the future years.

Grounding line positions
We use interferometric synthetic aperture radar (InSAR) data from the ESA's ERS radar satellites acquired in 1992 and 2011 with a 3 d revisit time in 1995/1996 from the ERS-1 and ERS-2 tandem mission and in 2011 from a 3 d revisit time ERS-2 mission ending. Data are downloaded via the ESA Online Dissemination Service as Single Look Complex (SLC) scenes and processed using the GAMMA Software (Werner et al., 2000). We measure the tide-induced vertical motion of ice using a quadruple differential SAR interferometry approach (QDInSAR) . In order to maintain good phase coherence in fast-flowing regions, we first coregister SLC data by estimating the pixel offset using classical speckle tracking technique (Michel and Rignot, 1999;Scheuchl et al., 2016). Interferograms are then assembled by calculating the phase difference between the two coregistered SLCs. To obtain the grounding line position, we combine two interferograms spanning the same time interval, correct for topography, and differentiate them . We use the GIMP v1 DEM time-tagged in 2007 to remove the topographic signal (Howat et al., 2014).
For the time period 2014 to 2021, we use observations collected by S1 with a repeat cycle of 6 to 12 d in Interferometric Wide (IW) swath mode, a Terrain Observation with Progressive Scans SAR (TOPSAR) mode. In order to avoid phase jumps at burst boundaries, interferograms from this sensor are processed using a precise TOPS coregistration methods, which also registers the Doppler history of the slave data . We use the GIMP DEM v2 timetagged in 2014 to correct for the topographic phase for S1 . For both ERS-1/2 and S1, the differential interferograms are formed after geocoding the interferograms in geographic coordinates (north polar stereographic projection) at 25 m posting. As in Rignot et al. (2014), we map the inward limit of detection of vertical motion, where the glacier first lifts off its bed.
To complement the S1 data in 2020 and 2021, we use highresolution observation of the grounding line from COSMO-SkyMed ® (CSK). CSK grounding line measurements use two satellites (CSK2 and CKS4) each with a 16 d repeat cycle. The temporal baseline between CSK2 and CSK4 is 1 d. In order to avoid changes in the horizontal velocity of the glacier, we combine a double-difference interferogram using two 1 d interferograms acquired 16 d apart. CSK acquires in Stripmap mode three consecutive frames to cover the entire glacier and its ice shelf. We stitch all the frames in order to combine a single SLC covering a 40 × 120 km swath. Following the approach described in Milillo et al. (2017) and Brancato et al. (2020), we apply eight looks in both range and azimuth to improve InSAR phase coherence. The final geocoded product has a resolution of 25 × 25 m. As for the S1 case, we use the GIMP v2 DEM to remove the topo- graphic signal. Finally, the grounding line is characterized using the same QDInSAR approach, by differencing two interferograms, in the vicinity of the flexure zone, as for ERS-1/2 and S1 data (see above). We document in Supplement Table S1 the complete list of interferograms used to map the grounding line.

Ice velocity
In addition to the migration of the grounding line, we monitor the evolution of glacier velocity by calculating the surface displacement from four different satellite sensors using images collected between 2013 and 2021. Three of them, ESA's Sentinel-2 (S2) and NASA's Landsat-7 (L7) and Landsat-8 (L8), are optical imagers and one, ESA's S1, is a synthetic aperture radar. We use persistent surface features or speckle to map ice displacements between two consecutive images. We calculate the normalized cross correlations between the reference and slave image chips using repeat cycles shorter than 30 d for L7/8 and S2 and 12 d for S1 Millan et al., 2019Millan et al., , 2022a. Between 1999 and 2012 we supplement our L7 ice velocity record with repeat cycles ranging from 336 to 400 d. For L7/L8, S2, and S1, sub-images of 32 × 32, 32 × 32, and 192 × 48 pixels are used, respectively. We calibrate our displacement maps by taking advantage of the ice velocity products from prior surveys in Mouginot et al. (2017). The final calibrated maps are resampled to 150 m posting in the north polar stereographic projection (EPSG:3413). The time series established is completed by historical measurements made from ERS-1/2, RADARSAT-1, ALOS PALSAR, Envisat ASAR, Landsat-4 to 7, and TerraSAR-X Derkacheva et al., 2020;Joughin et al., 2018). The satellite-derived measurements between 2015 and 2021 are post-processed with locally weighted polynomial regression (Derkacheva et al., 2020) to improve the signal-to-noise ratio and data redundancy. Relative changes in velocity are measured from the 2014-2015 winter mean velocity (January to March).
In order to monitor changes in rates of ice deformation, we derive the evolution of the shear strain rate for 2000 and 2019 using annual ice velocity mosaics . The choice of these two periods is made according to the quality of the surface flow velocity product and to cover time periods preceding and following large events of grounding line migrations (see Sect. 3.1). Strain rates were retrieved us-ing the same methodology as described in Alley et al. (2018), where the shear strain rate is defined aṡ ε shear = ε y −ε x cos α sin α +ε xy (cos 2 α − sin 2 α), (1) whereε x ,ε y , andε xy are the components of the strain rate tensor, calculated from the two components of the ice velocity field (Nye, 1959), and α is the flow angle, defined counterclockwise from the x axis (positive in the x direction) (Alley et al., 2018).

Grounding line evolution
We formed more than 800 quadruple difference interferograms, with 90 % from S1 while the rest was acquired from ERS and CSK, which allows grounding line monitoring at a high temporal frequency. Overall, we manually digitized 161 of the total number of interferograms (see Supplement  Table S1). In summer, decorrelation due to surface melting prevents us from mapping the grounding line on most interferograms. Differential interferograms and grounding lines are shown in Fig. 2 for 1992 , 1996, 2011, 2015, 2016, and every year from 2018 to 2021. For every single year, we display all the grounding line locations, which allows us to document its interannual spatial variability caused, for example, by oceanic tides. The high number of ERS and S1 double-difference interferograms allows us to provide a rough quantification of the grounding zone width of the glacier, which averaged 550 m in 1992 and 630 m in 2015, before the glacier started to retreat substantially (Fig. 3a). Considering an error in grounding line delineation of 110 m (Mohajerani et al., 2021), these estimates are consistent with those of Hogg et al. (2016) with 470 m in 1992-2011. In 1992, the grounding line position of Petermann Glacier was located 70 km from the ice front and remained relatively stable until 2011, meaning that it did not exceed variations induced by tidal modulations (Figs. 2, 3). Between the period 1992-2011 and 2016, the western (T1) and eastern (T4) margins of the grounding line retreated by 3 km. During the same time frame, we observe a retreat of 2 and 1.1 km along transects T2 and T3, respectively (Figs. 2 and 3). In 2015-2017, we were only able to map the central and eastern portion of the grounding line.
Between 2017 and 2019, the grounding line position along the central transect T2 retreated by more than 4 km, while the rest of the grounding line remained relatively stable. Between 2019 and 2021, the retreat continued along T2 and the glacier retreated by another kilometer (Figs. 2 and 3), for a total retreat of 7 km compared to 1992. Within the same timeframe, the eastern (T4) and western (T1) sections of the grounding line retreated by 1 km, for a total retreat 4 km each, since 1992.

Ice shelf shear fracturing
In addition to the grounding line signature in the differential interferograms, we note the formation between 1992 and 1996 of a decorrelation structure in the eastern side of the shelf parallel to the flow direction (Fig. 2b), which seems associated with the development of fractures. The fracture zones are manifest in the interferograms as regions of abrupt discontinuity in phase, indicating that the two sides of the fractures are not flexing exactly at the same pace, hence suggesting deep crevasses rather than surface cracks. In 2011, another parallel fracture zone is visible in the western end of the ice shelf (Fig. 2d-e). Between 2016 and present, these features became more pronounced and now extend on both sides of the grounding line section that retreated (Fig. 2ei). The surface manifestation of these fractures is also visible using optical imagery from L7/8 (Fig. 4c-d), with two distinct and large crevasses developing on the shelf. Using insights from calculated strain rates, it is clearly visible that the development of these fractures corresponds to regions of particularly high shear (Fig. 4a, b). Indeed, while the shear along the eastern fracture averaged −0.02 yr −1 in 2000, it doubled to −0.04 yr −1 in 2019. On the other hand, the fracture on the western side of the grounding line is a region of particularly high shear with a measured shear strain rate of up to −0.08 yr −1 in 2019 (Fig. 4b).

Time series of ice velocity
In Fig. 3b, we show the evolution of Petermann Glacier's velocity at a point located 1 km upstream of its 2021 grounding line for the time period 1992 to 2021. After 2014, we display in Fig. 3e the relative surface velocity change (see Sect. 2.2) across a 75 km long profile (C-C ) that coincides with transect T2 (Fig. 1). We find that between 1992 and 2021 the average speed increased by 14 %, or 150 m yr −1 (Fig. 3b-e). Before 2014, no multi-year change in ice velocity is observed, with winter speed staying near 1050 m yr −1 . A significant speedup started after the summer of 2015, with a winter speed that increased from 1075 m yr −1 to a winter speed of 1200 m yr −1 in 2018. Superimposed on this interannual trend, the frequent velocity measurements available after 2013 reveal a strong summer acceleration of 10 % to 15 % compared to the winter speed (Fig. 3b-e). The summer acceleration occurs uniformly along the ice shelf and propagates 40 km upstream of the grounding line (Fig. 3e). No delay between the acceleration at the grounding line and 40 km upstream is detectable on the weekly time series (Fig. 3e).
We also assembled a time series of surface flow velocity across two gates along the glacier width (Fig. 1). The profile at the grounding displays a pronounced asymmetry in glacier flow, with an ice velocity of 1000 m yr −1 to the west (0 to 10 km in Fig. 3c) versus an ice speed of more than 1300 m yr −1 on the east side of the grounding line (> 10 km of Fig. 3c). The velocity profiles across the ice shelf (Fig. 3d) consistently show two steep transitions (5 and 12 km of Fig. 3d), where the ice velocity abruptly increases from 1100 m yr −1 to more than 1300 m yr −1 over a distance of less than 3 km. These velocity discontinuities are detected more than 30 km downstream of the grounding line (Fig. 4a, b) where the fracture zones with high shear are located as previously described. Figure 2d also shows that these discontinuities evolve over time, with a velocity differential at the fracture location of 100 m yr −1 before 2016 to over 200 m yr −1 afterwards (Fig. 3d).

Ocean thermal forcing
We reconstruct the history of ocean thermal forcing by compiling conductivity, temperature, and depth (CTD) measurements from the Hadley Centre (https://bodc.ac.uk, last access: February 2022) within Petermann fjord and Nares Strait, spanning from 1960 to 2019, combined with CTD from NASA's Ocean Melting Greenland Earth Venture Suborbital mission from 2016-2021 . A warming signal was detected at depth from the 1970s to the 2000s, with a warming of roughly 0.1 • C. An even stronger signal has taken place in the last decade (Fig. 5), with a temperature that increased from 0.1 • C in 2000 to 0.3 • C in 2020. Overall ocean temperature in the fjord at 350-450 m depth (the maximum grounding line depth is ∼ 500 m) is > 0.3 • C warmer in the 2020 than in the 1970s-1980s. This warming signal has been documented elsewhere with a change in temperature of 0.23 • C between 2003and 2009(Washam et al., 2018.

Discussion
Our results indicate that the surface velocity has remained relatively constant between 1986 and 2010 (Fig. 4a). The significant glacier speedup observed after summer 2015 coin-cides with the development of the two along-flow shear fractures within the same time period (Figs. 2e-i and 3a). It is worth noting that the pronounced abrupt transitions in ice flow velocity (Fig. 3c-d) are spatially consistent with the development of these breakup zones along the ice shelf length (Fig. 4), showing a partially decoupled flow regime, meaning that the three different parts of the ice shelf flow at a slightly different pace. The development of these large fractures along the shelf may also be responsible for the difference in flow regime between the floating and the grounded ice, where the abrupt transitions are not present (Fig. 3c and  d). The loss of ice shelf resistance to flow, following the 2012 calving event, may have triggered the later increase in ice flow velocity observed in 2015, causing further ice shelf thinning and consequently a grounding line retreat that happened later on between 2017 and 2018. This result is consistent with the potential reduction in ice shelf buttressing after 2012 as it was calculated by Rückamp et al. (2019). After 2018, the surface velocity of Petermann remained stable, at 1200 m yr −1 (Fig. 3b), while the grounding line continued to retreat rapidly on a retrograde bed slope (Fig. 3a). This is consistent with the hypothesis that changes in ice dynamics may be the cause of the glacier retreat and that the recession of the grounding line did not affect the flow of Petermann Glacier yet.
A comparison of our 30-year-long time series of grounding line migration with BedMachine v3  shows that its retreat is correlated with bedrock geometry. Indeed, the recent grounding line retreat (> 5 km between 2017-2021) occurred across a section of retrograde bed topography that deepens from 470 to 517 m below sea level (Fig. 1b, green line). This pattern of retreat is similar to what has been observed for the Thwaites Glacier in the Amundsen Sea Embayment, West Antarctica, or Humboldt Glacier in North Greenland, where the retreat proceeded faster along topographic depressions of the bedrock (Milillo et al., 2019;Carr et al., 2015). In contrast, the eastern and western portions of the grounding zone have migrated slower than the central part, probably because of slightly prograde bedrock at these locations (Fig. 1b). The eastern portion of the grounding line has remained stable since 2016, on a high rise in the bedrock grounded 490 m below sea level. However, we note that the bedrock deepens over the next 8 km on a retrograde slope down to 540 m depth (Fig. 1c). Since the stability of Petermann Glacier is primarily determined by the level of buttressing at the grounding, which depends on the bedrock geometry and ice rheology, the observed down-sloping bedrock geometry along with the increased ice shelf fracturing may promote the glacier retreat in the coming years (Schoof, 2007;Gudmundsson et al., 2013).
Glacier thinning due to flow acceleration (Thomas et al., 2004;Flament and Rémy, 2012) contributes to grounding line retreat since a thinner glacier reaches flotation sooner. Based on the simple geometrical relationship from Rignot (1998), we calculate the grounding line retreat rate as a function of thickness change, grounding line depth, bed and surface slopes. With a dynamic thinning at about 1 m yr −1 (Smith et al., 2020) on a −0.11 % bed slope and a surface slope of 0.8 % measured along the profile C-C , we calculate an expected retreat of about 130 m yr −1 Hogg et al., 2016), which is below the satellite observations.
A warming ocean signal, however, is observed since about 2000 (Figs. 3-5). A warmer ocean will indeed erode the ice shelf faster by enhancing the submarine melt rate at the grounding line, hence with the potential of promoting glacier retreat An et al., 2021). With an increase in thermal forcing of about 0.3 • C (Fig. 5) and neglecting the role of subglacial runoff, the model of Rignot et al. (2016), which expresses submarine melt as a function of thermal forcing, suggests that the melt would increase by 13 m yr −1 . A thinning of 13 m yr −1 would correspond to a retreat of 1.6 km yr −1 with the same surface and bed slopes consid-ered previously, which is within 200 m yr −1 of the observed grounding line migration in 2018-2021 (Fig. 3b). Although these calculations must be interpreted with caution, as they rely on simple geometric relationships, they seem to indicate that the recent retreat of the grounding line is mainly caused by rising ocean temperatures rather than dynamic thinning.
Finally, the ice flow of Petermann has not yet increased as much as other glaciers in Greenland such as Jakobshavn Isbrae (Motyka et al., 2011) or Zachariae Isstrøm (Mouginot et al., 2015), which are experiencing increased rates of mass losses. However, the 7 km grounding line retreat and the recent fracturing suggest that Petermann Glacier may be entering a phase of destabilization, which has the potential to change the dynamic of this sector of the Greenland Ice Sheet in the coming years.

Conclusions
In this study, we report on recent developments concerning the stability of Petermann Glacier. In 2015, the glacier surface flow velocity increased markedly for the first time in the last few decades before stabilizing in 2018. Markedly after 2016, large fractures formed and split up the ice shelf in three sections with high strain rates and a partially decoupled flowing regime. Since 2017, the grounding line has begun to retreat rapidly, while it was considered stable in the previous 25 years (Hogg et al., 2016). The central section of the grounding line is now more than 7 km upstream of the 1992 position. This recession initiated after 2017 and is proceeding along topographic depressions of the bedrock. Large sections of the grounding line sitting on down-sloping beds could retreat further in the coming years. We posit that these changes are the consequences of the rapid rise in ocean temperature observed in Nares Strait.
Author contributions. RM and JM conceived and designed the research. RM, JM, and EC processed and analyzed data. AD processed ice velocity time series. ER processed and analyzed ocean data. All authors participated in the writing of the manuscript.
Competing interests. The contact author has declared that none of the authors has any competing interests.
Disclaimer. Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.