Automatic delineation of cracks with Sentinel-1 interferometry for monitoring ice shelf damage and calving

. Monitoring the evolution of ice shelf damage such as crevasses and rifts is important for a better understanding of the mechanisms controlling the breakup of ice shelves and for improving predictions about iceberg calving and ice shelf disintegration. Nowadays, the previously existing observational gap has been reduced by the Copernicus Sentinel-1 synthetic aperture radar (SAR) mission that provides a continuous coverage of the Antarctic margins with a 6 or 12 d repeat period. The unprecedented coverage and temporal sampling enables, for the ﬁrst time, a year-round systematic monitoring of ice shelf fracturing and iceberg calving, as well as the detection of precursor signs of calving events. In this paper, a novel method based on SAR inter-ferometry is presented for an automatic detection and delineation of active cracks on ice shelves. Propagating cracks cause phase discontinuities that are extracted automatically by applying a Canny edge detection procedure to the spatial phase gradient derived from a SAR interferogram. The potential of the proposed method is demonstrated in the case of Brunt Ice Shelf, Antarctica, using a stack of 6 d repeat-pass Sentinel-1 interferograms acquired between September 2020 and March 2021. The full life cycle of the North Rift is monitored, including the rift detection, its propagation at rates varying between 0.25 and 1.30 km d − 1 , and the ﬁnal calving


Introduction 25
Because of their buttressing effect that regulates the upstream flow of the grounded ice sheet, ice shelves play a key role in the mass balance of the Antarctic ice sheet. Ice shelf calving, especially for ice shelves that originate from large tributary glaciers, constitutes one of the main contributions to the mass loss in Antarctica (the IMBIE team, 2018; Rignot et al., 2019).
However, predictions about this contribution remain uncertain due to the poor understanding of the mechanisms controlling https://doi.org/10.5194/tc-2021-296 Preprint. Discussion started: 12 October 2021 c Author(s) 2021. CC BY 4.0 License.
ice shelf break up and the difficulty to model them (Sun et al., 2017;Cook et al., 2018), partly due to a lack of 30 comprehensive observational data. Ice shelves can be structurally weakened by processes such as ice shelf thinning, that leads to ice flow speedup and increased shearing, or such as ice shelf retreat that results into unpinning. These processes may trigger a feedback response, thereby enhancing damages such as deeply crevassed areas and open fractures, increasing the ice velocity gradient and further weakening the ice shelf structure (Lhermitte et al., 2020). The complex response of an ice shelf to rifting, the difficulty to predict ice shelf disintegration and the resulting uncertainties in mass balance models 35 highlight the need for systematic monitoring of the damage evolution (Pattyn et al., 2017).
On-site measurements of ice shelves and active rifts from e.g. ground penetrating radar, time-lapse camera or GPS provide valuable insights for monitoring damages and better understanding the mechanisms leading to rift propagation (Banwell et al., 2017;King et al., 2018;De Rydt et al., 2019). However, field missions are expensive, necessitate heavy logistics and only focus on a specific area (usually close to a base station) for limited periods in time. Therefore, despite their 40 unquestionable value, they provide no feasible solution for continuous long-term and large-scale monitoring of ice shelf rifting systems.
Nowadays, the majority of the Antarctic ice shelves is routinely monitored with optical and radar satellites, providing dense image time series that enable the continuous observation of fracture opening, propagation, widening and iceberg calving in near real-time. For damage monitoring, Synthetic Aperture Radar (SAR) sensors constitute a good alternative to optical 45 satellite imaging thanks to their all-day/all-weather/year-round observing capability. Furthermore, SAR signal penetrates through dry snow that can fill and obscure cracks in optical images. In particular, the acquisition strategy of Sentinel-1 provides a continuous coverage of almost the entire ice sheet margin of Antarctica with 6-and 12-day repeat intervals, which enables the systematic surveillance of ice shelf fracturing with radar imaging for the first time (Torres et al., 2016).
Previous surveillance of cracks was performed through visual analyses and only few studies investigated automatic methods 50 for mapping the fracturing of ice shelves. Moctezuma-Flores and Parmiggiani (2016) proposed the use of a morphological filter for reducing speckle noise in SAR backscatter data, followed by a stochastic segmentation for mapping the pre-collapse fractured area of Nansen Ice Shelf. However, this approach was applied on a subset of the SAR image focusing on a widely opened fracture. In practice, SAR backscatter images have a limited sensitivity to narrow cracks in the early stage of their development because of the sensor resolution and the variable contrast dependent on the viewing geometry. Edge detection 55 performed on wide SAR images often misses thin cracks and may confuse fractures with other topographic features, e.g. calving fronts, if no contextual information is used.
In addition to SAR backscatter imaging, a few studies reported on the potential of SAR interferometry (InSAR) for mapping rifting activity Larour et al., 2004;Hogg and Gudmundsson, 2017;De Rydt et al., 2018).
These studies showed that, in an interferogram, opened fractures correspond to well-defined and visually identifiable phase 60 discontinuities. Rignot and MacAyeal (1998) identified rifts as branch-cut discontinuities and interpreted the fringe patterns over downstream ice shelf fragments as due to a rigid-body rotation about an axis perpendicular to the ice shelf surface and located at the tip of the rift. The analysis of double difference interferograms and the modelling efforts presented in the https://doi.org/10.5194/tc-2021-296 Preprint. Discussion started: 12 October 2021 c Author(s) 2021. CC BY 4.0 License. complementary paper by MacAyeal et al. (1998) support the hypothesis that this rigid-body rotation originates from gravitational ice flow. At the time of writing, to our knowledge, no study proposed a method for extracting the crack location 65 automatically from the complex phase information supplied by an interferogram.
In this paper, we present an automatic method for delineating ice shelf fractures using Sentinel-1 Interferometric Wide (IW) SAR interferometry (Yague-Martinez et al., 2016). The proposed method exploits the changes in the ice stress field and the shearing of the ice flow caused by the rifting activity, which both translate into a discontinuous fringe pattern in an interferogram. We show that an active crack separates an ice shelf into distinct regions, characterized by fringe patterns with 70 different orientations and different fringe rates that can be quantitatively derived by calculating the phase gradient, and that active cracks correspond to spatial phase discontinuities that can be mapped with an edge detection procedure. Because Brunt Ice Shelf (BIS) showed significant rifting activity in the past years (De Rydt et al., 2019), we select it as test site. The performance of the method is demonstrated using a set of 6-day repeat-pass Sentinel-1 interferograms acquired over BIS between September 2020 and March 2021. In particular, we track the activation and the propagation of a new rift, the North 75 Rift, that led to the calving of iceberg A74 on 26 February 2021. Based on this study case, we demonstrate that SAR interferometry is sensitive to the dynamical response of an ice shelf to rifting activity and has potential to provide early indications of fracturing, not yet visible in SAR backscatter or optical satellite images.
In Section 2, the BIS test site is briefly described and, in Section 3, the ability of InSAR at capturing rifting activity is introduced. The delineation method and the processing line are described in Section 4, along with illustrating examples of 80 intermediate processing steps. The capabilities and limitations of the method are also discussed. The Sentinel-1 dataset used for tracking the North Rift propagation on BIS and the processing parameters are detailed in Section 5. In Section 6, we analyze the results and present the evolution of the North Rift extent as captured by the InSAR-based delineation.
Furthermore, we illustrate the gain of information provided by SAR interferometry versus SAR backscatter images and we attempt to identify the stress field variation using double difference interferograms. Finally, Section 7 summarizes the 85 potential, the benefits and the limitations of the proposed method.

