Automated Detection of Marine Glacier Calving Fronts Using the 2-D Wavelet Transform Modulus Maxima Segmentation Method

Changes in the calving front position of marine-terminating glaciers strongly influence the mass balance of glaciers, ice caps, and ice sheets. At present, quantification of frontal position change primarily relies on time-consuming and subjective manual mapping techniques, limiting our ability to understand changes to glacier calving fronts. Here we describe a newly developed automated method of mapping glacier calving fronts in satellite imagery using observations from a representative sample of Greenland’s peripheral marine-terminating glaciers. Our method is adapted from the 2-D wavelet transform modulus maxima (WTMM) segmentation method, which has been used previously for image segmentation in biomedical and other applied science fields. The gradient-based method places edge detection lines along regions with the greatest intensity gradient in the image, such as the contrast between glacier ice and water or glacier ice and sea ice. The lines corresponding to the calving front are identified using thresholds for length, average gradient value, and orientation that minimize the misfit with respect to a manual validation data set. We demonstrate that the method is capable of mapping glacier calving fronts over a wide range of image conditions (light to intermediate cloud cover, dim or bright, mélange presence, etc.). With these time series, we are able to resolve subseasonal to multiyear temporal patterns as well as regional patterns in glacier frontal position change.

ice caps, and ice sheet contributes ∼43% to the contemporary global sea level rise [30]. The Greenland ice sheet, the largest contributor to sea level rise in the twenty-first century [18], was responsible for an average of ∼0.47 ± 0.23 mm · yr −1 of sea level rise over 1991-2015 [38].
Changes to the terminus position (i.e., calving front position) of marine-terminating glaciers influence glacier mass balance through the direct loss of mass at the calving front. In addition, loss of mass at the glacier terminus modulates the forces governing the glacier ice flow [10], [12]. Loss of the resistive stress generated at the terminus can result in acceleration and thinning of the ice, contributing further to "dynamic" glacier mass loss [3]- [6], [12], [25], [27]. Prior analyses of glacier terminus position change and its influence on the dynamics and mass loss relied on manual mapping of glacier calving fronts (see [5], [21], [27]), which is time intensive and dependent on human interpretation. The accuracy of these manual delineations may drift over time for each analyst and may introduce biases when multiple analysts contribute to data sets. Automated methods for delineating glacier calving fronts are more efficient, repeatable, and objective.
Previously published automated glacier terminus detection methods were based on image classification using various multispectral image bands (see [29], [31], [35]) and simple edge detection methods (see [34]). The application of multispectral band ratios for image classification, such as the normalized difference snow index (NDSI) or normalized difference water index (NDWI), can be applied with reasonable success to the satellite images acquired in months of the year when the snow cover does not extend beyond the glacier margins [31], [35]. However, the seasonal variability in snow cover directly affects the mapped glacier area and may result in large errors in delineation of the ice margin when applied at other times of the year, which limits the temporal resolution of these studies [31]. Furthermore, these image classification techniques require additional processing (e.g., thresholding and vectorization) to identify the continuous regions or boundaries within the image, which can introduce additional errors [32]. The use of gradient-based edge detection algorithms mitigates this processing step. Seale et al. [34] developed an automated method to identify the glacier terminus positions using a simple Canny edge detector, which identified largest gradients in image brightness along a glacier flow line. While this method This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ was successfully applied to the glaciers in the study, the method has not been adopted for other glacier terminus change studies and does not measure positions at other locations along the width of the glacier terminus (i.e., other than along the glacier centerline).
In order to efficiently and objectively quantify changes in the length of marine-terminating glaciers in a more robust manner, we adapt an automated image segmentation technique used previously in a wide variety of applied science fields called the 2-D wavelet transform modulus maxima (WTMM) segmentation method [17]. This edge detection method calculates 2-D gradients in the image brightness and is applied at a variety of spatial scales. The 2-D WTMM segmentation method places 1-pixel thick lines representing contours in intensity change for each spatial scale of analysis. We apply the method to analyze panchromatic Landsat satellite images, chosen for their high image acquisition frequency and radiometric resolutions, of marine-terminating glaciers. This gradient-based and multiscale automated delineation method allows for the detailed analysis of glacier terminus position changes at hundreds to thousands of time points. With the delineations from this method, we generate time series of terminus positions for several glaciers peripheral to the Greenland Ice Sheet. The glaciers along Greenland's periphery are one of the largest contributors to global glacier ice loss [36], yet their contributions to sea level rise through terminus position and associated dynamic changes have not yet been quantified.

