Slope estimation influences on ice thickness inversion models: a case study for Monte Tronador glaciers, North Patagonian Andes

Abstract

Glacier ice thickness is crucial to quantifying water resources in mountain regions, and is an essential input for ice-flow models. Using a surface velocity inversion method, we combine ice thickness measurements with detailed surface elevation and velocity data, and derive ice thickness and volume estimates for the Monte Tronador glaciers, North Patagonian Andes. We test the dependence of the inversion model on surface slope by resampling glacier slopes using variable smoothing filter sizes of 16–720 m. While total glacier volumes do not differ considerably, ice thickness estimates show higher variability depending on filter size. Smaller (larger) smoothing scales give thinner (thicker) ice and higher (lower) noise in ice thickness distribution. A filter size of 300 m, equivalent to four times the mean ice thickness, produces a noise-free thickness distribution with an accuracy of 35 m. We estimate the volume of the Monte Tronador glaciers at 4.8 ± 2 km³ with a mean ice thickness of 75 m. Comparison of our results with earlier regional and global assessments shows that the quality of glacier inventories is a significant source of discrepancy. We show that including surface slope as an input parameter increases the accuracy of ice thickness distribution estimates.

Introduction

Glacier ice volume is key to quantifying water resources in mountain regions and their possible contribution to sea-level rise. It is also crucial for linking surface and subglacial topographies, a prerequisite for ice-flow modeling (Farinotti and others, 2017). There are several methods for inferring the total volume of glaciers: volume–area scaling approaches (Bahr and others, 1997), parameterization schemes (Haeberli and Hoelzle, 1995) and physical models based on ice-flow dynamics and mass conservation (Farinotti and others, 2009; Morlighem and others, 2011; Gantayat and others, 2014). Recently, there has been an increase in the number of studies using various numerical inversion approaches to recover the distribution of ice thickness from surface measurements and glacier characteristics (Farinotti and others, 2017). Gantayat and others (2014) proposed a numerical model based on an inversion of the parallel flow approximation, using only glacier surface slope and surface velocity maps as input data.

The parallel flow approximation assumes that glaciers deform only by simple shear; with Glen's flow law exponent n equal to 3, the flux depends largely on ice thickness and surface slope (U _s ~ H ⁴α³). Thus, small changes in slope or glacier geometry can significantly alter the shear flow (Cuffey and Paterson, 2010). Longitudinal stress gradients (compression and extension) produce short-wave variations in surface topography which influence ice thickness distribution. Therefore, to meet the parallel flow approximation, it is important to smooth surface topography (Oerlemans, 2001; Gantayat and others, 2014), otherwise, ice thickness distribution errors could propagate exponentially (Bahr and others, 2014).

Using ground penetrating radar (GPR) ice thickness measurements, a detailed digital elevation model (DEM) and ice surface velocity maps, we implement a parallel flow inversion model to derive ice thickness distribution maps and total ice volumes for the Monte Tronador glaciers. We analyze the sensitivity of these models to surface topography (Gantayat and others, 2014), specifically to the size of the smoothing filter (i.e. the distance), used to resample surface slope. We then compare our results with previous ice thickness and volume estimates for the Monte Tronador glaciers, derived from non-calibrated ice thickness distribution models (Carrivick and others, 2016; Farinotti and others, 2019), and volume–area scaling approaches (Bahr and others, 1997).

Study area

Monte Tronador (41.15^°S, 71.88^°W, 3475 m a.s.l.; meters above sea level) is an extinct stratovolcano located in the North Patagonian Andes along the Argentina-Chile border (Fig. 1). Its upper slopes host one of the most extensive contiguous ice covers in the region (~57 km² in 2012; Ruiz and others, 2015). Based on morphological characteristics, the 13 Monte Tronador glaciers can be grouped into nine mountain glaciers, all located above 1400–1500 m a.s.l. (Alerce, Castaño Overa, Condor, Frías, Norte, Peulla, Mistral, Parra and Vuriloches), and four valley glaciers with debris-covered tongues descending to lower elevations (Verde, Casa Pangue, Manso and Blanco) (Ruiz and others, 2017).

