Initial outlines for glacier limits and water bodies over the NPI area were obtained from online datasets. The Global Land Ice Measurements from Space (GLIMS) database and the Randolph Glacier Inventory version 5.0 (RGI 5.0) offer access to glacier outlines from previous studies and inventories over the Southern Andes region (Rivera et al., 2007; Davies and Glasser, 2012; Pfeffer et al., 2014; Mouginot and Rignot, 2015). Following visual assessment, RGI5.0 was selected as the best existing outline, corresponding to the ice extent of year 2001 (RGI6.0 was not available at the time of our study). RGI5.0 was not used “as is” but was improved as explained in section Improvement of the Glacier Outlines. Water body outlines were acquired from the USGS Earth Explorer SRTM water body data (earthexplorer.usgs.gov), and the ice divides derived from velocity fields (Mouginot and Rignot, 2015) were used to split the icefield into individual glaciers.

The RGI 5.0 glacier outline was corrected manually using late-summer and cloud-free LANDSAT, ASTER, and SPOT5-HRS images as references to adjust the position of the ice for the years 2000 and 2012. RGI 5.0 ice divides were replaced by the ice divides from Mouginot and Rignot (2015). A visual comparison between our final outlines and best existing previous ones can be observed in Supplementary Figure S2. Significant errors and gaps were found in the SRTM water body data, so manual corrections were also performed. Furthermore, a cloud mask was drawn manually to exclude the cloudy pixels over the SPOT5-HRS DEM. No cloud mask was created for ASTER DEMs because the lack of correlation over the clouds leads directly to data voids in the DEMs.

As explained above, two different methods were applied to derive the map of the rate of elevation changes over the NPI: SRTM and SPOT5 DEM differencing (SPOT-SRTM hereafter) and linear trend analysis of multiple ASTER DEMs (ASTER_trend hereafter).

The maps of elevation change rates (dh/dt) are obtained differently for the two methods. For SPOT-SRTM, a simple DEM difference is computed and divided by their time separation (12.1 yr) to calculate dh/dt_SPOT−SRTM. In the case of ASTER_trend, for each pixel, a first linear trend is fit through all the available elevation values (from all ASTER DEMs and the SRTM DEM) during 2000–2012. The final dh/dt_{ASTER_trend} is calculated from a second fit considering only the ASTER DEMs within the 95% confidence interval of the first regression line and excluding the SRTM DEM to avoid possible bias due to signal penetration. For the majority of the NPI area, the first valid elevation measurements used in this second linear fit come from a DEM acquired in April 2000, and the last measurements from two DEMs acquired in March 2012 and April 2012 (Supplementary Figure S3). We did not observe any obvious temporal concentration of ASTER DEMs. Thus, the time-stamp of the dh/dt maps from the two methods is similar, i.e., 2000–2012.

In the conversion from dh/dt to mass balance, only reliable pixels from glaciated terrain are considered. For ASTER_trend, we further exclude all pixels for which the uncertainty of the second linear fit is larger than 3 m yr⁻¹ (at the 95% confidence level).

NPI individual glaciers have shown different volume change rates during recent years, with differential responses within glaciers and a strong contrast between the eastern and western margins (Rivera et al., 2007; Willis et al., 2012a). To take into account these local effects, separate calculations were made for each glacier, and then the overall volume change rate and mass balance rate of the NPI were calculated as the area-weighted sum of every glacier. Each individual glacier is divided into 50 m altitude bands using the void-filled version of the SRTM DEM. For each altitude band, we compute the mean dh/dt after excluding all values lying outside of ±3^*NMAD (normalized median absolute deviation).This is an efficient way to exclude outliers (Berthier et al., 2004; Gardelle et al., 2013). The band average dh/dt is then weighted by its area to calculate the glacier-wide average dh/dt. In the rare cases where no measurement is available for an altitude band, the rate of elevation change is assumed to be 0.

Following Fischer et al. (2015), we considered the varying areas of each individual glacier between 2000 and 2012 to convert volume change rate to mass balance rate. We use a value of 850 ± 60 kg m⁻³ as the volume-to-mass conversion factor, which is reasonable regarding the duration of our study (>5 years) (Huss, 2013).

The ASTER_trend method allows rough estimation of the SRTM penetration depth (Wang and Kääb, 2015; Berthier et al., 2016). By extrapolating dh/dt to the date of acquisition of the SRTM DEM, a reconstructed SRTM DEM (called hereafter SRTM_rec) can be obtained. The differencing of SRTM_rec-SRTM maps the penetration depth of the C-band radar signal.

Formal uncertainties for mass balance estimates were calculated taking into account six main sources of error: errors on dh/dt, errors in glacier area, error on the density conversion factor, errors due to data voids, errors due to seasonal cycle sampling irregularities and departure from a linear trend.

