Co-Registration Methods and Error Analysis for Four Decades (1979–2018) of Glacier Elevation Changes in the Southern Patagonian Icefield

Abstract

The main goal of this paper is to compare two co-registration methods for geodetic mass balance (GMB) calculation in 28 glaciers making up the Upper Santa Cruz River basin, Southern Patagonian Icefield (SPI), from 1979 to 2018. For this purpose, geospatial data have been used as primary sources: Hexagon KH-9, ASTER, and LANDSAT optical images; SRTM digital radar elevation model; and ICESat elevation profiles. After the analyses, the two co-registration methods, namely M1, based on horizontal displacements and 3D shift vectors, and M2, based on three-dimensional transformations, turned out to be similar. The errors in the GMB were analyzed through a k index that considers, among other variables, the error in elevation change by testing four interpolation methods for filling gaps. We found that, in 63% of the cases, the relative error in elevation change contributes 90% or more to k index. The GMB throughout our study area reported that a loss value of −1.44 ± 0.15 m w. e. a⁻¹ (−3.0 Gt a⁻¹) and an ice thinning median of −1.38 ± 0.11 m a⁻¹ occurred within the study period. The glaciers that showed the most negative GMB values were Upsala, with an annual elevation change median of −2.07 ± 0.18 m w. e. a⁻¹, and Ameghino, with −2.31 ± 0.22 m w. e. a⁻¹.

1. Introduction

In general, studies of mass balance and ice sheet volume changes are relevant in the context of global warming and the consequent rise in mean sea level [1]. A sustained trend of rapid glacier retreat has been observed in most mountainous and cold regions of the world [2], including the Patagonian Andes [3,4], among others. The Southern Patagonian Icefield (SPI) region has increased its rate of frontal thinning and retraction in recent years and is in a phase of unstable retreat.

Therefore, given the sensitivity to climate change in the glacierized areas, it is necessary to quantify the changes that occur over time. In that regard, remote sensing provides effective tools, and thus, has been widely used over the last decades to extend and improve analysis at a regional and at a global scale. This is due to the fact that systems based on satellite missions continue improving imaging performance, including better spatial and spectral resolution [5]. Therefore, satellite platforms provide a wide range of geospatial data at different scales of analysis for environmental monitoring. Consequently, products and data such as digital elevation models (DEMs) of glaciated terrain are commonly used to measure changes in glacier geometry and volume [6]. Satellite-based data have been an excellent choice for analyzing glaciers in inaccessible areas, whereas in situ measurements present challenges in terms of logistics and the ensuing high economic cost and can be applied only in reduced zones. It should be noted that the use of unnamed aerial vehicles (UAV) in glacier monitoring has advantages due to their good geometric resolution but is subject to limitations in remote areas at a regional scale.

Over the last decades, the geodetic mass balance (GMB) has been widely used to estimate volumetric changes in glaciers and can be generated for entire mountainous ranges [7]. The geodetic method which consists of multitemporal satellite DEMs differentiation has proved to be an effective technique to be implemented in glaciated zones all over the world. The applicability of geodetic mass balance records depends on the accuracy of such records and DEMs resolution [8]. The impact on the accuracy of the generated products is also related, for example, to the quality of the DEMs used and to spatial resolution differences [9]. The distortions introduced during photogrammetric DEM generation affect the final quality of any glaciological results obtained. Thus, it is important to quantify the distortion and noise introduced in the models, as well as to determine the degree of uncertainty in the estimation of ice volume changes [10]. To maximize the value of geodetic mass balance records, realistic uncertainty assessments are necessary [8].

Subtracting DEMs requires analyzing multiple factors [9,11,12] to be able to detect the uncertainties or restrictions involved in each case. Differencing two DEMs can lead to relative planimetric and rotational biases [13]. These biases can spread and affect later analyses such as glacier volume changes and retreat rates [14], and for that reason these systematic errors should be calculated and corrected before DEM differencing. The process that minimizes this error is called co-registration or surface matching and consists of adjusting a digital model to another model [15]. The bias induced by not co-registering may lead to erroneous estimates of glacier volume changes or to false behavior detections [16]. Accordingly, co-registering two surfaces is key to estimate the elevation changes arising from DEM differencing [16,17].

The various existing co-registration methods available can be firstly grouped, according to [18], into those leading to solutions based on planimetric displacements and, secondly, into those whose solutions are obtained through 3D transformations. In addition, these methods can be grouped according to their reliance on the selection of common points (CP) or not [19]. It should be noted that a CP is a three-dimensional point showing a univocal correspondence between DEMs.

We establish two categories for the estimation of transformation parameters, the first one consisting of methods that provide solutions based on horizontal displacements [20] and 3D shift vectors [16], which are intended to minimize altimetry biases between both DEMs. In the field of glaciology, Kääb et al. [21], Berthier et al. [20], and Nuth and Kääb [16], among others, have applied this type of co-registration. Mass balance studies in the Central Andes and Southern Patagonia have used this method (e.g., [22,23,24,25]). We highlight that, within this group, there are studies that use more complex algorithms to calculate the transformation parameters (e.g., [26,27]).

In the second category, there are methods that lead to more complex numerical solutions based on three-dimensional transformations (e.g., [13,14,28]). There are very few applications of the second type in the field of glaciology (e.g., [13]). In the Central Andes, Lenzano [10] and Lenzano et al. [29] achieved co-registration by means of three-dimensional transformations. It is worth mentioning that, in Southern Patagonia, there are no applications of this method.

Nevertheless, it is important to estimate and evaluate co-registration accuracy. As a result of differenced DEM misalignment, a bias occurs in derived elevation changes which is directly related to the combined distribution of the slope and aspect of a glacier [30]. An indicator of DEM differencing uncertainty is the existence of terrain elevation differences where there should be no changes (stable terrain). In this respect, thanks to their accuracy and global coverage ICESat/GLAS elevations are a widely used source of data (e.g., [31,32,33]). Nuth and Kääb [16] and Paul et al. [9] have provided some examples of ICESat/GLAS data utilization for the quantitative evaluation of co-registration uncertainty. This makes it possible to estimate the influence of co-registration uncertainty on derived product error estimation (e.g., [34,35,36]), as is the case with GMB. Finally, we highlight the importance of performing glaciological studies dealing with GMB error impact propagation for each component involved in the mentioned balance calculation.

The main goal of this study is to apply two co-registration methods, called M1 and M2, and to assess their influence on GMB errors in the SPI during the 1979–2018 period. Note that M1 was proposed by Nuth and Kääb [16] and has been commonly used in glaciological applications. On the contrary, M2 is based on Li et al. [18] and has not been utilized in the field of glaciology. Images such as Hexagon (1979), ASTER (2018), Landsat 7 (2000), and SRTM elevation model (2000) were used. For this purpose, the following actions were undertaken: (1) Generation of DEMs by means of the photogrammetric method for years 1979 and 2018 by using Hexagon and ASTER images; (2) application and evaluation of two co-registration methods (M1 and M2) for the study periods by means of three parameters (three translations) and seven parameters (three rotations, three translations, and one scale factor), respectively; (3) calculation of area evolution and GMB using different interpolation methods for the study period for gap filling; (4) error impact evaluation of variables involved in balance calculation; (5) presentation of glaciological results. Thus, a contribution is made to optimizing the choice of co-registration methods between two DEMs and the impact of errors on glaciological results.

Study Area

The SPI represents, along with the Northern Patagonian Icefield (NPI), the largest glaciated region of the Southern hemisphere [37]. The glacial cover of the whole Argentinean territory had an area of 5741 km² in 2016 and almost 60% corresponded to the Argentinean portion of the SPI [38]. The SPI is located between latitudes 48°20′ and 51°30′ S and extends 350 km along the Andes, showing a width ranging between 30 and 40 km [39]. It occupied an area of around 13,000 km² in 1986 [40] and of 12,573 km² for the 2000–2015/16 period [24]. The ice surface of the SPI has a mean elevation of 1500–2000 m a. s. l. [41] and an average height of the outlet glaciers of 1355 m a. s. l. [39]. The slopes in the valleys are covered by native Nothofagus forests surrounding the glacier fronts [42].

The Santa Cruz River basin is the second largest hydrographic system fully situated within the national territory whose source is mainly glacial. It is formed by Viedma Lake subwatershed (1087 km²), Argentino Lake Northern Branch subwatershed (1211 km²), and Argentino Lake Southern Branch subwatershed (606 km²). Within it are situated the most important calving glaciers on the east SPI side in terms of size and dynamics which feed freshwater into Argentino and Viedma lakes.

According to the limited coverage of the Hexagon KH- 9 images, we decided to study 28 glaciers and part of the plateau on the Eastern watershed covering an area of 2087 km² which represent a major portion of the Upper Santa Cruz River basin (Figure 1). These glaciers are included in Randolph (RGI) and GLIMS (Global Land Ice Measurements from Space) international glacier inventories.

2. Input Data

2.1. Imagery, DEMs, and Elevation Profiles Used

Table 1. Optical images employed, their characteristics and uses.

Sensor	Product ID		Acquisition Date mm/dd/yyyy	Bands	Uses
Hexagon KH-9	DZB1215-	500030L001001	03/20/1979	Pan	Photogrammetric DEM
		500030L002001
