A Systematic Study of the Fracturing of Ronne-Filchner Ice Shelf , Antarctica , Using Multisource Satellite Data from 2001 to 2016

We propose a new framework of systematic fracture mapping and major calving event prediction for the large ice shelves in Antarctica using multisource satellite data, including optical imagery, SAR imagery, altimetric data, and stereo mapping imagery. The new framework is implemented and 10 applied for a comprehensive study of the fracturing of Ronne-Filchner Ice Shelf (RFIS), the second largest ice shelf in Antarctica, using a long time dataset dating back to 1957. New remote sensing data that have been made available in the past decade, including Landsat 8, WV-2, ZY-3 and others, greatly enhance our abilities to detect new fractures and monitor large rifts in three dimensions. Two large rifts, Rifts 1 and 2, were newly detected and are comparable to the Grand Chasm that caused a major calving event in 15 the region in 1986. Three-dimensional rift models generated from quasi real-time stereo ZY-3 images revealed important topographic information about the large rifts that can be used to improve the reliability of ice shelf modeling and support enhanced analyses of ice shelf stability. Based on the results of the 2D and 3D fracture mapping, the spatial and temporal analyses of the overall fracture changes and large rift evolutions, i.e., the level of fracturing in RFIS, were slightly increased, particularly at the front of the ice 20 sheet. The overall fracture observations do not seem to suggest immediate significant impacts on the stability of the shelf. However, the most active regional fracturing activities occurred at the front of Filchner Ice Shelf (FIS). A potential upcoming major calving event of FIS is estimated to occur in 2051. The stability of the ice shelf, particularly with regard to the developments of Rifts 1 and 2, should be closely monitored. 25


Introduction
Ice shelves play a very important role in buttressing the Antarctic ice sheet (AIS).They are also the the MOA (Mosaic of Antarctica) data for various ice shelf studies (Hulbe et al., 2010;Albrecht and Levermann, 2014;Bassis and Ma, 2015).Furthermore, SAR (Synthetic Aperture Radar) images, such as those from RADARSAT and ERS (European Remote Sensing) satellites, work regardless of clouds and brightness conditions and have also been frequently used in fracture studies (Fricker et al., 2005b;Li et al., 2016).The depth retrieval of large rifts has been mainly restricted to satellites and airborne laser altimetry, which only resolve rifts along sparse profiles or airborne stereo photogrammetric pairs (Fricker et al., 2005 a;Howat et al., 2012;Massom et al., 2015).An extensive survey of 78 rifts in 13 Antarctic ice shelves from 2002 to 2012 was performed by Walker et al. (2013).The location, time, rift length and change rate information were derived from satellite imagery.
In recent years, the availability of high-resolution stereo mapping satellites, such as Worldview-2 (WV-2) and Chinese ZY-3, have opened up the possibility of accurate three-dimensional (3D) rift measurements from the along track images in addition to cross-track stereo mapping (Sefercik et al., 2013;Tang et al., 2015;Tong et al., 2015).The along track stereo mapping technique using a ZY-3 images has the advantages of quasi-real-time capturing of the 3D ground movement, in addition to the increased image matching quality and strengthened geopositioning capability (Li et al., 2009).The stereo images used in this study presents a peak into the rifts using multiple viewing angles and maps them three dimensionally.
This paper presents the results of a systematic study of the fracturing of Filchner -Ronne Ice Shelf, Antarctica, by taking advantage of the data acquired by multiple satellite sensors from 2001 to 2016, including optical imagery, SAR imagery, altimetric data and others.Fractures of the entire ice shelf were mapped using the recent Landsat-8 OLI (Operational Land Imager) images from 2014/2015, with an improved resolution of 15 m.Compared with a fracture map of 2003/2004 derived from the MOA data and the two sets of ice flow speed maps of 1997/2000 and 2013-2016, systematic spatial and temporal analyses of fracture changes and trends were performed.The front of FIS was found to be the most active fracturing region, where two large rifts, Rift 1 and Rift 2, were detected, which may have the potential to trigger a major calving event, such as the one that occurred in 1986.The long-term behavior of these rifts were observed and analyzed using a series of datasets, including optical and InSAR images and altimetric data from 2001 to 2016.The high-resolution 3D topography of Rift 1 was reconstructed and its changes The Cryosphere Discuss., https://doi.org/10.5194/tc-2017-178Manuscript under review for journal The Cryosphere Discussion started: 11 September 2017 c Author(s) 2017.CC BY 4.0 License.
were investigated using ZY-3 and WV-2 stereo images.The 3D model revealed a valuable insight into the rift and included the wall slope, floor mé lange shape, depth changes and others.Finally, an empirical experiment was carried out to estimate the timing of the next possible major calving event at FIS.

