Long-Term Monitoring and Change Analysis of Pine Island Ice Shelf Based on Multi-Source Satellite Observations during 1973–2020

Abstract

Pine Island Glacier (PIG) is one of the largest contributors to sea level rise in Antarctica. Continuous thinning and frequent calving imply significant destabilization of Pine Island Glacier Ice Shelf (PIGIS). To understand the mechanism of its accelerated disintegration and its future development, we conducted a long-term monitoring and comprehensive analysis of PIGIS, including ice flow velocity, ice shelf fronts, ocean water temperature, rifts, and surface strain rates, based on multi-source satellite observations during 1973–2020. The results reveal that: (1) ice flow velocities of PIGIS increased from 2.3 km/yr in 1973 to 4.5 km/yr in 2020, with two rapid acceleration periods of 1995–2009 and 2017–2020, and its change was highly correlated to the ocean water temperature variation. (2) At least 13 calving events occurred during 1973–2020, with four unprecedented successive retreats in 2015, 2017, 2018, and 2020. (3) The acceleration of ice shelf rifting and calving may correlate to the destruction of shear margins, while this damage was likely a response to the warming of bottom seawater. The weakening southern shear margin may continue to recede, indicating that the instability of PIGIS will continue.

1. Introduction

West Antarctica exhibited persistent ice loss between 1994 and 2012 [1], and it alone has contributed approximately 6.9 mm to sea level rise since 1979 [2]. Additionally, Pine Island Glacier (PIG) has been the largest contributor to sea level rise in Antarctica in recent years [2,3]. PIG is undergoing basal melting at a greater rate than previously [4,5,6]. Recent findings indicate that atmospheric and oceanic variability, especially the Circumpolar Deep Water (CDW), is a dominant driver of destabilization [7,8]. Targeted observations and ice sheet models suggest that an irreversible retreat of PIG is underway [9].

Ice flow velocity is one of the key indicators that reflect glacier dynamics and estimate the mass balance of the ice sheet. Improved ice motion measurement techniques applied to remote sensing data, including data from optical and synthetic aperture radar (SAR) sensors, have advanced the frequent and accurate observation of surface ice velocity [10,11,12]. In the 1990s, satellite radar interferometry data showed that the grounding line of PIG retreated several kilometers and the ice flow accelerated by nearly 18% in 8 years [13]. By late 2008, the glacier had accelerated to nearly 4 km/yr [14], although it remained relatively constant from 2009 to 2014 and resumed its velocity growth by early 2015 [3,11]. The ice flow of PIG is very fast and still accelerating. Tracking the evolution of ice flow velocity in a long time series can reveal dynamic changes.

The fractured structure on the ice shelf is another indicator of instability. Pine Island Glacier Ice Shelf (PIGIS) has undergone multiple disintegrations as a result of interactive driving mechanisms [8,15]. Satellite radar interferometry observations revealed the retreat of the PIG hinge-line, as well as the thinning of ice by ~3.5 m/yr between 1992 and 1996 [16]. Bindschadler et al. [5] pointed out that the ice flow acceleration can form basal crevasses through a longitudinal extension increase. Jeong et al. [17] suggested that the disintegration of ice mélange around basal crevasses can reduce the support of the ice shelf and cause rifting. Arndt et al. [18] proposed that the formation of rifts after 2006 at PIGIS was favored by back stress fluctuations. Moreover, Lhermitte et al. [19] found that the damage feedback processes heavily influenced the instability of the future ice shelf. However, the mechanism of continuously formed rifts and rapid retreat remains poorly known.

According to Liu et al. [20], subbasin systems with negative mass balance, including the PIG, indicate that net mass loss caused by iceberg calving is comparable in magnitude to net mass loss caused by basal melt, which indicates that the role of iceberg calving may be overlooked in shrinking ice shelves. Considering frequent calving events have occurred at PIGIS in recent decades, we wanted to probe the ice shelf calving at PIG by long-term monitoring. In this study, we first processed available multi-source satellite data since 1973 to obtain the surface ice velocity of PIGIS and analyzed the velocity changes over a long time series. Then, by extracting the ice shelf fronts from 1973 to 2020, we explored the evolution of PIGIS over the past 47 years. We used the latest version of Estimating the Circulation and Climate of the Ocean (ECCO) project to simulate the ocean water temperature near PIG. In addition, we tracked the development of four rifts that led to four recent calving events and linked the surface strain rates to investigate the driving factors. Through a comprehensive analysis of dynamic changes of PIGIS, this study helps to further estimate ice sheet mass balance and provides significant guidance for the quantitative analysis of global sea level rise in the background of global warming.