ASTER	AST_L1A_	00310222018143801_20200309150631_11245		VNIR_3N VNIR_3B	Photogrammetric DEM
		00310222018143810_20200309150551_7618	10/22/2018
		00310222018143819_20200309150530_3445
	AST14DMO_	00310222018143801_20200308212920_11514		V3N V3B	Glacier outlines
		00310222018143810_20200308212930_11579	10/22/2018
		00310222018143819_20200308212940_11770
Landsat 7 *	LE07_L1TP_	231094_20000113_20170215_01_T1	01/13/2000	B1 B2 B3	Viedma glacier and its southern limits
		231095_20000113_20170215_01_T1	01/13/2000
		231094_20000402_20170212_01_T1	04/02/2000		Viedma glacier, Upsala glacier, and its Eastern tributaries
		231095_20000402_20170212_01_T1	04/02/2000		Upsala Eastern limits
		231094_19991126_20170216_01_T1	11/26/1999
		231095_19991126_20170216_01_T1	11/26/1999		Upsala Western limits, glaciers at the Western zone of Upsala glacier
		231095_20010320_20170206_01_T1	03/20/2001		Plateau zone between Upsala and Perito Moreno glaciers, glaciers at the western zone of Upsala and Perito Moreno
SRTM2000 (Global Version 3 1 arc-sec)	s50	_w073_1arc_v3	02/11/2000	C-band	SRTM2000 co-registration to KH91979 for error assessment, contour map digitalization for void filling
		_w074_1arc_v3
	s51	_w073_1arc_v3
		_w074_1arc_v3

* The Landsat 7 images are not listed in chronological order. Instead, they are in order of highest to lowest priority as used to digitize glacier outlines according to the date of the SRTM model.

summarizes the main optical products and accessories described below and used in this study.

2.1.1. Hexagon KH-9 Model