Test site
Brunt Ice Shelf is located along the Caird Coast in East Antarctica, in the eastern sector of the Weddell Sea (Figure 1(a)). It is connected to a larger ice shelf that is made up of the Stancomb-Wills Ice Tongue (SWIT) and the Riiser-Larsen Ice Shelf.
The border with the SWIT is defined at the northeast of BIS by the Brunt-Stancomb Chasm (Anderson et al., 2014). 90 BIS presents a rich rifting system made of old and new fractures experiencing variable propagation rates, that are described in Thomas (1973) and more recently in De . The McDonald Ice Rumples (MIR), on the northeast part of the ice shelf tip, originate from the only pinning point constraining the ice flow, and play therefore a key role in the rifting system De Rydt et al., 2019). On the southwest part of BIS, the Chasm 1 rift reactivated in November 2012 and progressively propagated towards the MIR to reach its current extent, after having been dormant for three decades. In October 2016, a new rift named the "Halloween crack" appeared close to the MIR and expanded with a variable speed in the east direction . More resilient than first expected, to date neither Chasm 1 nor the Halloween crack have reached their rupture point yet, although Chasm 1 remains only connected by a short ice bridge of a few kilometers in length.   (Thomas, 1973). In the following, we describe the InSAR-based method for delineating crack and illustrate it with an example on BIS. We then apply the method to a time series of repeat-pass interferograms and demonstrate its efficiency by tracking the activation and propagation of the North Rift. 110