2. Study Area

PIGIS is located at 75° S, 101° W, flowing into Pine Island Bay, Amundsen Sea, West Antarctica. Figure 1 shows the extent of PIGIS bounded by the grounding line of 2011 [21,22,23]. The whole ice shelf is divided into three parts: the central ice shelf, the northern ice shelf, and the southern ice shelf. The central ice shelf moves much faster than the northern and southern ice shelves [24]. Our research focused on the central ice shelf and its surrounding shear margins.

Hughes [25] proposed that Pine Island Bay may be the weak underbelly of the West Antarctic ice sheet, and it began to retreat in the 1940s and has been one of the fastest-moving glaciers in Antarctica [26]. Multi-source observations during 1947–2000 have revealed the weakening of PIGIS with retreating calving front and grounding line [24]. Ice discharge of PIGIS increased from ~78 Gt/yr in 1973 to ~132 Gt/yr in 2013, by ~69% in 41 years [27]. From 1992 to 2011, the grounding line of PIGIS retreated 31 km at the glacier center compared to 7 to 10 km on its two sides, indicating a faster thinning rate at the center [22].

3. Methodology

3.1. Experimental Data

In this study, the Landsat optical series, European Remote Sensing Satellite ERS-1/2 SAR image series, and Sentinel-1A SAR images were collected for spatial and temporal intermittency. The medium-resolution Landsat images have been widely used owing to their abundant resources and free access. We chose Landsat-1/2/3 Multispectral Scanner (MSS) Visible red images with a 60-m resolution and Landsat-4/5 Thematic Mapper (TM) Visible red images with a 30-m resolution. For the Landsat-7 Enhanced Thematic Mapper Plus (ETM+) and Landsat-8 Operational Land Imager (OLI), we chose panchromatic images with a 15-m resolution. However, this optical satellite cannot obtain valid data during the winter in Antarctica, especially from May to August. Thus, ERS-1/2 synthetic aperture radar (SAR) Level-1 Single Look Complex (SLC) image products and Sentinel-1A Level-1 SLC Interferometric Wide-swath (IW) products were used to cover the shortage of optical images.

Table 1. Summary of the datasets used in this study.

	Landsat Series	Sentinel-1A	ERS-1/2
Time	1973–2020	2014–2020	1992–1999
Repeat cycle	18-day (Landsat-1/2/3) 16-day (Landsat-4/5/7/8)	12-day	35-day
Spatial resolution	60 m (Landsat-1/2/3) 30 m (Landsat-4/5) 15 m (Landsat-7/8)	20 m	25 m
Quantity	61 scenes	6 scenes	9 scenes
Usage	Extraction of ice flow velocity, ice shelf fronts, and rifts.	Extraction of rifts.	Extraction of ice flow velocity, and ice shelf fronts.

shows a summary of the datasets used in this study.

3.2. Pre-Processing

The Landsat series of Earth Observation satellites have continuously acquired images since 1972. Landsat-1/2/3/4/5 Level-1 products were corrected radiometrically and geometrically. The Radarsat Antarctic Mapping Project Version 2 Digital Elevation Model [28] was used as the elevation reference for ortho-rectification. Landsat-7/8 Systematic Terrain Correction (L1GT) products have already been well-calibrated. However, the Enhanced Thematic Mapper (ETM+) sensor of Landsat-7 has acquired images with gaps caused by the Scan Line Corrector failure since June 2003. We selected multi-temporal images for gap filling. Considering the position change caused by ice movement, the pixels at gap areas could be restored by the corresponding area in another image through the Least Squares Matching method. The gap-filled images were used for ice flow velocity calculation and surface fracture extraction.

ERS-1/2 satellites supplemented the data for the 1990s, and Sentinel-1A data were obtained from 2014. Basic processing, including multi-looking, geocoding, and radiometric calibration, was conducted. The final resolutions of Sentinel-1A and ERS-1/2 intensity images were multi-looking processed to the scale of 20 m/pixel and 25 m/pixel, respectively.

