An Early Warning Sign of Critical Transition in The Antarctic Ice Sheet. A New Data Driven Tool for Spatiotemporal Tipping Point

In this paper, we introduce a new tool for data-driven discovery of early warning signs of critical transitions in ice shelves, from remote sensing data. Our approach adopts principles of directed spectral clustering methodology considering an asymmetric affinity matrix and the associated directed graph Laplacian. We applied our approach generally to both reprocessed the ice velocity data, and also remote sensing satellite images of the Larsen C ice shelf. Our results allow us to (post-cast) predict fault lines responsible for the critical transitions leading to the break up of the Larsen C ice shelf crack, which resulted in the A68 iceberg, and we are able to do so, months earlier before the actual occurrence, and also much earlier than any other previously available methodology, in particular those based on interferometry.


Introduction
Warming associated with global climate change causes global sea level to rise [18]. Three major reasons for this are oceans expansion [17], ice sheets lose ice faster than it forms from snowfall, and also glaciers at higher altitudes melt. During the 20 th century, sea level rise has been dominated by the retreat of glaciers. Still, this contribution starts to change in the 21 st century because of the ice shelves cracks. Ice sheets store most of the land ice (99.5%) [18], with a sea-level equivalent (SLE) of 7.4m for Greenland and 58.3m for Antarctica. Ice sheets form in areas where the snow that falls in winter does not melt entirely over the summer. Over thousands of years of this effect, the layers grow thicker and denser as the weight of new snow and ice layers compresses the older layers.
Ice sheets are always in motion, slowly flowing downhill under their own weight. Near the coast, most of the ice moves through relatively fast-moving outlets called ice streams, glaciers, and ice shelves. When a marine ice sheet accumulates mass of snow and ice at the same rate as it loses mass to the sea, it remains stable. Most of Antarctica has yet to see dramatic warming. However, the Antarctic Peninsula, which juts out into relatively warmer waters north of Antarctica, has warmed 2.5 degrees Celsius (4.5 degrees Fahrenheit) since 1950 [28]. A large area of the Western Antarctic Ice Sheet is also losing mass, probably due to warmer water up-welling from deeper ocean near the Antarctic coast. In Eastern Antarctica, no clear trend has emerged, although some stations report slight cooling. Overall, scientists believe that Antarctica is starting to lose ice [28], but so far, the process is not considered comparably fast as the widespread changes attribute in Greenland [28].
Geologic record of the icing of Antarctica reveals beginnings in the middle Eocene epoch, about 45.5 million years ago [16], and escalated during the EoceneOligocene extinction era, event about 34 million years ago. However, the Western Antarctic ice sheet declined somewhat during the warm early Pliocene epoch, Figure 1: A-68 iceberg. The fractured berg and shelf are visible in these images, acquired on July 21, 2017, by the Thermal Infrared Sensor (TIRS) on the Landsat 8 satellite. Credit: NASA Earth Observatory images by Jesse Allen, using Landsat data from the U.S. Geological Survey. approximately 5 to 3 million years ago. During this time, the Ross Sea opened, [16,24], but there was no significant decline in the land-based Eastern Antarctic ice sheet.
Since 1957, modern record of the continent-wide average reveals a surface temperature trend of Antarctica that has been positive and significant at > 0.05 • C/decade [23,8]. Western Antarctica has warmed by more than 0.1 • C/decade in the last 50 years, and this warming is most active during the winter and spring. Although this is partly offset by autumn cooling in Eastern Antarctica, this effect is prevalent to the 1980s and 1990s [23].
Of particular interest to us in this presentation, the Larsen Ice Shelf extends as a ribbon of ice shelf, down from the East Coast of the Antarctic Peninsula, from James Ross Island to the Ronne Ice Shelf. See Fig. 9. In fact, it consists of several distinct ice shelves, separated by headlands. The major Larsen C ice crack was already noted to have started in 2010 [13], but it was initially very slowly evolving, and there were no signs of radical changes according to Interferometry processing of the remote sensing imagry [11]. However, since October 2015, the major ice crack of Larsen C has been growing much faster, until the point more recently that it finally failed, resulting in calving the massive A68 iceberg. See Fig. 1; this is the largest known iceberg, with an area of more than 2,000 square miles, or nearly the size of Delaware. In summary, A68 detached from one of the largest floating ice shelves in Antarctica, and floated off in the Weddell Sea. Interestingly, two and a half years later, it remains largely intact, and has finally drifted from the near Antarctica seas into the more turbulent open Artic Ocean where it is expected to break apart more quickly.
In [9], the authors presented a structural glaciological description of the system, and subsequent analysis of surface morphological features of the Larsen C ice shelf as seen from satellite images spanning the period 19632007. The results and conclusions of the research stated that: Surface velocity data integrated from the grounding line to the calving front along a central flow line of the ice shelf indicate that the residence time of ice (ignoring basal melt and surface accumulation) is 560 years. Based on the distribution of ice-shelf structures and their change over time, we infer that the ice shelf is likely to be a relatively stable feature and that it has existed in its present configuration for at least this length of time..
In [12], the authors modeled the flow of the Larsen C and northernmost Larsen D ice shelves using a model of continuum mechanics of the ice flow. They applied a fracture criterion to the simulated velocities to investigate the ice shelfs stability. The conclusion of that analysis shows that the Larsen C ice shelf is Figure 2: Directed partitioning method. We see the image sequence to the left, and to the right, we reshape each image as a single column vector. Following the resultant trajectories, we see that the pairwise distance between the two matrices will result in an asymmetric matrix.
inferred to be stable in its current dynamic regime. This work published in 2010, and in the same year that the Larsen C ice crack already existed, but for its slow-growing rate according to analytic studies. There were no expectations at that time for the fast-growing and collapse that happened for Larsen C.
Interferometry has traditionally been a main technique to analyze and predict ice cracks based on remote sensing. Interferometry [2,15], is based on a family of techniques in which waves, usually electromagnetic waves, are superimposed, causing the phenomenon of interference patterns, which in turn are used to extract information concerning the underlying viewed materials. In fact, interferometers are widely used across science and industry for the measurement of small displacements, refractive index changes, and surface irregularities, and so it is considerd a robust and fimilar tool that is successful in the macro-scale application of monitoring structural health of the ice shelves. So it is our job here to contrast our methodology to interferometry. Here we will take a data driven approach directly fromt he remote sensing imagry, to infer structural changes of the impending tipping point to the critical transition of the break of the Larcen C. Fig. 10 shows the interferometry image as of April 20, 2017, and although it clearly shows the crack that already existed, but may provide no information or forecasting powers indicating what can happen next. In fact, just a couple of weeks after the image shown in Fig. 10, the Larsen C ice crack changed significantly, and took a different dynamic that quickly thereafter divided into two branches, as shown in Fig. 11. Our methods as we will show achieve a much more successful and earlier data driven indicator of this important outcome.