Study area
RFIS consists of two individual ice shelves, Ronne Ice Shelf (RIS) and Filchner Ice Shelf (FIS), connecting the dynamically-thinning West Antarctic Ice Sheet (WAIS) and relatively stable East Antarctic Ice Sheet (EAIS), as indicated in Fig. 1.It has a total ice volume of 443140 km 3 .During the ICESat (Ice, Cloud and land Elevation Satellite) observational period of 2003-2009, the ice steams of RFIS drained the catchment basins from both WAIS and EAIS at a rate of 241±42 Gt a -1 (Moholdt et al., 2014); the calving process accounted for a discharge of -237±15 Gt a -1 and the basal mass budget was estimated to be -124±66 Gt a -1 .All the contributions, considering the net mass budget, was -35±60 Gt a -1 .
Both RIS and FIS have experienced recent calving cycles that resulted in immense calving areas.The most notable disintegrations in RIS occurred in 1998 and 2000, with a total area of loss of ~16300 km 2 , which removed the ice front advance accumulated since the late 1940s and early 1950s (Ferrigno et al., 2005).On the other hand, the last major calving event in FIS occurred in mid-April of 1986 and onward, and resulted an area loss of ~11500 km 2 , which was broken into three large tabular icebergs.The volume of the mass loss was estimated as 4850 km 3 , approximately 40 times the average annual discharge along the front of FIS (Neuburg et al., 1959;Ferrigno and Gould, 1987;Ferrigno et al., 2005).
In general, rapid developments of large rifts precede major calving events.A large rift in FIS, Rift 1 in Fig. 1, was detected by Walker et al. (2013) based on the MODIS images of 2011 and was reported, among rifts in other ice shelves, as a rift feature called "F6", without parameters such as length and rate of change.Using a set of images from 4 satellites over a period of 15 years in this study, Rift 1 was found to have expanded significantly and another rift, Rift 2, was also detected and monitored.Rifts 1 and 2 were approximately 72 km from the FIS front in 2013.and ice streams (IS) labeled on the background MOA mosaic (Scambos et al., 2007).

Data
The data used in this study are summarized in Table 1.There are four datasets with different satellite sensors covering the period of 2001-2016.The first dataset is used for ice-shelf-wide 2D-fracture detection and propagation speed estimation.It contains 41 Landsat-8 OLI L1GT cloud-free scenes of the Austral summer from October 3, 2014, to February 28, 2015, which seamlessly cover the RFIS region (north of 82.4 °S ).The scenes were orthorectified and made available by the United States Geological Survey (USGS) with a ground resolution of 15 m (panchromatic) and a geolocation accuracy of 18.1 m (Storey et al., 2014).To estimate the fracture propagation, the speeds of the fracture features of RFIS extracted by Hulbe et al. (2010)  125 m and a geolocation accuracy better than one pixel (Scambos et al., 2007).
The second dataset was used for a detailed three-dimensional mapping of Rift 1 and contains three scenes of stereo ZY-3 images (forward, nadir and aft looking) and covers almost the entire Rift 1 area (Fig. 2).ZY-3 is the first of the Chinese surveying and mapping satellite series, launched on January 9, 2012.Its sensor system is equipped with a three-line camera to build along track and quasi real-time (~30 seconds maximum) stereo pairs by combinations of forward, nadir and aft looking images.These stereo images can avoid radiometric variations in images caused by temporal changes and differing solar illuminations (Tong et al., 2015).The forward and aft cameras are inclined ±23.5° from the nadir one, forming a base-to-height ratio of 0.87 (Tang et al., 2015).An accuracy of 10.8 m (horizontal) and 6.1 m (vertical) can be achieved without GCPs (ground control points) (Tang et al., 2015).The satellite was first time programmed on February 2014 to adjust its sensor parameters to fit the Antarctic setting and collected a set of Antarctic stereo images for a number of selected sites.Two pairs of WV-2 panchromatic stereo images were provided by the DigitalGlobe, which is derived from the Ortho Ready Standard (OR2A) product to which radiometric, sensor and geometric corrections were applied.These images were taken in 2012 and 2016 and cover the eastern part of Rift 1 (Fig. 2).Both pairs are along track stereo images taken at a ground resolution of 0.5 m.The nominal geolocation accuracy without using GCPs was reported to be 2.3 m (horizontal) and 2.2 m (vertical), which are converted from the CE90 and LE90, respectively (Toutin et al., 2012).Two tracks of ICESat GLAS (Geoscience Laser Altimeter System) points cut across Rift 1 and were applied to detect the rift depths with a longer span from 2003 to 2009, in addition to the above stereo satellite imagery.The data set used is the GLA12 R633 product with an overall vertical accuracy of 14 cm (5 cm in flat areas) (Schutz et al., 2005;Shuman et al., 2006).The data have been applied for the global and regional mass balance estimations of AIS and Greenland Ice Sheet (GIS) (Ewert et al., 2012;Zwally et al., 2015;Xie et al., 2016).    of December 6, 2013, superimposed with a color-coded bottom marine ice thickness (Lambrecht et al., 2007).This is an enlarged area of the white box in Fig. 1.First, a 2D mapping of fractures was mainly carried out using optical and SAR images.The fractures are extracted from as many images as possible to cover the entire ice shelf and build a time series of observations.One of the key techniques is to identify and follow the same fractures so that changes in their rates can be estimated, in addition to detecting new fractures and trends.Second, high-resolution stereo satellite images are processed to build 3D digital elevation models (DEMs) of large rifts that may evolve rapidly and cause major calving events.Although the 3D satellite mapping data have been available only within the past decade in this area, the multiple sets of DEMs and rift change parameters allow us to gain an unprecedented insight into the temporal 3D interior changes of a rift.Finally, the results of the ice shelf wide 2D fracture mapping and detailed 3D large rift parameters can be used to perform a fracture propagation analysis.They can also be applied in many aspects of Antarctic research, such as ice shelf stability, iceberg prediction, ice shelf model support, mass balance estimation and other studies.In the following sections, a few important techniques and methods of the above framework are introduced and presented.