Fig. 1. Glaciers and available GPR measurements at Monte Tronador (* indicates unofficial names given by Ruiz and others, 2017). Background image: false-color pan-sharpened Pléiades satellite image, 7 March 2012, PGO, CNES-Airbus D & S (Ruiz and others, 2015). Individual glacier limits are indicated by thin green lines. The thick, graduated blue lines show the location of the GPR profiles discussed in the text, including two transversal transects (A^′-A, B^′-B) over the debris-covered tongue of Manso glacier (Supplementary Fig. SM1). The star shows the location of the ice thickness observations by Rivera and others (2001).

Aside from recent studies on the geodetic mass balance and surface velocities of the Monte Tronador glaciers (Ruiz and others, 2015, 2017), little is known about other glaciological characteristics, and especially ice thickness. The only ground-based ice thickness measurement in the area was acquired in 2000 on the debris-covered tongue of Casa Pangue glacier (Fig. 1). Using a low-frequency impulse radar system, Rivera and others (2001) determined a mean ice thickness of 170 ± 10 m at 860 m a.s.l. More recently, in 2012, the Chilean Water Cadaster, Dirección General de Aguas, conducted a series of radar surveys onboard helicopters on the Chilean side of the study area, providing more comprehensive coverage of the ice thickness on Monte Tronador (Dirección General de Aguas, 2014).

Previous regional and global assessments of ice thickness distribution include the Monte Tronador glaciers. Carrivick and others (2016) modeled the volume and ice thickness for Patagonian glaciers ($41 {-}55^{\circ }{\rm S}$) using the perfect plasticity approximation. For Monte Tronador, they reported a mean ice thickness of 40 m with a maximum of 143 m, and a total ice volume of 2.6 km³. In their global assessment of ice thickness and glacier volume, Farinotti and others (2019) used different inversion models to recover the ice thickness distribution of all the glaciers in the Randolph Glacier Inventory (RGI) 6.0 (RGI Consortium, 2017). Their results for the Monte Tronador glaciers report a mean ice thickness of 61 m with a maximum of 262 m, and a total ice volume of 4.3 km³. None of the models in either study were calibrated with in situ ice thickness measurements, nor did they use ice surface velocity as input data.

Data and methods

Ice thickness observations

The most extensive ice thickness dataset was derived from a series of airborne surveys with a low-frequency GPR, operating at a central frequency of 25 MHz (Dirección General de Aguas, 2014). The survey was conducted during the summer of 2012 (Fig. 1), and the measurements add up to 23.8 km of GPR profiles on nine glaciers (Vuriloches, Verde, Peulla, Parra, Norte, Mistral, Frías, Casa Pangue and Blanco). The GPR system was controlled and operated from the helicopter cabin. The antennas were mounted on an aluminum structure hanging 40–50 m below the helicopter, and connected to the control unit via a fiber optic cable. The elevation of the antennas was measured in real time with a Differential Global Positioning System (DGPS) and laser altimeter. Dirección General de Aguas (2014) provide a detailed description of the methods used to retrieve ice thickness estimates from the raw radargrams.

For this study, the 2012 helicopter-borne measurements were complemented by 3.2 km of GPR surveys acquired on the debris-covered tongue of Manso glacier during the autumn of 2018 (Fig. 1; Supplementary Fig. SM1). We used a lightweight, low-energy consumption, GPR system, designed for on foot glacier surveys in hard-to-reach areas (Oberreuter and others, 2014). A 5 MHz center frequency antenna coupled to a DGPS was used to accurately determine the location of the radar profiles. We used Radan 6 software and a wave propagation velocity of 0.168 m ns⁻¹ to recover ice thickness from the raw radiograms.

Although Dirección General de Aguas (2014) carefully controlled the crossing of GPR profiles to validate their interpretations, they did not provide uncertainty estimates for their measurements. The theoretical vertical resolution of GPR data is typically assumed to be between a quarter and half of the electromagnetic wavelength (Sheriff and Geldart, 2006). Here, we take half the wavelength as an estimate, and calculate an uncertainty of 4 m for the airborne GPR measurements in the Dirección General de Aguas (2014) dataset, and 17 m for the ground-based measurements acquired on Manso glacier for this study.