Phase discontinuities across rifts
The efficiency of repeat-pass SAR interferometry for mapping ice motion and detecting ice surface deformation and stress (e.g. ice drift, grounding lines, grounding of pinning points, etc.) has long been established (e.g. Rignot et al., 1995Rignot et al., , 2000Mouginot et al., 2019). SAR interferometry is able to capture displacements and deformations of an ice shelf as small as a fraction of the sensor wavelength, making it highly sensitive to changes in comparison with SAR backscatter imagery 115 that can only provide information about cracks at the scale of the pixel resolution. Since 2014, the Sentinel-1 constellation offers interferometric capabilities at C-band with systematic acquisition, wide coverage, medium resolution and a reduced revisit cycle compared to former missions, therefore providing potential for a regular monitoring of ice shelf rifting activity with interferometry for the first time. The main acquisition mode of Sentinel-1 is the Interferometric Wide swath (IW) mode, based on Terrain Observation with Progressive Scans SAR (TOPSAR) (De Zan and Monti Guarnieri, 2006). In IW mode, 120 SAR images are divided in three swaths, each swath being made of Single Look Complex (SLC) tiles called bursts (Torres et al., 2012). Interferometric processing of Sentinel-1 IW acquisitions requires some specific processing steps, such as deramping or burst stitching, that are explained in Yague-Martinez et al. (2016).
Let us consider a repeat-pass interferogram computed with a master image acquired at epoch and a slave image acquired at epoch . After subtraction of the flat earth and topographic phase components, the phase in the interferogram can be 125 expressed as a sum of the following phase components: where is the ice flow motion phase component, is the tidal phase component, is the random phase noise and accounts for any deformation or displacement that does not originate from ice flow or tidal deformations, e.g.
iceberg drift, mechanical stress or structural deformation. In the following, for the sake of simplicity, we assume that the 130 phase noise is negligible. While the tidal component is determined by the change of the vertical position ⃗ ⃗ of the ice shelf between epochs and resulting from the balance between oceanic tides and atmospheric pressure (Padman et al., 2003), where is the local incidence angle.

145
In practice, interferograms over an ice shelf exhibit a clear segmentation of the fringe pattern, with discontinuities corresponding to the rifting system. This has already been observed in several cases for different ice shelvese.g. Brunt, Larsen-C, or Ronne ice shelves Larour et al., 2004;Hogg and Gudmundsson, 2017;De Rydt et al., 2018) and it is further illustrated in this paper for Brunt Ice Shelf. In Fig. 1, the rifting system of BIS is pictured and compared to a repeat-pass interferogram generated from Sentinel-1 images acquired 6 days apart in September. We observe 150 that active cracks, rifts and chasms correspond to phase discontinuities in the fringe pattern, dividing the ice shelf into regions with different fringe rates and fringe orientations. As highlighted by Eq (1), the segmented fringe pattern suggests https://doi.org/10.5194/tc-2021-296 Preprint. Discussion started: 12 October 2021 c Author(s) 2021. CC BY 4.0 License.
that SAR interferometry captures a spatially discontinuous ice flow velocity, tidal response and mechanical stress field created by rifting activity.
In the region of BIS, the vertical position of the floating ice, determined by the amplitude of oceanic tides and the variation 155 of atmospheric pressure, may vary by more than 1 m between acquisitions 6 days apart. However, away from the grounding line, the change of vertical position for a given pair of acquisitions shows smooth spatial variations of only a few centimeters. In comparison, the ice flow can reach velocities up to 2.5 m d -1 (ENVEO CryoPortal), partially captured by the LOS orientation. As a consequence, the tidal phase component corresponds to largely spaced fringes, and the phase signal is dominated by the ice flow motion contribution. 160 The segmentation of the interferogram into regions with distinct phase ramps related to the rifting activity constitutes the basis of the method presented in the following for delineating cracks automatically. Because the response of ice shelves to rifting and the mechanisms driving the crack propagation are still poorly understood, we assume that all phase contributions can have a spatially discontinuous behaviour and the crack detection is performed on interferograms containing the information of all three phase components. Let us note that, in the following, we will refer to repeat-pass interferograms 165 whose flat earth and topographic phase components have been subtracted as differential interferograms, or simply repeatpass interferograms. The terminology double difference interferogram will be used for interferograms computed as the difference of two differential interferograms.