Extraction of two dimensional fractures
Fractures of RFIS can be extracted from the scenes of the Landsat-8 panchromatic band data with a resolution of 15 m.Based on the sampling theorem, fractures with an opening of 30 m (2 pixels) or wider can be detected.Since the extracted fractures are compared with those extracted in previous studies, such as the MOA mosaic, validation of the geometric alignment of the two datasets is necessary.In this study, the datasets involved are all on the same coordinate system, i.e. the Polar Stereographic South projection.
A set of 36 check points that are outcrop points in the region and have not moved during the time of two image sets were selected.A comparison between the check points measured in the two datasets resulted in an average difference of 71 m, with a standard deviation of 29 m.Thus, the average geometric difference is within 5 pixels of the Landsat-8 images and 1 pixel of the MOA images.No further geometric transformations of either of the datasets were made.
Then, a median filter is applied to reduce the salt and pepper noises.After that, a set of image enhancement techniques, such as adaptive histogram equalization, histogram matching, haze removal, Wiener filtering, and Wallis adaptive filtering can be used so that the images used are of uniform quality while maintaining a high level of visible image texture in the fracture areas (Fahnestock and Schowengerdt, 1983;Lim, 1990;Zuiderveld, 1994).Some areas may need separate special local treatments and a combination of the preprocessing techniques.
The automation of Antarctic surface feature extraction from remote sensing imagery is still a challenging research topic.For example, an experiment using the automatic extraction of tidal cracks in Antarctic fast ices using a Canny edge detector and Landsat-8 images was reported to have a relative success rate of approximately 51%; success is defined as the ratio of total lengths of all detected crack lines to the manually detected ones (Hui et al., 2016).An automation experiment for crevasses detection in the Shaune Garang glacier, a small 4 km 2 Himalayan valley glacier, using Landsat-8 images achieved an accuracy of approximately 6 m for the extracted crevasse lengths (Bhardwaj et al., 2016).Another experiment of automatic surface crevasse extraction mapped in the Antarctic interior region using a support vector machine technique and ASTER images, resulted in a classification accuracy of approximately 86%, given a large training sample size of 30% (Xu et al. (2011).Considering the relatively small experimental areas and the large numbers of training samples used in the above experiments, a manual extraction of the fracture features is appropriate in this study to achieve the desired accuracy and to consider the complexity of the fracture system of RFIS.The manual extraction of the fracture features was performed according to the guidelines in Glasser and Scambos (2008) and Holt et al. (2013).Special attentions were made to distinguish the low-contrast undulations in some areas caused by blowing-snow events from the actual fracture features (Scarchilli et al., 2010).

Fracture propagation speed
When two sets of fracture features from different times are extracted, their correspondences can be others, but represent different parts of the facture.There are two principles that we applied in selection of the corresponding points.First, the two end points of a fracture should not be selected because their locations may change as the fracture propagates over the time.Second, intersections between a fracture and a long-standing ice flow line can serve as corresponding points.In addition, corner points along a fracture that are visible in both images and have sustained the fracture propagation over the time period can also be selected as corresponding points.A verification using a 3D anaglyph check is often a necessity.
The resulting propagation speed of a fracture is the average calculated from all the corresponding points of the fracture.
A validation of the estimated fracture propagation speed can be carried out by a comparison with existing gridded speed maps.To avoid the influences of the errors in the speed maps, a moving window smoothing using a 5 by 5 grid window is performed.After this, the fracture propagation speed data can be used for trend and stability analyses.

Three-dimensional reconstruction of rift geometry
Three overlapping ZY-3 images (forward, nadir and aft) are used for an exterior orientation by a bundle adjustment (BA) procedure, and then the estimated orientation parameters are employed to photogrammetrically triangulate the ground points (Li, 1998;Grodechi and Dial, 2003;Tang et al., 2015).
Because the area of the rifts is dynamic, with a maximum ice flow speed of ~1200 m a -1 , stable GCPs are not available.Thus, a relative orientation of the images is performed using only the tie points selected as, e.g., fracture features and other relatively distinct surface features.The maximum displacement (forward and aft) of a tie point between the images due to ice flow is limited to within 2 mm and is negligible.The tie points are selected to be uniformly distributed in the area for a better BA performance.A further process is performed, which uses an affine model to reduce the BA residuals of the tie points (Tong et al., 2010).The absolute accuracies of the ground points should thus be ~10 m (horizontal) and ~6 m (vertical), which are the same as the case of the geolocation without GCPs in Tang et al. (2015) and Wang et al. (2014).However, the relative accuracy should be generally higher (Wolf et al., 2000;Li et al., 2011 and2017).After that a stereo pair can be formed using a combination of either forward-nadir or nadir-aft bundle adjusted images with a smaller intersection angle of 23.5 o to ensure that the rift mé lange floor is as largely covered as possible.This is also accompanied by the cost of a lower base-height ratio of 0.4 and a slightly lower ground accuracy.The above workflow was implemented using the ERDAS 2013 Leica Photogrammetric Suite.The adaptive Automatic Terrain Extraction (aATE) module is used for the DEM derivation, which allows for a coarse-to-fine hierarchical matching process (Leica Geosystems, 2006;Li et al., 2011 and2017).The normalized correlation coefficient matching window was set as 15 by 15 pixels and the initial correlation coefficient threshold was as 0.6.The aATE method relies on the assumption that the threshold decreases if the points have correlation coefficients lower than the threshold and increases otherwise (Leica Geosystems, 2006).The produced DEM has a grid spacing of 5 m.
The photogrammetric processing and DEM reconstruction from the WV-2 images follow a similar procedure as that of the ZY-3 images.The DEM grid spacing is set as 1 m.The generated DEMs can be used to estimate a set of rift parameters, including the depth, width, shape of the walls, mé lange floor topography and others, such as for ice shelf stability analysis and the improvement of physical models of ice shelves.
ICESat tracks are generally perpendicular to the rifts of the Antarctic ice shelves.A number of points along an ICESat track may fall in a rift opening and provide repeated depth information in addition to the DEMs from the stereo image pairs.The ICESat/GLAS product of GLA12 R633 was used.No data filtering based on product flags were performed in order to retain as many ICESat points as possible inside the rift (Fricker et al., 2005 a).The data were then corrected for saturation and Gaussian-Centroid (G-C) bias using correction information provided by the National Snow and Ice Data Center (NSIDC, 2015;Xie et al., 2016).These ICESat points were then projected from the TOPEX ellipsoid to the WGS84 ellipsoid, and then to the Universal Polar Stereographic coordinate system.A registration of each DEM to the ICESat data was carried out to compensate for residual height biases and tidal effects.Rifts with visible openings (blue lines in Fig. 4) are mostly located in ice shelf front region and some of them occur earlier on the way to the front.Among the 118 detected rifts, there were 11 longitudinal rifts along the frontline; there are four large rifts, Rifts 1 and 2 in FIS, and Rifts 3 and 4 in RIS.
The fractures extracted from the MOA mosaic of 2003/2004 (Hulbe et al., 2010) are listed as a dataset in Table 1 and   while propagating with the other fractures together toward the front, Rift 3 almost kept its length of ~56.0 km and its width has increased from 6.0 km to 6.4 km; Rift 4 expanded by 8.7 km to reach about the length of 36.7 km and its opening was measured 5.6 km, widened by 2.5 km.However, Rift 1 (52 km) are compared with the ice flow speed map in Fig. 5, derived from Landsat-8 images of 2013-2016 (Fahnestock et al., 2016;Scambos, 2017).An average difference of 5 m a -1 with a standard deviation of 27 m a -1 between the two speed datasets shows a relatively high level of consistency.
The ice shelf frontline was mapped using Landsat-8 band 8 images from 2014/2015 (orange line in Fig. 4) and the MOA mosaic of 2003/2004 (red line).The frontline propagation rates are, on average, 1270 m a -1 in RIS and 1200 in FIS, which are at a similar level as those of the front fractures.

Temporal 2D mapping of Rift 1 and Rift 2 in Filchner Ice Shelf
Rifts 1 and 2 in FIS were newly detected in the satellite images in 2001.They are close to the front and presented unstable patterns, unlike Rifts 3 and 4 in RIS.Thus, they may have the potential to cause a major calving event like the one that occurred in 1986 (Ferrigno and Gould, 1987).100 scenes of different 10 satellite images acquired from 2001 to 2016, including ASTER, Envisat ASAR, Landsat-7 ETM+, and Landsat-8 OLI (Table 1), were used to measure the temporal changes of the width and length of the two rifts (Fig. 6).Length was measured between the two ends and along a rift.Width was tracked at the same location of a rift, for example, at the intersection of the rift and a flow line as it propagated toward the by 37 km, and width by 520 m.The change process may not have taken more than 10 months because our observation window was limited by the data availability.Thereafter, the width is kept at a steady rate of change of ~100 m a -1 and so did the length at ~1.3 km a -1 (or 3.5 m d -1 ).
The maximum width measured are over 1900 m for Rift 1 and 1800 m for Rift 2 during the observation period.Compared to the status before their abrupt expansion events, both rifts were in a steady process of accelerated expansion in both their lengths and width.

3D mapping of Rift 1 in Filchner Ice Shelf using multisource satellite data
A set of stereo ZY-3 images covering Rift 1 (Table 1 and Fig. 2) were collected on February 28, 2014.
After processes of a BA and dense matching a DEM of Rift 1 was generated using forward and nadir looking images.There are areas blocked by cloud coverage, shadows, and the forward-looking angle (details see the section accuracy analysis).After trimming the blocked areas, the resulting DEM has a size of 9200 × 1000 grid points with a grid spacing of 5 m.A color coded 3D view and an enlarged view inside the rift are illustrated in Fig. 7.A video clip of the fly through of the 3D rift model is made available online: https://youtu.be/a4dSJ_z6qNc.The DEM reconstructed 3D topography of Rift 1 including its fine longitudinal flow lines, nearly upright rift walls, surface of mé lange floor, and others.The east end of Rift 1 is not covered because of the data coverage limit.
The total length of the rift measured along the centerline (Fig. 7a) of the DEM is 44463 m.Fig. 8 shows a set of key parameter measurements along the centerline, starting from the west end.A specific measurement at a location is calculated via an average of the corresponding elevations over a transit across the centerline.A moving average with a window of 500 m was performed to reduce the noise in the measurements.There is an easement transition of the depth or bottom elevation from the western tip for ~ 2.3 km (Fig. 8).The width of the rift reaches its maximum of ~1400 m in the middle of the rift (approximately 23 km from western end) and decreases toward both ends, with a minimum of 40 m at the west end point.The rift surface is relatively flat, with an average elevation of 72 m.Immediately after the transition from the western tip, there is a cavity that may not have been filled by ice mé lange initially.Thereafter, the rift floor is relatively flat, with an average depth of 47 m.The rift is generally deeper in the eastern part.The bottom roughness was estimated using the standard deviations of the surface slopes using a neighborhood of 5 × 5 grid points (Bamber et al., 2009).Furthermore, Profiles 1 to 3 were located roughly evenly along the rift (Fig. 7a) and are shown in  approximately stable (Fig. 10b).However, the rift widened at a rate of ~190 m a -1 .Therefore, the abrupt changes in length and width of Rift 1 and Rift 2 occurred at different times.
Furthermore, there is also not a strong correlation between the abrupt changes in the length and width and that of the depth.However, it should be noted that the change of depth has an approximately opposite pattern of length and width.Namely, during the period of relatively small changes of rift length and width (Fig. 6a), the rift depth increased.On the other hand, the depth remained at almost the same level while the length and width significantly increased.

Accuracy analysis
The reference system used in this study is a coordinate system of a Polar Stereographic South projection with a central meridian at 0° and standard parallel at 71 °S .The vertical datum is a WGS 84 ellipsoid, which is widely used in Antarctic studies including, e.g., the LIMA (Landsat Image Mosaic of Antarctica) product.LIMA has a resolution of 15 m and an accuracy of 54 m (Lee et al., 2004;Bindschadler et al., 2008).Other existing data sets in Table 1 were also provided in the same reference system to facilitate the analysis of a unified reference system.The data sets, such as 3D rift models that were produced in this study, were generated in the same reference system.
We carried out a validation of the Landsat-8 OLI L1GT images against the LIMA mosaic in the study area using a total of 27 control points that are outcrops in FRIS.The experiment achieved a standard The various images used to evaluate the length and width changes of Rifts 1 and 2 in Fig. 6 and their geolocation accuracies are listed in Table 1.Using a simplified error propagation method, the accuracies of the estimated lengths and widths derived from Image K with geolocation error k (Table 1) are 1.4 k, which ranges from 25 m to 76 m.
As indicated in the section of the ZY-3 image BA, the absolute accuracy of the ZY-3 ground points is ~10 m (horizontal) and ~6 m (vertical), which depicts the overall location accuracy of the ZY-3 DEM in the ground coordinate system.The relative accuracies of the parameters, such as distance and speed, measured from the images depend on the qualities of the matched points.We adopted two metrics from Li et al. (2017) and Li (1998).The first metric is correlation coefficients calculated using a window of 25 by 25 pixels.The result shows that the correlation coefficients of the points matched inside of Rift 1 (floor and walls) have a mean value of 0.72, which is significantly higher than 0.46 of the ice shelf surface because of the higher level of image texture of rift mé lange on the floor.The second metric is the image location differences (residuals) between the automatic and manually matched points.A grid is defined and overlaid on the matched points.In each grid cell, one matched point is selected as a check point and is compared with a manually matched point.A total of 135 such check points were used to perform the validation.The relatively flat ice shelf surfaces on both sides of Rift 1 have mean residuals of 0.70 pixels and that inside the rift is 0.77 pixels.the width of the combined affected area across the rift.Therefore, the area in the resultant DEM affected by these two factors is a strip with an 82 m width along the entire rift, which is interpolated when visualized.A related issue is the base-height ratio that is an important parameter for accurate stereo mapping.It is recommended that this parameter is managed to be within the range of 0.5~1 whenever possible (Li, 1998;Hasegawa et al., 2000).Both WV-2 pairs possessed good ratios.However, the ZY-3 pair was formed using the nadir and daft images to avoid large blocked area and thus has a relatively low base-height ratio of 0.44, which is close to the lower bond of the recommended range.However, the image matching technique benefited from the improved similarities between the two images with similar viewing angles.

Most active fracturing region in Ronne-Filchner Ice Shelf
Over the period from 2003/2004 and 2014/2015, the detected fractures in RFIS increased from 1168 to 1562, which is partly attributed to the improved resolution of the Landsat-8 satellite.Although the two sets of detected fractures appeared in about the same regions, there were other fractures that occurred in new areas close to the Korff Ice Rise and on Berkner Island (Fig. 4).On the other hand, as the existing fractures propagated forward in the ice shelf front region, some of them were found to have been expanded and merged laterally across suture lines to form longer crevasses.Along the ice shelf frontline there were three additional new large longitudinal rifts detected, including LR 1 (32 km) in RIS and LR 3 (20 km) and LR 4 (5 km) in FIS.Furthermore, two new large traverse rifts, Rift 1 (52 km) and Rift 2 (53 km), were detected in FIS and experienced rapid expansions during the observational period (Fig. 6).The fractures formed fracture bands along the ice flow directions, which were separated by suture lines (Fig. 4).The origins of the fracture bands can be traced back to two types of sources, ice streams and rupture such as islands, ice rises, and rumples.Crevasses originated from ruptures in low ice flow speed areas (e.g., Henry Ice Rise and Berkner Island, Figs. 4 and 5) appeared to remain local.Others initiated from higher speed ruptures and (e.g., Kershaw Ice Rumples, Korff Ice Rise, Doake Ice Rumples) have the potential to reach the ice shelf front.Furthermore, most crevasses that originated from ISs appeared to traverse the entire journey to the sea.
A further observation is that the ice flow speed of the RFIS, more specifically that at its front, changed its pattern such that the speed decreased by less than 40 m a -1 in the western part of the RIS, as illustrated in the ice speed change map in Fig. 11, which was derived using an InSAR dataset from 1997/2000 (Joughin and Padman, 2003) and a Landsat-8 dataset from 2013-2016 (Fahnestock et al., 2016;Scambos, 2017); a speed increase of less than 40 m a -1 occurred in the eastern part of RIS; and then there was a significant increase of over 100 m a -1 at the front of FIS, which was also reported by Scheuchl et al. (2012) in an analysis of ice velocity changes in the RFIS during the shorter period of 1997-2009.Correspondingly, the fractures in the decreased speed area in the western part of the RIS did not appear to show dramatic changes as they propagated forward over the observational period from 2003/2004 to 2014/2015.For example, Rift 3 maintained its length but had a slight growth of 0.4 km in its width.The most noticeable change in the FIS is the occurrences of Rifts 1 and 2, which had opposite speed change trends across the rifts, suggesting uneven longitudinal motions (Fig. 11).A profile across Rift 1 (Fig. 11) reveals that the ice flow speed of 1997/2000 (blue line in Fig. 12) had a rather continuous curve over the entire profile, increasing gradually from upstream to the front.However, the speed during 2013-2016 (red line) stayed slightly below the line of that of 1997/2000 before it reached Rift 1; at Rift 1, the speed jumped by approximately 220 m a -1 , although it is unclear when this difference started.In addition, the thickness distribution of the bottom marine ice layer (Fig. 2) seems to be related to the ice flow speed (Fig. 11) where the ice shelf had thicker marine ice and more stable fractures in low speed areas and vice versa.This is consistent with the role of marine ice in rift termination and suture zone properties (Kulessa et al., 2014;McGrath et al., 2014).Finally, the average rate of melting for FRIS is relatively low (Depoorter at al., 2013;Rignot et al., 2013), but the total melt in the frontal region can be substantial (Moholdt et al.,  The latest major calving event in the FIS occurred in 1986.The earliest documented description of the rift that led to that event was published by Neuburg et al. (1959), depicting the rift as the "Grand Chasm," consisting of two nearly equal parts with a slight offset at the center.A set of parameters of the Grand Chasm were measured from aerial photographs, including rift length, and minimum and maximum widths and depths (Table 2).To better track the change trend, the minimum width is replaced by the width, which was measured consistently along the same ice flow line, as indicated in Fig. 6c, where the flow line and width of Rift 1 in 2015 are illustrated.The shape and location of the Grand Chasm in 1986 is illustrated in Fig. 2, from which the distance from each position of the rift was measured.Additional measurement information and computational details are presented in the Supplement.Table 2. Temporal measurements of Grand Chasm and Rifts 1 and 2 in FIS The set of rifts 1 and 2 and the Grand Chasm share two major commonalities: a) they were both initially two traverse rifts offset from each other at the central part of the ice shelf; and b) they were both located approximately 70 km from the ice shelf frontline.The FIS is now in another cyclic period in which it will 5 eventually advance beyond its confining embayment or pinning points and subsequently retreat by calving into tabular bergs, just as the Grand Chasm did in 1986 (Joughin and MacAyeal, 2005;Kim et al., 2007).
The rift size (lengths and widths) measurements of Rifts 1 and 2 are generally smaller than those of the Grand Chasm, as are their rates of change (Table 2).Thus, Rifts 1 and 2 are in an earlier stage of the cycle.2011-2016 0.05km a -1 0.11km a -1 0.08km a -1 1.05km a -1 59.0m a -1 0.44 m a -1 The Cryosphere Discuss., https://doi.org/10.5194/tc-2017-178Manuscript under review for journal The Cryosphere Discussion started: 11 September 2017 c Author(s) 2017.CC BY 4.0 License.
complex approach and is beyond the scope of this study.