Directed Partitioning
In our previous work [1], we developed the method of Directed Affinity Segmentation (DAS), for which we showed high performance in successsfully detecting coherent structures in fluidic systems, observed from "movie data" and without the need for the intermediate stage of finding the vector field responsible for underlying advection.
Two of the most commonly used and successful image segmentation methods are based on 1) the kmeans [14], and 2) spectral segmentation [20], respectively. However, wile these were developed successfully for static images, these methods need major adjustments for successful application to sequences of image, for the spatiotemporal probleem of motion segmentation associated with coherence, despite that traditionally they are considered well suited to static images [22]. The key difference is what underlies a notion of coherent observations, that we must also understand directionality associated with the arrow of time.
An affinity measure is the phrasing for a comparison, or cost, between states, and as such a loss function of some kind of often the starting point for many algorithms in machine learning. However, when there is an underlying arrow of time, the loss functions that most naturally arise when tracking coherence are inherently not symmetric. Correspondingly, affinity matrices associate the affinity measure for each pairwise comparison across a finite data set. Also it is useful to consider the undirected graph associated with the affinity matrix, where there is an edge between each state for which there is a nonzero affinity, and generally in the symmetric case these graphs are undirected. Now consider that if the affinity matrices are not symmetric, then these are associated with directed graphs. This is a theoretical complication to standard methology since much of the theoretical underpinnings of standard spectral partitioning assumes a symmetric matrix corresponding to an undirected graph and then considers the spectrum of its corresponding symmetric Laplacian matrix that follows. This can all be accomodated by methods considering the spectral theory of a graph Laplacian for a weighted directed graphs, built upon the theoretical work of F. Chung [5], and as we built upon in [1].
Before proceeding with our directed partitioning method, we formulate the (movie) imagry data set as the following matrices; where each X i is the i th image (or the image at i th time step) reformed as a column vector, See Fig. 2, τ is the time delay, X 0 and X τ are the images sequences stacked as column vectors with a time delay at the current and future times respectively. Choosing the value of the time delay τ , can results in significant differences in the segmentation process. Consider that in the case of a relatively slowly evolvoing dynamical system, where the change between two consecutive images is not significantly distinguishable, then choosing a large value for τ may be a better suited. Note that the rows of X 0 , X τ ∈ R d×T −τ represent the the change of the color of the pixel at a fixed spatial location z i . Then, we introduced [1] an affinity matrix in terms of a pairwise distance function between the pixels i and j as, where S : R 2 → R is the spatial distance between z i and z j , and C : R T −τ × R T −τ × R → R is a distance function describing "color distance" the i th and the j th color channels. The parameter α ≥ 0 regularizes balancing these two effects. In this work, we choose the functions S and C each to be L 2 -distances, and We see that the spatial distance matrix S is symmetric, however, the color distance matrix C is asymmetric for all τ > 0. Then, while the matrix generated by C(X 0 i , X τ j , 0) refers to the symmetric case of spectral clustering approaches, we see that the matrix given by C(X 0 i , X τ j , τ ), τ > 0 implies an asymmetric cost naturally due to the directionality of the arrow of time. Thus we require an asymmetric clustering approach should be adopted.
First we define our affinity matrix from Eq. 3 as, This has the effect that both spatial and measured (color) effects have "almost" Markov properties, as far field effects are almost "forgotten" in the sense that they are almost zero, and near field values are largest. Notice we have suppressed including all the parameters in writing W i,j , and that besides time parameter τ that serve as sampling and history parameters, together the parameters α and σ serve to balance spatial scale and resolution of color histories. We proceed to cluster the spatiotemporal regions of the system, in terms of the directed affinity W by interpreting the problem as random walks through the weighted directed graph, G = (V, E) designed by W as a weighted adjacency matrix. Let, where is the degree matrix, and P is a row stochastic matrix representing probabilities of a Markov chain through the directed graph G. Note that P is row stochastic implies that it row sums to one. This is equivalently stated that the right eigenvector is the ones vector, P1 = 1, but the left eigenvector corresponding to left eigenvalue 1 represents the steady state row vector of the long term distribution, which for example if P is irreducible, then u = (u 1 , u 2 , ..., u pq ) has all positive entries, u j > 0 for all j, or say for simplicity u > 0. Let Π be the corresponding diagonal matrix, and likewise, which is well defined for either ± sign branch when u > 0. Then, we may cluster the directed graph by concepts of spectral graph theory for directed graphs, following the weighted directed graph Laplacian described by Fan Chung [4], and a similar computation has been used for transfer operators in [7,10] and as reviewed [3]. The Laplacian of the directed graph G is defined, [4], The the first smallest eigenvalue larger than zero, λ 2 > 0 such that, allows a bi-partition, by, by sign structure. Analogously to the Ng-Jordan-Weiss symmetric spectral image partition method [20], the first k eigenvalues larger than zero, and their eigenvectors, can used to associate a multi-part partition, by assistance of k-means clustering these eigenvectors.