Method
To achieve an automatic delineation of active cracks on ice shelves, we propose a method that aims at detecting 170 discontinuities in the phase image. The method consists of four major steps: 1) the generation of a repeat-pass differential interferogram, 2) the derivation of the phase gradient map, 3) a Canny edge detection applied to the image of the phase gradient magnitude, and 4) the vectorization and cleaning of the edge detection results. The processing line is presented in Fig. 3. The flowchart describes the input data, the auxiliary data, the processing steps and the output product of each step. In the following, each processing step is described in detail and exemplified using the Sentinel-1 interferogram of BIS shown in 175 First, Sentinel-1 repeat-pass interferograms are generated from IW SLC acquisitions according to the method presented in Andersen et al. (2020), which is optimized for ice velocity measurements with TOPSAR interferometry. Due to the steering of the antenna from the aft to the fore during each burst acquisition, hence introducing different viewing angles at the 180 overlap of one burst and the next, Sentinel-1 repeat-pass TOPSAR interferograms may suffer from phase jumps caused e.g.
by along-track ice motion (De Zan et al., 2014). Furthermore, ice flow motion, especially in fast-flowing regions, shifts phase centers of the slave image with respect to their location in the master image, leading to potential decorrelation. In addition to account for the precise state vectors, the coregistration procedure compensates for the local shifts between the https://doi.org/10.5194/tc-2021-296 Preprint. Discussion started: 12 October 2021 c Author(s) 2021. CC BY 4.0 License. master and slave burst SLCs caused by the along-track and across-track ice flow components, derived e.g. from offset-185 tracking, in order to reduce phase jumps and improve the coherence. After coregistration of the master and slave burst SLCs, the interferogram is generated at the burst level, the flat earth and topographic phase components are subtracted and the burst interferograms are stitched together to form an area-wide differential interferogram.

