12 Years of Area Variation by the Drygalski Ice Tongue as Measured With COSMO-SkyMed

One important coastal polynya around Antarctica is Terra Nova Bay (TNB) polynya. Its formation and persistence are due to the combined effect of katabatic winds, regional ice conditions, and the Drygalski Ice Tongue (DIT). The combined effect of these elements arranges a delicate balance in TNB. To evaluate the interacting elements, we focused our attention on the DIT, one of the largest ice tongues in Antarctica. The DIT is a 70-km long tongue, which juts out like a pier from the icy land into northern McMurdo Sound of Antarctica's Ross Dependency. The analysis presented in this article was carried out using Cosmo-SkyMed (CSK) images. The methodology developed for the analysis consists of three steps: 1) DIT shape enhancement, 2) fuzzy Bayesian segmentation which associates the analyzed set of pixels with a binary label field, and 3) parametric representation of the shape of the DIT to compare the area of the ice tongue as it was 12 years ago and as it is today. The computed mean value of the DIT surface velocity was 703 m y<inline-formula><tex-math notation="LaTeX">$^{-1}$</tex-math></inline-formula>. To estimate the closeness of the derived parameter to the true value, using a statistical approach, the 95% confidence interval for the true mean value is 703<inline-formula><tex-math notation="LaTeX">$\pm$</tex-math></inline-formula>4.236 m y<inline-formula><tex-math notation="LaTeX">$^{-1}$</tex-math></inline-formula>. Our analysis confirms a trend of advance. That is, after the significant calving events of the tongue in 2005 and 2006, the DIT has resumed its normal growth, and its length has increased by 12% in the last 12 years.