Results
Now we apply the above directed affinity segmentation to satellite images of Larsen C ice shelf and ice surface velocity data as follows. Here w show that the directed affinity segmentation of spatiotemporal changes can work as an early warning sign tool for critical transition in marine ice sheets. We will apply our post-casting experiments on images of Larsen C before the splitting of the A68 iceberg, and then we will compare our forecasting based on segmentation results to the actual unfolding of event. In Fig. 3 we see different snapshots of the ice surface velocity data set [29,21,19], which is part of the NASA Making Earth System Data Records for Use in Research Environments (MEaSUREs) Program, and it provides the first comprehensive [29], high-resolution, digital mosaics of ice motion in Antarctica assembled from multiple satellite interferometric synthetic-aperture radar systems. We apply the directed affinity partitioning to the available data set, and the results are shown as a labeled image in Fig. 4.
As shown on Fig. 4 we note the following: • The data collected from eight different sources [29,26], with different coverage and different error range, and interpolating the data from different sources explains the smooth curves in segmentation around the region on interest. Result of the directed partitioning is shown in Fig. 4. Source of data: [29]. to the region of interest shows large variations of ice surface velocity within a small area, to give a clearer focused view of the differences in speeds. In Appendix, Fig. 12 shows the surface plot for the same result. Figure 5: Two coherent sets dynamic. As the inner set contact the boundaries of the outer one, than give the chance for a new reactions that "may" cause critical transition.
• The directed partitioning shows the Larsen C ice shelf as a nested set of coherent structures that contain successively within each other.
• The zoom picture highlight shown in the right of Fig. 4 shows the region where the Larsen C ice crack starts. Furthermore, we see that within narrow spatial distance (4 miles) there is a large change of velocity. More precisely, the outer boundaries of the different coherent sets become spatially very close (considering the margin of error in the measurements [26]. We conclude a high probability that they contact).
Directed partitioning gives us informative clustering, meaning that each cluster has homogeneous properties, such as the magnitude and the direction of the velocity. In general, for the coherent sets are nested, A 1 ⊂ A 2 ⊂ ... ⊂ A n , physically, as each set A i−1 keeps its coherence within A i because of a set of properties (i.e., chemical or mechanical properties) that rules the interaction between them. However, observe that the contact between the boundaries of the sets A i−1 and A i , see Fig. 5, means a direct interaction between A i−1 and A i+1 .
In the case of including the ice velocity when partitioning, we collect these observations and discuss them here. However, since the sets boundaries are not completely contacted, and the direction of the velocities reveal no critical changes; we believe this results implicitly from the nature of the data preprocessing that includes interpolation and smoothing of the measurements. We state nothing more than such close interaction between coherent sets boundaries can be an early warning sign that should be considered and investigated by applying "what if" assumptions and analyzing the consequences from any change or any error in the measured data.
As a matter of declaring success of our approach over standard methodology, observe that our directed partitioning method achieves better results using the remote sensing satellite images [27] as in contrast to the standard interferometry concept. Further, to reduce the obscuration effects of noise (clouds and image variable intensity), we used the averaged images, over one month, as a single snapshot for the directed affinity constructions. Fig. 6, the directed affinity partitioning for two time windows starts from December 2015. Notice that the directed partitioning begins to detect the significant change in the Larsen C ice shelf on July 2016. In Fig. 7 we see that by September 2016, we detect a structure very close in shape to the eventual and actual iceberg A-68, which calved from Larsen C on July 2017. Moreover, by November 2016, see Fig. 13, the boundaries of the detected partitions match the crack dividing into two branches that happened in later in May 2017 and shown in Fig. 11.