190
The phase signal in the interferogram is a sum of the ice motion component, the tidal component and the random phase noise, that results in distinct phase ramps throughout the ice shelf corresponding to separate regions of the rifting system ( Figure 1). The wrapped differential interferogram is geocoded and subsequent processing steps are applied to the geocoded fringes.
In the second step, the phase gradient is calculated pixelwise for each geocoded differential interferogram. Given , the 195 value of the wrapped phase for a pixel with coordinates ( , ) in the interferogram, and being the discretization indices in the x-and y-directions respectively, the phase gradient is written as ) . ( The temporal indices and are here omitted for the sake of readability. The x-and y-directions refer to the axes of the map projection, which is the Antarctic Polar Stereographic projection (EPSG 3031) in this case. The discrete phase derivatives 200 are computed by averaging the phase differences between adjacent pixels along the x-and y-directions over a square window with an odd size of pixels. The averaging is performed using the complex representation of the phase, as expressed by where ∠ represents the argument of the complex exponential. Equations (6) and (7) provide a phase variation per pixel, 205 assuming that pixels are square. If the aspect ratio of the pixel is different than one, then a scaling factor should be applied.  The calculation of the phase gradient translates the complex information provided by the spatially variable fringe pattern into a two layer real image (i.e. the x-and y-components corresponding to the projection axes in the cartesian case, the gradient magnitude and angle in the polar case). Moreover, by computing the phase gradient directly from the wrapped phase (expressed as a complex number), the tedious step of phase unwrapping is avoided, as well as the phase artifacts that it may introduce. 215 In practice, the phase gradient is converted to a polar vector, whose magnitude holds the information about the local fringe rate and whose angle indicates the direction of the phase ramp. For the demonstration case, the images of phase gradient magnitude and angle are shown respectively in Fig. 4(a) and 4(b). In both the magnitude and direction images, the location of the phase discontinuities, i.e. the active cracks, is enhanced and has become easy to identify (e.g. North Rift, Chasm 1 or Brunt-Stancomb Chasm). For most of the identified rift structures, the edges in the magnitude and direction images indicate 220 similar locations. However, for wide open chasms with a complex structure (e.g. the widest part of Chasm 1 or the Brunt-Stancomb Chasm), the magnitude and direction of the phase gradient may picture a slightly different fractured area. Given that the crack locations mapped by both indicators are mostly similar and that the angles are wrapped by nature, which makes them difficult to manipulate, we neglect the phase gradient direction and focus on the information held by the phase gradient magnitude. 225 In the phase gradient magnitude image, active cracks correspond generally to a well-defined step-edge with variable contrast.
Crack delineation is hence performed by applying a Canny edge detection to this image (Canny, 1986). In order to reduce the noise while preserving the edges, a median filter is applied beforehand. In practice, the Canny edge detection consists of computing the intensity gradient of the input image and applying a double threshold for mapping the edges: the upper threshold discriminates the strong edges; the lower threshold is used for selecting the weak edges, meant to connect the 230 strong edges present in their neighborhood. An additional Gaussian filtering, with tunable standard deviation, is performed as part of the Canny edge detection procedure for reducing the noise while computing the intensity gradient.
We focus on the rifting activity on the ice shelf and therefore mask the areas of grounded ice. The masking is performed with a simple thresholding of the TanDEM-X global digital elevation model at a 50 m height, that shows a rough agreement with the grounding line location. Patchy decorrelation can also cause erroneous edge detection and areas with low coherence (< 235 0.12) are therefore also excluded.
As shown by the blue lines in Fig. 5, applying the Canny edge detection to the gradient magnitude efficiently maps the cracks present on BIS. However, it also maps small "dangles" (loose curvy lines) caused by phase artifacts, topography or crevasses. As a final step, the raster mask generated by the Canny edge detection module is thinned, vectorized and cleaned using GRASS GIS tools. The cleaning consists of removing small dangles, as they are less likely to correspond to a major 240 propagating crack. The cleaning step is critical and the result should be carefully evaluated, because poor thresholds can lead to a substantial loss of detectable cracks. Though necessary for removing noise and obtaining readable information, the cleaning process sets a bound on the minimum size of detectable cracks, depending on the dangle size threshold.
As shown by the red lines in Fig. 5, the noisy aspect is reduced after the cleaning step. Some residual errors remain, especially in the region near the grounding line. This area is highly crevassed due to the tidal bending and the rapid height 245 change at the transition between floating and grounded ice. These residual dangles thus correspond to a damaged area, though crevasses are not the damage that we aim at mapping. Future studies might investigate the density of detected edges in the vicinity of the grounding line as an indication of the degree of crevasse damages.
Although the InSAR-based crack delineation performs well under conditions that preserve the phase coherence, its applicability is primarily limited by the quality of the SAR interferogram. In case of fast-flowing ice, snowfall, surface melt 250 or snow drift caused by katabatic winds, the InSAR signal decorrelates and the method cannot be applied. Even so, the regularity of Sentinel-1 acquisitions offers an increased likelihood of coherent interferometric pairs. Another limitation directly originates from TOPSAR interferometric processing, which is strongly affected by coregistration errors, uncorrected ionospheric delays and uncompensated tidal displacements, that may leave residual phase discontinuities at the burst overlap (e.g. see the western tip of the ice shelf in Fig.1(b)), hence causing potential false detections (De Zan et al., 2014). 255 Let us note that, in the illustrating example, the Halloween crack is only partially mapped even though it is visible in the brightness images. The fringe pattern is similar on both sides of the crack and the corresponding gradient discontinuity appears faint in the gradient image, which might indicate that the rift was inactive during the investigated time period.

Dataset and processing
To capture the propagation of the North Rift and the calving of A74, the InSAR-based method for crack delineation is tested on a dataset of Sentinel-1 SLC images acquired every 6 days between 1 September 2020 and 6 March 2021, along track 50 265 ( Figure 6). We selected the frames covering BIS and generated all available 6-day repeat-pass interferograms. Overall, 32 interferograms were generated, out of which 9 could not be used for crack delineation because of signal decorrelation. In particular, the interferograms spanning the periods directly before and after the calving event (i.e. between 10 February 2021 and 6 March 2021) could not be used for mapping the North Rift with InSAR. Because the amount of fringes caused by ice motion increases with the temporal baseline, most of the Sentinel-1 interferometric pairs acquired along track 164, with 12-270 day repeat pass, were decorrelated. Given their limited quality and temporal resolution, these pairs are not used for crack delineation, but only to support the interpretation and analysis of the temporal evolution of fracturing.

275
In order to reduce the phase noise, the interferograms are computed with a multilooking factor of 9 x 3 pixels in the slant range and azimuth directions, respectively, and an adaptive Goldstein filtering is applied for further noise reduction (Goldstein and Werner, 1998;Baran et al., 2003). For the coregistration, the average ice flow motion is compensated using a multiannual ice velocity map with 200 m pixel spacing generated within the ESA Antarctic Ice Sheet CCI project (ESA Climate Office: Ice Sheets Antarctic). The Antarctic velocity map is calculated using offset-tracking applied to all Sentinel-1 280 6-and 12-day repeat pairs available for the mission lifetime (2014-today). The offset-tracking processing is described in Nagler et al. (2015), and an additional correction of the vertical displacement induced on floating ice by differential tides (CATS2008; Erofeeva et al., 2019) and atmospheric pressure (ERA-5; C3S, 2017) is applied to offset-tracking results. The topographic phase is estimated and subtracted from the repeat-pass interferograms using the TanDEM-X polar DEM with 90 https://doi.org/10.5194/tc-2021-296 Preprint. Discussion started: 12 October 2021 c Author(s) 2021. CC BY 4.0 License. m grid resolution, extended to cover ice shelves (Wessel et al., 2021). Finally, the interferograms are geocoded on a grid 285 with a 40 m pixel spacing in the Antarctic Polar Stereographic reference system, that matches approximately the pixel size in the radar geometry.
For each geocoded repeat-pass interferogram, the phase gradient is calculated over a 9 x 9 pixel window and the image of the phase gradient magnitude is further filtered using a median kernel of 9 x 9 pixels. The Gaussian kernel has a standard deviation equal to 5, and the lower and upper thresholds for edge detection are set respectively to 0.15 and 0.21. 290 Furthermore, all areas above 50 m height or having a coherence lower than 0.12 are masked out for the edge detection procedure. The raster output of the edge detection is thinned, vectorized and cleaned by applying a threshold of 2000 m on the dangle size.
6 Results and discussion

Interferogram time series 295
The sensitivity of Sentinel-1 repeat-pass interferograms to the North Rift propagation is illustrated by a time series of interferograms presented in Fig. 7. These differential interferograms are all acquired along track 50 and computed as described in the previous section. The black dashed line represents the manually delineated calving front location, after the calving event of A74, and it is used to picture the North Rift's maximum extent. In this figure, it is visible that the discontinuity line expands along the crack location from one interferogram to the other. The magnitude and direction of the 300 phase gradient is highly variable and provides a first qualitative indication of the propagation timeline of the North Rift: in November 2020, the rift propagates as a straight line in a given direction and the phase ramp on the northern side of the rift varies gently; at the end of December, after a complete rotation of the fringe pattern, the crack changed direction to propagate toward the Brunt-Stancomb Chasm; at the end of January, when part of the northern plate is decorrelated and the fringe pattern is very tight, it is visible that the crack almost reached the chasm. 305

Automatic delineation
The proposed InSAR-based method for crack delineation is applied as described in the previous section to produce a time series of vector files mapping the BIS crack system with a 6-day resolution, only interrupted whenever the interferometric phase is decorrelated. Though a rough estimate could already be obtained from the visual analysis of the interferogram time series, the automatically delineated cracks enable a refined and more objective determination of the North Rift propagation 310 from the start of the crack propagation until the calving of A74. In Fig. 8, some selected repeat-pass interferograms are displayed at a larger scale together with the phase gradient vector field and the automatically delineated cracks. It is seen that the North Rift started to open up a few kilometers between 18 and 24 November 2020. Indeed, the phase discontinuity indicated by the yellow line expands between 12-18 November 2020 and 18-24 November 2020. Around the expanding end of the crack, we observe a rotation of the phase gradient that corresponds to a change in the stress field. Compared to the 315 https://doi.org/10.5194/tc-2021-296 Preprint. Discussion started: 12 October 2021 c Author(s) 2021. CC BY 4.0 License.  The propagation history is summarized in Fig. 9, where automatically delineated cracks representing the major steps of the 330 North Rift advance are shown. This figure results from a selection of the relevant dates and the coherent interferograms in the processed dataset. The background is a Sentinel-1 brightness image acquired on 6 March 2021, that shows a good agreement between the delineated cracks and the calving front after the iceberg A74 formed. As already stated above, the crack started to open sometime between the 18 and 24 November 2020, and continued to rapidly propagate in early December. As the rift opened, both sides of the crack slowly drifted away from each other, leading to a massive coherence 335 loss on the northern side of the rift and some erroneous detections after mid-December. Despite the phase noise, the crack delineation method was able to map the easternmost part of the North Rift and the results agree well with the shape of the manually derived calving front. From the interferogram of 17-23 January 2021, about a month before the iceberg broke off, the quasi-full extent of the calving could already be mapped.

Propagation rates
In addition to the delineation, the North Rift length is also estimated for the selected dates and the evolution of the rift length 345 over time is plotted in Fig. 10. The crack length is estimated as the cumulative length of the delineated segments corresponding to the North Rift. In the case of incomplete or interrupted segments, e.g. due to decorrelation, the missing https://doi.org/10.5194/tc-2021-296 Preprint. Discussion started: 12 October 2021 c Author(s) 2021. CC BY 4.0 License. segment length is estimated manually. As expected, the estimated length consistently grows over time, except for the pair of 12-18 November 2020, whose crack extent is estimated to be smaller by about a hundred meters compared to the earlier pair of 7-13 September 2020. This inconsistency stresses the inaccuracy of the estimation method, that is used for analysis 350 purposes rather than accurate quantitative estimates. Based on the estimated crack extent, the propagation rates are derived as the segment slope between two consecutive points (see the annotated values in Fig. 10). The propagation rate is the largest (about 1.3 km d -1 ) when the North Rift activates at the end of November and reaches a similar value at the latest stage of propagation in January. In between, the progress of the North Rift extent is in the order of several hundred meters per day.

Comparison with SAR backscatter imagery
In order to stress the added value of SAR interferometry for mapping active rifts, we present in Fig. 11 the evolution of the North Rift as seen with a time series of Sentinel-1 backscatter images. While SAR interferometry is sensitive to subcentimeter deformations between two epochs, SAR backscatter imagery characterizes the static situation at a given time. 360 Furthermore, its application to crack monitoring is limited by the image resolution and the lack of backscatter contrast over thin cracks. In the brightness images, the North Rift is visible in November, but shows a relatively small extent compared to that mapped with interferometry (see interferograms in Fig. 7). The well-advanced breach is visible in the brightness images only a few days before the calving of the iceberg. Comparatively, the interferogram of 17-23 January 2021, fully captures the extent of the crack that becomes visible only a month later in the brightness image acquired on 22 February 2021. 365

Stress field variations
The discontinuity of the fringes constitutes the assumption behind mapping active cracks with SAR interferometry and the variability of the phase gradient magnitude and orientation has been illustrated in Fig. 7 and Fig. 8. However, the origin of https://doi.org/10.5194/tc-2021-296 Preprint. Discussion started: 12 October 2021 c Author(s) 2021. CC BY 4.0 License.
these temporal variations remains difficult to determine from differential interferograms, as the differential phase is a combination of ice flow, ice deformation and intrinsically variable tidal displacements (see Eq. (1)). Temporally variable 370 components of the phase signal can be isolated by computing differences between repeat-pass interferograms.  To highlight the variable phase contributions, Equations (2)-(4) can be reformulated in terms of the variation with respect to 385 the bulk contributions of the ice flow velocity , tidal vertical displacement and deformation : = 4 ( + Δ ) cos , and (9) A double difference interferogram computed from two repeat-pass interferograms spanning respectively epochs , and 390 , has a phase Δ , written as: Thanks to the almost exact repeat orbit of Sentinel-1, the line-of-sight direction and the incidence angle can be assumed to be constant for all interferograms along a given track. The phase difference can therefore be calculated easily using repeat-pass interferograms geocoded on a common grid. If the repeat-pass interferograms have the same temporal baseline, i.e. Δ = 395 Δ = Δ , inserting Eqs. (8)-(10) into Eq. (11), the phase of the double difference interferogram can be formulated as: Equation (12) shows that the differentiation of two interferograms with the same temporal baseline removes the bulk contribution of the ice flow velocity, the oceanic tides and the mechanical stress field. In the absence of ice flow speedup and other deformations than the tidal ones, the double difference interferogram contains only the change of tidal bending 400 between 4 epochs (3 epochs in case of a common date) . Nevertheless, if the ice shelf undergoes changes in flow velocity or creep deformation, such differentiation would remove the bulk contributions of the ice flow and of the tidal deformation, leaving only the natural variations caused by the change of tidal amplitude and the deformation between 4 (or 3) epochs, caused e.g. by rifting activity.
For isolating the temporally variable component of the phase signal introduced by rifting activity, differences between 405 consecutive 6-day repeat-pass interferograms along track 050 are computed. These double difference interferograms are shown in Fig. 12. The first one is the difference between the interferograms of 12-18 November and 18-24 November 2020.
It shows curved fringes on the expanding tip of the crack, likely caused by the mechanical stress of the diverging ice plates as the crack starts to propagate. In the second double difference interferogram (18-24 November and 24-30 November 2020), the stress field becomes almost evenly distributed on both sides of the rift (14 fringes on the north plate and 12 on the 410 south plate). This number of fringes corresponds to a displacement of about 35 cm in the direction of the line-of-sight, which is almost parallel to the fringe orientation. Finally, in December 2020, when the crack propagation changes direction, the fringe pattern becomes very tight and remains almost evenly balanced on the north and south plates. Without the information https://doi.org/10.5194/tc-2021-296 Preprint. Discussion started: 12 October 2021 c Author(s) 2021. CC BY 4.0 License. from another viewing direction, the observed phase ramp could either correspond to a vertical deformation or to a horizontal displacement. For solving this ambiguity, we also display the difference between 12-day interferograms acquired in 415 November along the opposite direction (track 164). Interferograms along this orbit have a lower temporal resolution and quality, and they are therefore only used here to support the analysis. The double difference results in this case into a homogeneous fringe pattern with fringes oriented nearly parallel to the crack. For both tracks, the double difference interferograms show fringes nearly parallel to the line-of-sight direction around the North Rift. Peltzer et al. (1994) simulated the fringe pattern due to a rigid-body rotation around an axis that is perpendicular 425 to a horizontal surface and demonstrated that it would create such fringes parallel to the viewing direction. Later on, Rignot and MacAyeal (1998) also observed fringes parallel to the LOS direction on Ronne Ice Shelf and interpreted it as rigid-body rotation around an axis perpendicular to the ice shelf surface and located at the tip of the crack. Because they observed this pattern in the differential interferograms but not in the double difference ones, they attributed this rigid-body rotation to a velocity difference between both sides of the rift, not to a transient phenomenon. In the case of the North Rift, fringes 430 parallel to the line-of-sight are observed with two different viewing directions and a rigid-body rotation around the tip of the crack seems therefore likely. Since this pattern is present in the double difference interferograms, it would be in this case associated to a time-varying response to rift propagation that dominates the variation of vertical tidal displacements.
The hypothesis of the rigid body rotation is further strengthened by ice velocity measurements performed with offsettracking and shown in Fig. 13. The offset-tracking measurements show an acceleration of the ice flow around the southwest 435 part of the North Rift from December 2020 onwards. Simultaneously with the North Rift propagation in November-December, a speedup of about 0.5 m d -1 is measured on the tip of the northern plate. In January, the measured velocity further increases to reach a speedup of about 2 m d -1 on the tip. The amplitude of the speedup varies from the southwest to the northeast, similarly to the stress field observed in the double difference interferograms, and this acceleration comes together with a rotation of the velocity vector field towards the northeast and southwest directions respectively north and 440 south of the rift. The change in orientation of the measured velocity field fits indeed with a clockwise rotation around the tip of the North Rift. Although this could originate from ice flow speedup due to the loss of constraint at the MIR, it could also be caused by a drift/rotation of the whole iceberg that would be misinterpreted as ice flow acceleration by offset-tracking.

Conclusion
The InSAR-based method proposed for automatic crack delineation has been successfully tested and qualitatively validated 445 on BIS with a dataset of Sentinel-1 6-day repeat-pass interferograms spanning a period of 6 months. The applicability of the method has been demonstrated by tracing the propagation history of the North Rift, from the rift activation up to the calving of the iceberg A74.
For the North Rift, the shape of the delineated crack agrees well with the calving line location after the iceberg calving, thereby demonstrating the suitability of the approach. In general, phase artifacts in the interferograms may introduce noise in 450 the delineated cracks, but the InSAR-based method is still less impeded by topography and structural heterogeneity of the ice shelf than SAR backscatter imaging. A limiting factor of the method is the decorrelation caused by snow drift, snow melt or fast-flowing ice, as it prevents any InSAR measurements. https://doi.org/10.5194/tc-2021-296 Preprint. Discussion started: 12 October 2021 c Author(s) 2021. CC BY 4.0 License.
In addition to the propagation history, the temporally variable phase contribution could be isolated and interpreted as rigidbody rotation about the expanding tip of the North Rift in response to the rifting activity. Without further information, it is 455 not possible to determine whether the rotation origin is an ice flow speedup or ice drift.
Combined with the continuous 6-day coverage of the Antarctic margins by Sentinel-1, the InSAR-based crack delineation opens the possibility for operational monitoring of natural or climate change-driven damages and rifting activity over most ice shelves, as well as the detection of ongoing breakoff processes and precursor signs of calving events. Thanks to the high sensitivity of InSAR to dynamic changes in the ice shelf stress field, the rifting activity can be captured well before it is 460 visible in SAR backscatter images. With SAR interferometry, the quasi-full propagation of the North Rift could be mapped from image pairs acquired about a month before the iceberg broke off. This is a major advantage for predicting future calving events and to improve modelling of the response of ice shelves to damages.
Future work should focus on testing the method on other ice shelves, with various structures and geometries, and on improving the post-processing for reducing the errors in the detected cracks. 465 Data availability. Sentinel-1 satellite data are freely available on the ESA Open Access Hub (https://scihub.copernicus.eu).
Interferometric products and delineated crack generated in this study can be made available on request.
Author contributions. LL developed the algorithm, performed the processing, led the data analysis and drafted the manuscript. JW and TN contributed to the data interpretation and the discussion of the results, and revised the manuscript.
Competing interests. The authors declare that they have no conflict of interest. 470