Images may not be well aligned, even after being ortho-rectified owing to orientation errors. The stable rock outcrops around were considered as zero-motion areas, which could be used for geometric control and velocity accuracy evaluation. We co-registered all of the ortho-images to the Landsat Image Mosaic of Antarctica (LIMA) [29] through translation transformation. The resulting co-registration error of the imagery in reference to LIMA was approximately 100 m for Landsat-1/2/3, 60 m for Landsat-4, 30 m for Landsat-5, 25 m for ERS-1/2, 20 m for Sentinel-1A, and better than 15 m for Landsat-7/8.

3.3. Ice Flow Velocity Extraction

We extracted ice flow velocity by using the Co-registration of Optically Sensed Images and Correlation (COSI-Corr) software, which allows for precise sub-pixel measurement of ground surface deformation [30]. The frequency correlator of COSI-Corr estimates the phase difference in the Fourier domain and uses a bell-shaped function to avoid edge effects, and is less sensitive to outliers [31]. Fang et al. [32] suggested that phase correlation in COSI-Corr with larger matching windows can achieve the best and most robust glacier motion estimation when applied to optical Landsat images. For fast-moving ice, Landsat image pairs with short repeat cycles (between 32 and 80 days) and ERS-1/2 image pairs with a 35-day repeat cycle are preferred. In the post-processing step, the pixels with signal-to-noise ratio values below 0.98 were removed, and a low-pass median filter (5 × 5) was further applied. We evaluated the velocity extraction accuracy for different satellites by using the stable rock outcrops as check points. We calculated the mean velocity value of the check points. The average error was ±354.2 m/yr for Landsat-3, ±190.6 m/yr for Landsat-4, ±37.8 m/yr for Landsat-5, ±53.3 m/yr for Landsat-7, ±50.1 m/yr for Landsat-8, and ±111.7 m/yr for ERS-1/2. Compared with the ice flow velocity of several kilometers per year, the measurement errors were acceptable.

3.4. Surface Strain Rates Retrieval

Strain rates can help understand ice shelf behavior, such as crevasse and rift development and ice shelf calving. The ice velocity field can be used to calculate surface strain rates. The flow-oriented strain rates [33,34] are given by:

$${\dot{\epsilon }}_{longitudinal}={\dot{\epsilon }}_x{\cos }^2\alpha +2{\dot{\epsilon }}_{xy}\sin \alpha \cos \alpha +{\dot{\epsilon }}_y{\sin }^2\alpha ,$$

(1)

$${\dot{\epsilon }}_{transverse}={\dot{\epsilon }}_x{\sin }^2\alpha -2{\dot{\epsilon }}_{xy}\sin \alpha \cos \alpha +{\dot{\epsilon }}_y{\cos }^2\alpha ,$$

(2)

$${\dot{\epsilon }}_{shear}=\left({\dot{\epsilon }}_y-{\dot{\epsilon }}_x\right)\sin \alpha \cos \alpha +{\dot{\epsilon }}_{xy}\left({\cos }^2\alpha -{\sin }^2\alpha \right),$$

(3)

where, ${\dot{\epsilon }}_{longitudinal}$ is the longitudinal strain rate, which is parallel to the local flow direction; ${\dot{\epsilon }}_{transverse}$ is the transverse strain rate, which is perpendicular to the local flow direction; ${\dot{\epsilon }}_{shear}$ is the shear strain rate; and $\alpha$ is the local flow direction, which is measured counter-clockwise from the image x-axis. Additionally, the strain rates ${\dot{\epsilon }}_x$, ${\dot{\epsilon }}_y$, and ${\dot{\epsilon }}_{xy}$, which center at point ($x_i,y_j$) and correspond to the image axes, are calculated as:

$${\dot{\epsilon }}_x=0.25\left(b+d-a\right)+0.75c,$$

(4)

$${\dot{\epsilon }}_y=0.25\left(b+d-c\right)+0.75a,$$

(5)

$${\dot{\epsilon }}_{xy}=0.5b-0.5d,$$

(6)

where, the strain rates a, b, c, d correspond to the four directions of 0°, 45°, 90°, and 135°, respectively:

$$a=0.5\left(l_{0,1:0,0}+l_{0,-1:0,0}\right),$$