To avoid autocorrelation between GPR measurements, the radar data were resampled to match the 16 m grid cells of the final model outputs. When a cell contained more than one GPR measurement, the mean value was used. After resampling the 27 km of GPR transects to a 16 m resolution, we recovered 1292 measurements over 10 of the 13 main glaciers on Monte Tronador. The mean ice thickness measured for each glacier ranges between 76 and 150 m, with an overall mean of 88 m. The maximum ice thickness varies between 140 and 210 m, with values over 150 m for Manso, Verde, Peulla and Blanco glaciers. The maximum ice thickness is typically found at an altitude of around 2000 m a.s.l., although on the Manso glacier tongue, the maximum of 210 m is found at lower elevation, around 1000 m a.s.l. (Fig. 1).

Ice thickness model description

The parallel ice-flow approximation (Eqn (1)) assumes that glaciers deform by simple shear, and therefore flow lines are parallel (Cuffey and Paterson, 2010):

(1)$$U_{{\rm s}} = U_{{\rm b}} + {2A\over n + 1}\tau_{{\rm b}}^nH\comma\; $$

where U _s and U _b are the surface and basal velocities, n is Glen's flow law exponent, A is the creep parameter (a measure of the smoothness of the ice), H is the ice thickness and τ_b is the basal stress (τ_b = fρgHsinα), where ρ is the ice density, g the acceleration due to gravity, f the shape factor and α the surface slope.

The distribution of ice thickness for the Monte Tronador glaciers was estimated by re-arranging the terms of Eqn (1) into:

(2)$$H = \left({\lpar U_{{\rm s}}-U_{{\rm b}}\rpar \lpar n + 1\rpar \over 2A \lpar fg\rho\sin{\alpha}\rpar ^n} \right)^{1/\lpar n + 1\rpar }\comma\; $$

where we use n = 3, A = 2.4 × 10⁻²⁴ Pa⁻³ s⁻¹ (glacier temperate ice; Cuffey and Paterson, 2010), a constant value for ρ = 917 kg m⁻³, g = 9.8 m s⁻² and f = 0.8 (Haeberli and Hoelzle, 1995), and assume basal velocity to be proportional to surface velocity (U _b = a × U _s), where a is the constant of proportionality.

Data input

Ruiz and others (2015) derived ice surface velocities for the Monte Tronador glaciers from three pairs of Pléiades satellite images acquired between March and May 2012. Gaps under 0.15 km² were filled using a bilinear interpolation method, resulting in a full coverage, 16 m resolution map, of mean annual surface velocities for all glaciers on Monte Tronador (Fig. 2b). Ruiz and others (2017) manually digitized glacier outlines (Fig. 1) from a panchromatic Pléiades ortho-image from 21 April 2012. Ice divides were then corrected using surface displacement vectors from the 2012 surface velocity map (Ruiz and others, 2015). Here, we estimate surface slope from a Pléiades DEM (PLEI DEM) generated from a triplet (back, nadir and front) of Pléiades images acquired on 21 April 2012. The derived slopes were resampled to a 16 m grid to stay consistent with the 2012 surface velocity map (Ruiz and others, 2015). Glacier outlines (Ruiz and others, 2017) were then used to mask ice-free terrain (Fig. 2a).

Fig. 2. Input data used for the model. (a) Ice surface elevation from the PLEI DEM. (b) Surface slope from a DEM resampled using the SSL approach at a spatial scale of 18×. (c) Surface velocity for 2012 (Ruiz and others, 2015). (d) Ice thickness estimate for a surface slope resampled at a spatial scale of 18×.

Experimental setup

We apply a resampling filter on the PLEI DEM to assess the influence of surface slope on ice thickness calculations. First, the original 16 m resolution PLEI DEM was resampled using a weighted average smoothing filter at different scales, from 1× (unfiltered, 16 m resolution) to 45× (equivalent to a cell size of 720 m). The resampled surface grid was then downscaled to the original 16 m cell size using spline interpolation. Finally, for each smoothed DEM, the surface slope was obtained using the Zevenbergen and Thorne (1987) algorithm implemented in SAGA GIS software (3.2 or higher).