The NPI had a total glacierized area of 3,856 ± 211 km² in 2000, and 3,740 ± 200 km² in 2012. The total area loss for the period 2000–2012 is 116 km² (3%), at a rate of −9.7 km² yr⁻¹.

Areal changes are heterogeneous and do not follow any obvious spatial pattern. Largest areal changes are observed on glaciers U4 (19%) located in the south-east, and glaciers Circo (12%) and Grosse (10%) located in the north of the NPI. The largest frontal retreats were observed in the tongues of two glaciers in the south-west, Steffens (3.3 km frontal retreat) and HPN1 (2.5 km retreat over its main tongue and even 3.1 km retreat over one of its tributaries), and San Quintin in the north-west (2.5 km frontal retreat). Glaciers Gualas, Grosse and Reicher in the north and glacier Acodado in the south also showed significant frontal retreats of the order of 1 and 2 km between 2000 and 2012.

Our area estimate for 2000 is in good agreement with previous studies, especially with Jaber et al. (2016), suggesting that our error bars are conservative (Table 2). Davies and Glasser (2012) obtained a slightly larger value for the NPI area in 2001 because they considered some of the rock nunataks as part of the icefield area. Similarly, Rivera et al. (2007) included adjacent small icy areas as part of the icefield. To our knowledge, we provide the first estimate of NPI area for 2012.

Glacier areal changes are not homogeneous through the NPI. In general the largest frontal retreats were observed for glacier tongues situated in the south-west of the icefield (glaciers Steffens and HPN1), also in agreement with previous studies (Rivera et al., 2007; Lopez et al., 2010; Davies and Glasser, 2012; Willis et al., 2012a). No convincing relationship was observed between glacier-specific area loss and mass balance. This conclusion is in agreement with earlier findings for Alaskan glaciers (Arendt et al., 2002).

We accounted for these area changes in our mass balance rate estimates. At the scale of the entire icefield, the effect is minor as the NPI mass balance rate is changed by 0.01 m w.e. yr⁻¹ if only the 2000 inventory is available. Yet, the mass balance rate can vary by up to 10% for the fastest retreating glaciers.

Some positive dh values were found above 2,900 m a.s.l. with the SPOT-SRTM method. These positive values are anomalous as, to our knowledge, no study documented thickening in the upper part of the NPI since 2000. Figure 7 shows the distribution with altitude of our estimate of the SRTM penetration depth (SRTMrec-SRTM), a side product of the ASTER_trend method. Positive values are observed above 2,900 m a.s.l. This is consistent with the apparent thickening observed above 2,900 m a.s.l. in the SPOT-SRTM dh/dt and thus suggest possible penetration of the SRTM C-band signal on the NPI snow located above 2,900 m a.s.l.

The SRTM_rec-SRTM and the SPOT-SRTM curves in Figure 7 do not match perfectly over this high-altitude zone. Likewise, the SRTM_rec-SRTM curve shows positive values between 1,000 and 2,900 m a.s.l. However, these differences remain within the uncertainties of the SRTM penetration depth, assumed to be roughly ±3.1 m in Berthier et al. (2016), and cannot be further interpreted. We further confirmed that this ±3.1 m is conservative by extrapolating dh/dt_{ASTER_trend} to reconstruct two DEMs at the dates of acquisition of the SPOT5 2012 DEM and of an additional SPOT5 DEM acquired 18 May 2005 (not used in the study for measuring elevation changes but with similar specifications as the SPOT5 DEM 2012). Differencing between SPOT5_{rec_2012}-SPOT5₂₀₁₂ and SPOT5_{rec_2005}-SPOT5₂₀₀₅ yielded mean dh on the NPI of 1.03 m and 0.14 m, respectively, both values well below the 3.1 m assumed error.

If we consider the region above 2,900 m a.s.l., the mean penetration depth is 10 m, a value similar to the 10 m C-band penetration depth found on the dry and cold firn of Greenland and Alaska (Rignot et al., 2001). Averaged over the whole icefield, the mean penetration depth is negligible (0.02 m), in agreement with a recent study in New Zealand (also in the southern hemisphere), where no penetration was found (Wang and Kääb, 2015).

Up to now, SRTM penetration depth has been considered negligible in Patagonia due to almost constant wet conditions at the glacier surface in February 2000 (Jaber et al., 2013). Here we show that penetration is actually occurring at the highest elevations, where temperatures are lower and firn is probably dry even in summer. Nevertheless, only 0.75% of the NPI area (29 km²) is found above 2,900 m a.s.l., a fraction small enough not to affect the overall mass balance rate of the icefield (the corresponding bias is <0.01 m w.e. yr⁻¹). Therefore, SRTM penetration can be assumed negligible over the NPI. Still, future geodetic mass balance studies in Patagonia using SRTM C-band DEM, should carefully consider the penetration depth, especially when a significant fraction of the glacier lies at high elevations.

