Figure 1 . Study area. a . Larsen Ice Shelf (A, B, C, and D) from north to south. The white area is the ice shelf extent in 2015. The color-coded area shows the surface topography of the grounded region. We divided Larsen C into six flow units (C1, C2, C3, C4, C5, and C6). b . The ice front locations in 1963, 2000, 2010, and 2017 with the ice shelf surface elevations from REMA DEM (Howat et al., 2019), twelve flow inlets, and suture zones flowing from peninsulas and islands.
3 Data and Methods
3.1 Mapping ice front, rifts, and flow velocities
We used multi-source satellite images acquired between 1963 and 2020 to map the ice front, identify rifts, and derive flow velocities in different time periods. The images include declassified intelligence satellite photographs (DISP) acquired by the ARGON mission, optical images acquired by the Landsat-4/5 TM, Landsat-7 ETM+, and Landsat-8 OLI sensors, and synthetic aperture radar (SAR) images acquired by the ERS-1, ERS-2, Radarsat-1, Envisat, and ALOS satellites. We performed image orthorectification, image co-registration, noise reduction, and image enhancement on the above set of images (Text S1).
We mapped the ice front locations of Larsen C for the period 1963 to 2020 (Figure S1) and delineated two rifts in the vicinity of Gipps Ice Rise that were closely related to the 2017 calving event. We derived flow velocities over the ice shelf for multiple periods using an image-matching-based feature tracking method (e.g., Heid & Kääb, 2012; Liu et al., 2012; Scambos et al., 1992). The feature tracking method tracks the spatial displacements of moving and identical surface features from two sequential images acquired at different times. The derived flow velocity represents the average velocity between the acquisition dates of sequential images. The image-matching algorithm automatically searches for matching points between two images by maximizing the cross-correlation coefficient between conjugate areas in the two images. The cross-correlation-based method can achieve an accuracy of better than 0.5 pixel (Scambos et al., 1992). For our study, the image space was examined every 300 m to search for matching points, with a matching window of 64x64 pixels. Automated directional and local statistical filters (Liu et al., 2012) plus human inspection were employed to remove erroneous matches and verify the matching results. The accuracy of the derived velocity fields is determined by the image co-registration accuracy, image matching accuracy, and time separation of sequential images, which was estimated to be 3.7~25 m/year, depending on the image pairs used (Table S1).
3.2 Ice shelf modeling experiments
We used the Ice-sheet and Sea-level System Model (ISSM) (Larour et al., 2012) to conduct the ice shelf modeling experiments. We used the shallow shelf approximation (SSA) model (MacAyeal 1989; Weis et al. 1999; Morland 1987) to examine how ice shelf flow velocities and stress fields vary with front retreat and rift development (Text S2). For the model configuration, we used the anisotropic mesh adaptation method to generate the numerical mesh net nodes. The boundary condition was specified as kinematic for the interior ice shelf and dynamic for the ice shelf front. The flow law exponent n was set to 3. The ice rigidity \(\overset{\overline{}}{B}\) was initialized as 1.5×108 Pa s1/3. For the inversion of \(\overset{\overline{}}{B}\), we used the M1QN3 optimization algorithm in ISSM (Morlighem et al., 2015) by minimizing a cost function defined as the weighted sum of the difference between the modeled and observed flow velocity fields and the gradient of the optimized\(\overset{\overline{}}{B}\) at each iteration. We further calculated the deviatoric stresses, including longitudinal deviatoric stresses and principal stresses (Text S1). We also calculated the backstress across the ice shelf. Backstress describes the stress transmitted upstream that opposes spreading of an ice shelf (Cuffey and Paterson, 2010), arising from side drag by lateral walls and slower flowing ice, and from basal resistance from grounded ice rises or ice rumples. These stress fields are tightly related to rift growth and ice shelf instability development.
4 Observations derived from satellite data
4.1 Retreat history of Larsen C during 1963–2020
During the period from August 1963 to March 2020 (Figure S1), the area of Larsen C decreased from 55,511 km2 to 43,182 km2, a reduction of ~22%. In the six decades, Larsen C underwent two large retreat cycles, during 1976–1986, and during 2004–2017. Each cycle began with northern retreat near the Bawden Ice Rise (C1–C3 units), followed by southern retreat near the Gipps Ice Rise (C4–C6 units). For both cycles, the time interval between the northern and southern retreat was approximately ten years. In the first cycle, the northern retreat occurred at the C2–C3 units in 1976 due to a collision between an iceberg and the ice shelf (Ferrigno & Gould, 1987); afterwards, two large icebergs calved off from the C4–C6 units in 1986 (Skvarca, 1994), causing a retreat of ~100 km in the south. The ice shelf slowly advanced between 1986 and 2004. In the second cycle, the northern front detached from the C2–C3 units in 2004/2005, with an area loss of ~1500 km2. The major retreat of the southern front happened on 12 July 2017, when a ~5800-km2 tabular iceberg detached from the C3–C6 units.
4.2 Rifting related to the 2017 calving event
The 2017 calving event resulted from the propagation of two rifts (labeled R1 and R2, Figure 2). Long and transverse fractures (Figure 2a) tend to form between Kenyon Peninsula and Gipps Ice Rise due to a stress condition change, where ice shelf flow becomes less constrained by the lateral shear margins after passing by Kenyon Peninsula. The majority of the fractures terminate near the Joerg suture zone (Figure 2a). Suture zones can inhibit fracture growth because of the material properties of marine ice underneath (Kulessa et al., 2014, 2019; McGrath et al., 2014). The rift R1, a transverse rift propagating perpendicular to the flow direction, developed from a preexisting fracture (L1 in Figure 2a) that is identifiable on the 1963 DISP imagery. The lateral extent of R1 changed little until 2008, when it began to propagate towards R2. By 2010, the southern rift tip of R1 had traveled ~14 km from its 2006 location (from A’ to A, Figure 2c), to a point less than 1 km from R2. At the northern end, R1 developed a new rift trajectory from point B’ to B (Figure 2c), and subsequently propagated to the north along this trajectory. R2 was an oblique rift that formed between Gipps Ice Rise and a transverse fracture (L2 in Figure 2a). R2 is not present on the 1963 DISP imagery. Its formation was very likely due to the front geometry change caused by the 1986 retreat. R2 widened over time, while its rift tips exhibited little change until 2013. Then, it began developing towards the south, and the rift tip almost reached Gipps Ice Rise (Figure 2e). Simultaneously, R1 penetrated through the Joerg suture zone and started rapid propagation to the north (Figure 2d). R1 cut through the C5 unit in 2014 and had propagated across the Tonkin suture zone by January 2016 (Figure 2d). By March 2017, the length of R1 dramatically increased by ~65 km, and its northern rift tip (‘B’) reached the Cole Peninsula suture zone (connecting the C2 and C3 units). R1 continued penetrating through the Cole Peninsula suture zone, but did not extend into the C2 unit. Instead, it was developing along the C2 unit boundary towards the ice front, resulting in the complete detachment of the downstream portion in July 2017.
The propagation of R1 appears to have had a strong association with the development of R2. The extension of R2 to Gipps Ice Rise was immediately followed by the rapid growth of R1 towards the north. This suggests that rifting can promote the growth of another rift by weakening the lateral shear margins near an ice rise. Ice rises are locally grounded topographic highs in ice shelves. They often act as pinning points, buttress upstream ice flows and play an active role in ice shelf stability (Favier & Pattyn, 2015; Matsuoka et al., 2015). Weakening around ice rises could reduce their buttressing effect, leading to enhanced longitudinal stretching and creating favorable conditions for rift propagation. We further explore this hypothesis through modeling experiments in Section 5.