The resulting smoothed slope (henceforth abbreviated SSL) grids are resampled at different spatial scales (1× to 45×) with elevation values depending on the local morphology of the cells, but have the same spatial resolution as the surface velocity map (16 m) (Fig. 2c and Supplementary Fig. SM2). Finally, following Gantayat and others (2014), we eliminated abrupt breaks in the ice thickness grids using a 3 × 3 cells smoothing filter. To test the influence of slope calculations, the model (Eqn (2)) was evaluated for each of the slope grids.

Since there are no basal velocity measurements for the Monte Tronador glaciers, we estimate basal velocities as a proportion of surface velocities (a = U _b/U _s). The parameter a for each glacier was obtained by iteratively testing the model (Eqn (2)) in 10 000 runs, with a between 0 and 1, and H values randomly selected from half of ice thickness measurements available for each glacier and each of the surface slope grids. The rest of the thickness measurements were reserved for model validation. We selected the best value of a for each glacier based on the smallest difference between model and ice thickness measurements in terms of the bias and the root mean square error (RMSE) of the residuals (Table 1}). For glaciers without ice thickness measurements (Parra, Alerce, Castaño Overa, Casa Pangue and Condor), we used a values from glaciers with similarities in size and morphology, as well as in ice thickness and surface velocities (Ruiz and others, 2015, 2017).

Table 1. Residuals of a with a resampling scale of 18×

Minimum and maximum values of a are for the whole set of spatial scales. ^aDue to similar morphological characteristics, the same a values as those derived for Mistral glacier were used for Parra and Alerce glaciers. ^bDue to similar morphological characteristics, a values for Peulla glacier were used for Castaño Overa and Condor. ^cDue to similar morphological characteristics, a values for Manso glacier were used for Casa Pangue glacier.

Ending in a proglacial lake, the front of Manso glacier has experienced a recent increase in surface velocity likely associated with enhanced basal sliding due to higher subglacial water pressure (Ruiz and others, 2015). Similar to Vieli and others (2000) and Stuefer and others (2007), to take this effect into account in our derivation of the parameter a, we allowed a to increase exponentially with decreasing elevation in the lower tongue of Manso glacier. In the accumulation area and the upper part of the glacier tongue, a was calculated iteratively as for the rest of the glaciers, but as it approaches the terminus, a increases exponentially, reaching a maximum value of 0.997 at the ice front.

The uncertainty of the ice thickness model was quantified by taking the partial derivatives of each variable in Eqn (2) (see Supplementary Eqn (SM1)). The accuracy of the model outputs for each slope grid (see Supplementary materials) was evaluated in terms of the differences between the modeled ice thickness (after calibration of the parameter a) and the remaining 50% (n = 646) of ice thickness measurements. The bias (BIAS), median (MED), RMSE, interquartile range (IQR) and confidence interval (CI) statistics were used to quantify the error distribution. Following Farinotti and others (2017), we also expressed the accuracy of the selected final model outputs relative to the mean of the ice thickness measurements (Table 2).

Table 2. Glacier volume and ice thickness of the Monte Tronador glaciers derived with the SSL approach at a resampling scale of 18 ×

Also shown are the values from Farinotti and others (2019) and Carrivick and others (2016). RMSE, root mean square error in meters; RMSE (%), RMSE relative to the mean of ice thickness measurements; IQR, interquartile range; CI, confidence interval for the percentile deviations from ice thickness measurements.

Results

Ice thickness distribution