(7)

$$b=0.5\left(l_{0,1:1,0}+l_{0,-1:-1,0}\right),$$

(8)

$$c=0.5\left(l_{-1,0:0,0}+l_{1,0:0,0}\right),$$

(9)

$$d=0.5\left(l_{-1,0:0,1}+l_{0,-1:1,0}\right),$$

(10)

Additionally,

$$l_{k,l:m,n}=\frac{1}{t}\ln \frac{L_f}{L_i},$$

(11)

is the strain rate along the line between ($x_{i+k},y_{j+l}$) and ($x_{i+m},y_{j+n}$), and t is the time interval from the initial length $L_i$ to the final length $L_f$.

3.5. Surface Features Extraction

The surface features of the ice shelf include crevasse, rift, the ice shelf front, ice flowline, etc. [35]. Generally, surface crevasses are small in scale and are often covered by snow or ice debris. They are not easy to identify in remote sensing images. Crevasses can develop into rifts under the influence of various forces, such as atmospheric, oceanic, and glaciological factors [35]. Rifts are obviously visible fracture features, which can penetrate the entire ice shelf and finally cause disintegration. Long-term mapping and analysis of the ice shelf fronts and rifts can be used to explore the ice shelf disintegration mechanism [20,36].

Surface features, including rifts and fronts of PIGIS, were extracted on the Landsat image series, ERS-1/2, and Sentinel-1A intensity images at the same scale. We delineated the position of ice shelf fronts from 1973 to 2020 by manual digitization at 1:20,000 scale. The rifts in two-dimension were also delineated at a 1:20,000 scale. Manual digitizing was chosen following the guidelines of Glasser and Scambos [35] and Holt et al. [37] owing to the challenging nature of automatic feature extraction.

3.6. Ocean Water Temperature Estimation

The ECCO project Version 4 release 4 (V4r4) can derive the latest ocean state estimates from 1992 to 2017 [38,39]. It uses the Massachusetts Institute of Technology general circulation model (MITgcm) to solve the time-varying non-linear inverse problem [40]. We achieved the ECCO-V4r4 to simulate the average ocean water temperature near PIG area (74–76° S, 97–106° W). In consideration of the depth of the CDW in this region [4,8], the water temperature deeper than 300 m was calculated.

4. Results

4.1. Interannual Change of Ice Surface Velocities

The annual ice flow velocity of PIGIS from 1973 to 2020 were obtained. Among them, velocities of 1973–1975, 1986–1988, 1996, 2002, and 2006 were retrieved via the National Snow and Ice Data Center (NSIDC) due to limited data [27,41]. We extracted a longitudinal profile F-F’ along the direction of the ice flow (red line in Figure 1), and took sample points every 200 m from F to F’. Variations of ice flow velocity of PIGIS along F-F’ are shown in Figure 2. For intuitive comparison, we took the same points within 56 km from F to calculate the annual mean velocity (the blue dotted line in Figure 3).

The velocity value doubled over the 47-year study period. From 1973 to 2000, it increased by ~0.6 km/yr but did not exceed 3 km/yr before 2000. However, it accelerated greatly from 2000 to 2009, by ~1.2 km/yr, and exceeded 4 km/yr in 2009. The velocity increment in the 2000s was twice that of the previous 27 years. In 2010, the velocity started to decrease slowly. Milillo et al. [42] pointed out that the slowdown in velocity during 2011–2013 was due to colder ocean conditions. After 2013, the velocity gradually increased, exceeding 4.5 km/yr by 2020. We also noticed that when a rift on the ice shelf was severely ruptured, the velocity value from the rift to the front edge of the ice shelf increased sharply (e.g., profiles of 2012 and 2020 in Figure 2). Overall, the velocity of the whole PIGIS has risen substantially since 1973, which indicates that PIG has been undergoing accelerating mass loss into the ocean.

4.2. Ice Shelf Evolution from 1973 to 2020

We manually traced the ice shelf fronts on the co-registered images at the same scale. The selected images are either at the beginning of every year or at the end of last year to ensure that the time interval is approximately the same. The evolution of PIGIS fronts from 1973 to 2020 are shown in Figure 4a–e. For clarity, we divided into 5 parts from the 1970s to the 2010s and marked the fronts with the same color while the ice shelf continued to develop. We also delineated the ice shelf fronts of the most recent period after the ice shelf collapsed (Figure 4f).