Conclusions
This paper presents the results of a comprehensive study of the fracturing in Ronne-Filchner Ice Shelf, Antarctica, using a new framework of systematic ice shelf fracture mapping and a long record of data dating back to 1957.The new remote sensing capabilities made available in the past decade, including Landsat-8, WV-2, ZY-3 and other datasets, allow us to detect new fractures and others that were not detected previously.Among the newly detected fracture features were Rifts 1 and 2, which are comparable to the Grand Chasm that caused the major calving event of the FIS in 1986.Three-dimensional rift models generated from quasi real-time stereo ZY-3 images revealed insightful topographic information about the large rifts, which can be used to improve the reliability of ice shelf modeling and as support for further analyses of ice shelf stability.
In summary, based on the analysis using data involving fracture changes, propagation speeds, ice shelf melting, and bottom marine ice thickness, the level of fracturing in the RFIS was slightly increased, particularly at its front, from 2003/2004 to 2014/2015.The overall fracture observations do not seem to suggest an immediate significant impact on the stability of the shelf.However, with the detected rapid changes and the 3D measurements of Rifts 1 and 2, the most active fracturing activities occurred at the front of the FIS from 2001 to 2016.A potential upcoming major calving event in the FIS is estimated to occur in 2051.The stability of the ice shelf, particularly the developments of Rifts 1 and 2, should be closely monitored.