Ice thickness distribution maps derived from the different SSL grids (1× to 45×) show the same general pattern (Fig. 3; Supplementary Fig. SM3), with mean total ice volume estimates of 4.6 km³, and a range of 0.35 km³ (Table 2; Supplementary material SM4). Smaller resampling scales (1× to 15×) show higher noise, with overall thinner glacier ice (mean thickness between 70 and 72 m), but higher maximum ice thickness estimates (> 360 m). On the other hand, larger resampling scales (30× to 45×) show smoother, thicker glacier ice (mean thickness between 75 and 80 m), and a maximum ice thickness between 240 and 264 m. Model runs with resampling filter sizes between 18× and 25× also show smoother ice thickness distributions, with a mean between 71 and 73 m, and a maximum between 254 and 282 m. For smaller spatial filter sizes (1× to 15×), areas of thicker ice (> 150 m) are found around small concave areas (accumulation spots) on almost all glaciers (Figs 3a, b). For larger spatial sizes (> 18×), thicker ice areas are located on the Norte and Frías mountain glaciers, and on the debris-covered tongues of Verde and Manso glaciers (Figs 3c, d). Although results are noisier at smaller resampling scales, the SSL approach yields a smoother ice thickness distribution without the presence of sharp boundaries. Accuracy estimates show no strong variations among the different calculated slopes, with BIAS between −1 and −4 m, RMSE between 29 and 49 m, and CI between 98 and 178 m (Table 2; Supplementary Figs SM6 and SM7).

Fig. 3. Ice thickness distributions derived from SSL resampled at different spatial scales (a) 1×: no resampling, (b) 10×: resampling scale of 160 m, (c) 18×: resampling scale of 290 m (this spatial size was selected as the most appropriate), and (d) 45×: resampling scale of 720 m.

Ice volume estimates

Total ice volume estimates for all Monte Tronador glaciers combined, as well as those obtained for each individual glacier, are relatively close at all resampling scales (Table 3). For simulations with different filter sizes (1× to 45×), total volumes for individual glaciers vary between 4 and 27% from the mean of all simulations combined. The largest relative differences are found on the smallest glaciers; Condor (0.5 km²), Parra (2.5 km²) and Alerce (2.5 km²), with relative differences of 27, 26 and 20%, respectively.

Table 3. Glacier volume and maximum ice thickness modeled for individual glaciers

Discussion

The numerical inversion model (Eqn (2)), based on the parallel flow approximation proposed by Gantayat and others (2014) to calculate glacier ice thickness from surface velocity and slope, has been shown to perform in a similar way to other models requiring more input data (Farinotti and others, 2017). The parallel flow approximation is valid when ice deformation is driven by simple shear stress alone (i.e. no compressional, extensional or other shear stress in a direction other than the xz plane along flow lines). Due to the exponential relationship between ice deformation and shear stress (n = 3), glacier flow is highly sensitive to slope changes (U _s ~ H ⁴α³; Cuffey and Paterson, 2010). To test the model sensitivity to the slope parameter, we use the incrementally smoothed PLEI DEM to calculate surface slope.

Uncertainty and sensitivity of the model

Gantayat and others (2014) noted that the uncertainty of the model (Eqn (SM1)) depends on the uncertainty in (1) the velocity estimates U _s and U _b, (2) the creep parameter A, (3) the shape factor f, (4) the density of ice ρ and (5) the surface slope estimate α.

Using the null test over motionless ice-free areas, Ruiz and others (2015) estimate that the uncertainty of the ice surface velocity data derived over three separate periods in March–May 2012 varies between 11 and 22 m a⁻¹. However, this estimate is highly conservative, as it includes forested terrain where cast shadow patterns introduce larger displacement errors. When considering non-forested areas, errors vary between 5 and 10 m a⁻¹. We decided to evaluate the uncertainty of the model using both these higher and lower error ranges. The sum of the independent errors for the three separate displacement grids gives an overall uncertainty for the 2012 mean annual surface velocities, of 26 or 12 m a⁻¹, whether forested areas are taken into account or not.

To quantify the uncertainty in U _b, we propagated the uncertainty in a, taking the maximum standard deviation for each glacier (0.2), and in U _s, and obtained an uncertainty of 50 or 23 m a⁻¹. Due to the high accuracy of the PLEI DEM (Berthier and others, 2014), the uncertainty in surface slope is supposed to be small. Nevertheless, since slope uncertainty depends on local topography, we assumed α uncertainty to be between 1^° and 3^°. Finally, following Gantayat and others (2014), we set the uncertainties in A, f and ρ to 8 × 10⁻²⁵ Pa⁻³ s⁻¹, 0.1 and 90 kg m⁻³, respectively.