Discussion
Here we have presented a new approach for predicting possible critical transitions in spatiotemporal systems, specifically marine ice sheets, based on remote sensing satellite imagery. Our approach shows reliability in detecting coherent structures, and when the object of concern is a rigid body such as ice sheets. The main idea  Fig. 7. We see that during 2016, there was no significant change in Larsen C crack at the beginning of the year. In July 2016, the directed affinity segmentation propose a large change in the crack dynamics, and this change keeps going faster as Fig. 7 shows.  is that observing a significant and perhaps topological form change of a coherent structure, over time, may be indicative of an important underlying critical structural change of the ice. The computational approach is based on spectral graph theory in terms of the directed graph directed graph Laplacian. In the case of the Larsen C ice shelf this is born out and we successfully observe the calving of the A68 iceberg months before the primary competing method based on interferometry. This transition of the coherent structure can indicate a possible fracture along the edges of directed affinity partitioning. We see that the directed affinity partitioning can be a useful early warning sign that indicates the possibility of critical spatiotemporal transitions, and it can help to bring the attention for specific regions to investigate different possible scenarios in the analytic study, whether computational, or possibly even supporting further field studies and deployed aerial remote sensing missions.
In our future work, we plan to pursue the idea of connecting our data driven approach of computing boundaries by directed partitioning, with the computational science approach in terms of stress/strain analysis of rigid bodies and an understanding of hte underlying physics. In addition to expressing the risk of the possible critical transitions of multiple coherent structures that surround each other, see Fig. 5, in terms of Lyapunov exponent of the minimum distance between two evolving shape coherence curvatures that surround each other.

Acknowledgments
This work was funded in part by the Army Research Office, the Naval Research Office, and also DARPA. Figure 9: Location of ice shelves on the Antarctic Peninsula. Source [6]. We can measure the iceberg crack propagation much more accurately when using the precise surface deformation information from an interferogram like this, rather than the amplitude (or black and white image) alone where the crack may not always be visible. Source [25].  Figure 12: Directed affinity partitions with the mean velocity (speed) of the partition assigned for each label entries. The spatial distance between the arrows tips is less than two miles, while the difference in the speed is more than 200 m/year.