Figure 4a shows the position of fronts in the 1970s. Because of the scarcity of image resources, we obtained only two years’ Landsat images. The red arrow in Figure 4a points out a clear crack at the front of PIGIS. According to Figure 4f, at least one calving event occurred between 1975 and 1982 because the position of the 1982 front is behind the 1975 front. Figure 4b shows the position of fronts in the 1980s. At least one calving event occurred between 1982 and 1987 and one calving event occurred between 1987 and 1989 due to the changes in the position of the ice shelf front. In addition, three cracks (marked by red arrows in Figure 4b) indicate the possible calving of PIGIS. Figure 4c shows the position of fronts in the 1990s. There were at least three calving events in the 1990s, i.e., between 1990 and 1991, between 1992 and 1993, and between 1995 and 1996. Figure 4d shows the position of fronts in the 2000s. At least two calving events occurred during observation, i.e., between 2001 and 2002, and between 2007 and 2008. A big crack shown on the background image, which was acquired on 04 January 2001, supports the calving event during 2001–2002 (marked by a red arrow in Figure 4d). Figure 4e shows the position of fronts in the 2010s. The damage in this decade was extremely serious. A big crack (marked by a red arrow in Figure 4e) propagated almost the entire ice shelf. Shortly thereafter, PIGIS disintegrated in November 2013. There were at least five calving events in the 2010s, i.e., between 2013 and 2014, between 2015 and 2016, between 2017 and 2018, between 2018 and 2019, and between 2019 and 2020. Figure 4f clearly reveals the inland retreat of PIGIS over a long period of time. From 1973–2013, the calving positions were concentrated between 101°15′ W and 101°30′ W, which means that the ice shelf maintained a relatively stable terminus position. Although the subsequent calving events in 2015 and 2017 took place at almost the same location, the northern shear margin receded ~22 km compared with the disintegration position of 2013. Moreover, in November 2018, PIGIS retreated by ~5.5 km, and in February 2020, it retreated by ~7.6 km (the average vertical distance between the disintegration lines). Hence, the instability of PIGIS is obvious, with frequent disintegrations and reductions of the ice shelf.

In addition, we measured the distance from the 2011GL to the ice shelf fronts along the profile F-F’ (the red dotted line in Figure 3). According to the drops of the line chart, at least 13 calving events occurred in PIGIS from 1973 to 2020. It experienced periodic disintegration every 6 years between 1996 and 2013. However, this regularity was short-lived, and it was not uncommon for the disintegration frequency to be 1–2 years. Between 1987 and 1995, although the ice shelf calved four times, the calving position did not change much. However, since 2015, not only did the ice shelf calve four times, but the calving position continued to migrate inland. Compared with the longest distance to the 2011GL (over 90 km) in 2013, the ice shelf front was only ~50 km away from the 2011GL in 2020.

4.3. Rifts Development since 2014

An irreversible retreat of PIG is underway. We noticed that the calving position retreated considerably since the 2015 calving event. To understand the future development of PIGIS, we looked at the development of four rifts that led to recent four major calving events. We wanted to investigate whether this frequent instability will continue and what the drivers are. We traced four rifts, all of which initiated at the center of the central ice shelf, underwent substantial changes, and caused huge collapses. Here, rift R1 was first clearly visible on 11 January 2014 and collapsed between 25 July 2015 and 6 August 2015. Rift R2 was first detected on 18 November 2014 and collapsed between 21 September 2017 and 23 September 2017. Rifts R3 and R4 were both detected on 24 December 2017, while R3 collapsed between 25 October 2018 and 06 November 2018, and R4 collapsed between 05 February 2020 and 17 February 2020. All of the rifts were first detected on the Landsat-8 OLI panchromatic images. We depicted each rift and kept the sampling interval as equal as possible. The images are displayed in a fixed position at the same scale (Figure 5). R1 and R2 almost calved at the same location, about 61 km away from the 2011GL. R4 was closer inland than other rifts. We also compared the time series changes of the length and width of each rift (Figure 6). The distance between the two ends of the rift was measured as the length of the rift. Along profile F-F’ (red line in Figure 1), the distance between the two boundaries of the rift was measured as the width of the rift.