Using maximum and minimum uncertainty values for the different components, we estimated that the mean uncertainties in the modeled ice thickness range between 107 and 57 m. Although, for low lying areas (where both slope and surface velocity are close to zero), relative ice thickness uncertainties could be much more considerable. The large range of uncertainties shows that the model is highly sensitive to surface velocity and slope inputs. For example, a 1 m a⁻¹ increase in the uncertainty of surface velocity represents an increase of 6% in ice thickness uncertainty. Meanwhile, a 0.1^° increase in the uncertainty of α results in a 1% increase in ice thickness uncertainty.

Ice thickness distribution and accuracy

Through the analysis of ice thickness distribution maps (Fig. 3; Supplementary Fig. SM3), together with the statistical evaluation of ice volume estimates (Table 2 Supplementary Figs SM5 and SM6), we evaluated the trade-offs between the different smoothing filter sizes. Overall, the accuracy of ice thickness estimates is similar for SSL grids smoothed at different spatial scales, with a bias of −0.9 to −3.5 m, and RMSE of 29–47 m, which gives a relative precision for the mean ice thickness estimates of 30–50%. While average thickness is also similar between filter sizes, as mentioned earlier, maximum thickness values and where the ice is thickest vary considerably depending on resampling scale. These variations in ice thickness distribution reflect the sensitivity of surface velocity inversion methods to surface slope estimates. The similarity in the accuracy estimates is best explained by the fact that most of the ice thickness measurements were obtained in areas of gentle slope, which do not change considerably between SSL grids at different spatial scales. Since error statistics do not show substantial differences in method performance, we suggest that care be taken when analyzing model results based on error statistics alone.

Our results indicate that larger spatial filter sizes (> 15×) perform better in recovering a smoother and more realistic ice thickness distribution. However, the spatial scale used to smooth the slope should be evaluated before the final results are retrieved. To account for the influence of longitudinal stress gradients and fulfill the parallel flow approximation, the slope should be smoothed over a distance between one to four times the ice thickness (Kamb and Echelmeyer, 1986) or, as suggested by Cuffey and Paterson (2010), more than three to four times. This rule of thumb could be used to search for the most appropriate distance over which to resample the slope (Gantayat and others, 2014). In the case of the Monte Tronador glaciers, considering a distance four times the mean ice thickness derived from our simulations (73 m), or the mean thickness calculated from GPR measurements (90 m), would suggest that a spatial filter size between 240 and 360 m (15× to 22×) would be enough to avoid the effect of longitudinal stress gradients. The 290 m (18×) smoothing filter is equivalent to four times the mean ice thickness, and yields a noise-free thickness distribution with the same pattern as larger filter sizes. Since the accuracy does not vary considerably among simulations at different scales, we suggest that the SSL at the 18 × scale is representative of the ice thickness distribution on Monte Tronador. The SSL method, or other smoothing approaches which conserve local topography, appear to provide a straightforward solution to calculate glacier slopes for inversion methods that require surface velocity as input data.

In 2000, Rivera and others (2001) conducted the first GPR ice thickness measurement on Monte Tronador, on the debris-covered lower tongue of Casa Pangue glacier (Fig. 1), and measured an ice thickness of 170 ± 10 m. Casa Pangue has experienced considerable mass loss between 2000 and 2012, over 80 m in 12 years (Ruiz and others, 2017). Subtracting the ice thickness measured by Rivera and others (2001) from the ice surface elevation from the 2000 SRTM (SRTM band-X; Ruiz and others, 2017) puts the bedrock elevation at this location at 690 ± 20 m (by quadratic sum of the SRTM and GPR measurement errors). In comparison, the bedrock elevation calculated from the ice thickness from the SSL approach (with the 18× smoothing filter) and the 2012 PLEI DEM is 703 ± 31 m. Although based on a single measurement, this simple analysis shows that the model performed within the estimated accuracy even for an area not calibrated through measurements.

Comparison with existing reconstructions