Figure 1 .
Figure 1.Study area of Ronne-Filchner Ice Shelf in Antarctica: ice shelf boundary and features (red lines) adopted from Depoorter et al. (2013), Rift 1 and Rift 2 (blue lines) in Filchner Ice Shelf, and ice rises (IR) from MOA were adopted, which was developed based on MODIS (Moderate Resolution Imaging Spectroradiometer) images taken between 2003 and 2004 at a resolution The Cryosphere Discuss., https://doi.org/10.5194/tc-2017-178Manuscript under review for journal The Cryosphere Discussion started: 11 September 2017 c Author(s) 2017.CC BY 4.0 License.
The third dataset is a group of multisource satellite images from 2001-2016 used for the temporal 2D mapping of the length and width changes of Rifts 1 and 2. It includes ASTER (Advanced Spaceborne Thermal Emission and Reflection Radiometer) L1T images, Landsat-7 ETM+ (Enhanced Thematic Mapper Plus) L1GT images, and Landsat-8 OLI L1GT images.Two different types of SAR images are also used as complements of the visible images to address the difficulties caused by clouds and austral 5 winter, including the ASAR IMP (Image Mode Precision) images and Sentinel-1 A EW (Extra Wide Swath) Ground Range Detected Medium C-band images.

Figure 2 .
Figure 2. Details of the data coverage and features of FIS, including Rifts 1 and 2, frontlines, Grand Chasm and longitudinal rifts of the 1986 major calving event.The background is a Landsat-8 OLI scene