Figure 5a is the evolution of R1. It took about one and a half years from the first detection to the disintegration. Both the width and length grew with time. However, in the last four months before disintegration (17 March 2015–25 July 2015), the length and width changed dramatically (Figure 6a). In particular, the length increased by 14 km and quickly expanded to the shear margins. Figure 5b shows the evolution of R2. It moved seaward, gradually increasing in width and length. In the last two months before disintegration (14 July 2017–21 September 2017), both the length and width changed considerably (Figure 6b). R2 reached the position where R1 collapsed but had not yet expanded to both shear margins. The reason for calving may be the weight of the front ice and the weak margins with ice mélange. Figure 5c is the evolution of R3 and rift R3s. The ice shelf collapsed during observation for less than a year. We found that the rapid growth of R3s contributed substantially. We first detected R3s on 17 September 2018, and it was a small rift opening at the southern shear margin. However, only half a month later, it expanded to the middle of the central ice shelf. During this period, the width and length of R3 grew rapidly (Figure 6c). On 25 October 2018, R3 and R3s connected. Soon, the ice shelf disintegrated. Similarly, R4 collapsed under the influence of rift R4s, which formed to the south of R4 and behind it (Figure 5d). Before R4s was detected, the length and width of R4 had not changed much (Figure 6d). R4s was first visible on 13 December 2018 and grew to the middle of central ice shelf within 14 months. On 5 February 2020, several cracks occurred between R4 and R4s. Soon after, an approximately 350 km² iceberg calved from PIGIS.

The development of these four rifts showed similarities and differences. They all developed slowly in the early stages but accelerated in the later stages. However, unlike R1 and R2, R3 and R4 both encountered other rifts and then changed their original development pattern. This indicates that fracturing shear margins may have a destructive effect.

4.4. Rifts Development since 2014

We investigated the deformation of shear margins by comparing surface strain rates between 2014 and 2020 (Figure 7), which were calculated from the velocity fields of the corresponding years. Areas with a positive longitudinal strain rate indicate that the ice extended along the ice flow. The shear margins and the grounding line areas had higher positive values (Figure 7a,d), indicating that extending deformations, such as heavy crevassing or thinning, were happening. We also noticed that the high longitudinal strain rate at the ice front of the central ice shelf in Figure 7d rightly corresponded to the rifts there. In 2014, the highest longitudinal strain rate occurred at the northern shear margin near the ice shelf front, exceeding 1.0 yr⁻¹ (Figure 7a), which showed evidence of crevasse extension. Thus far, the northern margin has retreated ~22 km. In 2020, a high longitudinal strain rate, exceeding 1.2 yr⁻¹, extended from the southern shear margin edge to inland (Figure 7d), indicating that the weak margin may continue to recede. In addition, the positive longitudinal strain rate gradually increased along slow-flow margins near the 2011GL regions (Figure 7a,d), where Bamber and Dawson [3] reported that thinning rates were at their highest since 2010.

Transverse strain occurs under the action of longitudinal stress. A positive transverse strain rate indicates the divergence of the ice flow, whereas a negative value indicates the convergence of the ice flow. The shear margin areas with negative transverse strain rate indicate additional ice inflowed into the central ice shelf (Figure 7b,e). At the northern shear margin parts, positive transverse strain rate in 2014 (Figure 7b) gradually changed to a negative rate in 2020 (Figure 7e), which showed more ice inflow, and was consistent with the increase in ice flow velocity at that margin area. The shear strain rate maps have both positive and negative values, which refer to ice deformations with right-lateral movement that have a positive shear strain value [34]. A high shear strain rate magnitude illustrates that velocity changed rapidly across the shear zones. In Figure 7c,f, the magnitudes of the shear strain rate along the shear zones were over 0.5 yr⁻¹, whereas the values of the central ice shelf were less than 0.01 yr⁻¹, indicating that lateral drag along the shear margins can form crevasses.

5. Discussion

5.1. Effects of CDW on Ice Flow Velocity and Ice Shelf Disintegration