The recent availability of glacier inventories and advances in inversion methods have resulted in a large number of studies mapping the thickness distribution of glacier ice in different mountainous regions around the world (Farinotti and others, 2009, 2019; Huss and Farinotti, 2012), including the southern Andes (Carrivick and others, 2016; Millan and others, 2019). Among these studies, Carrivick and others (2016) and Farinotti and others (2019) reported the ice thickness distribution on Monte Tronador in sufficient detail to allow comparison with our results (Figs 4a, b).

Fig. 4. Comparison of ice thickness distribution for the glaciers on Monte Tronador in relation to glacier outlines from Ruiz and others (2017). Ice thickness distribution following (a) Carrivick and others (2016), (b) Farinotti and others (2019) and (c) using the SSL approach. (d) and (e) Close-up of the ice divide between Frías and Casa Pangue glaciers, with ice thickness derived with (d) the SSL approach, and (e) the Farinotti and others (2019) consensus ice thickness distribution model. Note the absence of rock outcrops and internal outcrops in (e), as well as the offset between the ice divides on the glacier outlines versus the thickness map.

Using a perfect plasticity approach, Carrivick and others (2016) modeled the ice thickness of the Monte Tronador glaciers (referred to as Vicente Pérez Rosales National Park glaciers), and reported a mean thickness of 40 m, a maximum thickness of 138 m, and a total ice volume of 2.6 km³, corresponding to almost half the volume according to our best estimates (SSL 18× 4.8 ± 2 km³).

Farinotti and others (2019) analyzed the results of three different inversion methods (Huss and Farinotti, 2012; Frey and others, 2014; Maussion and others, 2019) to recover the ice thickness distribution of glaciers included in the RGI 6.0 (RGI Consortium, 2017), and used the Glacier Thickness Database (GlaThida) 2.0 (GlaThiDa Consortium, 2019) as a source of ice thickness measurements to calibrate and validate model results. Although it is possible to analyze the results of the different models independently, we refer to the consensus composite result, a weighted mean of the different models used to estimate the global volume of glaciers. For the Monte Tronador glaciers, Farinotti and others (2019) reported a composite solution with a mean ice thickness of 55 m, a maximum thickness of 266 m, and a total ice volume of 4.3 km³, which is closer to our best estimates.

The statistical analysis of the differences in ice thickness between our SSL results and those of Carrivick and others (2016) reveals a larger negative bias and almost twice the error in the latter study (− 38 m bias and 50 m RMSE). The ice thickness distribution from Farinotti and others (2019) shows a better agreement with our results (RMSE of 31 m), nevertheless, it still exhibits a considerable negative bias (−16 m). Neither study used ice thickness measurements for model calibration, which could explain the differences between their results and those presented in this study. However, a detailed analysis of their ice thickness distribution maps reveals other sources of discrepancy worth mentioning (Fig. 4).

First, glacier representation (e.g. the glacier inventory used as input data) appears to be a major source of discrepancies. Carrivick and others (2016) used the inventory from Davies and Glasser (2012), which mapped the glaciers on Monte Tronador as a single contiguous ice mass. Although their inventory did not include the lower elevation debris-covered tongues of the Manso, Verde and Casa Pangue glaciers, they reported an ice-covered area slightly larger (by 65 km²) than ours (Fig. 4a). Based on this data, Carrivick and others (2016) suggested that all the glaciers on Monte Tronador, and their correspondent ice volumes, were above 1500 m a.s.l., while our results, consistent with Farinotti and others (2019), show that 12–14% of the ice volume is located below this elevation. The distribution of ice thickness with elevation has important implications for future projections of ice volume change.