Figure 3 .
Figure 3. Framework of systematic fracture mapping and analysis of large Antarctic ice shelves, such as RFIS
established and their propagation patterns and speeds can be derived.During implementation, the identification of two matched fracture features is performed manually because automated algorithms can be confused by repeated similar fracture shapes of candidate fractures among other factors.Subsequently, the fracture propagation speed is estimated by the distances measured between the corresponding points of the two fractures.A high level of effort should be made to avoid the selection of corresponding points which may only appear similar because of shadows, illumination conditions, short term textures and The Cryosphere Discuss., https://doi.org/10.5194/tc-2017-178Manuscript under review for journal The Cryosphere Discussion started: 11 September 2017 c Author(s) 2017.CC BY 4.0 License.
of ice shelf wide fractures in 2014/2015 in Ronne-Filchner Ice Shelf A total of 41 scenes of Landsat-8 (band 8) images of 2014/2015 were used to extract the fracture features in RFIS.A total of 1562 surface fractures were detected and shown in Fig. 4, including 1444 crevasses covered by snow (illustrated in green lines) and 118 rifts with openings visible in the images The Cryosphere Discuss., https://doi.org/10.5194/tc-2017-178Manuscript under review for journal The Cryosphere Discussion started: 11 September 2017 c Author(s) 2017.CC BY 4.0 License.(bluelines).In terms of their spatial distribution, there is a concentration of crevasses in the western half of RIS, which appear to originate from the upstream region and a complex setting of peninsulas, IS, Carlson inlet, and Korff IR; the crevasses then propagated toward the ice shelf front with an average spacing of ~4 km in the upstream area and ~9 km in the ice shelf front area.The crevasses at this time were generally separated by the suture zones (white lines) that originate from and run downstream in peninsulas.The suture zones are capable of binding together neighboring ice flows and are typically composed of ice mé lange consisting of snow, sea ice, and ice blocks(McGrath et al., 2014).In the eastern half of RIS, there are two groups of fractures, one in the shear margin of Berkner Island and the other downstream of Henry Ice Rise.The smaller ice shelf, FIS, has a concentrated set of crevasses in its eastern half.In addition, there are crevasses on the east side of Berkner Island, where two small local streams are located.
are shown as yellow lines in Fig. 4. The dataset contains a total of 1168 fracture features.Thus, the total detected fracture features increase from 1168 in 2003/2004 to 1562 in 2014/2015.This increase may be attributed to the improvement of fracture detection by moving from the MOA mosaic (resolution 125 m) to Landsat-8 (resolution 15 m) images; on the other hand, the ice shelf itself is gradually approaching its next major calving event, as it has a cycle of decades long.Some longitudinal rifts seen in 2003/2004 disappeared in the regular calving process, but other new ones occurred in 2014/2015.The large longitudinal rift in RIS, LR 1, was 32 km long in the Landsat-8 image, but appeared only as a set of disconnecting crevasses in the MOA mosaic.LR 2 consisted of two connecting rifts cut inland into the ice shelf by approximately 34 km, with an increase of 19 km.Two new longitudinal rifts, LR 3 (20 km) and LR 4 (5 km) in FIS, are detected in the 2014/2015 dataset.The Cryosphere Discuss., https://doi.org/10.5194/tc-2017-178Manuscript under review for journal The Cryosphere Discussion started: 11 September 2017 c Author(s) 2017.CC BY 4.0 License.