Early on, scientists pointed out that relatively warm CDW from the southern Pacific Ocean had caused basal melting of PIGIS [43,44]. To investigate the relationship between ice flow velocity, ice shelf disintegration, and ocean water temperature, we simulated the ocean water temperature near PIG from 1992 to 2017 by using the model ECCO-V4r4 (the orange line in Figure 3). Figure 3 shows that changes in ice flow velocity were well linked to the ocean water temperature, and there was a good correlation between them. When the seawater warmed up, it increased the basal melting and thinned the ice shelf [14,45]; thus, the ice flow accelerated (e.g., during 1995–2010). When the seawater temperature decreased, the ice flow slowed down (e.g., during 2010–2014).

Even though the acceleration of the ice flow could speed up the movement of rifts towards the ice shelf front, there was no clear correlation between ice flow velocity and disintegration of the ice shelf. During the 2000s, the ice flow increased rapidly, whereas disintegration regularly followed a 6-year cycle and maintained a relatively stable calving position. However, during 2010–2014, the ice flow decreased by the ocean water temperature, but a calving event still occurred in 2013. Nowadays, the disintegration of PIGIS no longer follows a regular pattern. The seawater has kept warming up since 2014, and four successive calving events in 2015, 2017, 2018, and 2020 occurred. Jeong et al. [17] have pointed out that a combination of ice thinning, basal melting, and increased deviatoric stresses can accelerate rifting. Thus, we attribute ocean force to an important driving factor of the disintegration of PIGIS.

5.2. The Mechanism of Continuously Formed Rifts and Rapid Retreat

The development of four rifts are compared and discussed in Section 4.3. We found an important clue that fracturing shear margins contributed destructive effect on the rapid disintegration of the ice shelf. However, this damage is irreversible. Surface strain rates shown in Figure 7 revealed the changes in the shear margins of PIGIS between 2014 and 2020. The severely damaged northern shear margin edge retreated to the 2011GL after the collapse of R1 and R2. Subsequently, the southern shear margin edge was fractured more prominently than before, contributing to a destructive effect, and even accelerating the calving of R3 and R4 by forming new rifts and extending to the central ice shelf. The shortening of the disintegration period may be related to the damage of the shear margins on both sides, while this damage was a response to the warming of the bottom seawater. Ocean-induced basal melting can thin the ice shelf and therefore provide a potential positive feedback between fracturing shear margins and acceleration of retreat [36,46].

Typical ice shelf calving happened as marginal rift formed at the heavily fractured shear margins and propagated across the ice shelf [47,48,49]. However, recent rifts shown in Figure 5 initiated from the center of the ice shelf and propagated toward the sides. If the central rift got the assistance of marginal rift, the frequency of ice shelf disintegration could greatly shorten. When we look at the surface strain rates in 2020 (Figure 7), especially in the southern part, it shows a clear inland increase. Bamber and Dawson [3] pointed out that the highest thinning rates in PIG are now along the shear margins. Owing to the continuous action of CDW [9] and the accelerated ice output of PIG, the future of PIGIS is grim.

6. Conclusions

We combined Landsat 1–8 series, ERS 1–2 series, and sentinel-1A images to study PIGIS from 1973 to 2020. Time series of annual ice flow velocity were calculated. We extracted long-term ice shelf fronts to analyze the ice shelf disintegration and the frequency of disintegration. Ocean water temperature near PIG area was simulated by the latest ECCO-V4r4 model. In addition, we traced the development of four rifts that led to four recent calving events and linked the surface strain rates to further investigate the factors influencing rifts’ development and the ice shelf disintegration. We drew the following major conclusions: (i) The ice flow velocity accelerated with the warming of the ocean water. It increased from ~2.3 km/yr to ~4.5 km/yr by 2020 over the 47-year study period. (ii) PIGIS has been in nearly 50 years of instability, and at least 13 calving events occurred at PIGIS from 1973 to 2020. Especially since 2015, continuous retreats in 2015, 2017, 2018, and 2020 were unprecedented. Compared with the furthest position of 2013, the ice shelf had retreated inland over 40 km. Moreover, (iii) an investigation of the development of four rifts, as well as the comparison of surface strain rates between 2014 and 2020, showed that the fractured shear margins accelerated rifting and calving. Our results suggest that the shortening of the disintegration period may correlate to the damage of the shear margins on both sides, while this damage is likely a response to the warming of the bottom seawater. Furthermore, the weakening southern shear margin with high longitudinal strain rates exceeding 1.2 yr⁻¹ extended from the southern edge to inland, indicating that the instability of PIGIS will continue.