Farinotti and others (2019) used the RGI 6.0 (RGI Consortium, 2017), which shows a better agreement with the glacier outlines from Ruiz and others (2017) used in this study (Figs 4b, c). Although no information on the date of the base image used to map the RGI glacier outlines is reported, the extent of the Manso and Casa Pangue tongues suggests a date prior to 2012. In addition to differences in the ice extent of the debris-covered glacier tongues, we identified other factors impacting ice thickness distribution maps. Since the RGI 6.0 does not include internal outcrops, Farinotti and others (2019) modeled ice thicknesses of up to 70 m in areas without current ice cover. The automatic methods for dividing ice-covered areas into individual glaciers described in Kienholz and others (2013), and implemented in the RGI 6.0 (RGI Consortium, 2017), generate a greater number of glaciers (24 vs 14) with their respective ice divides. Since the models used by Farinotti and others (2019) use glacier outlines and ice divides as an indication of zero ice thickness, these artificial divisions are located in sectors of ice up to 130 m thick (Figs 4d, e). The greater overall extent of glaciers and the inclusion of internal outcrops as ice-covered areas seem to compensate for the underestimation of ice thickness by the consensus model. However, thinner and smaller glaciers suggested by Farinotti and others (2019) will have implications for estimating the future impact of climate changes on these glaciers. Similar inconsistencies between the consensus model of Farinotti and others (2019) and ice thickness observations have been recorded by Millan and others (2019) along the northern and southern Patagonian Ice Fields.

Volume–area scaling for the Monte Tronador glaciers

One of the most widely used techniques for retrieving glacier volume is the so-called volume–area scaling method (Bahr and others, 1997). This empirical method needs to be calibrated with known glacier volume data (Bahr and others, 2015), however, most of the time it has been used without sufficient calibration data (Radić and Hock, 2010; Grinsted, 2013). To test the applicability of the volume–area scaling method in our study area, we performed two simple exercises.

First, we used the coefficients k (0.027) and γ (1.36) initially proposed by Bahr and others (1997) for the area–volume power law based on data for 144 glaciers. We then derived these coefficients empirically from the relationship between area and volume based on our results using the SSL approach at the 18× scale (Table 3). Plotting the area and volume for the glaciers in our study reveals values for the parameters k and γ of 0.0567 and 1.16, respectively. Using this method, the total ice volume on Monte Tronador is underestimated by 31% (3.25 km³) relative to our best results (Fig. 5).

Fig. 5. Volume–area relationship for the Monte Tronador glaciers (red dots). The coefficient of determination (R ²) and equation are based on data from our study.

Conclusions

Following the methodology of Gantayat and others (2014), we applied the parallel flow law of ice dynamics to estimate the ice thickness distribution for the Monte Tronador glaciers, using surface slope and ice surface velocity maps at a spatial resolution of 16 m. In order to test the sensitivity of the inversion model to surface slope, we analyzed different sizes of smoothing filters. While we found that the accuracy of ice thickness maps derived using different smoothing scales was similar (bias ~0 and RMSE ~31 m), the spatial distribution of ice thickness shows higher variability depending on filter size. The spatial distribution of ice thickness is greatly improved through the SSL approach, particularly using smoothing filters with sizes similar to the theoretical distance range suggested to account for the effect of longitudinal stress gradients.

Our best results suggest a total ice volume for the Monte Tronador glaciers of 4.8 ± 2 km³, with average and maximum ice thicknesses of 75 and 300 m, respectively. A comparison with recently published ice thickness distributions for these glaciers (Carrivick and others, 2016; Farinotti and others, 2019) shows that, even without considering the different modeling approaches, the accuracy of glacier inventories has a high impact on ice thickness distribution. This has considerable implications for ice-flow modeling studies, total glacier volume estimates and future quantification of water resources in mountain regions. The use of low-pass spatial filters to smooth glacier topography and surface velocity data could significantly improve inversion models, allowing better detection of glacier boundaries and more accurate reconstructions of ice thickness distribution.

We also tested a volume–area relationship approach to estimate total ice volume. Results showed that, using the empirically established coefficients proposed by Bahr and others (1997), the total volume could be underestimated by at least 31%. Therefore, care must be taken when using simple, uncalibrated approaches to recover glacier volume in areas without ice thickness measurements.

We conclude that SSL provides a significant improvement in feeding numerical inversion models when it comes to recovering the distribution of ice thickness from surface measurements or glacier characteristics. The use of different modeling approaches and a larger number of ice thickness measurements for calibration and validation is required to better quantify water storage in glaciers and estimate the impacts of climate change across the Andes.