Figure 4 .
Figure 4. Fractures of Ronne-Filchner Ice Shelf: crevasses (green lines), rifts (blue lines), and ice shelf front lines (orange line) mapped from the Landsat-8 band 8 images of 2014/2015 in this study: crevasses and rifts (yellow lines), suture zone lines (thin white lines), and ice shelf front line (red line) mapped from the MOA mosaic of 2003/2004 (Hulbe et al., 2010); and MOA mosaic as the background

10and
Rift 2 (53 km) were not in the MOA fracture dataset.As indicated in the later sections, they occurred recently and experienced rapid changes.They may potentially lead to major calving events.The results of the temporal and three-dimensional sensing and propagation analyses of Rifts 1 and 2 are presented The Cryosphere Discuss., https://doi.org/10.5194/tc-2017-178Manuscript under review for journal The Cryosphere Discussion started: 11 September 2017 c Author(s) 2017.CC BY 4.0 License.later in the paper.4.2 Fracture propagation speed from 2003/2004 to 2014/2015 in Ronne-Filchner Ice ShelfThe fracture propagation speed of RFIS was estimated using the fracture features from the Landsat-8 images of 2014/2015 and the MOA fractures of 2003/2004(Hulbe et al., 2010).A total of 651 corresponding points between the two sets of fracture features were identified and measured.Their corresponding vectors represent the average advection speeds and direction information and are illustrated as arrows in Fig.5.The speed in RIS increases from ~700 m a -1 in the upstream region to ~1000 m a -1 in the central area, and further to ~1400 m a -1 at the front of the ice shelf.In FIS, the fracture propagation speed increases from ~600 m a -1 in the central part to ~1200 m a -1 in the front.The estimated propagation speeds of the detected and matched fractures in the entire ice shelf ofRFIS from 2003RFIS from  /2004RFIS from   to 2014RFIS from  /2015