Between the 1960s and the 1980s, a series of intelligence missions were conducted by the U.S. government to map land areas. These missions started under the name of Corona and ended up being called Gambit and Hexagon [43]. In all, more than 900,000 images acquired through these programs and originally secret were gradually declassified [44] in 1995, 2002, and 2012 [45,46,47]. In this study, two images that were 14 µm digitized photograms from the last Hexagon mission (

Table 2. Hexagon KH-9 sensor characteristics.

Mission Parameters	Value
Mission number	1215-5
Orbit (km)	177.0 × 256.0
Film (In)	9 × 18
Camera resolution (ft)	30–35
Scanning resolution	14 µm
Film type Focal length Lens type	SO-315
	304.8 mm
	f/6 Biogon

) were used to generate KH9₁₉₇₉ DEM. These images are available through the United States Geological Survey (USGS, https://earthexplorer.usgs.gov/ (accessed on 8 May 2017)).

2.1.2. ASTER Model

ASTER (Advanced Spaceborne Thermal Emission and Reflection Radiometer) system is an image acquisition instrument mounted on the U.S. Terra satellite (https://asterweb.jpl.nasa.gov/ (accessed on 21 July 2020)), launched in December 1999. The ASTER sensor captures images on 14 bands of the electromagnetic spectrum obtained through three different subsystems: VNIR (visible and near infrared), SWIR (shortwave infrared) and TIR (thermal infrared). The stereoscopic images (Band 3 Nadir and Band 3 Backward) contained in “AST_L1A” product (https://search.earthdata.nasa.gov/search (accessed on 20 March 2020)) (see Table 1) enabled the generation of AST₂₀₁₈ DEM for GMB calculation. Orthorectified images were also generated for glacier outline delineation in 2018.

2.1.3. SRTM Model

SRTM (Shuttle Radar Topography Mission) was a joint project between USA, German, and Italian agencies whose vehicle was launched on 11 February 2000 [48]. Different versions of SRTM global coverage elevation model are currently available through USGS (https://earthexplorer.usgs.gov (accessed on 19 March 2019)). This study used a mosaic based on four SRTM Global Version 3 files with 1 arc-sec resolution (Digital Object Identifier (DOI) number: /10.5066/F7PR7TFT). This product, thereafter, called SRTM₂₀₀₀, was used for co-registration between DEMs and for their later evaluation.

2.1.4. Landsat 7 Images

The images were used for glacier outline delineation (Table 1) and were selected according to the proximity of the SRTM₂₀₀₀ model acquisition date.

2.1.5. ICESat Profiles

The Ice Cloud and Land Elevation Satellite (ICESat) was designed to accurately measure thickness changes in large polar ice layers in the Antarctic and Greenland [49]. Fifteen products (GLAH01 to GLAH15) derive from Geoscience Laser Altimeter System (GLAS) sensor, out of which GLAH12 to GLAH14 correspond to elevation products of different areas [50]. The GLAS sensor products acquired were downloaded from NASA official website (https://search.earthdata.nasa.gov/ (accessed on 12 July 2020)) (

Table 3. ICESat profiles used in this study.

File Name	Acquisition Date	Laser Campaign	Footprint Size
			Major Axis Mean ± St. Dev.	Eccentricity Mean ± St. Dev.
GLAH14_634_2107_002_0043_0_01_0001	2004-February-25	L2B	89.52 ± 4.93	0.822 ± 0.045
2107_003_0043_0_01_0001	2004-May-26	L2C	88.37 ± 19.12	0.892 ± 0.044
2109_002_0043_0_01_0001	2004-October-12	L3A	55.79 ± 0.43	0.567 ± 0.043
2111_002_0043_0_01_0001	2005-February-27	L3B	79.53 ± 11.55	0.753 ± 0.051
2111_003_0043_0_01_0001	2005-May-29	L3C	55.41 ± 1.84	0.633 ± 0.034
2113_002_0043_0_01_0001	2005-October-30	L3D	52.04 ± 1.06	0.523 ± 0.010
2115_002_0043_0_01_0001	2006-March-02	L3E	52.31 ± 1.60	0.483 ± 0.040
2115_003_0043_0_01_0001	2006-June-01	L3F	51.20 ± 1.63	0.480 ± 0.023
2117_002_0043_0_01_0001	2006-November-02	L3G	53.41 ± 1.51	0.510 ± 0.037
2119_002_0043_0_01_0001	2007-March-20	L3H	55.61 ± 0.48	0.521 ± 0.019
2121_002_0043_0_01_0001	2007-October-11	L3I	57.28 ± 0.57	0.590 ± 0.013
2123_002_0043_0_01_0001	2008-February-25	L3J	58.66 ± 1.52	0.575 ± 0.036
2125_002_0043_0_01_0001	2008-October-12	L3K	51.99 ± 1.12	0.611 ± 0.036
2129_002_0043_0_01_0001	2009-March-17	L2E	NA	NA
2131_002_0043_0_01_0001	2009-October-09	L2F	NA	NA

NA, not available. This table was constructed using the metadata included in the HDF5 files and the official metadata table provided by the ICESat Science Investigator-led Processing System (I-SIPS), 2014.

). This study used the latest version available at the time of processing of GLAH14 product from the GLAS/ICESat L2 Global Land Surface Altimetry Data, version 34, HDF5 format [51]. This was used over stable (off-glacier) terrain to evaluate co-registration accuracy.

2.2. Accessory Data Sources

2.2.1. Glacier Outlines

Open-source GLIMS vectorial layers were used for glacier outline delineation (http://www.glims.org/maps/glims (accessed on 12 July 2020)). Only the outlines dated in 1986, 2000, and 2016 were selected. Glacier limits were manually modified following Hexagon KH9, ASTER, and Landsat data sources (Table 1). The outlines corresponding to 2016 were updated to 2018 by using ASTER orthorectified images. In the case of 2000 outlines, they were modified employing a series of Landsat 7 images. Next, the outlines of 1986 were modified to 1979 through KH9 orthophoto. Finally, these outlines were applied to the glaciers area evolution between 1979 and 2018.

2.2.2. Ground Control Points

The fifteen ground control points (GCPs, Figure 1) used to generate the KH9₁₉₇₉ and AST₂₀₁₈ DEMs were measured in situ by Differential Global Positioning System (DGPS). For this purpose, we performed campaigns in the study area where a GPS network with strategically positioned points was installed. A Trimble double-frequency GPS receiver was used. The data were processed with Bernese 5.0 software [52]. GNSS carrier phase data were processed with precise ephemeris from the International GPS Service (IGS) and yielded millimetric accuracy. The GPS network was linked to the global system WGS84 through the IGS continuously operating reference stations (CORS).

3. Methodology

The methodology presented in this study follows the scheme in Figure 2: (1) The input data from images and elevation products were prepared; (2) DEM extraction based on Hexagon and ASTER images was carried out through photogrammetric processing; (3) elevation products were co-registered by using methods M1 and M2; (4) co-registration procedures were analyzed and assessed; (5) GMB calculation from AST₂₀₁₈ and KH9₁₉₇₉ models was co-registered by M1; (6) differencing voids were filled through four combined interpolation methods called Global 1 (G1), Global 2 (G2), Local 1 (L1), and Local 2 (L2); (7) GMB errors were estimated and the impact was analyzed for each component involved.

3.1. Dataset Preparation

3.1.1. SRTM₂₀₀₀ Model

Although the product used was provided with void filling by the National Geospatial-Intelligence Agency (NGA), the study area still showed some gaps which were interpolated by the inverse distance weighting (IDW) algorithm. A mask was also created in these areas to be able to identify the location of the interpolated data in future analyses. Product elevations were originally referenced to EGM96 geoid [48,53]. In order to keep consistency between the datasets, elevation values were converted to WGS84 ellipsoid with EGM96 geopotential model [54,55]. Lastly, the final product was projected according to UTM system, zone 18S, with a resolution of 30 m.

3.1.2. ICESat Elevation Data

The fifteen ICESat/GLAS products contain, in all, more than 2000 parameters [50], among which there are several “flags” that can be used to filter and correct elevations [56]. Some inconsistencies were found as regards the names of flag parameters and their uses. A thorough bibliographic search was performed which made it possible to choose which parameters to apply.

ICESat Filtering

The quality parameters used were: (a) elev_use_flg = 0 index, which shows that the elevation calculated in the product is valid [57,58]; (b) sigma_att_flg = 0, which shows that the altitude calculated in the precision attitude determination (PAD) is good quality [58]; and (c) sat_corr_flg ≤ 1 saturation index, which shows whether the signal received by the sensor is not saturated (=0) or the saturation correction is negligible (=1) [57,59].

As regards waveform indicators, we used the i_numPk parameter which corresponds to the number of peaks in the wave received by the sensor. Thus, we chose to accept the elevation value when i_numPk ≤ 5, which made it possible to eliminate most of the points corresponding to forest reflections [60]. For cloud filtering, ICESat elevation values differing by more than 50 m from SRTM₂₀₀₀ model were excluded [61]. For comparison purposes, the d_DEM_elv geophysical parameter available was used, which corresponds to SRTM model elevation for each footprint [58] and which is also referenced to TOPEX/Poseidon (T/P) system as well as ICESat/GLAS elevations. In addition, as a general rule, we took into consideration that, if the quality of ICESat/GLAS data is reasonable, quality indicators are set at zero [56].

ICESat Elevations Correction

In our case, saturation correction (d_satElevCorr) was not applied. According to Sun et al. [62], it is not advisable to use it when terrain slope is higher than 7°. ICESat/GLAS elevations had to be converted to WGS84 system as they were originally referenced to T/P system. For this purpose, d_deltaEllip geophysical parameter was used, defined as the difference between elevations above T/P ellipsoid and WGS84 ellipsoid, respectively [63]. Thus defined, this parameter was subtracted from each elevation measurement. To obtain the final profiles, the elevations were filtered and corrected by using sample codes provided by HDF-EOS Tools and Information Center [64].

3.2. DEM Processing and Extraction

3.2.1. KH9₁₉₇₉ Image Correction and Model Generation

Each KH-9 mission photogram is scanned by USGS and provided in two files with an overlap of around 1 cm [65,66]. Hence, a mosaic is generated by searching for common features within the overlap area. The classified nature of Corona/Hexagon missions does not allow the adoption of a traditional photogrammetric approach due to the difficulty of finding information related to camera type, flight altitude, and orbital data [67,68,69]. Although the images contain fiducial marks, the lack of camera calibration details hinders the application of a classical internal orientation step to solve scene distortions. Nevertheless, the photograms have a reseau grid with a spacing of 1 cm [70,71], which makes it possible to restore the geometry at the moment of photograph acquisition [66]. In this study, we made use of this approach to correct the geometry of the images.

For image geometric correction, the steps were as follows: (1) Automatic detection of crosses; (2) calculation of original positions; (3) calculation of distortion vectors; (4) removal of crosses with a Laplacian filter (this step is required not to hinder later automatic correlation processes for common feature detection [66,72]); (5) affine transformation of the image utilizing the vectors obtained in step (3). We implemented this transformation due to its geometric advantages for shape preservation [73] and the distortion pattern observed in Figure 3. The first results exhibited a rotational pattern, as shown in [65]. This rotational component was removed because it masked the pattern shown in Figure 3. Both images, the left-side and the right-side ones, showed similar patterns. The resulting distortion vectors were on average equal to 8.0 ± 3.1 and 7.7 ± 3.1 pixels, respectively. Finally, these quantities were corrected by affine transformation.

After that, we proceeded with a traditional photogrammetric processing by using Photomod 4.2 software [10,74] in order to generate the KH9₁₉₇₉ model. In the first place, for internal orientation, we worked with a focal distance of 304.8 mm, considering the principal point coordinates that matched the central cross of the reseau with a pixel size equal to 13.9 µm. Then, for the Aero Triangulation (AT) step, twelve GCPs measured on the terrain were arranged (Figure 1). Following that, for relative orientation 72 tie points were added to cover all Von Gruber areas. Finally, the block was adjusted by the independent stereopair method, obtaining subpixel results. The mean absolute residuals of GCP were X = 3.22 m, Y = 2.74 m, and Z = 2.56 m.

The automatic correlation for DEM generation shows some complications due to low-contrast white zones (e.g., [65,75]). In addition, most difficulties appeared in high glacier zones where a poor terrain morphometric definition was observed. Consequently, local artifacts had to be implemented and data gaps arose in accumulation regions. For automatic point cloud extraction, we used a high correlation threshold (0.9), and the final point cloud was filtered due to the gross errors involved. For this process, the decision was made to compare each point against its corresponding SRTM₂₀₀₀ elevation, excluding those points with an absolute difference larger than the mean plus three standard deviations. Then, the bilinear interpolator was used to create the model to fill the gaps in the high zones. It should be noted that the bilinear method was chosen due to the simplicity of its application [76,77] and the scarce spatial distribution of the points in these zones to digitize the Triangulated Irregular Network (TIN). Nevertheless, the terrain detail was reinforced by manually incorporating topographic breaklines. Finally, KH9₁₉₇₉ was generated from the point cloud filtered and interpolated from TIN with a pixel size of 30 m.

3.2.2. AST₂₀₁₈ Model Generation

Three stereoscopic models were extracted by means of digital photogrammetric software Photomod 4.2 [78] from the selected ASTER_L1A images (Table 1). The GCPs used in the generation of the KH9₁₉₇₉ model were also applied to solve the AT stage of AST₂₀₁₈. Then, between 50 and 60 tie points were manually added during relative orientation. Finally, the numerical adjustment was made and the mean absolute residuals of GCP were X = 5.41 m, Y = 3.99 m and Z = 10.84 m.

Once generated, the AST₂₀₁₈ model was filtered by applying different criteria. For this purpose, there are authors who choose to use global statistical thresholds either based on the mean and the standard deviation [79] or on any existing antecedents (e.g., [79,80,81]). In this study, the first filtering we applied consisted of calculating a minimum and maximum global threshold of −560 m and 190 m, respectively (Equation (1)), based on the elevation change information for the 2000–2016 period provided by Malz et al. [24]:

$$\mathsf{\Delta }{h}_{A-S} \in \left[-560,190\right]$$

(1)

where $\mathsf{\Delta }{h}_{A-S}$ is the height difference of AST₂₀₁₈ against SRTM₂₀₀₀. Subject to this condition, 1.8% of the elevations were classified as outliers and eliminated. Even after filtering, several anomalous values such as pits or bumps similar to the ones specified by Arefi and Reinartz [82] were observed, which showed that establishing global thresholds does not necessarily solve anomalies occurring in local environments.

Then, a second filtering was carried out, for the whole model, by applying a moving window criteria, similar to the one proposed by Lovell et al. [83], but using the median and the normalized median absolute deviation (NMAD) [84]. The size of the window had to cover at least twice the size of the anomaly to be filtered. On this basis, the decision was made to use a kernel size of 41 × 41 pixels. Finally, AST₂₀₁₈ was generated with a resolution of 30 m.

3.3. Digital Elevation Models Co-Registration

The co-registration process enabled the mitigation of systematic errors generated by the use of models of different origins. In this study, we tested two co-registration methods: M1 proposed by Nuth and Kääb [16] and M2 adapted from Li et al. [18]. For purposes of transformation parameters calculation, co-registration methods always require the selection and use of points belonging to both DEMs, which we will call co-registering points (CRP). The choice of CRPs is key for co-registration, as it requires identifying points in areas that have remained invariant between the dates of the DEMs to be matched [13]. M1 and M2 methods employed CRPs selected within areas defined as stable (off-glacier) terrain. In this study, the selection procedure of these points was different for each method (see details in below sections). Then, the parameters were obtained and applied to the DEM considered to be slave in order to align it with respect to the master DEM.

3.3.1. Method 1 Application

This method uses the elevation differences (dh) and relates them to the slopes ($\alpha$) and their orientations ($\psi$) through Equation (2) [16]:

$$dh=a\times \cos \left(b-\psi \right)\times \tan \left(\alpha \right)+\overline{dh}$$

(2)

where $a$ is the degree of horizontal displacement and $b$ is its direction.

With this method, displacements are obtained in East, North, and Up directions to be applied to the slave DEM [16]. For CRP_M1 selection, a random search was performed on a lithological information layer and outside the glacier limits according to Randolph Glacier Inventory 6.0 (RGI6.0) classification. Then, elevation differences were calculated as against the SRTM and CRP_M1 points were filtered.

The mean and the root mean square error (RMSE) are highly sensitive to gross errors in input observations [85]. Thus, we decided to use a robust statistical indicator such as the median ($\widetilde{\Delta h}$) as central measurement and NMAD as dispersion measurement. The NMAD results from multiplying the median absolute deviation (MAD, [86]) by 1.4826. Finally, the observations ($\Delta h$) of stable zones used in M1 were filtered by applying the condition stated under (3):

$$\Delta h \in \left[\widetilde{\Delta h}\pm g\times NMAD\right]$$

(3)

where $g$ factor is the Hampel identifier [87,88]. To filter $\mathsf{\Delta }h$ observations in this study, g values of 2.5 [89] and 3 were sequentially tested. The obtained shifts were similar, and as a result we decided to use $g$ = 3.

In our experience, RMSE convergence in the iterative process for DEM co-registration with M1 is not always clear, while the shift usually tends to a stable minimum. For this reason, we chose to follow the recommendations of Nuth and Kääb [16] to stop the process when RMSE drops less than 2% or the total shift is under 0.5 m.

3.3.2. Method 2 Application

This corresponds to an adaptation of the proposal by Li et al. [18] that consists of calculating seven transformation parameters among the reference systems of each DEM by using the CRP_M2 points. The particular feature of this method lies in selecting the CRP_M2 points since the use of centroids is proposed in subwatersheds delineated within continuous areas defined as stable without using the lithological layer. In this case, the CRP_M2 points are internally selected in the matching process. Each subwatershed with its associated centroid in the master DEM was matched or not with another one and its centroid in the slave DEM. For this matching, enhanced invariant moments [90] were calculated for each subwatershed from which the centroids were finally extracted. Apart from using the invariant moments, the spatial distance between centroids was added. In addition, the input DEMs resolution was decreased to 90 m to reduce the noise generated by the delineation. Hence, this allowed the matching of nearby subwatersheds. As an example, Figure 4 shows a workflow with the settings for how subwatershed delineation carried out in our study for a section of SRTM₂₀₀₀.

3.4. Analysis and Assessment of Co-Registration Procedures

The sources of elevation data and their participation in the methodological steps of the study are listed in

Table 4. Data sources and their utilization in the co-registration procedure and error assessment of M1 and M2.

Methods			Data
			KH91979	SRTM2000	AST2018	ICESat
Co-registration	M1		x	x	x	x
	M2		x	x
Error assessment	M1	Triangulation sum	x	x	x	x
		GMB co-registration error	x		x	x
	M2	Cross evaluation with M1	x	x

. First, all the DEMs and ICESat profiles used here were involved in the co-registration stage with M1. In this sense,

Table 5. DEMs and ICESat in co-registration role process. Note that SRTM₂₀₀₀→KH9₁₉₇₉ is the only co-registration transformation that could be reached with both M1 and M2 methods.

Co-Registration Role			Vector
Slave		Master
SRTM2000	$\to$	KH91979	SK
AST2018	$\to$	KH91979	AK
KH91979	$\to$	ICESat	KI
SRTM2000	$\to$	ICESat	SI
AST2018	$\to$	ICESat	AI
AST2018	$\to$	SRTM2000	AS

shows all possible combinations of data sources and their respective roles (slave or master). In contrast, for the M2 application, we selected the KH9₁₉₇₉ and SRTM₂₀₀₀ DEMs for co-registration (Table 5). The reason for this selection was that the SRTM₂₀₀₀ and KH9₁₉₇₉ DEMs were the only data sources with sufficient coverage and continuity in off-glacier areas for watershed delineation, which was indispensable for the application of M2. Next, in M1 Error Assessment all data sources and their resulting co-registration vectors were involved in the triangulation sum procedure. Then, the triangulation sum corresponding to KH9₁₉₇₉ - SRTM₂₀₀₀ - ICESat was selected for GMB co-registration error calculation. Finally, the error assessment of M2 was performed by comparing the KH9₁₉₇₉ and SRTM₂₀₀₀ DEMs co-registered with M1 in off-glacier terrain.

3.4.1. Method 1 Evaluation

ICESat elevation profiles and SRTM₂₀₀₀ model were used to evaluate M1. The triangulation method described by Nuth and Kääb [16] and Paul et al. [17] was used, which consists of co-registering combinations of three different databases from the ones available (KH9_1979, SRTM_2000, AST₂₀₁₈ and ICESat), thus, obtaining three co-registration vectors with components $\Delta \mathrm{X},\Delta \mathrm{Y}$ and $\Delta \mathrm{Z}$ for each combination. Therefore, in theory, the sum of the three resulting translation vectors should be equal to zero [17]. However, in practice, this is not usually the case and the total error calculated with Equation (4) is considered to be an estimator of co-registration uncertainty.

$${\mathsf{\epsilon }}_T=\sqrt{{\mathsf{\epsilon }}_x^2+{\mathsf{\epsilon }}_y^2+{\mathsf{\epsilon }}_z^2}$$

(4)

where ${\mathsf{\epsilon }}_x$ is the sum of the three $\Delta \mathrm{X}$ obtained from the combinations, ${\mathsf{\epsilon }}_y$ is the sum of the three $\Delta \mathrm{Y}$, and ${\mathsf{\epsilon }}_z$ is the sum of the three $\Delta \mathrm{Z}$.

DEM co-registrations with ICESat profiles were performed by using Table 3 data. DEM elevations for each ICESat footprint were calculated by bilinear interpolation [91,92,93].

Calculating ${\mathsf{\epsilon }}_T$ value in the case of AST₂₀₁₈–KH9₁₉₇₉–ICESat triangulation is key to analyze errors in elevation and GMB changes. Since our study period extends from 1979 to 2018, the ${\mathsf{\epsilon }}_T$ value in this triangulation represents the co-registration error ${\sigma }_{\mathsf{\Delta }{h}_{co-reg}}$ finally used in error propagation.

3.4.2. Method 2 Evaluation

M2 accuracy was evaluated by cross evaluation with M1. For this purpose, SRTM₂₀₀₀ was co-registered to KH9₁₉₇₉ by means of M2. Next, cross-sectional profiles were then generated including stable (off-glacier) terrain. Then, average elevation differences were calculated between SRTM₂₀₀₀ and KH9₁₉₇₉, making a comparison before and after co-registration with both methods. This evaluation methodology is an alternative to ICESat data triangulation, because watersheds can only be delineated on DEMs and not on elevation profiles.

3.5. Glacier Elevation Change and Geodetic Mass Balance

Based on the results, the decision was made to use M1 for GMB calculation. Then, considering ICESat data triangulation, the co-registration error was calculated and included in GMB error equation.

3.5.1. Calculation of ∆h and Use of Interpolation Methods

On the basis of AST₂₀₁₈ and KH9₁₉₇₉ models co-registered with M1, elevation changes were calculated in the glaciated zones of the study area. This calculation was made by creating a common point grid where, for each point, the values of both DEMs were extracted. Notwithstanding this, the initial model had to have its voids filled due to the low correlation of both optical DEMs, mainly in high zones above 1100 m a.s.l. In addition, we used SRTM₂₀₀₀ to digitize a contour map for the voids filling procedure.

There are different methods to fill these voids, such as constant, spatial, and hypsometric methods (McNabb et al. [34]). Among the spatial methods, there is bilinear interpolation, which consists of filling gaps in a DEM or in their differentiation [34] by independent linearization in the two perpendicular directions, i.e., latitude and longitude [94]. Examples of glaciological applications include those of Zheng et al. [95], Brun et al. [96], Podgórski et al. [97], and others. In addition, hypsometric methods are classified as local or global depending on whether the dataset is used over large areas (global) or data from a particular glacier are considered (local) [34]. There are widely used applications in glaciological studies (e.g., [2,98]). Hypsometric filling consists of taking a reference DEM, dividing the area into elevation bins, and assigning to each band a representative elevation change value, estimated from the parametrized dataset [99]. These bins are usually established every 50 m (e.g., [23,100]) or, in other cases, every 100 m (e.g., [4,101]). Then, on the basis of the data contained in each bin, the mean (e.g., [24,99]) or the median (e.g., [102,103]) is calculated for each one of them. Thus, each missing value is replaced with the mean or median of the bin to which it belongs.

Some advantages of one method over the other are in some cases related to the extent of the area to fill in the gaps. In this sense, the bias generated by the bilinear method turns out to be very small [104]. In addition, hypsometry is widely used although the performance can be unstable when applied to a wider dataset [104], but according to Mcnabb et al. [34] its use on a regional scale is well-founded. However, when voids are large, as might be the case with glacier upper zones, hypsometry can be a possible alternative due to the fact that there are no meaningful data to constrain other methods [104].

In our study, we decided to test the following four void filling methods: G1 (global hypsometric); G2 (global hypsometric combined with bilinear); L1 (local hypsometric); and L2 (local hypsometric combined with bilinear). We calculated the glacier hypsometry by means of 50 m elevation bins considering the final scale. For statistical purposes, we chose the median, calculated at a global level for G1 and G2. For L1 and L2 fillings, the median was calculated for each separate glacier. In G2 and L2 options, we combined the hypsometric filling with bilinear interpolation. The bilinear method was applied in front glaciers zones up to the equilibrium-line altitude (ELA) of each glacier, using the data obtained from [105]. The decision to combine bilinear and hypsometric methods was due to the fact that, in lower zones, there was higher data density and bilinear interpolation preserved variability better than the hypsometric method [35]; however, in the SPI, plateau data were scarce. Hence, by using these combinations we were able to decide which of the performances best fits our study.

3.5.2. Geodetic Mass Balance

Once the different $\Delta h$ values were obtained, specific values for mass balances $b$ were calculated following Equation (5), adapted from Thibert et al. [106].

$$b=\frac{\rho }{A_m}{\displaystyle \sum }_{j=1}^n\Delta h_j{A}_{grid}=\frac{\rho \times r^2}{A_m}{\displaystyle \sum }_{j=1}^n\Delta h_j$$

(5)

where $\rho$ is the density factor and $A_m$ is the mean glacier area calculated according to Equation (6)

$$A_m=\frac{A_{1979}+A_{2018}}{2}$$

(6)

where $r$ is pixel size, $\Delta h_j$ is the elevation change in “j” pixel, and $n$ is the number of pixels of the glacier area.

Mass balance values b were divided by the period of 39.58 years to obtain the results in m. w. e. a⁻¹. In addition, in the study area different values of density factor $\rho$ were used. For example, Dussaillant et al. [98] and Farías-Barahona et al. [101] used a conversion factor $\rho =850$ kg m⁻³ [107], whereas Jaber et al. [108] and Richter et al. [109] chose to use $\rho =900$ kg m⁻³ based on Cogley [110]. In our case and in order to compare our results with those of other authors, mass balances were calculated with both conversion factors stated as $b_{850}$ and $b_{900}$, respectively.

3.6. Error Estimation

3.6.1. Error in Glacier Areas

Different criteria exist to estimate ${\sigma }_{A_{1979}}$ and ${\sigma }_{A_{2018}}$ errors when calculating areas, which were estimated by calculating ±1/2 pixel buffer around each glacier. Some authors have considered that the glacier outline delineation error is ±2 pix [20] or ±1 pix [111] and, on this basis, have calculated a buffer area around the perimeter of the glacier. This error is also estimated at 5–10% of the glacier area [9,112,113]. We decided to estimate the area error by using ±1/2 pixel buffer, which probably entailed an error overestimation.

3.6.2. Error in the Mean Area

Mean area error ${\sigma }_{A_m}$ was calculated according to Equation (7) for each glacier. Mean area $A_m$ is stated as a function of the area covered by the glaciers in 1979 and 2018 (Equation (6)). Therefore, and assuming their independence, their error is determined as a function of ${\sigma }_{A_{1979}}$ and ${\sigma }_{A_{2018}}$. Consequently, the error stated by Equation (7) derives from error propagation in Equation (6) as follows:

$${\sigma }_{A_m}=\frac{1}{2}\sqrt{{\sigma }_{A_{1979}}^2+{\sigma }_{A_{2018}}^2}$$

(7)

3.6.3. Error in $\Delta \mathrm{h}$

The error in median elevation change ${\sigma }_{\widetilde{\Delta h}}$ was estimated by Equation (8) adapted from McNabb et al. [34]:

$${\sigma }_{\widetilde{\Delta h}}=\sqrt{{\sigma }_{\Delta h_{co-reg}}^2+{\sigma }_{\Delta h_{rand}}^2}$$

(8)

where ${\sigma }_{\mathsf{\Delta }{h}_{co-reg}}$ is the co-registration error obtained by ICESat data triangulation calculated with Equation (4) and ${\sigma }_{\Delta h_{rand}}$ is the random component of the elevation change error.

To calculate the random component of the elevation change error, the co-registered elevation differences in stable zones are used. This component is calculated by dividing some dispersion estimator of elevation differences, for example, RMSE [34] or STDV [22,23], by the square root of the number of independent observations (Equation (9)). In our case, NMAD was used as a dispersion estimator as follows:

$${\sigma }_{\Delta h_{rand}}=\frac{NMAD}{\sqrt{n}}$$

(9)

It should be noted that $n$ represents the number of independent observations that could be equal to the total number of pixels of each glacier. However, we should take into consideration that these observations are correlated [114]. In consequence, the number of independent observations results from incorporating this correlation as shown in Equation (10).

$$n=\frac{N}{\pi \times {\left(\frac{L}{r}\right)}^2}$$

(10)

where $N$ is the number of pixels of the glacier under study; $L$ is the distance of the auto-correlation, which is considered to be equal to 500 m [80]; and $r$ is the pixel size.

3.6.4. GMB Error Estimates

The GMB error was estimated by starting with the formulation of the mass balance calculation (Equation (5)) and following the same error propagation logic used for mean area error, thus obtaining the error estimate (Equation (11)). Equation (11) was an adapted formulation of the ones used by the glaciological studies of Braun et al. [4], Malz et al. [24], and Huber et al. [35] as follows:

$${\sigma }_b=\left|b\right|\sqrt{{\left(\frac{{\sigma }_{\rho }}{\rho }\right)}^2+{\left(\frac{{\sigma }_{A_m}}{A_m}\right)}^2+{\left(\frac{{\sigma }_{\overline{\mathsf{\Delta }h}}}{\overline{\mathsf{\Delta }h}}\right)}^2}$$

(11)

where ${\sigma }_{\rho }$ is the density factor error, here considered as being equal to 60 kg m⁻³ [107] for the $\rho$ values used (850 kg m⁻³ and 900 kg m⁻³), and ${\sigma }_{A_m}$ is the mean area error calculated in the previous section (Equation (7)).

3.6.5. Impact of Error Component Analysis

The process of calculating the GMB results implied estimating the errors for each glacier. There are few records of error impact analysis for each component involved in this mass balance calculation. For this reason, we decided to propose a methodology taking that into consideration, where Equation (11) is restated as Equation (12) as follows:

$${\sigma }_b=\left|b\right|\times k$$

(12)

where $\left|b\right|$ is the absolute mass balance value of a given glacier and $k$ is defined as shown in Equation (13):

$$k=\sqrt{a_1+a_2+a_3}$$

(13)

where the terms $a_1$, $a_2$, and $a_3$ correspond to ratios as a function of the errors in density ${\left(\frac{{\sigma }_{\rho }}{\rho }\right)}^2$, mean area ${\left(\frac{{\sigma }_{A_m}}{A_m}\right)}^2$, and elevation change ${\left(\frac{{\sigma }_{\overline{\Delta h}}}{\overline{\Delta h}}\right)}^2$, respectively.

In Equation (12) we can see that the mass balance error (${\sigma }_b$) is proportional to the absolute mass balance value $b$ multiplied by k, which can be defined as an error or proportionality index. If k is higher than 1, ${\sigma }_b$ error will be larger than mass balance. On the contrary, if k is between 0 and 1, ${\sigma }_b$ error will be smaller than mass balance.

This allowed us to define $k$ coefficients for each glacier and for each interpolation method (G1, G2, L1, and L2), and then, to compare them and establish their relationship with the control line k = 1.

The terms $a_1$, $a_2$, and $a_3$ also represent the contributions to k total error index. Thus, each specific contribution can be written in percentages as follows:

$$a_{\mathrm{i}\left(\%\mathrm{k}\right)}=\frac{a_{\mathrm{i}}}{k^2}\times 100, i=1,2,3$$

(14)

Percentage values $a_{1\left(\%\mathrm{k}\right)},a_{2\left(\%\mathrm{k}\right)},\mathrm{and} a_{3\left(\%\mathrm{k}\right)}$, thus, state the proportion to the square of k index and represent relative density, mean area, and elevation change errors, respectively.

4. Results

4.1. Co-Registration Evaluation

The results of the M1 evaluation based on ICEsat and SRTM data triangulation reported smaller-than-pixel values. Starting from the four possible combinations of the four data sources available, co-registration residual was calculated based on the scheme of Nuth and Kääb [16].

Table 6. Estimated uncertainty of co-registration Method 1. ${\epsilon }_X,{\epsilon }_Y, \mathrm{and} {\epsilon }_Z$ are expressed in meters. Each error shift vector component, i.e., ${\epsilon }_X,{\epsilon }_Y, \mathrm{and} {\epsilon }_Z$ is obtained following its corresponding error vector equation.

Triangulation	Error Vector Equation	${\mathit{\epsilon }}_{\mathit{x}}$	${\mathit{\epsilon }}_{\mathit{y}}$	${\mathit{\epsilon }}_{\mathit{z}}$	${\mathit{\epsilon }}_{\mathit{T}}$$/{\mathit{\sigma }}_{\mathit{\Delta }{\mathit{h}}_{\mathit{c}\mathit{o}-\mathit{r}\mathit{e}\mathit{g}}}*\left(\mathbf{m}\right)$
AST2018–SRTM2000–ICESat	$AS+SI-AI$	−14.0	−3.8	−3.1	14.8
SRTM2000–KH91979–ICESat	$SK+KI-SI$	22.1	7.4	4.4	23.7
AST2018–KH91979–ICESat	$AK+KI-AI$	5.6	0.2	1.7	5.9 *
AST2018–SRTM2000–KH91979	$AS+SK-AK$	2.5	3.4	−0.4	4.2

shows final error values where the last column was calculated through the formulation of ${\mathcal{R}}^3$ vector norm (Equation (4)) and corresponds to the method uncertainty. The comparison of ICEsat profiles reported higher values due to the smaller amount of available data and the scarce spatial distribution of profiles as compared with digital models. At the same time, it should be noted that the slopes of the terrain where areas were hit by the ICESat sensor resulted in <20° for 80% of the data. Finally, the co-registration error ${\sigma }_{\mathsf{\Delta }{h}_{co-reg}}$ used in GMB error calculation resulted from AST₂₀₁₈–KH9₁₉₇₉–ICESat triangulation, which was chosen because it is the one that connects ICESat data with the two models (AST₂₀₁₈ and KH9₁₉₇₉) involved in this calculation.

The M2 performance was evaluated by means of cross validations with M1. This was carried out by calculating MAD using the data of the profiles generated on the Upsala and Viedma glaciers (Table 7). Table 7 shows the numerical results of these profiles both for stable terrain and ice zones. The first columns (stable terrain) show the performance of both methods where differences should be close to zero. However, elevation differences in stable terrain dropped with both methods but were not equal to zero. For example, in the Upsala glacier, they decreased by 37% (M2, Upsala 3-3 profile) and the values in the Viedma glacier fell by up to 64% (M1, Viedma 4-4 profile). This remaining altimetry error between the SRTM₂₀₀₀ and KH9₁₉₇₉ DEMs can be attributed to their different origins (radar versus optical). Added to these are the differences in area quality and definition, with SRTM₂₀₀₀ being a more accurate and continuous model than KH9₁₉₇₉. The last columns (ice zone) enable the evaluation of the impact of applying, or not, a co-registration method prior to calculating ice thickness change. In this respect, in glaciated zones, height differences are initially in the order of 90 to 160 m. After co-registration, these differences were reduced at least by 17% (M1, Upsala 1-1 profile) and at most by 56% (M2, Viedma 3-3 profile).

Table 7. Glacier elevation change statistics obtained from profiles before and after applying co-registration methods. The MAD between SRTM₂₀₀₀–KH9₁₉₇₉ elevations is presented for each profile in Figure 5, separating stable from ice zone.

Glacier	Profile	Stable Terrain—MAD (m)			Ice Zone—MAD (m)
		Before	After		Before	After
			Method 1	Method 2		Method 1	Method 2
Viedma	2-2	42.1	19.1	15.8	93.6	66.0	43.4
	3-3	38.0	14.9	16.6	92.1	64.4	40.6
	4-4	39.3	14.0	14.4	103.0	75.4	50.8
	1-1	-	-	-	91.8	64.2	41.2
Upsala	3-3	26.6	14.0	16.8	154.8	127.0	120.7
	4-4	29.4	17.1	17.5	137.7	110.5	104.8
	5-5	34.4	15.6	16.9	103.6	76.0	64.3
	1-1	-	-	-	157.7	130.2	120.9
	2-2	-	-	-	136.0	108.5	98.8

Table 7. Glacier elevation change statistics obtained from profiles before and after applying co-registration methods. The MAD between SRTM₂₀₀₀–KH9₁₉₇₉ elevations is presented for each profile in Figure 5, separating stable from ice zone.

Glacier	Profile	Stable Terrain—MAD (m)			Ice Zone—MAD (m)
		Before	After		Before	After
			Method 1	Method 2		Method 1	Method 2
Viedma	2-2	42.1	19.1	15.8	93.6	66.0	43.4
	3-3	38.0	14.9	16.6	92.1	64.4	40.6
	4-4	39.3	14.0	14.4	103.0	75.4	50.8
	1-1	-	-	-	91.8	64.2	41.2
Upsala	3-3	26.6	14.0	16.8	154.8	127.0	120.7
	4-4	29.4	17.1	17.5	137.7	110.5	104.8
	5-5	34.4	15.6	16.9	103.6	76.0	64.3
	1-1	-	-	-	157.7	130.2	120.9
	2-2	-	-	-	136.0	108.5	98.8

4.2. GMB Error and Interpolation Methods Analysis

Applying different global and local interpolation methods (G1, G2, L1, and L2) had a direct influence on the mean elevation change result of each glacier, and thus, on the k index. Therefore, a k index associated with each interpolation method was calculated for each glacier. For this reason, analyzing k index variability between different interpolation methods and glaciers provides us with valuable information for decision-making purposes. Figure 6 shows the four k error indexes (kG1, kG2, kL1, and kL2) calculated for each glacier. To enable their comparison, they were ordered increasingly with respect to the mean glacier area and were represented by means of their identifier (glacier ID). In addition, the dotted line in the graph represents the control line ($k=1$), above which GMB calculation error is higher than the GMB value itself. In the case of global interpolation methods G1 and G2, more than 40% of the 28 glaciers under study are above this line.

As regards local interpolation methods L1 and L2, only 10% of the glaciers under study are above control line $k$. In addition, towards the left side of Figure 6 we can observe a slight increase in the number of $k>1$ indexes. This situation could indicate that the smaller the glacier area, the higher the error indexes, and could be related to the lower data density of high zones (plateau), which has a larger impact in the case of smaller glaciers. An example of this is the Delgado Bertrand glacier (ID = 5), which has a mean area smaller than 5 km² and $k>1$ considering any kind of interpolation method. Furthermore, this glacier shows the highest k value of all glaciers (k = 36, with G1). It should be noted that the k index of this glacier has remained outside the vertical scale of Figure 6. Consequently, Figure 6 has been escalated on the vertical axis up to k = 15 because that is where most of the glaciers under study are concentrated. On the contrary, larger glaciers such as the Upsala (ID = 28) and Viedma (ID = 27) glaciers have all of their k error indexes below the control line. Finally, due to the performance observed regarding the k index, the decision was made to present the GMB results obtained by local interpolation Method L2. Although local Methods L1 and L2 yielded similar results, Method L2 was chosen considering that it uses a combination of the bilinear interpolation and the local hypsometric methods for void filling, while Method L1 only utilizes the hypsometric method for filling.

Table 8. Percentual contributions from $a_1$. $a_2$, and $a_3$ to k index considering each interpolation method.

ID	G1			G2			L1			L2
	${\mathit{a}}_{1\left(\mathit{\%}\mathit{k}\right)}$	${\mathit{a}}_{2\left(\mathit{\%}\mathit{k}\right)}$	${\mathit{a}}_{3\left(\mathit{\%}\mathit{k}\right)}$	${\mathit{a}}_{1\left(\mathit{\%}\mathit{k}\right)}$	${\mathit{a}}_{2\left(\mathit{\%}\mathit{k}\right)}$	${\mathit{a}}_{3\left(\mathit{\%}\mathit{k}\right)}$	${\mathit{a}}_{1\left(\mathit{\%}\mathit{k}\right)}$	${\mathit{a}}_{2\left(\mathit{\%}\mathit{k}\right)}$	${\mathit{a}}_{3\left(\mathit{\%}\mathit{k}\right)}$	${\mathit{a}}_{1\left(\mathit{\%}\mathit{k}\right)}$	${\mathit{a}}_{2\left(\mathit{\%}\mathit{k}\right)}$	${\mathit{a}}_{3\left(\mathit{\%}\mathit{k}\right)}$
1	0.0	0.0	100.0	0.0	0.0	100.0	5.1	9.8	85.1	4.8	9.3	85.9
2	0.0	0.0	100.0	0.0	0.0	99.9	12.8	3.9	83.3	12.1	3.7	84.2
3	0.1	0.1	99.8	0.1	0.1	99.8	0.5	0.7	98.8	0.7	0.9	98.5
4	0.6	0.2	99.1	0.9	0.3	98.8	6.1	2.3	91.5	5.7	2.2	92.1
5	0.0	0.0	100.0	0.0	0.0	100.0	0.1	0.0	99.9	0.1	0.0	99.9
6	0.0	0.0	100.0	0.0	0.0	100.0	12.5	4.8	82.7	12.5	4.8	82.7
7	0.0	0.0	100.0	0.1	0.0	99.9	7.4	1.9	90.7	8.1	2.1	89.8
8	0.1	0.0	99.9	0.1	0.0	99.9	14.3	6.3	79.4	14.3	6.3	79.3
9	5.3	1.8	92.9	1.3	0.4	98.3	0.0	0.0	100.0	0.0	0.0	100.0
10	0.1	0.0	99.9	0.1	0.0	99.9	10.5	3.6	85.9	10.5	3.6	85.9
11	0.3	0.1	99.6	0.6	0.3	99.2	10.3	5.0	84.6	10.5	5.1	84.5
12	1.3	0.3	98.4	1.3	0.3	98.4	1.2	0.2	98.6	1.0	0.2	98.7
13	3.6	0.7	95.7	2.0	0.4	97.6	1.5	0.3	98.3	1.5	0.3	98.3
14	0.1	0.0	99.8	0.1	0.0	99.8	1.9	0.5	97.5	1.9	0.5	97.5
15	7.5	1.9	90.6	3.3	0.8	95.9	2.8	0.7	96.5	2.8	0.7	96.5
16	19.2	3.1	77.7	16.7	2.7	80.6	15.3	2.5	82.3	15.0	2.4	82.6
17	0.2	0.0	99.8	0.2	0.0	99.8	0.2	0.0	99.8	0.1	0.0	99.9
18	0.0	0.0	100.0	0.0	0.0	100.0	9.5	0.7	89.9	9.5	0.7	89.9
19	11.0	0.8	88.2	11.0	0.8	88.2	10.3	0.8	88.9	12.1	0.9	87.0
20	4.3	0.3	95.4	0.0	0.0	100.0	1.2	0.1	98.7	1.0	0.1	98.9
21	28.0	1.5	70.5	20.0	1.0	79.0	38.9	2.0	59.1	39.1	2.0	58.8
22	0.0	0.0	100.0	0.0	0.0	100.0	0.1	0.0	99.9	0.1	0.0	99.9
23	68.5	5.4	26.1	68.0	5.3	26.6	67.3	5.3	27.4	64.6	5.1	30.3
24	3.2	0.1	96.7	0.3	0.0	99.7	8.2	0.2	91.6	5.6	0.2	94.3
25	13.0	1.0	86.0	5.0	0.4	94.6	2.3	0.2	97.5	3.0	0.2	96.8
26	0.2	0.0	99.8	0.0	0.0	100.0	2.6	0.1	97.3	2.6	0.1	97.3
27	47.0	0.3	52.7	45.0	0.3	54.7	51.6	0.3	48.1	51.6	0.3	48.1
28	33.6	0.3	66.1	33.6	0.3	66.1	66.3	0.6	33.1	66.3	0.6	33.1

shows the results of Equation (14) for the contribution of the errors involved ($a_1$,$a_2$, and $a_3$) divided by the square of k error index according to each interpolation method. These contributions were calculated after differencing KH9₁₉₇₉ and AST₂₀₁₈ models co-registered with M1. Having said that $a_3={\left(\frac{{\sigma }_{\overline{\Delta h}}}{\overline{\Delta h}}\right)}^2$ represents the relative error in elevation change, we found that in 63% of the cases $a_3$ contributes 90% or more to k² value. It should be noted that co-registration error contribution is included in $a_3$.

4.3. Glaciological Results

Table 9. Median glacier elevation change and geodetic mass balance 1979–2018 period for various SPI glaciers. For a better understanding and comparison, all glaciers, except the Viedma glacier, were also named considering the RGI 6.0 ID and the names in GLIMS database.

ID (Current Study)	RGIId (RGI6.0)	glac_name (GLIMS)	Area 1979 (km2)			Area 2018 (km2)			$$\begin{gathered} \widetilde{\mathit{\Delta }\mathit{h}} \\ \left(\mathbf{m}\right) \end{gathered}$$			$$\begin{gathered} \frac{\widetilde{\mathit{\Delta }\mathit{h}}}{\mathit{\Delta }\mathit{t}} \\ \left(\mathbf{m}{\mathbf{a}}^{-1}\right) \end{gathered}$$			GMB850 (m w. e. a−1)			GMB900 (m w. e. a−1)
28	17.00172; 17.00395; 17.01268	Upsala	933.14	±	7.68	808.68	±	8.88	−85.46	±	4.26	−2.16	±	0.11	−2.07	±	0.18	−2.19	±	0.18
27	17.05076	Viedma *	484.65	±	3.08	443.02	±	3.98	−62.84	±	4.29	−1.59	±	0.11	−1.50	±	0.15	−1.59	±	0.15
26	17.04928	Spegazzini	119.42	±	2.14	116.64	±	2.84	−10.21	±	4.43	−0.26	±	0.11	−0.29	±	0.13	−0.31	±	0.14
25	17.04939; 17.04935	Bolados & Onelli	82.92	±	1.88	68.90	±	2.24	−11.20	±	4.51	−0.28	±	0.11	−0.81	±	0.33	−0.86	±	0.35
24	17.04903	SPI-121	69.56	±	1.04	66.08	±	1.19	−15.74	±	4.56	−0.40	±	0.12	−0.41	±	0.13	−0.43	±	0.13
23	17.04902	Ameghino	67.22	±	1.57	56.87	±	1.89	−94.62	±	4.57	−2.39	±	0.12	−2.31	±	0.22	−2.45	±	0.23
22	17.04951;	Agassiz	53.42	±	1.19	51.24	±	1.32	−2.27	±	4.65	−0.06	±	0.12	−0.13	±	0.26	−0.13	±	0.28
21	17.04908	Mayo	44.69	±	0.95	41.42	±	1.02	−54.65	±	4.73	−1.38	±	0.12	−1.14	±	0.14	−1.21	±	0.14
20	17.04915	SPI-44	41.02	±	0.98	38.42	±	1.17	−6.77	±	4.77	−0.17	±	0.12	−0.06	±	0.04	−0.06	±	0.04
19	17.00368	Cerro de Mayo Norte	26.62	±	0.60	25.15	±	0.80	−26.61	±	5.04	−0.67	±	0.13	−0.64	±	0.14	−0.67	±	0.15
18	17.04916	Aguilera	24.57	±	0.59	23.81	±	0.70	−23.43	±	5.10	−0.59	±	0.13	−0.40	±	0.10	−0.42	±	0.10
17	17.04959	Ante-cumbre Bertrand Sur	13.15	±	0.46	11.70	±	0.47	−2.60	±	5.74	−0.07	±	0.15	0.07	±	0.17	0.07	±	0.18
16	17.04906	-	11.05	±	0.43	10.23	±	0.43	−36.07	±	5.99	−0.91	±	0.15	−0.64	±	0.13	−0.68	±	0.14
15	17.04932	Heim	9.90	±	0.42	8.59	±	0.50	−14.94	±	6.16	−0.38	±	0.16	−0.34	±	0.16	−0.36	±	0.17

* The calculated area of the Viedma glacier was limited to the incomplete coverage of this glacier in 1979 images.

shows the glaciological results for glaciers with an approximate area >10 km² in 1979. It should be noted that the nomenclature of some glaciers contains more than one RGI code. This is because some glaciers have become independent of their tributaries by 2018, due to retreat since 1979 (e.g., Upsala glacier). The table shows the area changes between both dates, the elevation change, and the GMB, along with the associated errors.

Area calculations for 2018 were performed considering the area change undergone by some tributary glaciers as independent ice bodies. An example of this is the Upsala glacier which, considering its tributaries (RGI6.0-17.00395 and RGI6.0-17.01268) in 1979, showed an area reduction of 13.3 % and a thickness loss of −1.97 ± 0.11 m a⁻¹, being the glacier with the highest loss rate after the Ameghino glacier. In fact, according to area reduction ranking the Ameghino glacier shows a 15.4% loss with a thinning rate of −2.39 ± 0.12 m a⁻¹. The Viedma glacier also had its area reduced by 8.6%. We stress that, in the case of this glacier, its change was only partially calculated because KH9₁₉₇₉ DEM did not cover all the glacier area. In addition, GMB values are shown according to the two density values chosen. It can be observed that the Ameghino glacier is the one that shows the highest negative balance (−2.31 ± 0.22 m w. e. a⁻¹). The Upsala glacier, the largest among the ones studied here, reported a GMB₈₅₀ of −1.94 ± 0.18 m. w. e. a⁻¹. By contrast, SPI-44 glacier reached an annual loss of −0.06 ± 0.04 m w. e. a⁻¹, the lowest one that was detected.

Figure 7 shows the spatial distribution of elevation changes during the 1979–2018 period for the 28 glaciers under study. The elevation change map results from applying M1 co-registration and L2 interpolation, which yielded the results appearing in Table 9. It is worth stressing that the glacier outlines used to show the spatial distribution of $\frac{\Delta h}{\Delta t}$ signal in Figure 7 correspond to the year 1979. The highest loss signals in reddish tones denote greater ice thinning intensity. In addition, the positive thickness change values observed in the high zones (plateau) may correspond to anomalous signals. These values probably denote noise due to white low-contrast zones, interpolation, and lack of data. The median elevation change of the glaciated area is −1.38 ± 0.11 m a⁻¹, with a GMB₈₅₀ of −1.44 ± 0.15 m w. e. a⁻¹.

5. Discussion

5.1. Application of Co-Registration Methods for an Optimized GMB

The implementation of co-registration methods in the field of glaciology for GMB estimation is limited to a few alternatives. Among the methods that provide solutions based on 3D shift vectors, it is usual to apply the method of Nuth and Kääb [16], as illustrated by studies in the Andes (e.g, [4,22,25,115]). By contrast, the method that applies horizontal displacements [20] has been implemented also in the Andes [23,25]. Co-registration results between [20] and [16] applications show negligible differences according to Dussaillant et al. [25] and Falaschi et al. [23]. However, when dealing with large data volumes, the Nuth and Kääb [16] method presents a computationally faster performance [17]. Therefore, the application of these methods is simple, and the codes are available to be used. Instead, the application of three-dimensional transformation methods usually requires more effort because they involve closed codes built into the differencing process [10,29]. In our case, on the basis of the methodological description presented by Li et al. [18] we adapted the method to the case under study for GMB estimation.

It is essential to evaluate the accuracy of co-registration methods. This is usually done by testing stable (off-glacier) terrain [4,16,108]. However, the co-registration method error is not always calculated and taken into consideration in the GMB error equation (e.g., [116]). In general, the co-registration error is evaluated by using a statistical measure (mean or median) before and after co-registration, and then, a comparison is made. Note that our results using MAD such as the best error estimator are in line with those of other authors (e.g., [24,117]).

DEM co-registration errors derived from Hexagon images present an elevation error in the order of the tens of meters [75]. This is in turn related to the generation of DEMs by means of optical sensors which find some difficulty generating them in glaciated zones [118]. In turn, this is evidenced by the automatic estimation of elevations in low-contrast white zones [106,114]. Difficulties also arise with mountainous reliefs due to the shadows they generate and the geometrical differences between images. All of this hinders correlation calculation procedures [119], and therefore, jeopardizes the quality of the DEMs generated. In general, working with DEMs showing higher spatial resolution and data density leads to differences closer to zero in stable terrain. It is frequent to calculate co-registration accuracy by triangulation with other data sources such as ICESat or SRTM (e.g., [9,16,17]). In our case, using ICESat resulted in subpixel values (14.8 m, 23.7 m, and 5.9 m) higher than those arising from triangulation between DEMs (AST₂₀₁₈–SRTM₂₀₀₀–KH9₁₉₇₉ triangulation, Table 4). This is similar to the situation presented by Nuth and Kääb [16] regarding triangulation results with ICESat data and can be explained by the smaller number and aligned distribution of points. Furthermore, the values obtained can be accounted for by the fact that the accuracy of ICESat/GLAS data is highly susceptible to the influence of ground slope [120] and drops quickly in steep areas [121].

5.2. Relationships with Other Glaciological Studies

Over the last decades, the general trend of glacier mass loss due to climate change has accelerated all over the world [122]. The SPI is no exception as far as these changes are concerned and total loss values in the order of 10.7 Gt a⁻¹ have been reported during the 2000–2019 period [115]. According to White and Copland [123], the eastern SPI showed an area reduction of 269.35 km² at an average rate of 4.68% during the period 1976/79–2008/10. In the specific case of some glacier units, White and Copland [123] found, for example, that during the same period, the Upsala glacier lost 7.64% of its area just as the Ameghino glacier lost 7.24%, while the Viedma glacier showed a reduction of 1.73%, and the Onelli and Bolados glaciers were the ones that shrank the most (11.07%).

During the 2000–2012 period, the SPI showed a total area reduction of 1.17% [108]. According to these authors, between 2000 and 2012, the following ice loss percentage is estimated for these glaciers: Viedma (−0.77%), Upsala (−2.99%), Ameghino (−0.84%), Agassiz (−0.18%), Onelli (−2.62%), while the Spegazzini and Mayo glaciers did not show any area variations.

In our case, glaciers exhibited different area percentage reductions during the 1979–2018 period, for example, some of the variations included: Viedma (−8.59%), Upsala (−13.34%), Agassiz (−4.08%), Bolados and Onelli (−16.9%), Spegazzini and Heim (−3.16%), Mayo (−7.32%), and Ameghino (−15.39%). Although our study period does not coincide exactly with that of other authors, some consistency is observed in the values obtained, with the Ameghino and Upsala glaciers being the ones with the largest area loss percentage. In relation to annual elevation changes, these area losses correspond to the maximum thinning signals detected in the front zones of the Ameghino glacier (−13 m a⁻¹) and the Upsala glacier (−17 m a⁻¹) over the 1979–2018 period. This is in line with Malz et al. [24], who reported a maximum value of −19 m a⁻¹ for the Upsala glacier during the 2000–2015/16 period.

In the SPI, ice mass loss is induced by front ablation, which is dominated by the process of calving [115]. For the 1968/75–2000 period, Rignot et al. [124] found a loss of −13.5 ± 0.8 km³ a⁻¹ by extrapolating their data to the whole SPI. Jaber et al. [108] detected a total ice loss in the SPI of −13.386 ± 0.712 Gt a⁻¹ for the 2000−2012 period and of −10.674 ± 1.839 Gt a⁻¹ for the 2012–2016 period. In our study, the ice thinning signals observed (Figure 7) showed patterns consistent with those reported by Malz et al. [24] and Jaber et al. [108].

The average GMB of the glaciers under study within the Santa Cruz basin was −1.44 m w. e. a⁻¹ (−3.0 Gt a⁻¹) for an area of 2087 km² (year 1979), while for the whole SPI it was −0.954 m w. e. a⁻¹ according to Malz et al. [24] for the 2000–2015/16 period. In the specific case of the glaciers under study, they showed dissimilar values. Among the ones with more loss, there is Ameghino, with −2.31 ± 0.22 m w. e. a⁻¹, and then Upsala, with −2.07 ± 0.18 m w. e. a⁻¹. On the one hand, despite having undergone a sharp change in ice retraction and ablation starting from 2010, the Viedma glacier shows higher percentage values than Upsala, but during the study period it reported values of −1.50 ± 0.15 m w. e. a⁻¹. On the other hand, Agassiz showed one of the lowest values (−0.13 ± 0.26 m w. e. a⁻¹). Historical data collected by Aniya et al. [125] accounted for volumetric changes due to thinning in the front zone of the Upsala glacier of −3.6 m a⁻¹ between 1968 and 1990. In addition, Aniya and Sato [126] found values close to −2.3 m a⁻¹ on Ameghino front during the 1949–1993 period. In this case, an ice loss of −2.6 ± 0.3 m a⁻¹ and a retreat of 55 m a⁻¹ were reported [127]. According to the values estimated in our study, glacier behavior, as regards ice loss, is highly variable. In the case of the Ameghino glacier and according to previous studies, thinning rates have been steady, without much variation. Starting in 2000 and between 2016 and 2019, the Upsala glacier has shown an average ice loss rate of around −3.19 m a⁻¹ [24,108,115], different from our findings of −2.16 m a⁻¹. This difference could be explained, in the first place, by the fact that we considered a longer period (39.58 years) and, in the second place, because it is well-known that the SPI glaciers have shown a significant retreat and ablation acceleration over the last decades [109]. Hence, any studies covering recent decades could estimate higher rates.

6. Conclusions

In this study, we have presented a co-registration methodology comparison belonging to a group of methods based on horizontal displacements and 3D shift vectors (M1) and three-dimensional transformations (M2) adapted and calibrated to be used in a portion of the Upper Santa Cruz River basin, eastern SPI side. Applying M2 entailed some subwatershed outline delineation limitations due to the high cloudiness in area images. In general, the methods are considered to be robust; however, M2 generates false matchings, so we developed specific filters to enhance performance. By contrast, M1, which is open-source and easy to apply, does not entail major configurations and inter-software migrations. Although M1 also requires filtering, its application is easy and does not require programming. We used the median as a robust estimator because the average did not solve gross errors due to the presence of clouds. Despite their dissimilar application, both M1 and M2 provide similar results for DEM matching. However, if we consider the impossibility to calculate M2 error and the ease of M1 use, we suggest that the application of M1 is preferred.

Based on the results obtained, we stress the importance of evaluating and applying, whenever necessary, some co-registration method for surface matching. Using the statistical data of elevation profiles in ice zones, it was observed that failure to perform co-registration would affect elevation change results in a range of 27 to 56% for the Viedma Glacier and in a range of 17 to 38% for the Upsala Glacier. This shows the impact that the lack of a co-registration method could have on GMB results. It also evidences the need to evaluate the accuracy of the method to be used and incorporate it in GMB error calculation.

Considering the contrast difficulty generated in data provided by optical sensors, it was vital to use interpolators to enhance data density in high glaciated zones (plateau) with pixel voids. Out of the four methodologies implemented, we decided to use interpolation method L2 due to its performance in relation to GMB error impact.

In addition, the analysis and impact of errors made it possible to decide the use of the co-registration and interpolation method. In this respect, the term $a_3$ related to co-registration impact proved to be a decisive contributor to k index, which, in turn, represents the overall GMB error. By analyzing k index, a possible association was found in that the smaller the glacier area, the higher the error index. It should be remembered that co-registration error is common to all glaciers and is calculated for each glacier by dividing the elevation change by the median. The random error is, in turn, divided by the number of independent observations. A possible explanation is that the glaciers that change the most tend to be the largest ones, and this reduces the error. By contrast, in smaller glaciers, the same error turns out to be significant in the GMB result because of their lesser area and the fewer changes undergone.

Finally, the methods applied and analyzed made it possible to contribute highly reliable methodological instruments to obtain an accurate GMB estimation that reflects the glaciological process in the study area. The results obtained proved to be consistent with previous studies in the zone. In addition, for the first time, we estimated almost four decades of ice volumetric changes in the main calving glaciers of the Santa Cruz River basin.