I. INTRODUCTION
S ATELLITE imagery is a practical alternative for tracking changes in glaciers over time. Chronologically, the first remote observation of changes in the East Antarctica glaciers occurred by Landsat 1 data acquired in late 1972, a few months after the launch of this satellite [1]. In this early remote sensing era, in 1976, Landsat imagery allowed mapping of the Victoria Land coast [2], and the velocity estimation of DIT using sequences of Landsat images from the 1970s and 1980s [3]. Preceding the civilian satellite missions, Argon was a satellite-borne camera system for reconnaissance of earth resources [4], [5]. To collect spaceborne photographs over the Antarctic continent, three missions were completed between 1962 and 1963 [6]. This strategic collection was released in 1995 and is known as the Declassified Intelligence Satellite Photography (DISP) dataset. Novel Antarctic products of this historical material are a mosaic to delineate the 1963 coastal margin [7] and combining Synthetic Aperture Radar (SAR) imagery of 1997, a study of the variability over time of the ice shelves of the Atlantic coast [8].
For the observation of the David Glacier (DG)/Drygalski Ice Tongue (DIT), a variety of data have been used, mainly optical and SAR imagery, interferometry data, in-situ measurements of flow and stratification and complementary ancillary information. Long term studies of glaciers flowing into the western Ross Sea have also incorporated digital scans of historical maps and a collection of antique aerial photographs and contemporary remote sensing imagery. Involving outlet glaciers, a review of remote sensing of ice motion in Antarctica is provided in [9].
To detect changes in the ice tongues, a natural option is to identify and track distinctive features. In [10], the ice velocity at the terminus of Columbia Glacier, Alaska, was estimated by the joint use of sequential photography, known landmark points, vertical photogrammetry techniques, topographic maps, and aerial photography. Using sequential Landsat images, the displacement of common features, such as crevasses, was the basis for calculating the average surface velocity at the Thwaites Glacier (TG)/iceberg tongue [11] and the Smith Glacier ice tongue [12]. Feature-tracking technique on ERS amplitude images and SAR interferograms was accomplished in [13], and the derived velocity variations of TG tongue were modelled by rigid-body rotation and spatially inhomogeneous divergence influence. Global Positioning System (GPS) allows continuous time monitoring. In [14], GPS stations placed around the Mertz Glacier ice tongue were used to investigate the glacier dynamics and its potential implication in rift propagation. In [15], GPS were employed to monitoring flow velocity and direction in the Shirase Glacier tongue, East Antarctica. Cross-correlation matching is an alternative for tracking features between pairs of images. It was applied to Sentinel-1 SAR imagery and aerial photographs to derive the surface velocities of Parker Ice Tongue [16], and using Radarsat-1 ScanSAR images, to estimate the motion field of the Mertz Glacier Tongue ice surface [17]. In another study, direct ground observations of a research expedition on Dalk Glacier allowed a precise estimate of ice flow velocity [18]. By Sentinel-2 multiband imagery at 10 m spatial resolution, a two-phase deep learning method [19] was adapted to delineate marine-terminating outlet glaciers in Greenland.
The dynamical behavior of glacier termini and ice-shelf fronts are linked to climate and ocean interactions, internal dynamic This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ forces, the bed topography, and the inherent cyclic behavior of calving [20]. The ice tongues are also affected on its entire surface by atmospheric and ocean forces, ice thickness, and the melting at the grounding zone. Iceberg collisions are also an important factor of calving. The crevasses can also indicate a pattern of eventual calving events.
The DIT, which is a floating extension of the land-based DG, is a 70-km long tongue which juts out like a pier from the coastline into western Ross Sea. DG is the largest outlet glacier in Victoria Land. Ice flowing from DG comes from East Antarctica. Ice tongues, like the DIT, normally change in size and shape as waves and storms weaken their ends and sides, breaking off pieces which float as icebergs. Regarding the calving phenomena, in 1957, as a result of a severe storm, a major calving event was observed and the DIT lost about 675 km 2 of its extent. In February 2005, two pieces from the north side of the snout calved off. By Advanced Synthetic Aperture Radar (ASAR), the estimated area of icebergs released was about 74 km 2 [21]. In April 2005, the drifting iceberg B-15 A collided with the ice tongue, reducing the tip extent by about 64.5 km 2 . In 2000, calving of the northern ice front of the Ross Ice Shelf produced the C-16 iceberg. In March 2006, strong katabatic winds blowing down from DG drove the C-16 iceberg to hits the tip of the DIT and caused the calving of iceberg C-25 [22]. The computed area of C-25 iceberg was about 143 km 2 .
The occurrence and persistence of TNB polynya are assumed to be due to the combined effect of katabatic winds and regional ice conditions driven by fluctuations in both air temperature and wind speed [23]. Polynyas are able to induce regional ocean-atmosphere heat exchanges in winter. Some modeling studies [24] have linked the significance of polynyas for preserving oceanic thermohaline circulation and the likely future impact to climate variability [25]. Interacting with local coastal currents, the DIT is an important oceanographic element [26]. Besides, it is a central regulator of the Terra Nova Bay (TNB) polynya. In 1957, after an important calving event, the average size of the polynya was reduced of about 33%, and this highly impacted the sea-ice production in the Ross Sea. It is hence important to examine the progression and stability of the ice tongue. The dynamics of TNB winter polynya are also influenced by the DIT, which obstructs northward ice drift into TNB [27], [28], [29].
In an attempt to evaluate the interacting elements in TNB, we focus our attention on observing the area evolution of the DIT, one of the largest ice tongue in Antarctica.

A. Image Data Collection
Thanks to its all-weather capability, day/night acquisitions and observation even in the winter darkness, SAR imagery is well-suited for observing polar regions from space. The study was carried out using Cosmo-SkyMed (CSK) HugeScan images over TNB: three of them from the catalogue of the Italian Space Agency (ASI) and one from 2020. Four CSK images, two from 2008, one from 2012, one from 2020 and one Sentinel-1 SAR image from 2016 were selected over the area of interest. The Sentinel-1 image was retrieved from the data hub of the European Space Agency's Copernicus Programme. The five images were remapped onto a regular polar stereographic projection reference (pixel size = 50 m) using the commercial software TeraScan by SeaSpace Corp.
Three of the major outlet glaciers of Victoria Land, David, Reeves, and Priestley, flow into TNB. The Nansen Ice Sheet is the seaward extension of Reeves and Priestley Glaciers. Fig. 1 shows a geo-referenced CSK HugeScan image of TNB. To the north, the scene is delimited by the Cape of Washington, to the west by the Prince Albert Mountains and the Nansen Ice Sheet, and to the south by the DIT. DG flows through the Prince Albert Mountains.
Meteorological and geophysical variables influence tongue morphology. Ice strain is a significant element in the initial formation of fractures, and ice thickness regulates fracture propagation trend [30]. Ice is thinnest at the tip nose, and this region reacts more rapidly than thicker shelves to a temperature change [31]. Fast-ice extent plays a role in the marineterminating glacier stability. An inverse correlation between fast ice extent and ice tongue velocity was reported in [16]. Some important references on ice dynamics, glacier variations, and ice shelf environments are [32], [33], [34], [35], [36], but a study of these multidisciplinary issues is outside the scope of this article. The goal of our study is to follow the area evolution of DIT.
In the analyzed scenes, land/ice, and open-sea are contrasted regions, which means that a noncomplex segmentation is required. The first step of the analysis is shape enhancement in the selected DIT windows. This is accomplished using a nonlinear filtering technique. The second step is the segmentation of the DIT windows using fuzzy Bayesian labeling. In the final step, the ice shape and extent are computed.