II. METHODOLOGY
In order to efficiently analyze hundreds of images for each marine-terminating glacier, we developed an automated image processing technique to map the calving front positions. Images were analyzed over the Landsat 8 record (2013-2020) available from the US Geological Survey (USGS) on the Amazon Web Services (AWS) cloud. Ten glaciers were chosen to represent the distribution in size, morphology, fjord geometry, and sea ice conditions at the terminus for the 641 marineterminating glaciers that are peripheral to the Greenland ice sheet (Fig. 1). The peripheral marine-terminating glaciers were identified based on their Randolph Glacier Inventory (RGI) classifications. Although we focus on Greenland, these conditions are broadly representative of conditions across Arctic glaciers and glaciers along the Antarctic Peninsula.
Aside from the Landsat images, the only inputs required for our automated workflow are (1) a bounding box for the glacier terminus; (2) the Landsat scene boundaries available through the USGS; and (3) an ice velocity map. The glacier bounding boxes (i.e., terminus boxes) were manually drawn to span the narrowest portion of the glacier terminus, in line with the traditional box method for mapping termini [12], [20], [27]. The lengths of the boxes encompass the glaciers' recent terminus positions and extend several kilometers inland in order to account for potential retreat. The Landsat scene boundaries file was downloaded from the USGS (https://www.usgs.gov/media/files/Landsat-wrs-2scene-boundaries-kml-file). For ice velocities, we used the MEaSUREs 1995-2015 Greenland Ice Sheet Velocity Mosaic available at 250-m pixel resolution [14], [15]. The use of a time-averaged velocity data set was preferred over velocities derived from a single-year or less to minimize the effects of random errors and data gaps on our processing pipeline, particularly in regions of slow-flowing ice. Fig. 2 summarizes the steps used to process these downloaded satellite image subsets. The Landsat images available through AWS are Collection 1 Level-1 Precision and Terrain Correction (L1TP) and Systematic Terrain Correction (L1GT) data products from the Landsat 8 Operational Land Imager (OLI)/Thermal Infrared Sensor (TIRS).
To reduce the computational effort, the Landsat images are subset using buffered terminus boxes. The use of the buffer ensures that the automated delineations are not influenced by edge effects from the 2-D WTMM analysis. The buffer size corresponds to the maximum dimension (either length or width) of the terminus box. The subset 15-m resolution panchromatic Landsat images that fully overlap each terminus box, excluding images with >20% cloud cover in the terminus box, are automatically downloaded from AWS [ Fig. 2(a)]. The cloud confidence classification provided in the Landsat Quality Assessment (QA) band is automatically used to determine the fraction of pixels within the terminus box that correspond to high cloud confidence. If greater than 20% of the pixels are high cloud confidence pixels, the image is not downloaded.  [34]. This rotation establishes a common frame of reference for examining terminus changes in which increased distance from the left margin of the box corresponds to terminus advance and decreased distance corresponds with retreat.
The rotated image subsets are then analyzed using the 2-D WTMM segmentation method described in detail in Section II-A. Briefly, the WTMM method identifies points in the image that represent the maximum intensity gradients at the size scale of analysis. These points, named maxima points, are connected in lines called maxima chains that correspond to the regions with greatest intensity gradients throughout the image [ Section II-A describes this method in further detail, Section II-B describes the filtering of the maxima chains, and Section II-C describes the construction of terminus position time series.

A. Adaptation of 2-D WTMM Segmentation Method
The satellite images over the glaciers are analyzed using the 2-D WTMM segmentation method developed to perform automatic image segmentation for a variety of images across scientific fields including biomedicine, solar physics, astrobiology, etc. [2], [9], [16], [17], [22], [23], [33]. The 2-D WTMM segmentation method is a multiscale, gradientbased method that identifies contours representing the locally maximal changes in the intensity of an image. This is achieved by using the first-order (Gaussian) derivative of the 2-D smoothing function [17] φ( The continuous wavelet transform of the image f is calculated with respect to the partial derivatives of the smoothing function φ with respect to x 1 and x 2 1 and ψ 2 ( which amounts to taking the gradient of the convolution of the image with φ where * represents convolution, b represents the parameter of position, a represents a scale parameter, and T ψ 1 , T ψ 2 are the two components of the wavelet transform and The wavelet transform (T ψ ) is a gradient vector which has a magnitude (i.e., wavelet transform modulus, M ψ ) corresponding to the gradient in intensity and a direction (i.e., argument, A ψ ) that points to the highest intensity regions, which are expressed in polar coordinates as Maxima points represent the regions in the image where the intensity gradients (i.e., moduli) are maximal. WTMM are automatically connected along the maxima chains that act as the edge detection lines for the change in the intensity [ Fig. 3(d)]. The algorithmic procedure leading to the calculation of these one-pixel thick maxima chains is outlined in [23,Appendix]. These maxima chains are generated at  50 spatial scales (a = 0, 1, 2, 3, …, 49) that range from the minimum scale required to resolve the wavelet, 7 * 2 0/10 = 7 pixels, to 7 * 2 49/10 = 209 pixels, corresponding to a range of 105-3135 m in the 15-m resolution Landsat 8 panchromatic images. For the best results, glaciers wider than 105 m should be analyzed with this method.
In the satellite images of the glaciers, the 2-D WTMM calculates intensity gradients using the raw digital number in the pixels and generates maxima chains along the regions with high-intensity contrasts, such as the contrast between glacier ice and sea ice, open water, and land. The multiscale analysis allows for delineation of small-scale and large-scale features (Fig. 4) that allow the algorithm to adaptively delineate glacier termini of variable sizes, geometries, and environmental conditions without a priori knowledge of the image conditions [17].

B. Maxima Chain Filtering
Once all the maxima chains outside of the terminus box are removed, we objectively select the final maxima chain corresponding to the glacier terminus using the maxima chains' attributes. The properties for each chain include shape (closed or open), length (L), average modulus value of all maxima points (m), and arguments of the maxima points, which represent the chain's orientation. The closed loops are eliminated immediately, since the glacier termini are open linear features in the terminus box. The remaining attributes are extracted at each scale and objective thresholds are applied sequentially to filter the maxima chains. Generally, the maxima chain delineating the glacier terminus will have a high average modulus value (i.e., large-intensity contrast) and will be longer than most other chains, especially those corresponding to noise but also including features such as crevasses. Orientation is also used to filter the maxima chains because the glacier terminus should be primarily oriented vertically in the rotated reference frame, such that maxima points will mostly have arguments pointing left or right. The threshold for argument is applied to the fraction of arguments that make up each chain and the threshold for length is normalized by the maximum chain length determined for the image, so that the same thresholds can be applied to wide or narrow glaciers. Similarly, the threshold for average modulus value is normalized by the maximum average modulus value determined for the image, so that the same threshold can be applied to images with highor low-intensity contrasts.
In order to objectively determine the filtering sequence and thresholds used for filtering, we constructed a data set consisting of 512 manual terminus delineations for 5 of the 10 sample glaciers chosen for their diverse morphologies (see Fig. 1 for additional information). For threshold optimization, we randomly allocate 90% of the manual delineations to a training data set and the remaining 10% to a test data set for cross validation. We perform the cross validation four times and report the average errors from the four iterations. Using the training data set, we constrain the length, modulus, and argument thresholds (C L , C m , and C A , respectively) using an optimization technique that minimizes a cost function () that represents the misfit between the automated delineations and the manual delineations. The misfit for each image (i.e., timepoint) of analysis is the average difference between the automatically and manually delineated terminus positions along three glacier flow lines at one-quarter, onehalf, and three-quarters of the width of the glacier terminus box [ Fig. 5(a)] expressed as the variable X diff where |X auto − X manual | represents the Euclidean distance between the automatically and manually delineated terminus positions for a given flow line. These flow lines are extrapolated from the quarter points of the left and right sides of the glacier terminus boxes, which are determined from the box vertices using several midpoint calculations [ Fig. 5(a)]. The use of three flow lines at varying distances across the glacier's width allows the method to extract terminus positions from delineations that do not span the entire width of the terminus. In addition, patterns in terminus position change are likely to vary closer to the fjord walls versus the middle of the glacier width as the terminus position tends to be more stable along the margins [12]. This strategy could be expanded by calculating the terminus positions along the additional flow lines. Calculating terminus position along three flow lines allows our method to resolve the variability in terminus position change along distinct longitudinal profiles, which is useful for analysis of differences in terminus dynamics due to glacier geometry, velocities, and other factors that vary along the glacier width. In contrast, the box method would evaluate terminus change in terms of a single area-averaged change value [20], [24], [27]. The misfit for each image X diff is averaged for a total of N = 460 images in the training data set and divided by F 3 , where F is defined as the ratio of the number of intersections with the glacier flow lines (i.e., terminus positions calculated) to the total possible number of intersections (3N), in the final cost function The filtering sequence using the thresholds C m , C L , and C A affects the automated delineations of terminus position and thus the misfit. To identify the optimal sequence, we calculated the cost using a range of 0-1 for each of these normalized thresholds, for each combination of filtering order. The sequence that yielded the lowest median cost value was chosen as the optimal order for filtering: average modulus value, length, and then argument. Using this optimal filtering sequence, we computed the costs for combinations of the thresholds C m , C L , and C A . For C m , the cost was minimal at C m = 0.7 for all combinations of C L and C A . Therefore, we set C m to 0.7 and computed costs for a grid of C A and C L values. With C m set to 0.7, cost value was minimal for C A = 0.1 and C L = 0.4 (Fig. 6).
The WTMM chains that do not satisfy the thresholds set above are eliminated. Chains that do not correspond to the glacier terminus (e.g., delineations of long, high-contrast, and vertical features such as shadow boundaries and sea ice margins) may still remain but may be removed when the time series is filtered (see Section II-C). These remaining chains across the 50 scales of analysis are aggregated and the five chains that are most likely to delineate the terminus We choose up to five chains from each image with the highest metric values. If fewer than five chains remain after the thresholding, then all the remaining chains will be identified as the top chains, which will be allowed to pass on to the time series filtering.

C. Constructing Terminus Position Time Series
Using the procedures described above, there are up to five maxima chains that delineate the terminus position in each Contour plot of the cost as a function of length and argument thresholds, C L and C A , respectively. C m is set as 0.7. The white contour is the contour boundary separating the lowest cost interval, located at C L = 0.4 and C A = 0.1. available satellite image. For each maxima chain, we extract the points of intersection with the three glacier flow lines. These points are then filtered iteratively using terminus change rates calculated as a forward difference using the terminus position of the current time point and the subsequent time point. If the rate of terminus advance calculated from the delineation is three times greater than the maximum glacier speed, this delineation is considered inaccurate and the point is eliminated. Although the glacier terminus cannot advance faster than the rate at which the ice is flowing, we use a conservative threshold of three times the maximum flow speed within the glacier terminus box in order to account for temporal variations in glacier velocities. The same threshold is applied to points that yield a high rate of terminus retreat if followed by the rate of advance that violates the flow speed condition. The filtering is repeated three times; and at each time, the terminus change rates between time points are recalculated. If multiple terminus position points remain for one time point, the point associated with the highest metric value defined in (8) is chosen to represent the terminus position at that time point. Terminus positions along each flow line are filtered separately.
The resulting filtered time series are dense, showing changes in terminus position at subseasonal timescales (Fig. 7). The time series shown in Fig. 7 allow for observation of seasonal glacier retreat from ∼April through ∼October, which is the typical timing of retreat for glaciers in Greenland [11], [26] as well as the overall trend in retreat from 2013 to 2020.

III. EVALUATION OF THE METHOD'S PERFORMANCE
We evaluate our automated method's accuracy in comparison to the data set of manual delineations, with consideration of the various environmental and image conditions shown in Fig. 8 (clear, bright and dim lighting, thin cloud presence, sea ice presence, shadow presence, etc.). We compare the automated-manual differences to manual-manual differences in delineation. The manual-manual differences in terminus positions were determined experimentally from two human analysts delineating the same n = 50 images over the five sample glaciers (see Section II-B), including images with each of the conditions listed above. Experiments yielded a standard deviation in the manual delineation of ±31.0 m or ∼2 pixels. The accuracy of the automated method on an imageby-image basis, referred to from now onward as the point uncertainty, was calculated using the misfit (X diff ). Table I shows the median misfits between the automated and manual delineations by image condition for 10% (51 out of the 512-see Section II-B) of manual delineations excluded from the threshold optimization as well as the full manual data set. The cross-validation values in Table I are averaged from four iterations, with different training and testing data sets each time.
Based on these results, the automated method is capable of delineating terminus positions at an accuracy similar to that of manual uncertainties in images that are clear, uniform in brightness (either dim or bright), or contain thin cloud cover (Table I). For clear images, where there is high contrast between the glacier ice and open water [e.g., Fig. 8(a)], the automated method effectively delineates terminus positions within 1-pixel uncertainty from the manual delineations. Due to the method's gradient-based algorithm, it also effectively delineates glacier terminus positions in images that are uniformly bright or dim [ Fig. 8(b) and (d)] and where thin clouds are present [ Fig. 8(e)], with better than 1-pixel uncertainty. The median misfit for these three image conditions was 12.5 ± 7.5 m or <1 pixel. The automated method's point uncertainty for these image conditions is within the interanalyst uncertainty in delineation of ∼2 pixels.
Images where the fjord walls cast shadows across the glacier terminus, producing a higher intensity gradient than the gradient along the terminus, delineations are often along the shadow boundary [ Fig. 8(c)]. In addition, the method may delineate the shadow cast from the ice cliff at the calving front onto the sea ice. For the latter condition, the distance between the glacier terminus and the delineated shadow from the ice cliff will be greater when the incident angle of the sunlight is greater and will depend on the glacier's aspect. North-facing glaciers are particularly susceptible to these inaccuracies. If the ice cliff shadow is delineated using the automated method, then the movement of the shadow with the gradual change in sun angle over time may result in a false decrease in terminus position prior to the start of the actual glacier retreat in April, as seen in the time series shown in Fig. 7(a). It may be possible to correct for the shadow offset using the time of acquisition of the image, the viewing angle of the satellite, and precise knowledge of the surrounding terrain; but we do not explore these more complex methods here. Terminus positions delineated during low sun angle months at high latitudes should be interpreted with caution, but the automated delineations are robust for delineations in less challenging lighting conditions. Sea ice breakup can complicate the delineation of the terminus boundary. The automated method can delineate the terminus boundary when there is sea ice that is cohesive and extends out of the terminus box [ Fig. 8(g)]. However, sea ice breakup within the terminus box can present a challenge to delineation, since the high contrast between the bright sea ice and open water is often higher than the contrast between glacier ice and sea ice. For glaciers where gradual seasonal sea ice or mélange breakup is prominent, the automated method may delineate boundaries that are further out than the glacier terminus. Due to our selection of top five maxima chains and our time series filtering using glacier velocities, these false advances are often filtered out. The standard glacier velocity threshold of three times the maximum flow speed could be made stricter (e.g., 1 or 2 times the maximum flow speed) for more stringent time series filtering. However, lowering the glacier velocity threshold results in a tradeoff between the improvement in time series for these glaciers and a loss of temporal resolution for glaciers with large uncertainties in velocity, which are overfiltered as a result. In the absence of accurate, high spatial resolution time series of ice velocity, we recommend using the large standard threshold of three times the maximum flow speed and performing additional time series filtering outside the algorithm.
Snow cover is another seasonal condition that may influence the method performance. Snow covering the boundary between glacier ice and sea ice may greatly reduce the intensity contrast across the terminus boundary. In addition, spatial variability in snow covering the glacier ice Fig. 5(a) may create intensity gradients within the glacier margins that are greater than the gradients across the terminus boundary. Delineation of nonglacial snow margins will likely be similar to those corresponding to noisy features (shorter in length and variable in orientation) which would be filtered out during the thresholding process. In cases where the method does not accurately resolve the terminus boundary, the calculated terminus positions may also be filtered out through the time series filtering steps discussed in Section II-C.
As the glaciers retreat and advance, the movement of the glacier terminus outside the terminus box and to new positions along the fjords with variable orientations is possible. The glacier terminus boxes were drawn to encompass the glaciers' terminus positions in 2000 and 2015, accounting for the retreat along the fjords. Therefore, few adjustments to the box length or orientation were necessary for our analysis over the Landsat 8 record from 2013 to 2020. However, outside the 2000-2015 time period, the terminus may retreat or advance to a position in the fjord outside of the box and potentially exhibit a different orientation. In these cases, the terminus boxes will require adjustment, which can be implemented automatically using ice velocity vectors to determine glacier orientations along the fjord as well as a window in which the terminus position movement is possible.
Despite the number of temporal variations and environmental conditions that may challenge the method's ability to delineate glacier front positions, our method is able to produce dense time series of frontal positions. These time series resolve glacier length changes at subseasonal timescales and at numerous locations across the glacier width, which will allow for the detailed analysis of the length and the terminus geometry change in conjunction with the time series of glacier dynamics

IV. CONCLUSION
The use of the adapted 2-D WTMM segmentation method for delineation of glacier calving fronts is promising. This automated method is capable of delineating marine glacier calving fronts within 1-pixel uncertainty in images with clear conditions, dim or bright lighting conditions, and thin cloud presence. Images where there are shadows cast by the fjord walls across the glacier terminus and where there are shadows from the ice cliff at the calving front remain challenging to delineate using this method. However, knowing that two humans do not produce identical delineation lines, our experiments show that the automated method has an uncertainty that is within the interhuman variability. This automated method can be applied to accurately and efficiently resolve subseasonal to multiyear patterns in the glacier terminus position change in a variety of image conditions.
The resulting time series of glacier terminus position generated from this method can be used to assess glacier calving front shape changes as well as spatiotemporal patterns in glacier length change. We have shown that the time series are able to capture a range of glacier behavior including gradual advances and retreats, calving, and changes in terminus shape. Immediate future work will focus on extending the workflow to include images from earlier Landsat missions such as Landsat 7 and Landsat 5. Future work should explore its application to other image types (satellite radar images and other optical images such as Sentinel-2, etc.), in other regions, and for delineations of other features such as coastlines or closed features such as icebergs and lakes. Adaptation of the method to analyze satellite radar images could be particularly useful for increasing the temporal resolution of these terminus position time series, as radar images are not hindered by cloud presence or darkness in polar night. The code used for all image analyses discussed is available in a GitHub repository (https://github.com/julialiu18/automated-glacier-terminus).