The lack of SRTM penetration, except at the highest elevations, is in line with a backscattering analysis using the SRTM swath image data (SRTM-IMGR) in Jaber et al. (2016). They found low backscattering values over the plateau, indicating the presence of wet snow and consequent insignificant penetration depth. Less negative backscattering values are also observed over the highest zones, suggesting dryer snow where the C-band signal can penetrate.

Both methods led to similar estimates of the NPI mass balance rate for the 2000–2012 period (Table 3).

The SPOT-SRTM method has already been used in the southern Andes (Falaschi et al., 2016) but this method is all the more relevant now that we have verified the negligible penetration of the SRTM C-band signal at this latitude.

Even though individual ASTER DEMs have a higher noise level than the SPOT5 DEM when subtracted from the SRTM DEM (mean standard deviation of the residual for ASTER DEMs: 14 vs. 8 m for SPOT5 DEMs), their larger number led to final errors on stable ground very similar to those obtained with the SPOT-STRM method.

Icefield-wide mass balance error differs slightly between methods (±0.21 m w.e. yr⁻¹ for SPOT-SRTM and ±0.15 m w.e yr⁻¹ for ASTER_trend). This difference is mainly related to the amount of data gaps, higher in SPOT-SRTM. The distribution of data gaps coincides over the NPI plateau (~1,200–1,700 m a.s.l.) for both methods. Yet, over some glaciers in the north-east margins of the NPI, specifically Grosse, Bayo, and Exploradores glaciers, less than one-third of their area is covered with the SPOT-SRTM method due to clouds in the SPOT5 images. Coverage is also restricted on small steep glaciers like glacier Pissis in the south, where a high percentage of data gaps exists over the SRTM DEM. These latter glaciers are those for which the mass balance rate differs most between the two methods. The redundancy of the acquisitions used in the ASTER_trend method overcomes this issue.

The 2000–2012 Ṁ are in agreement within error bars between the two methods for all glaciers with coverage larger than 66%. Specific cases where Ṁ do not agree within error bars correspond to glaciers poorly sampled in the SPOT-SRTM method due to clouds or SRTM data gaps. The three largest glaciers of the icefield (51% of the total area) present Ṁ differences smaller than 0.2 m w.e. yr⁻¹.

Availability of the data is also important to consider when comparing the two methods. A strong advantage of the ASTER_trend method is that it is based on freely and nearly globally available ASTER L1A stereo images. SPOT5 stereo images were only acquired over specific areas and are not yet freely available. Thus, the ASTER_trend has the potential to be extended to wider areas.

Many studies already estimated the geodetic mass balance rates of the NPI for diverse time periods (Rignot et al., 2003; Aniya, 2007; Rivera et al., 2007; Glasser et al., 2011; Willis et al., 2012a; Jaber et al., 2016). However, only two of these studies examined the same time period as us and can therefore be directly compared with our results (Table 3). NPI-wide surface mass balances have also been estimated from atmospheric models (Schaefer et al., 2013; Mernild and Wilson, 2016), but these numbers cannot be compared with ours as they do not account for the mass loss by calving which is an important component of the mass loss for the NPI.

Generally, good agreement is found between our two estimates and previous geodetic estimates (Table 3). Although they used ASTER stereo images like us, Willis et al. (2012a) found a smaller mass loss. The difference may result either from the use of different ASTER DEMs, different coregistration methods, the amount of DEMs considered to extract dh/dt (55 in Willis et al., 2012a, 118 in our study) or the slight difference between time periods. We generated our own ASTER DEM using ASP, and they exhibit less artifacts than the DEMs from the 14DMO product used in Willis et al. (2012a) (Supplementary Figure S4).

Our mass loss estimates agree with Jaber et al. (2016) who compared the SRTM DEM with a TanDEM-X DEM mosaic. These two DEMs, derived from radar interferometry, are not affected by clouds and provide good coverage of the accumulation zone. Nevertheless, their coverage of the icefield was not complete, as northernmost NPI glaciers were not covered by TanDEM-X DEMs.

Error estimates in Jaber et al. (2016) and Willis et al. (2012a) are three to four times smaller than ours. One reason for our larger errors could be the factor of five that we conservatively allocated to errors over data gaps. However, mass balance differences of typically 0.2 m w.e. yr⁻¹ between our two methods for well-sampled glaciers (Figure 6, Supplementary Table S1) suggest that our error bars are reasonable and not too conservative.

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.