B. DIT Shape Enhancement
The backscattered SAR signal of the earth surface elements is a complex random process. Before the application of the DIT shape detector, a first task is the reduction of the degradation caused by speckle. An important requirement for accomplishing SAR filtering is to avoid the degradation of the image attributes; in our case, the main attribute is the edge boundary distribution. To achieve this attribute, the speckle reducing anisotropic diffusion (SRAD) filter [37] was used. The SRAD filter makes use of a local gradient operator to drive the diffusion process. Nonlinear diffusion provides low-pass filtering on homogeneous regions. The diffusion is reduced in perpendicular directions of high gradient elements, and, as it is carried out all along a detected edge, it preserves the edge features.

C. Fuzzy Bayesian Segmentation
The analyzed images display a raw bimodal distribution of pixel gray level, and this is a favorable condition for labeling the scene elements. To achieve a more automated procedure, our choice was to apply an unsupervised clustering segmentation.
Given the shape of the DIT, the next step is to assign each SAR pixel X =x to a set of mutually exclusive label field {W = lk : k = 1, 2}, where k is the number of labels or clusters. The Bayes formulation takes into account the image process information for driving a labeling process. The posterior distribution is where P (X|W, θ) is assumed as the mixture probability law of the gray-level distribution, θ is the parameter set of W , P (W ) is the marginal probability of the label field, and P (X) is the probability of observing the input image, X, which remains invariant during the process. With a maximum-a-posteriori (MAP) estimate, equation (1) can be approached by P (W |X) ∝ max W {P (X|W, Θ)P (W )}. The data under analysis show a backscatter signal from land/ice elements much brighter than from open sea regions. Thus, a binary segmentation is required (k= 2), where the estimation of θ of (1) is an implicit step. The probability density low of the numerator of (1) must be estimated from the available data. The mixture probability law P (X|W, θ) is the probability occurrence of the pixel gray-level for each class. The analyzed scene is composed of two selected features: land/ice and open sea pixels, hence it is expected to approach P (X|W, θ) by a bimodal cluster function. The membership function P (W ) is also a two-class array, which is the prior knowledge of the input data (the probability terms must be known before starting the process). The histogram modes are connected, and a given pixel can fit in both binary clusters. Hence, for building P (X|W, θ), it is desirable to incorporate a level of uncertainty into the decision criteria. Our choice was to approach the likelihood distribution P (X|W, θ) by the fuzzy c-Means (FCM) [38], [39] algorithm. The FCM cost function is where d(x i , Θk) is a similarity metric (the Euclidean distance), q is the fuzzifier parameter, and u ik is the membership The membership values can now be quantified in terms of probability, P (X|W, θ)∝ u k . For performing the fuzzy clustering, the set of pixels is projected in the u k functions, and a value close to unity corresponds to a high probability of selecting the cluster k. In a graphical representation, u k simulates the relationship between the gray-level distribution of the scene objects and the classconditional probability term P (X|W, θ) of each labeled class l k ∈ W . To tackle the prior probability knowledge, a k-Means clustering procedure provides an estimation of the P (X) term. Thus, (1) is solved by P (W |X) = max{u k P (X)}.
To assess the accuracy of the algorithms, a comparison with two canonical segmentation methods was performed. Fig. 2 shows: (a) a raw input sub-window, (b) the application of a pixel-based segmentation, the k-Means algorithm, (c) a deterministic relaxation, the Iterated Conditional Modes algorithm (ICM) [40], which approaches a maximum likelihood estimator by sequentially updating local label configurations, and (d) the result obtained by the proposed scheme. ICM requires an initial label configuration, which was taken from the k-Means result. ICM provides a more homogeneous result than k-Means, but by the proposed scheme, the iceberg form is fully detected and can be extracted from the binary field.