Figure 5 .
Figure 5. Fracture propagation vectors of RFIS estimated using the fracture features from the Landsat-8 images of 2014/2015 and the MOA mosaic of 2003/2004 (Hulbe et al., 2010) with a background of an ice flow speed map derived from Landsat-8 images of 2013-2016(Fahnestock et al., 2016;Scambos, 2017)

TheFigure 7 .
Figure 7. 3D mapping of Rift 1 using ZY-3 stereo images: (a) color coded rift DEM with locations of profiles and two ICESat tracks, and (b) fly through screen inside Rift 1 (viewing from Profile 1 toward Profile 2, depth exaggeration = 10 times) Fig. 9 to show the cross sections of the rift.Based on their profiles, the rift walls are almost vertical; the average rift depth is approximately 50 m, and the ice mé lange rises approximately 10 m from the rift base.The Cryosphere Discuss., https://doi.org/10.5194/tc-2017-178Manuscript under review for journal The Cryosphere Discussion started: 11 September 2017 c Author(s) 2017.CC BY 4.0 License.

Figure 8 . 5 Figure 9 .
Figure 8. Key parameters of Rift 1 derived from the ZY-3 DEM (February 28, 2014): rift width, ice shelf surface elevation, rift bottom elevation, depth, and bottom roughness measured along the centerline starting from the west end marked with S in Fig. 7a, all of which are based on the WGS 84 ellipsoid datum.

Figure 10 .
Figure 10.Temporal changes of Rift 1 along Profile 4 defined in Fig. 7a: (a) ICESat points of Track 1303 from October 2003 to October 2008, (b) ZY-3 and WV-2 DEMs along Profile 4 from December 2012 to February 2016, and (c) ICESat points of Track 0085 from October 2003 to March 2009.

The
Cryosphere Discuss., https://doi.org/10.5194/tc-2017-178Manuscript under review for journal The Cryosphere Discussion started: 11 September 2017 c Author(s) 2017.CC BY 4.0 License.deviation of 31 m, which is about two pixels of the LIMA mosaic (15 m).The relative accuracies of the extracted fractures from the Landsat-8 OLI images are less than one pixel (15 m).Thus, the overall accuracy of the Landsat-8 OLI fractures is 34 m.The fracture changes presented in Fig.4were derived from the Landsat-8 OLI fractures of 2014/2015 and the MOA fractures of 2003/2004 that did not have a reported accuracy.We used the MOA geolocation accuracy of 125 in Table1and 34 m for the Landsat-8 OLI fractures.Thus, the average fracture propagation speed over the studied 11 years has an accuracy of 12 m a -1 .The fracture propagation speed at each fracture feature was compared with the Landsat-8 speed map of 2013-2016 in Fig.4.The result showed an average difference of 5 m a -1 and a standard deviation of 27 m a -1 .The comparison result of the two data sets presents the difference in their speeds, one with an average of approximately 11 years and the other of approximately 3 years, both of which showed a high level of consistency.

The
Cryosphere Discuss., https://doi.org/10.5194/tc-2017-178Manuscript under review for journal The Cryosphere Discussion started: 11 September 2017 c Author(s) 2017.CC BY 4.0 License.Finally, mapping the area inside the rift is affected by two additional factors: a) the sun shadow cast by the northern edge of the rift and b) the floor area on the forward-looking image blocked by the northern edge of the rift.Given the local time of 07:35:16 on February 28, 2014, an average rift depth of 47 m and a rift azimuth angle of 89.5°, the length of the shadow across the rift is calculated to be 82 m.The actual measured shadow length across the rift is 42~114 m.Furthermore, given the forward-looking angle of 23.5°, the blocked area across the rift floor should be 15 m long.The actual measured blocked area width is 12~19 m.The affected areas of these two factors overlap with each other; whichever is greater defines

Figure 12 .
Figure 12.Ice flow speed of 1997/2000 (blue line) and 2013-2016 (red line) and their differences along a profile A-A' across Rift 1 in the FIS as indicated in Fig. 11

Table 1 . Summary of datasets used in this study ata and type Acquisition time (UTC) Quantity (scene) Spatial resolution Base- height ratio Geolocation accuracy
The Cryosphere Discuss., https://doi.org/10.5194/tc-2017-178Manuscript under review for journal The Cryosphere Discussion started: 11 September 2017 c Author(s) 2017.CC BY 4.0 License.