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

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−1 (−3.0 Gt a−1) and an ice thinning median of −1.38 ± 0.11 m a−1 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−1, 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 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 2 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 2 in 1986 [40] and of 12,573 km 2 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 2 ), Argentino Lake Northern Branch subwatershed (1211 km 2 ), and Argentino Lake Southern Branch subwatershed (606 km 2 ). 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 2 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. Figure 1. Study area and ground control points (GCP) distribution. The glacier outlines presented here correspond to the glacier limits in 1979 according to Hexagon images. The hillshades on which the glacier outlines are mapped were made from the ALOS Global Digital Surface Model (AW3D30, https://www.eorc.jaxa.jp/ALOS/en/aw3d30/data/index.htm (accessed on 12 July 2020)). * All the IDs from RGI6.0 inventory start with the prefix "RGI60-", which has been omitted for space reasons. Table 1 summarizes the main optical products and accessories described below and used in this study.

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) were used to generate KH9 1979 DEM. These images are available through the United States Geological Survey (USGS, https: //earthexplorer.usgs.gov/ (accessed on 8 May 2017)).  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 2018 DEM for GMB calculation. Orthorectified images were also generated for glacier outline delineation in 2018.

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 2000 , was used for co-registration between DEMs and for their later evaluation.

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

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). 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. 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.

Ground Control Points
The fifteen ground control points (GCPs, Figure 1) used to generate the KH9 1979 and AST 2018 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).

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 2018 and KH9 1979 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. 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.

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 2000 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].

KH9 1979 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 1979 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 2000 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 1979 was generated from the point cloud filtered and interpolated from TIN with a pixel size of 30 m.

AST 2018 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 1979 model were also applied to solve the AT stage of AST 2018 . Then, between 50 and 60 tie points were manually added during relative orientation. Fi-nally, 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 2018 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]: where ∆h A−S is the height difference of AST 2018 against SRTM 2000 . 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 2018 was generated with a resolution of 30 m.

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.

Method 1 Application
This method uses the elevation differences (dh) and relates them to the slopes (α) and their orientations (ψ) through Equation (2) [16]: 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 ( ∆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 (∆h) of stable zones used in M1 were filtered by applying the condition stated under (3): where g factor is the Hampel identifier [87,88]. To filter ∆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.

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 2000 .

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. First, all the DEMs and ICESat profiles used here were involved in the co-registration stage with M1. In this sense, Table 5 shows all possible combinations of data sources and their respective roles (slave or master). In contrast, for the M2 application, we selected the KH9 1979 and SRTM 2000 DEMs for co-registration ( Table 5). The reason for this selection was that the SRTM 2000 and KH9 1979 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 1979 -SRTM 2000 -ICESat was selected for GMB co-registration error calculation. Finally, the error assessment of M2 was performed by comparing the KH9 1979 and SRTM 2000 DEMs co-registered with M1 in off-glacier terrain. Table 4. Data sources and their utilization in the co-registration procedure and error assessment of M1 and M2.

KH9 1979 SRTM 2000 AST 2018 ICESat
Co Cross evaluation with M1 x x ICESat elevation profiles and SRTM 2000 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 2018 and ICESat), thus, obtaining three co-registration vectors with components ∆X, ∆Y and ∆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.
where ε x is the sum of the three ∆X obtained from the combinations, ε y is the sum of the three ∆Y, and ε z is the sum of the three ∆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 ε T value in the case of AST 2018 -KH9 1979 -ICESat triangulation is key to analyze errors in elevation and GMB changes. Since our study period extends from 1979 to 2018, the ε T value in this triangulation represents the co-registration error σ ∆h co−reg finally used in error propagation.

Method 2 Evaluation
M2 accuracy was evaluated by cross evaluation with M1. For this purpose, SRTM 2000 was co-registered to KH9 1979 by means of M2. Next, cross-sectional profiles were then generated including stable (off-glacier) terrain. Then, average elevation differences were calculated between SRTM 2000 and KH9 1979 , making a comparison before and after coregistration 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.

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.

Calculation of ∆h and Use of Interpolation Methods
On the basis of AST 2018 and KH9 1979 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 2000 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.

Geodetic Mass Balance
Once the different ∆h values were obtained, specific values for mass balances b were calculated following Equation (5), adapted from Thibert et al. [106].
where ρ is the density factor and A m is the mean glacier area calculated according to Equation (6) where r is pixel size, ∆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 −1 . In addition, in the study area different values of density factor ρ were used. For example, Dussaillant et al. [98] and Farías-Barahona et al. [101] used a conversion factor ρ = 850 kg m −3 [107], whereas Jaber et al. [108] and Richter et al. [109] chose to use ρ = 900 kg m −3 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.

Error in Glacier Areas
Different criteria exist to estimate σ A 1979 and σ 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.

Error in the Mean Area
Mean area error σ 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 σ A 1979 and σ A 2018 . Consequently, the error stated by Equation (7) derives from error propagation in Equation (6) as follows:

Error in ∆h
The error in median elevation change σ ∆h was estimated by Equation (8) adapted from McNabb et al. [34]: where σ ∆h co−reg is the co-registration error obtained by ICESat data triangulation calculated with Equation (4) and σ ∆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: 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).
where N is the number of pixels of the glacier under study; L is the distance of the autocorrelation, which is considered to be equal to 500 m [80]; and r is the pixel size.

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: where σ ρ is the density factor error, here considered as being equal to 60 kg m −3 [107] for the ρ values used (850 kg m −3 and 900 kg m −3 ), and σ A m is the mean area error calculated in the previous section (Equation (7)).

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: where |b| is the absolute mass balance value of a given glacier and k is defined as shown in Equation (13): where the terms a 1 , a 2 , and a 3 correspond to ratios as a function of the errors in density , respectively. In Equation (12) we can see that the mass balance error (σ 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, σ b error will be larger than mass balance. On the contrary, if k is between 0 and 1, σ 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: Percentage values a 1(%k) , a 2(%k) , and a 3(%k) , thus, state the proportion to the square of k index and represent relative density, mean area, and elevation change errors, respectively.

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 shows final error values where the last column was calculated through the formulation of 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 σ ∆h co−reg used in GMB error calculation resulted from AST 2018 -KH9 1979 -ICESat triangulation, which was chosen because it is the one that connects ICESat data with the two models (AST 2018 and KH9 1979 ) involved in this calculation. Table 6. Estimated uncertainty of co-registration Method 1. ε X , ε Y , and ε Z are expressed in meters. Each error shift vector component, i.e., ε X , ε Y , and ε Z is obtained following its corresponding error vector equation.

Triangulation
Error Vector Equation ε x ε y ε z ε T /σ ∆h co−reg * (m) 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 2000 and KH9 1979 DEMs can be attributed to their different origins (radar versus optical). Added to these are the differences in area quality and definition, with SRTM 2000 being a more accurate and continuous model than KH9 1979 . 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 coregistration methods. The MAD between SRTM 2000 -KH9 1979 elevations is presented for each profile in Figure 5, separating stable from ice zone.

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 2 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 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 1979 and AST 2018 models co-registered with M1. Having said that a 3 = σ ∆h ∆h 2 represents the relative error in elevation change, we found that in 63% of the cases a 3 contributes 90% or more to k 2 value. It should be noted that co-registration error contribution is included in a 3 . Table 8. Percentual contributions from a 1 . a 2 , and a 3 to k index considering each interpolation method.     Table 9 shows the glaciological results for glaciers with an approximate area >10 km 2 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.

Glaciological Results
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 −1 , 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 −1 . 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 1979 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 −1 ). The Upsala glacier, the largest among the ones studied here, reported a GMB 850 of −1.94 ± 0.18 m. w. e. a −1 . By contrast, SPI-44 glacier reached an annual loss of −0.06 ± 0.04 m w. e. a −1 , 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 ∆h ∆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 −1 , with a GMB 850 of −1.44 ± 0.15 m w. e. a −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 2018 -SRTM 2000 -KH9 1979 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].

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 −1 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 2 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%).
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 −1 ) and the Upsala glacier (−17 m a −1 ) over the 1979-2018 period. This is in line with Malz et al. [24], who reported a maximum value of −19 m a −1 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 3 a −1 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 −1 for the 2000−2012 period and of −10.674 ± 1.839 Gt a −1 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 −1 (−3.0 Gt a −1 ) for an area of 2087 km 2 (year 1979), while for the whole SPI it was −0.954 m w. e. a −1 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 −1 , and then Upsala, with −2.07 ± 0.18 m w. e. a −1 . 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 −1 . On the other hand, Agassiz showed one of the lowest values (−0.13 ± 0.26 m w. e. a −1 ). 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 −1 between 1968 and 1990. In addition, Aniya and Sato [126] found values close to −2.3 m a −1 on Ameghino front during the 1949-1993 period. In this case, an ice loss of −2.6 ± 0.3 m a −1 and a retreat of 55 m a −1 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 −1 [24,108,115], different from our findings of −2.16 m a −1 . 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.

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 coregistration 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.
Author Contributions: For this research, P.V. and M.G.L. carried out the study conception, design, and writing of the manuscript; P.V. processed the data acquisition; P.V., M.G.L., A.V. and L.L. analyzed and interpreted the data and results; M.G.L., A.V., and L.L. reviewed the article for important intellectual content. All authors have read and agreed to the published version of the manuscript.

Data Availability Statement:
The study did not report any data.

Acknowledgments:
The authors are very grateful to Esteban Lannutti, Silvana Moragues, and Andrés Lo Vecchio for their collaboration as team members of the Andean Geomatics Lab, to Ana María Sfer for her support in statistics, and to Javier Carelli for his contributions in photogrammetry knowledge. We thank the Los Glaciares National Park for permits to explore the study area. P.V. has been awarded CONICET PhD fellowship #14120150102306CO. Note that this research is part of the PhD thesis of P.V. Special thanks go to the anonymous reviewers for their valuable comments and suggestions, which led to an improved final manuscript.

Conflicts of Interest:
The author(s) declared no potential conflict of interest with respect to the research, authorship, and/or publication of this article.