III. RESULTS
The diffusion process of the SRAD filter was performed in 40 iterations, and to overcome the rift contrast variation, the local gradient value was fixed to 20. The fuzzy binary segmentation was implemented with k = 2, q = 2 and d(x i , Θk) = x i − Θk 2 . Fuzzy membership arises with q > 1. A minimal distance clustering convergence is attained in less than 50 iterations, when no further changes are made in u k . Once the parameters are configured, the joint use of the SRAD filter and the Bayesian scheme define an unsupervised clustering segmentation.  From the full geo-referenced CSK HugeScan scenes, windows were extracted, 1000×2000 pixel size. The result of the segmentation for the first and the last of the DIT image sequence is shown in Fig. 3. The geographical DIT grounding line was not visible in the image of 26 November 2012. The partial view on the west side of this image imposed a delimiting line on the DIT inland side of the image set.
Four parameters were used to evaluate the temporal variation of the tongue: the area, the major/minor axis length, segments AB and CD, respectively (see Fig. 3), and the mean advancing rate. The results are shown in Table I. In a period of about 12 years, the area increased by 6.3% and the major axis length by about 12%. Computed on the direction of the major axis, the length increase was 7.9 km. Thus, assuming no calving, the advancing rate was 699.  Some invariant structural features were identified as control points (CP) to estimate the average velocity of the ice stream. Analysing the displacement of 45 CP (see Fig. 3), the minimum and maximum distances were 7.1 and 8.17 km, respectively. Points with displacement greater than 8 km are indicated by white circles. Green circles indicate displacements of less than 8 km. As can be seen, the highest advances occur after the middle of the DIT extension, toward the tongue terminus. In the set of CP, the mean distance was 7.94 km and the standard deviation was 0.164 km. In the elapsed time, the computed average velocity was 703 m y −1 . Fig. 4 shows the time progression of the shape of the DIT. Some variations in the estimated parameters can be introduced by drifting ice and multiyear sea-ice becoming attached to the DIT [41]. In fact, large multiyear sea-ice sheets tend to accumulate for many years in the north margin of the DIT. These ice sheets are driven by the DIT flow into TNB. Then, the attached multiyear sea-ice can increase the transverse side of the DIT. However, in the analyzed images, a slight decrease is observed in the minor axis length, see Table I. The north margin is outlined by fractures, small rifts, and attached ice that alter the length of the minor axis. The DIT south side seems to be well preserved, but the opposite side shows a profile affected by the progression of irregular features, as indicated by the north outlines in Fig. 3 and by inspection of the contour progression in Fig. 4. A 100 km long segment was taken to estimate the area variation in both margins, see dotted squares in Fig. 3(b). In the elapsed time, subtle area variations of around -4.4% and 3.4% were computed on the north and south margins, respectively. Then, a slight variation in the length of the minor axis is formed.

IV. DISCUSSION
One relevant parameter of ice tongues is the ice velocity. Some reported studies are: by U.S. Geological Survey maps and Landsat images, between 1960 and 1993 the average velocity of the ice front was approximately 0.8 km y −1 [41]. Using Landsat images, 14 identifiable features on the glacier (arranged from the grounding line to the floating terminus) were the basis for measuring the DIT displacement. The average velocity for the period 1972-1973 was given as 0.84 km y −1 [42]. In a similar study, using 73 crevasse patterns referenced on a pair of Landsat images of 1973 and 1978, the average velocity was 0.7 km y −1 [3].
On a study of the effects of ocean currents on the dynamics of floating glacier tongues, observing the displacement of a prominent notch on the DIT north edge, the estimated velocity was about 730 m per year [43]. Precise measurements were performed by GPS. Then, three GPS stations were positioned on DIT. In the period 1992-1993, the GPS stations at the terminus and about 3 km of the DIT terminus provided an average velocity of 714 m y −1 . Close to the DIT grounding line, the average velocity was 553 m y −1 [44]. Incorporating interferometry data, a more automated approach tracked crevasse patterns on SAR images. In a three-year interval (1997)(1998)(1999)(2000) the reported advance rate of the front of the DIT was 733 m y −1 [45]. By our analysis, the advance rate, computed on the direction of the major axis, was estimated as 699.7 m y −1 .
In an analysis of the changing extent of the glaciers along the western Ross Sea, based on digital scans of antique maps and remote sensing imagery, the surface speed was estimated by averaging the displacement of some features on the DIT surface, i.e., large crevasses. In the period 2008-2012, the registered average speed was 711 m y −1 [46]. In the period 2016-2021, using DIT rift and ice front information, in the range 615-680 m per year, an average surface velocity of 670 m per year was reported in [47]. In our study, using the reference of 45 CP, the average DIT surface velocity was about 703 m y −1 , in the range 629-724 m y −1 . In the collected velocity measures, the obtained standard deviation was around 14.5 m y −1 . To reduce the uncertainty of the calculation, we sought to incorporate the largest possible number of CP.
Concerning the margins of error, the gray-level transitions can introduce deviations in the contour of the segmented shapes. For the length of the major and minor axes, the estimated error is about ±1 pixel at each end, so the largest error is about ±2 pixels (±100 m). A similar case occurs in the estimation of the iceberg area since the contour can deviate by ±1 pixel. Dilation and erosion operations were applied to simulate contour deviation and using this criterion the estimated area diverged by approximately ±1%. To estimate the closeness of the derived parameters of the DIT surface velocity to the true values, the number of CP was used to compute the margin of error (ME). By a statistical approach [48], the ME states intervals on the distribution of the likely error of the mean and the standard deviation. Assuming a Gaussian distribution in the sample data, the confidence interval for the mean (CIM) is computed as CIM = m ± zσ/ √ n, where m and σ are, respectively, the mean and the standard deviation of the sample data, z is the confidence level, and n is the number of CP. The term zσ/ √ n is the M. The confidence interval for the standard deviation (CISD) is CISD = σ ± zσ/ √ 2n. Using a confidence level of 95% (z=1.96), the confidence interval for the true mean value of the DIT surface velocity is CIM= 703 ±4.236 m y −1 . In other words, the ME (the error of the estimate) does not exceed 4.236 m y −1 , and the probability that the true mean value of the DIT surface velocity lies between 698.76 and 707.24 m y −1 is 0.95. The computed CISD was 14.5 ±2.996 m y −1 .

V. CONCLUSION
With regard to the processing scheme, the aim was to define a straightforward methodology for automatic detection of the shape of the DIT. This was accomplished by the combined application of a nonlinear SAR filtering technique and a fuzzy inference of a MAP segmentation. Both k-Means clustering and FCM perform automated progressions, but concerning the involved parameters, the number of clusters needs to be specified (two for the binary segmentation), as well as the fuzzifier q, which is selected by heuristic criterion. Operating iteratively, the FCM procedure drives a fuzzy inference of the Bayesian segmentation. The segmented map provides the basis for defining the shape of the DIT, but the main disadvantage is the dependence on scenarios where land/ice and the open sea are contrasted regions. The analysis of low-contrast scenes is still a challenging problem. The proposed method is well adapted to our image set, but wind-roughened sea surface and sea-ice in different phases could modify the performance of the proposed method. A coming study could include a segmentation contextual algorithm to detect the DIT in sea-ice landscapes and a gradient-based function to improve the recognition of CP.
In a time-window of 12 years, an increase in both the area and the major axis length was observed. By reported studies, we find that apparent stability can be supposed. The average displacement of the DIT is almost similar through the periods 1973-1973, 1992-1993, 2008-2012, and 2008-2020. Some of our estimations are in agreement with published data. After the big calving event of 2005 and 2006, the floating platform has resumed its natural advance into the Ross Sea, as demonstrated in this study.