Heterogeneous spatial and temporal pattern of surface elevation change and mass balance of the Patagonian ice fields between 2000 and 2016

The northern and southern Patagonian ice fields (NPI and SPI) have been subject to accelerated retreat during the last decades, with considerable variability in magnitude and timing among individual glaciers. We derive spatially detailed maps of surface elevation change (SEC) of NPI and SPI from bistatic synthetic aperture radar (SAR) interferometry data of the Shuttle Radar Topography Mission (SRTM) and TerraSAR-X add-on for Digital Elevation Measurements (TanDEM-X) for two epochs, 2000–2012 and 2012–2016, and provide data on changes in surface elevation and ice volume for the individual glaciers and the ice fields at large. We apply advanced TanDEM-X processing techniques allowing us to cover 90 % and 95 % of the area of NPI and 97 % and 98 % of SPI for the two epochs, respectively. Particular attention is paid to precisely co-registering the digital elevation models (DEMs), accounting for possible effects of radar signal penetration through backscatter analysis and correcting for seasonality biases in case of deviations in repeat DEM coverage from full annual time spans. The results show a different temporal trend between the two ice fields and reveal a heterogeneous spatial pattern of SEC and mass balance caused by different sensitivities with respect to direct climatic forcing and ice flow dynamics of individual glaciers. The estimated volume change rates for NPI are $-\mathrm{4.26}\pm \mathrm{0.20}$ km³ a⁻¹ for epoch 1 and $-\mathrm{5.60}\pm \mathrm{0.74}$ km³ a⁻¹ for epoch 2, while for SPI these are $-\mathrm{14.87}\pm \mathrm{0.52}$ km³ a⁻¹ for epoch 1 and $-\mathrm{11.86}\pm \mathrm{1.99}$ km³ a⁻¹ for epoch 2. This corresponds for both ice fields to an eustatic sea level rise of 0.048±0.002 mm a⁻¹ for epoch 1 and 0.043±0.005 mm a⁻¹ for epoch 2. On SPI the spatial pattern of surface elevation change is more complex than on NPI and the temporal trend is less uniform. On terminus sections of the main calving glaciers of SPI, temporal variations in flow velocities are a main factor for differences in SEC between the two epochs. Striking differences are observed even on adjoining glaciers, such as Upsala Glacier, with decreasing mass losses associated with slowdown of flow velocity, contrasting with acceleration and increase in mass losses on Viedma Glacier.

1 Introduction

The northern and southern Patagonian ice fields (NPI and SPI) are the largest contiguous temperate ice bodies in mid-latitudes of the Southern Hemisphere. They stretch from 46.5 to 47.5^∘ S and 48.3 to 51.6^∘ S, respectively, along the main ridge of the southern Andes and cover areas of about 4000 and 13 000 km² (Davies and Glasser, 2012). The perturbation of the strong and consistent westerly flow caused by the Andes leads to one of the strongest precipitation gradients on Earth (Garreaud et al., 2013). Because the ice fields are located on the only significant land mass between 45^∘ S and Antarctica, they offer unique possibilities for studying the impact of changes in Southern Hemisphere westerly flow on glacier evolution and for inferring Holocene climate history from glacial evidence (Rasmussen et al., 2007; Lopez et al., 2010; Glasser et al., 2011; Davies and Glasser, 2012; Garreaud et al., 2013).

Precise, spatially detailed data on changes of glacier area and volume and on the mass balance are essential for establishing reliable relations between climate signals and glacier records in order to reconstruct the past climate and to develop accurate predictive tools of glacier response to climate change (Fernández and Mark, 2016; Marzeion et al., 2017). The dynamic adjustment of a glacier to changing external forcing does not happen instantaneously. In particular for calving glaciers the dynamic behaviour and mass balance may be largely decoupled from direct climate forcing (Benn et al., 2007; Åström et al., 2014). The main outlet glaciers of the Patagonian ice fields are tidewater or freshwater calving glaciers, showing heterogeneous patterns of changes in frontal position and hypsometry (Warren and Aniya, 1999). This stresses the need for spatially detailed geodetic repeat observations covering different epochs in order to resolve the complex pattern of glacial responses. High-resolution topographic satellite data from synthetic aperture radar (SAR) interferometry, as employed for the work reported in this paper, provide an excellent basis for handling these issues.

There has been a general retreat of SPI and NPI glaciers since the Little Ice Age (Davies and Glasser, 2012); however, there is considerable variability in magnitude and timing for individual glaciers. Only a few glaciers advanced intermittently during recent decades. The most striking case is the Pio XI Glacier showing a large cumulative frontal advance since 1945, recently including a general advancing period of its southern and northern branches starting in 2000 and 2005, respectively (Wilson et al., 2016).

Geodetic mass balance estimates of NPI and SPI have been derived from various sources. The first remote-sensing-based estimates of NPI and SPI volume change and mass balance were reported by Rignot et al. (2003), comparing Shuttle Radar Topography Mission (SRTM) elevation data of February 2000 with topographic maps of 1968 and 1975 to estimate the volume change of the 63 largest glaciers. Elevation changes measured at low elevations were fitted to a polynomial as a function of elevation in order to extrapolate the results to higher elevations where the topographic maps have large gaps. A similar analysis was performed for 20 glaciers of SPI, where 1995 cartography is available, and scaled to infer the volume loss of the entire ice field.

The SRTM digital elevation model (DEM) was used in several studies as baseline for deriving volume change of the ice fields during periods spanning the subsequent 10 to 15 years. Willis et al. (2012a, b) analysed the ice volume change of NPI between 2000 and 2011 and of SPI between 2000 and 2012 by comparing the SRTM DEM with time series composed of 55 DEMs of NPI and 156 DEMs of SPI, which were derived from data of the Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER), operating on the Terra satellite. The authors used all suitable ASTER scenes and applied a weighted linear-regression to the time series of elevations on a pixel-by-pixel basis in order to obtain the elevation change for the full time spans. Willis et al. (2012b) assumed a 2 m radar signal penetration bias for SRTM without taking into account the melting state of the surface. Introducing the same penetration bias, Willis et al. (2012b) also revised their previous volume loss estimate of NPI. This underlines the importance of correct treatment of radar signal penetration in case of interferometric DEMs. Dussaillant et al. (2018) determined the NPI volume change for 2000 to 2012 with two methods, on the one hand by differencing the SRTM DEM and SPOT-5 DEM from March 2012 and on the other hand by fitting pixel-based linear elevation trends over 118 DEMs calculated from ASTER stereo images acquired between 2000 and 2012. Ice-field-wide rates of volume change by both methods agree very well.

Foresta et al. (2018) exploited swath-processed CryoSat-2 interferometric data to produce maps of surface elevation change over the Patagonian ice fields and estimated the mass balance for 6 years between April 2011 and March 2017. The maps cover 46 % of the total area of NPI and 50 % of SPI, with large gaps on the termini of most glaciers. Relations between CryoSat-2 elevation change and surface elevation in the SRTM DEM were used to fill the gaps in the surface elevation change (SEC) maps.

The TanDEM-X (TerraSAR-X add-on for Digital Elevation Measurements) mission (hereafter TDM) (Krieger et al., 2007), composed of the two formation-flying radar satellites TerraSAR-X and TanDEM-X, has been operational since December 2010. The mission opened up excellent new capabilities for high-resolution topographic mapping of the global land surfaces including glaciers and ice sheets. Abdel Jaber (2016) generated DEMs based on comprehensive TanDEM-X datasets, over NPI mostly acquired in austral summer 2014 and over SPI between March 2011 and March 2012. DEM differencing versus the SRTM DEM yielded volume loss rates over NPI for the period 2000 to 2014 and over SPI for the period 2000 to 2011/2012. Abdel Jaber (2016) also presented detailed information on the methods used for DEM generation and calibration, on error assessment and on the analysis of radar backscatter signatures of SRTM and TDM drawing conclusions on radar signal penetration. For Jorge Montt Glacier, featuring the highest thinning rate on the ice field, the geodetic mass balance for 2011 to 2014 was also derived from TDM DEMs. Malz et al. (2018) performed DEM differencing of TDM DEMs of December 2015 versus the SRTM DEM over SPI. For the southern part of the SPI they further differentiated the SRTM elevation versus TDM DEMs from January to March 2012, autumn 2014 and winter 2016 and computed the TDM-based SEC for the epoch January–March 2012 to December 2015. Seasonal variations in surface elevation and of radar signal penetration were not taken into account for retrieving SEC rates.

The studies cited above use different datasets and methodological approaches and cover, at least to some extent, different epochs. This impairs the comparison of results, the evaluation of temporal trends and the analysis of commonalities and differences between the two ice fields. In order to tackle these issues, we derive volume change of NPI and SPI exclusively from bistatic SAR interferometry (InSAR) data (SRTM and TDM) for two multi-annual periods, spanning from 2000 to 2012 and from 2012 to 2016. Furthermore, for the main retreating glaciers of SPI an estimate of subaqueous ice loss is provided for the period 2000–2011/2012. The generation of surface topography products and the analysis of elevation change build upon methods developed by Abdel Jaber (2016) with various upgrades regarding TDM data selection and processing methods. The results presented here are completely novel, unless clearly stated. Particular attention is paid to the acquisition dates of the different DEMs, applying corrections for deviations of repeat DEM coverage from full annual time spans in order to avoid seasonality biases when deriving annual SEC rates. Significant effort is dedicated to the assessment of the wetness of the snow and firn surface, through a careful analysis of the backscatter of SRTM and TDM data, and to modelling and quantifying different sources of uncertainty.

The paper presents the first spatially detailed analysis of surface elevation change and the derived total net mass balance over the Patagonian ice fields for two different epochs based on the same observation technique, including a catalogue of volume change for the epoch 2000 to 2012 and 2012 to 2016 for all glaciers >2 km² of NPI and >9 km² of SPI. The results indicate a different temporal trend between the two ice fields and reveal the complex spatial pattern of SEC and mass balance as result of intricate interdependencies between direct climatic forcing and effects of ice flow dynamics.

2 Data

For the generation of surface elevation change rate (SECR) maps we rely exclusively on pairs of multitemporal bistatic InSAR DEMs. This technique provides wide-coverage surface elevation, overcoming issues affecting optical DEMs, such as lack of contrast on smooth snow or the presence of clouds, as well as the limited spatial coverage and resolution of altimeters. In this study we exploit data from TanDEM-X and SRTM, the sole Earth observation systems equipped with single-pass radar interferometers.

2.1 TanDEM-X

The primary objective of the TanDEM-X mission is the generation of a global, consistent DEM with high resolution and accuracy (Krieger et al., 2007). The main payload of the twin satellites is a SAR instrument operating at X-band (9.65 GHz), capable of a swath width of 30 km in the operational Stripmap single-pol (HH) mode. The global DEM product (DLR-EOC, 2018), whose performances are analysed in Rizzoli et al. (2017) and Wessel et al. (2018), is the result of the combination of 4 years of bistatic data acquisitions with different baselines and geometries. This product is hence not suitable for the derivation of surface elevation changes; nevertheless, we exploited the 0.4 arcsec (∼12 m) release as a reference DEM of the region for various processing aspects (Sect. 3), after proper editing of unreliable samples.

In this study we processed single, selected TDM bistatic raw datatakes into so-called Raw DEMs (Rossi et al., 2012; DLR-CAF, 2010) using the operational Integrated TanDEM-X Processor (ITP) (Breit et al., 2012; Fritz et al., 2011), in order to generate two elevation maps completely covering the ice fields in the years 2012 and 2015. The TDM data selection for each coverage was based on various criteria like the reduction of temporal span and of the number of datatakes, warm seasons to minimize SAR signal penetration, and small height of ambiguity (HoA), to reduce interferometric noise and similar imaging geometry. Ideally, data acquisitions should be at the end of the ablation season when the surface is at its lowest, but most importantly the two coverages should be acquired at the same time of year in order to minimize seasonal changes, which can be significant on the Patagonian ice fields. Since data availability restricted our ability to fulfil the last criterion, the residual temporal gap had to be compensated for.

The TDM acquisitions used to generate the two elevation maps are summarized in Tables S1 and S2 in the Supplement for NPI and SPI, respectively. The footprints of the individual Raw DEMs are shown in Fig. S1. The first elevation map is composed of descending acquisitions from austral summer 2012. An exception are termini of several western NPI glaciers where, due to data unavailability, we had to rely on an acquisition from May 2011. The second coverage is achieved with descending acquisitions from December 2015 (beginning of austral summer). On part of SPI we additionally processed three acquisitions from December 2011 acquired with the same geometry of the 2012 datatakes in order to measure seasonal elevation changes during summer. Three TDM datatakes from December 2015 (scenes 6 and 7 on NPI and 13 on SPI) feature a steep look angle (<27^∘) leading to increased layover.

For each Raw DEM the ITP provides additional geocoded rasters (Rossi et al., 2010; DLR-CAF, 2010) which were used in different phases of this study: height error map (HEM), uncalibrated SAR amplitude, backscattering coefficient, interferometric coherence and flag mask indicating critical areas.

2.2 SRTM

The SRTM (Farr et al., 2007; Rabus et al., 2003) was launched 11 February 2000 and produced over 9 d of acquisition a near-global DEM (60^∘ N–56^∘ S) with 1 arcsec (∼30 m) posting. The main payload was a bistatic C-band (5.36 GHz) SAR capable of a 225 km swath, achieved by applying the ScanSAR technique to four sub-swaths featuring horizontal and vertical co-polarization (HH, VV, VV, HH) and look angles between 30 to 56^∘. The large number of interwoven acquisitions at higher latitudes contributed to both absolute and relative accuracy as well as to reducing voids: the nine ascending and nine descending datatakes covering the Patagonian ice fields (Seal and Rogez, 2000) are listed in Table S3. The performance of SRTM, among others, was assessed by Rodriguez et al. (2005), Brown et al. (2005), Carabajal and Harding (2006), and Wendleder et al. (2016). The main issue is the presence of long wavelength height errors with magnitudes up to ∼20 m globally and spatial variation scales of hundreds to thousands of kilometres, mainly caused by residual roll errors due to the attitude adjustment manoeuvres of the space shuttle and by the applied absolute calibration of the sub-swaths.

The NASADEM (Crippen et al., 2016) is a new version of SRTM DEM, consisting of a complete reprocessing of the raw data, with improved phase unwrapping (significantly reducing voids) and an ICESat-based calibration, tackling issues such as limited absolute vertical accuracy and long wavelength height errors. In this study we used a provisional version of NASADEM (NASA JPL, 2018) as the elevation map from the year 2000 for both ice fields. The choice was done after comparing over a vast region surrounding the Patagonian ice fields, the NASADEM and the SRTM version 3 (SRTMGL1) (NASA JPL, 2013) to the TDM global DEM rescaled to 1 arcsec. The SRTMGL1 dataset, besides suffering from a vertical offset of ∼1 m against the reference (statistics are given in Table S4), displays a stronger presence of long wavelength elevation and geolocation biases (Δh images are shown in Fig. S2) and a higher root mean square (rms) when compared to the NASADEM. On the ice fields the differences between the two SRTM datasets are larger on NPI and in the south of SPI.

Furthermore, we retrieved the SRTM radar brightness images (SRTMIMGR) (NASA JPL, 2014) for the sub-swaths covering the ice fields (Table S3) with the purpose of assessing the melting state of the glacier surface. We also used the SRTM Water Body Data (SWBD) (Farr et al., 2007) for statistical and visualization purposes.

2.3 Glacier outlines

We relied on the Randolph Glacier Inventory (RGI) version 6 (RGI Consortium, 2017; Pfeffer et al., 2014), which contains improved basin divides of NPI by Rivera et al. (2007). We manually updated the RGI outlines at the glacier termini (including internal rocks) using the SAR amplitude, the DEM and optical images in order to reflect the exact extent of the glaciers at the time of acquisition of each elevation map (2000, 2012, 2015).

3 Method and error estimation

3.1 Methods for SEC and mass balance

3.1.1 Raw DEM processing

The use of ITP to process the Raw DEMs allows a great degree of flexibility with respect to processing parameters and algorithms. The beginning and end times of each scene were adapted (up to ∼30 s total length) in order to minimize the number of scenes and to include the widest possible ice-free terrain suitable for DEM co-registration (Sect. 3.1.2). The ruggedness of the topography of the study region with its steep mountains and intricate water bodies poses a significant difficulty for the ITP operational algorithms of phase unwrapping (Lachaise, 2015) and absolute height determination (Rossi et al., 2012). We therefore used an alternative algorithm of ITP (Lachaise and Fritz, 2016), which tackles both issues by exploiting an external reference DEM (Sect. 2.1).

The absolute phase simulated from the reference DEM is subtracted from the interferometric phase of the data. The fringe frequency of the differential phase is significantly lower and its unwrapping is unproblematic as long as elevation differences versus the reference DEM are not too large (a maximum of half of the HoA). The absolute phase of the data is then reconstructed by summing the phase simulated from the reference DEM to the unwrapped differential phase, thus removing any influence of the latter on the relative elevation in output. The output Raw DEM is finally obtained by geocoding in ITP the absolute phase of the data, implicitly determining an absolute phase offset (APO) value, on which the absolute height and the across-track position of the Raw DEM depend. ITP allows us to manually update the APO value and perform a new geocoding, allowing us to fine-tune the co-registration with a reference DEM, as described in Sect. 3.1.2.

3.1.2 DEM co-registration

The master and slave DEMs may be affected by vertical biases with respect to each other, these can be constant (offset), linear (tilt) or even varying with low frequency. They can furthermore be affected by horizontal shifts causing an additional slope- and aspect-dependent elevation bias in the SEC, which couples with the vertical bias, resulting in a systematic error with high potential impact on the volume change rate estimated over large areas. To obtain two consistent TDM elevation maps, co-registered to each other and to the SRTM DEM, we co-registered the single Raw DEMs and the SRTM DEMs of NPI and SPI to the reference DEM (Sect. 2.1).

An error in the APO of a TDM Raw DEM leads to a vertical height offset, an across-track horizontal shift and a tilt around the master flight trajectory (in order of impact, the latter being negligible in our Raw DEMs). These three effects are solved by fine-tuning the APO through an accurate estimation of the height offset versus the reference DEM and by repeating the geocoding with ITP. This method assures high precision by exploiting the geometrical parameters of the SAR acquisition and allows for avoiding critical aspects of the generic co-registration problem, accurately tackled by Nuth and Kääb (2011), such as estimation of horizontal shifts, interpolation, etc.

To estimate the height offset we manually selected a large number of calibration regions (CRs) over stable terrain around the ice fields relying on the SAR amplitude, the TDM slope and optical imagery. Tall vegetation was avoided because of physical changes and varying scattering phase centre at different incidence angles and radar frequencies. The CRs were chosen to be as flat as possible in order to isolate the actual vertical height offset. Layover and shadow regions were avoided as well as water pixels affected by low coherence. The footprints of the CRs are visualized in Fig. S1 and their features are summarized in Table S5.

From the elevation difference Δh between the reference DEM and the single Raw DEMs (or the SRTM DEMs of NPI and SPI), we computed on each CR with index r the mean μ_r, the standard deviation σ_r and the standard error of the mean:

$$\begin{array}{}\text{(1)}& {\mathrm{SE}}_r={\displaystyle \frac{{\mathit{\sigma }}_r}{\sqrt{N_r}}},\end{array}$$

where N_r indicates the number of spatially uncorrelated samples on CR r and was estimated through a semivariogram analysis as described in Sect. 3.3.2. A height offset estimate for each DEM was obtained through the weighted average

$$\begin{array}{}\text{(2)}& \mathit{\delta }{h}_{\mathrm{off}}={\displaystyle \frac{{\sum }_r\frac{{\mathit{\mu }}_r}{{\mathrm{SE}}_r^{\mathrm{2}}}}{{\sum }_r\frac{\mathrm{1}}{{\mathrm{SE}}_r^{\mathrm{2}}}}}.\end{array}$$

Values of δh_off ranged in magnitude between 0 and 1.8 m for the TDM Raw DEMs and were used to fine-tune the APO. For the NASADEM the height offset δh_off was subtracted, this was equal to 0.3 and 0.1 m in absolute value on NPI and SPI, respectively. Residual horizontal shifts with spatially varying magnitude might still be present in the NASADEM, their effect on the 0.4 arcsec SECR of the ice fields is nevertheless limited, given that 80 % of the ice field surface has a slope lower than 15 and 23^∘ on NPI and SPI, respectively.

Furthermore, range and azimuth tilts caused by baseline errors (Hueso González et al., 2010) were verified and found to be negligible for all TanDEM-X Raw DEMs. Height consistency between overlapping Raw DEMs was also checked in order to ensure a seamless elevation map.

The co-registration procedure partly compensates for the crustal uplift rates due to the glacial isostatic adjustment affecting the region, characterized by rates up to 40 mm a⁻¹ on the plateau of SPI and decreasing with distance, as reported by Dietrich et al. (2010).

3.1.3 Seasonal correction

Seasonal variations in SECR should be taken into account for deriving annual rates of surface elevation change if the time span of the repeat DEMs does not exactly match yearly intervals. Commonly, the mean daily SECR of the given time span is used for filling temporal gaps or for subtracting the contribution of excess days. This approach introduces a bias in annual SECR in case of seasonal variations. The magnitude of the bias depends on the percentage of missing (or excess) days and the amplitude of the seasonal cycle. The seasonal correction, as elaborated here, refers to the difference between mean annual SECR over epochs of 4 years (2012 to 2016) and 12 years (2000 to 2012), respectively, taking seasonal differences in SECR for missing days into account versus mean annual SECR without accounting for such differences.

The temporal mismatch versus exactly 4 and 12 years varies for the two epochs and for different sections of the ice fields. For 2012 to 2016 the impact of seasonal corrections is more important due to the shorter time span. For NPI the number of missing days corresponds to 3.6 % to 4.8 % of the 4-year time span for the two main tracks acquired in February 2012. A small section of NPI (covering the lower termini of San Rafael, San Quintín and Benito glaciers) was instead acquired on 28 May 2011 (NPI scene no. 1; Table S1, Fig. S1a). For the main sections of SPI the percentage of missing days ranges from 3.6 % to 5.0 %, except for a small sub-area where it is 7 %. For SPI the 2000 to 2012 mismatch in percentage of the full period ranges from 0.1 % to 1.0 % of the 12-year period. Nevertheless, we also applied seasonal corrections to this dataset. For the two tracks covering the main sections of NPI the gap corresponds to only 0.1% of the 12 years; given the limited surface covered by the third track, no correction was applied to this dataset.

Here we explain details of the seasonal correction for the epoch 2012–2016 because of their larger impact. Depending on the availability of additional TDM DEM data, the following procedures were applied for different sections of the ice fields.

• Three additional TDM acquisitions of December 2011 (Table S2) cover the southern, central and north-western sections of SPI (59.4 % of SPI, Fig. 8). These data were used to compute daily SECR over summer 2011/2012 by DEM differencing versus TDM data of March 2012 (99 d) covering the main part of SPI and versus TDM data of 31 January 2012 (33 d) covering a subsection in the south-east. The SECR maps of summer 2011/2012, scaled to the length of the missing period, were used for substituting the missing days in summer 2015/2016, which is a valid approach as the mean temperatures of the two summers agree within 0.1 ^∘C. Air temperatures of the European Centre for Medium-Range Weather Forecasts Interim Re-Analysis (ERA-Interim) (Dee et al., 2011; Berrisford et al., 2011) show, at the 850 hPa level for the grid point 47.25^∘ S, 73.55^∘ W (NPI), a mean summer temperature of 5.9 ^∘C in 2011/2012 and 6.0 ^∘C in 2015/2016; for point 50.25^∘ S, 73.55^∘ W (southern SPI), the values are 3.6 ^∘C in 2011/2012 and 3.5 ^∘C in 2015/2016.

• For SPI glaciers not covered by the summer 2011/2012 SEC map (except Pio XI and Jorge Montt glaciers), we used daily SECR, as a function of elevation derived from the 99 d SECR maps of summer 2011/2012 (blue curve in Fig. S3). For the section of Pio XI Glacier that is not covered by the summer SEC map we use the hypsometric SECR curve of low-loss glaciers (green curve in Fig. S3). On NPI, except on the termini covered by scene no. 1, we also used the hypsometric green curve of Fig. S3 since the majority of glaciers are not calving or the calving fluxes are a small component of total mass balance (Schaefer et al., 2013).

• For the three NPI termini covered by scene no. 1 and for Jorge Montt Glacier, which is not covered by the summer 2011/2012 SECR maps and subject to significant dynamic downwasting, we separate the SEC component in the ablation areas related to surface melt from SEC due to dynamic downwasting. For this purpose, estimates of the elevation dependence of the specific surface mass balance (SMB) during summer and the full year are needed. To our knowledge, until now the only multi-year time series of ablation measurements on any glacier of SPI and NPI, including the separation of summer and annual periods, has been performed on Perito Moreno Glacier (Stuefer et al., 2007). The applicability of the Perito Moreno mass balance elevation gradient has been checked by means of model output on SMB for NPI west coast glaciers (Schaefer et al., 2013) and mass balance data of Chico Glacier (Rivera, 2004), accounting for the west–east difference in temperature lapse rate across the ice field (Bravo et al., 2019). Further details are given in Sect. S4. The ratio between daily SMB-related SEC during summer versus the rest of the year, as a function of elevation, is used for estimating the increased SEC contribution due to surface melt during summer on the Jorge Montt terminus and for the reduced melt contribution during May to December for the three NPI termini. For the dynamic downwasting component we used the average SECR of the full period.

For the epoch 2000–2012 a seasonal correction was applied on SPI to the TDM elevations of 2012 to correct for the gap in the SRTM acquisition mid-date (17 February). The effect of the correction is small, being of main relevance for TDM scenes 6 and 7, acquired on 15 March 2012, covering ∼6000 km² with a temporal gap of 38 d. The two 99 d beams of the 2011/2012 summer SECR were used pixel-wise where available. The hypsometric mean (blue curve in Fig. S3) was used elsewhere, with a reduction by 20 % to account for the late summer season (mid-February to mid-March). This scaling factor for late summer is based on a time series of daily air temperature measurements from 1995 to 2003 near the front of Perito Moreno Glacier and ablation data (Stuefer et al., 2007).

Finally the correction rasters were obtained by scaling the daily correction rate pixel-wise by the temporal gaps in days (Fig. S4).

3.1.4 Derivation of SECR maps and estimation of mass balance

Two DEM mosaics were obtained for each ice field from the main TDM coverages by means of stacking, where the most reliable scene (evaluated through the height error, the look angle and the backscattering) was prioritized for overlapping regions. The WGS84 (EPSG: 4326) projection with posting 0.4 arcsec was enforced through cubic convolution on a common geographic frame. Corresponding mosaics of the additional geocoded rasters computed by ITP were also obtained, as well as the SRTM DEM and its error layer. For each ice field two SECR rasters including seasonal correction were obtained differencing the DEM mosaics for the epochs 2000–2012 (12 years) and 2012–2016 (4 years), with the end of summer as the reference start and end time of each epoch. To avoid biases of the mass balance we masked out artefacts due to phase unwrapping, layover, shadow, etc. (Sect. S3). The elevation of the water surface subject to frontal retreat, usually decorrelated, was manually edited in order to correctly capture the freeboard SEC.

By multiplying the average SECR with the corresponding glacier area over elevation intervals of 50 m, the altitude-dependent volume change rate (VCR) was computed. The reference elevation used for the hypsometry is the 2012 TDM DEM (small voids are filled with the global DEM), which is common to the two investigation epochs. The mass balance was computed over the entire ice fields, as well as on single glaciers defined by the updated RGI glaciers outlines. The maximum extent of each glacier, either at the beginning or at the end of the observation period, was used to spatially capture all changes.

We used a glacier-wide density of 900 kg m⁻³ for the conversion of the VCR to mass change rate. This value is commonly used for geodetic mass balance measurements and provides traceability for comparisons with other studies (Cogley, 2009). The main mass losses on the Patagonian ice fields refer to ice areas, and for the accumulation areas assumptions about changes of the vertical profiles of snow and ice density would be speculative.

3.2 Impact of radar penetration

A critical issue affecting InSAR-based elevation data is the penetration of the radar signal in dry snow and firn. In this case the scattering phase centre is situated below the surface, causing an elevation bias in the DEM (Dall, 2007), ranging from decimetres to metres at the C- and X-band. This represents an important source of local systematic error on the SEC and consequently on the resulting total net mass balance. The penetration depth depends on the microstructure and the dielectric properties of the snowpack, which are in turn strongly dependent on the liquid water content (LWC). Several models (Tiuri et al., 1984; Mätzler, 1987) show how at the C- and X-band the penetration depth already drops rapidly below 0.2 m with a LWC of approximately 0.5 %^vol. We used the backscattering coefficient σ⁰ as a proxy to assess the wetness status of the snow and firn (Ulaby et al., 2014; Mätzler, 1987). The C- and X-band radar return from the bare rough ice of the glacier termini is dominated by surface scattering so that penetration is not an issue here.

3.2.1 Assessment of TanDEM-X backscatter

The TanDEM-X sensor features an absolute and relative radiometric accuracy of 0.6 and 0.3 dB, respectively (DLR-CAF, 2013), allowing precise measurements of backscatter. For each Raw DEM, we processed with the ITP the geocoded backscatter image including the annotated noise contribution. This typically varies between −29 and −17 dB along the range direction and can thus have a significant impact on σ⁰ of weak scatterers such as smooth wet snow. The σ⁰ mosaics corresponding to the 2012, 2015 and 2011 DEMs are shown in Fig. S5. No masking of artefacts was applied.

TDM austral summer datatakes were chosen in order to increase the likelihood of imaging wet snow and firn. The mid-range look angle (θ_l) ranges between 35 and 45^∘, except for scenes 6 and 7 of NPI and scene 13 of SPI, which have steeper look angles (Tables S1 and S2). The satellite overpasses were at approximately 06:00 local time (UTC−4 h), which is generally the coldest time of the day, although the plateaus of NPI and SPI usually feature limited daily variations in air temperature due to the dense clouds and strong precipitation occurring most of the year (Garreaud et al., 2013; Schaefer et al., 2013).

On the plateau of NPI (covered by scenes 2, 5 and 6 in Table S1) ${\mathit{\sigma }}^{\mathrm{0}}<-\mathrm{18}$ dB dominates in our dataset up to approximately 2300 m in altitude in 2012 and 2015, confirming high LWC on most of the surface. Above this altitude ${\mathit{\sigma }}^{\mathrm{0}}>-\mathrm{10}$ dB can be found on limited areas (particularly in December 2015), implying the presence of dry snow. Some regions with $-\mathrm{15}<{\mathit{\sigma }}^{\mathrm{0}}<-\mathrm{11}$ dB are found in scene 2 (θ_l=38.4^∘) at altitudes below 2000 m. Given the season and time of day, these can possibly be explained by the formation of a refrozen crust layer on top of wet snow or firn, implying an offset of the scattering phase centre up to a few decimetres (Reber et al., 1987; Mätzler, 1987).

The backscattering of SPI is more heterogeneous compared to NPI. The 2015 coverage features ${\mathit{\sigma }}^{\mathrm{0}}\le -\mathrm{19}$ dB, revealing wet snow on large parts of the plateau (particularly on the western margin). The σ⁰ of the 2012 coverage is on average higher (especially in scenes 4 and 5 acquired at the end of March). Over the main parts of the plateau σ⁰ is still lower than −16 dB, an indication for wet snow, possibly covered locally by a thin frozen crust that would introduce only a small elevation bias. The December 2011 coverage displays values of ${\mathit{\sigma }}^{\mathrm{0}}<-\mathrm{18}$ dB, indicating wet snow, on most of the plateau. Some isolated regions with higher σ⁰ in the southern sector have been conservatively masked out in the 2011/2012 summer SECR prior to using this dataset for seasonal correction (Sect. 3.1.3).

Based on the analysis of the backscatter and of the SEC maps, we manually outlined regions on the plateau that we considered prone to signal penetration in each DEM mosaic (Fig. S5). The outlining was performed manually in order to avoid areas on rough ice surfaces and on radar fore-slopes, where high σ⁰ is not an indicator for signal penetration. We assigned a potential penetration height offset to each of these outlined regions according to their average σ⁰. The offsets are based on a relationship between σ⁰ and height offset between wet snow surfaces and dry snowpack derived from multi-seasonal TDM Raw DEMs of NPI (Sect. S5). The extent of the areas likely affected by penetration is rather modest because we selected SAR scenes with widespread surface melt. Depending on the ice field and date, these cover from 0.9 % to 2.3 % of the total accumulation area. The penetration height offsets are taken into account in the error budget (Sect. 3.3.3).

3.2.2 Assessment of SRTM backscatter

The SRTM absolute and relative radiometric accuracy nominal values are 3 and 1 dB, respectively (Farr et al., 2007). The SRTMIMGR product provides the radar brightness β⁰ at 1 arcsec, corrected for flat Earth, for all the sub-swaths acquired during the mission. Lacking the orbital parameters of each acquisition, we coarsely removed the flat Earth correction using the mid-look angle of each sub-swath, introducing an error up to ±0.6 dB, and computed the backscattering coefficient using the provided local incidence angle (θ_loc) mask as ${\mathit{\sigma }}^{\mathrm{0}}={\mathit{\beta }}^{\mathrm{0}}\cdot \sin {\mathit{\theta }}_{\mathrm{loc}}$.

Figure S6 shows the arithmetic mean (${\stackrel{\mathrm{‾}}{\mathit{\sigma }}}^{\mathrm{0}}$) and the standard deviation computed pixel-wise from the sub-swath σ⁰ images covering the ice fields (four to seven stacked pixels are usually found). The measure of spread supports the interpretation of ${\stackrel{\mathrm{‾}}{\mathit{\sigma }}}^{\mathrm{0}}$. While σ⁰ is similar for the HH and VV polarizations of the sub-swaths, variations of several decibels are induced by the wide range of look angles (30 to 56^∘) at parity with snow conditions (Ulaby et al., 2014). The ERA-Interim data show mean 850 hPa air temperatures between 3.3 and 3.8 ^∘C in February 2000 over the ice fields, though temporal variations in LWC due to changing meteorological conditions cannot be excluded during the 9 d of acquisition. On the other hand, variations due to the diurnal temperature cycle are unlikely given the time of the space shuttle overpasses (Table S3).

Values of ${\stackrel{\mathrm{‾}}{\mathit{\sigma }}}^{\mathrm{0}}<-\mathrm{22}$ dB denoting the presence of wet snow are found on large sections of the plateaus. In the north-western part of SPI, a west–east gradient is visible (Fig. S6a). Values of ${\stackrel{\mathrm{‾}}{\mathit{\sigma }}}^{\mathrm{0}}$ up to −18 dB are found in the 1800–1900 m range (the mean elevation of the plateaus) and up to −16 dB at elevations up to approximately 2300 m. These values may be attributed to wet snow with a rough surface (Nagler and Rott, 2000). Above 2300 m, ${\stackrel{\mathrm{‾}}{\mathit{\sigma }}}^{\mathrm{0}}$ reaches up to −12 dB (excluding steep fore-slopes). Here nocturnal freezing of the upper snow layer is more likely, implying a displacement of the scattering phase centre in the order of decimetres (Floricioiu and Rott, 2001; Reber et al., 1987). This analysis was previously presented by Abdel Jaber (2016).

Considering the analysis of ${\stackrel{\mathrm{‾}}{\mathit{\sigma }}}^{\mathrm{0}}$, it can be concluded that the SRTM elevations are not affected by a bias due to C-band radar signal penetration except for few areas at high elevations with higher likelihood of penetration. These have been outlined (Fig. S6a) and assigned an estimated penetration height offset (Sect. S5) taken into account in the error budget (Sect. 3.3.3).

The good agreement of our volume change rates from 2000 to 2012 over NPI with the results of Dussaillant et al. (2018), based on optical data, supports the validity of our approach and the conclusions regarding signal penetration. Regarding possible penetration issues of the SRTM DEM, the analysis of Dussaillant et al. (2018) also indicates lack of penetration of the C-band SRTM radar signal into snow and firn except for a region above 2900 m a.s.l.

3.3 Uncertainty of SECR and mass balance

This section reports on the estimation of the different error sources affecting the SECR maps and the mass balance computed with the geodetic method.

3.3.1 Random error

The random error of each SECR sample σ_SECR was computed pixel-wise as the quadrature sum of the random errors σ_h of the elevation samples of master and slave DEM, divided by the corresponding Δt in years:

$$\begin{array}{}\text{(3)}& {\mathit{\sigma }}_{\mathrm{SECR}}\left(x,y\right)=\sqrt{{\mathit{\sigma }}_{h_{\mathrm{m}}}{\left(x,y\right)}^{\mathrm{2}}+{\mathit{\sigma }}_{h_{\mathrm{s}}}{\left(x,y\right)}^{\mathrm{2}}}/\mathrm{\Delta }t\left(x,y\right).\end{array}$$

For TDM elevations the random error is given in the HEM raster, which expresses the interferometric standard error for each sample computed assuming a normally distributed error as (Rossi et al., 2010):

$$\begin{array}{}\text{(4)}& {\mathit{\sigma }}_h\left(x,y\right)={\mathit{\sigma }}_{\mathit{\varphi }}\left(x,y\right){\displaystyle \frac{h_{\mathrm{a}}}{\mathrm{2}\mathit{\pi }}},\end{array}$$

where h_a is the height of ambiguity and σ_ϕ(x,y) is the standard deviation of the interferometric phase, which depends on the coherence and on the number of looks (Lee et al., 1994). The HEM does not include any systematic error components (phase unwrapping errors, etc.), these are discussed in Sect. 3.3.3. Concerning SRTM, the NASADEM also comes with a corresponding height error map providing σ_h. Where performed (the section of SPI covered by the 2011/2012 summer SECR), the random error contribution of the pixel-wise seasonal correction (Sect. 3.1.3) was also added in quadrature. The resulting random error maps for the two epochs are shown in Fig. S8.

3.3.2 Spatial correlation and spatial averaging

The standard error (SE) of a spatial average of several SECR samples is given by $\mathrm{SE}=\mathit{\sigma }/\sqrt{N}$, where σ is the random error and N is the number of uncorrelated samples. To determine N, the spatial correlation of the SECR maps was estimated by means of semivariograms. Two different regions of interest (ROIs), both verifying the assumptions of first- and second-order stationarity, were selected on ice-free terrain. ROI 1 features a relatively flat topography similar to the one of the CRs (Sect. 3.1.2); ROI 2 features varied slope and aspect distribution, simulating the ice field topography. The empirical omnidirectional semivariograms obtained on the two ROIs for the TDM–SRTM and TDM–TDM SECR were furthermore fitted with an exponential model and are shown in Fig. S9. Among the model parameters reported in Table S7 the range of the semivariogram is an estimate of the correlation distance d_c of the SECR map, which was conservatively increased by ∼40 %, to account for possible higher slopes on the averaged regions, among other factors. For the TDM–SRTM and TDM–TDM SECR maps we used, respectively, d_c=120 m and d_c=60 m to compute the standard error of the mean height offset on each CR (Sect. 3.1.2) and, respectively, d_c=200 m and d_c=100 m to compute the standard error of the mean SECR on the elevation intervals. For the estimation of N, the theory of geostatistics was applied as in Rolstad et al. (2009) by integrating the exponential semivariogram model (they used a spherical model) in polar coordinates over a circular integration area A. The assumption of a negligible nugget (representing the uncorrelated component of the variance for the applied sampling interval) leads to the following expression for the number of uncorrelated samples N within A:

$$\begin{array}{}\text{(5)}& N={\left[-{\displaystyle \frac{\mathrm{2}}{\mathrm{9}}}{\displaystyle \frac{A_{\mathrm{c}}}{A}}\left(\mathrm{3}\sqrt{{\displaystyle \frac{A}{A_{\mathrm{c}}}}}{e}^{-\mathrm{3}\sqrt{\frac{A}{A_{\mathrm{c}}}}}+e^{-\mathrm{3}\sqrt{\frac{A}{A_{\mathrm{c}}}}}-\mathrm{1}\right)\right]}^{-\mathrm{1}},\end{array}$$

where $A_{\mathrm{c}}=\mathit{\pi }{d}_{\mathrm{c}}^{\mathrm{2}}$ is the correlation area. Equation (5) simplifies to $N=\frac{\mathrm{9}}{\mathrm{2}}\frac{A}{A_{\mathrm{c}}}$ for the common case where A≫A_c.

3.3.3 Systematic errors

Systematic errors are not reduced when spatial averaging is applied, they can hence have a significant impact on the mass balance of large areas. We defined four systematic error components.

1. An error linked to the co-registration to the reference DEM (Sect. 3.1.2) was defined for each Raw DEM and for the SRTM DEM of NPI and SPI as the interquartile range (IQR) of all μ_r values (mean of Δh on CR r) used to estimate the height offset of that DEM. This error ranges between 0.04 and 0.3 m. The corresponding systematic error component on the SECR ε_reg is obtained pixel-wise as the quadrature sum of the co-registration errors of the master and slave DEMs scaled by Δt in years, in a similar way to Eq. (3).

2. To account for signal penetration we used the penetration height offsets assigned to critical regions on each DEM mosaic (Sect. 3.2) as local systematic errors. Furthermore, a bulk systematic error of 0.1 m was assigned to all remaining pixels above 1000 m a.s.l. to account for undetected regions and for possible small offsets on refrozen upper layer of snow and firn (Sect. 3.2). The systematic error component on the SECR ε_pen was obtained analogously to ε_reg.

3. An additional bulk systematic error was assigned to all glacier samples to account for unmodelled sources (e.g. residual glacial isostatic adjustment effects, residual tilts, unmasked local errors due to phase unwrapping or layover, etc.). This source includes effects of the curvature-dependent SEC bias caused by the different resolution of the SRTM and TDM DEMs affecting small regions mostly at high elevation (Abdel Jaber, 2016). This additional error was set to 0.05 m for TDM, while for SRTM it was set to 0.2 m on SPI and 0.3 m on NPI to account for residual low-frequency elevation biases (Sect. 2.2). The systematic error component on the SECR ε_add was obtained analogously to ε_reg.

4. To compute the systematic error linked to the seasonal correction (Sect. 3.1.3), the previous three systematic error components (ε_reg, ε_pen and ε_add) were estimated separately for the summer 2011/2012 daily SECR and expressed in m d⁻¹. Here ε_add was increased by a factor of 1.5 to account for the different temporal coverage. All three components were summed in quadrature and conservatively further increased by a factor of 3.0 on extrapolated regions (north of SPI and NPI). A pixel-wise scaling by the number of corrected days and by the Δt in years was applied, leading to a fourth systematic error component on the SECR, ε_seas.

The total systematic error ε(x,y) of each SECR sample was obtained pixel-wise as (omitting (x,y) for compactness):

$$\begin{array}{}\text{(6)}& \mathit{\epsilon }\left(x,y\right)=\sqrt{{\mathit{\epsilon }}_{\mathrm{reg}}^{\mathrm{2}}+{\mathit{\epsilon }}_{\mathrm{pen}}^{\mathrm{2}}+{\mathit{\epsilon }}_{\mathrm{add}}^{\mathrm{2}}+{\mathit{\epsilon }}_{\mathrm{seas}}^{\mathrm{2}}}.\end{array}$$

The mean values of ε(x,y) and of its components for the four SECR maps are reported in Table S8.

3.3.4 Geodetic mass balance error

The geodetic method was applied to estimate the average SECR on separate elevation bins and the corresponding volume and mass change rates. The total error of the mean SECR on elevation bin b was computed by summing in quadrature the mean systematic error ${\stackrel{\mathrm{‾}}{\mathit{\epsilon }}}_b$ on bin b and the standard error SE_b of the spatial average on bin b (which is generally negligible compared to ${\stackrel{\mathrm{‾}}{\mathit{\epsilon }}}_b$), obtained as follows:

$$\begin{array}{}\text{(7)}& {\mathrm{SE}}_b=\sqrt{{\displaystyle \frac{\stackrel{\mathrm{‾}}{{\mathit{\sigma }}_{\mathrm{SECR}}^{\mathrm{2}}}}{N_b}}},\end{array}$$

where N_b is computed according to Eq. (5).

In the geodetic method the mean of the valid SECR samples of bin b is extrapolated to the unsurveyed area of the bin. On such gaps the total error was increased by a factor of 1.5 when computing the mass balance of a single glacier basin and a by factor of 3.0 when computing it on the entire ice field, to account for the across-basin variability of the SEC, particularly at lower elevations.

To calculate the volume change rate a 2 % error was assigned to the glacier area obtained from the updated outlines (Sect. 2.3). This value is higher than the RGI error suggested by Pfeffer et al. (2014) and in line with the empirical findings of Paul et al. (2013). The uncertainty of the density used for the volume to mass change rate conversion was set to ±36 kg m⁻³ (4 %). This number is based on an estimated uncertainty of ±17 kg m⁻³ for density in the ice areas and ±54 kg m⁻³ in the firn areas. In Sect. 4 the average SECR and VCR errors estimated for each bin are visualized graphically on the hypsometric plots, while the errors estimated for the entire ice fields and for the individual glaciers are reported in the results tables.

4 Results

The SECR maps of NPI and SPI after seasonal correction are shown in Fig. 1 for the two main epochs 2000–2012 (epoch 1) and 2012–2016 (epoch 2), along with the TDM DEM mosaic of 2012 used as hypsometric reference to analyse the elevation dependence. Unsurveyed areas in the SECR maps are relatively small and geographically evenly distributed, with the exception of the eastern margin of NPI in 2012–2016 because of layover caused by the steep incidence angle of scenes 6 and 7. In Table 1, the SECR, the volume change rate (VCR), the mass balance and the contribution to sea level rise are specified for the entire ice fields. Table 2 provides SECR and VCR for NPI glaciers larger than 2 km², Table 3 for SPI glaciers larger than 35 km² and Table S9 for SPI glaciers with area between 35 and 9 km². The tables also report the measured basin areas (based on the updated RGI glacier outlines) and the percentage of SECR coverage for the two epochs. The reference hypsometry and the distribution of unsurveyed areas are shown in Figs. S11 and S12 for NPI and SPI, respectively. The altitude dependence of SECR and VCR is shown in Fig. 9 for NPI and its main glaciers and in Fig. 10 for SPI and its main glaciers, while plots for additional glaciers are reported in Figs. S13 and S14. SECR and VCR are assembled in 50 m elevation bins using the surface of the 2012 TDM DEM as reference. The SECR averaged over each glacier basin is visualized in Fig. S10 together with the 2012 TDM DEM average surface elevation.

Table 1 Results over NPI and SPI for the two epochs. The reported area refers to the beginning of the epoch, the coverage of the SECR map is also reported. Subaqueous ice changes are not included.

Table 2 Average surface elevation change rate (SECR) and volume change rate (VCR) for NPI and its glaciers larger than 2 km² for the two epochs. The reported area refers to the beginning of the epoch; the coverage of the SECR map is also reported. Subaqueous ice changes are not included.

Table 3 Average surface elevation change rate (SECR) and volume change rate (VCR) for SPI and its glaciers larger than 35 km² for the two epochs. The reported area refers to the beginning of the epoch; the coverage of the SECR map is also reported. Subaqueous ice changes are not included. The list is continued for glaciers up to 9 km² in Table S9.

Figure 1 SECR maps of NPI and SPI for the two epochs (a) 2000–2012 and (b) 2012–2016. Unsurveyed areas are marked in yellow. (c) The TDM DEM of 2012 used as reference for the geodetic mass balance.

Figure 2 SECR of San Quintín: (a) 2000–2012, (b) 2012–2016.

Figure 3 SECR of San Rafael: (a) 2000–2012, (b) 2012–2016.

Figure 4 SECR of Upsala: (a) 2000–2012, (b) 2012–2016.

Figure 5 SECR of Jorge Montt: (a) 2000–2012, (b) 2012–2016.

Figure 6 SECR of Pio XI: (a) 2000–2012, (b) 2012–2016.

Figure 7 SECR of Grey: (a) 2000–2012, (b) 2012–2016.

Figure 8 Daily SECR of SPI during summer 2011/2012. Acquisition dates north of the green line: 18 December–26 March (time span: 99 d); between the green and blue lines: 7 December–15 March (99 d); south of the blue line: 29 December–31 January (33 d). Unsurveyed areas are marked in yellow (including regions of high backscattering which were masked out). The SECR was further filtered with a median and a smoothing filter of kernel size 9.

Figure 9 Surface elevation, volume and mass change rates (SECR, VCR, MCR) versus altitude in 50 m intervals for NPI and its main glaciers for epochs 2000–2012 (red) and 2012–2016 (blue). The hypsometric curve of 2011/2012 is shown in grey.

Figure 10 Surface elevation, volume and mass change rates (SECR, VCR, MCR) versus altitude in 50 m intervals for SPI and its main glaciers for epochs 2000–2012 (red) and 2012–2016 (blue). The hypsometric curve of 2012 is shown in grey.

NPI shows a similar pattern of elevation change during the two epochs, with the highest rates of thinning on the lowest sections of the glacier tongues, gradually decreasing up-glacier. Equilibrium state is reached on average at about 1800 m elevation (Fig. 9). On the south-western sector of the ice field and on San Quintin Glacier the thinning rates at elevations below 1200 m are slightly higher than in the northern and eastern sectors. All glaciers with an area larger than 20 km² show volume losses during both epochs except Leones Glacier, which reveals a modest increase in ice volume (Table 2). The volume loss rate of NPI increased from epoch 1 (VCR = −4.26 km³ a⁻¹) to epoch 2 (VCR = −5.60 km³ a⁻¹). The three largest glaciers (San Quintin, San Rafael, Steffen) account for 50 % of the NPI volume loss during epoch 1 and for 48 % during epoch 2. During both epochs the highest SECR at basin scale was observed on HPN 1 (VCR = −2.50 and −3.25 m a⁻¹, respectively). On all glaciers larger than 20 km², except Arco Glacier and Leones Glacier (with positive mass balance), the loss rates were higher during epoch 2. San Quintin Glacier (Fig. 2, Fig. 9) shows the highest increase in volume loss (VCR = −0.60 and −0.92 km³ a⁻¹). On San Rafael Glacier the loss rate increased slightly from epoch 1 to epoch 2 (VCR = −0.81 and −0.87 km³ a⁻¹) but the loss pattern changed (Figs. 3, 9). On the terminus below about 800 m a.s.l. the rate of surface lowering decreased, whereas in the upper reaches loss rates became larger.

On SPI the spatial pattern of surface elevation change is more complex and the temporal trend is less uniform. Contrary to NPI, the volume loss of SPI decreased from epoch 1 (VCR = −14.87 km³ a⁻¹) to epoch 2 (VCR = −11.86 km³ a⁻¹). The three glaciers Upsala, Jorge Montt and Viedma account for 45 % of the SPI volume loss in epoch 1 and for 58 % in epoch 2. On Upsala Glacier the rate of surface lowering decreased on the terminus from epoch 1 to epoch 2 (Fig. 4) associated with a slowdown of calving velocity. The losses increased on Jorge Montt Glacier (Fig. 5) and Viedma Glacier. Very high loss rates are observed on the lower sections of the Jorge Montt terminus, with SECR (for a 50 m elevation bin) of −16.7 m a⁻¹ during epoch 1 and −25.6 m a⁻¹ during epoch 2. The loss rates decrease gradually up-glacier, but the main sections of the accumulation area of these glaciers, up to elevations of 1800 to 2000 m, were affected by downwasting during both epochs (Fig. 10).

Other glaciers with volume loss rates >0.5 km³ a⁻¹ are located in the northern sector of SPI (O'Higgins, Bernardo, Greve, Tempano, Occidental), in the central–western sector of the ice field (HPS 12) and in the south-west (Tyndall Glacier). Average thinning rates up to 40 m a⁻¹ are observed on the terminus of HPS 12 Glacier during epoch 1. The HPS 12 terminus, flowing through a deep, narrow fjord, retreated by almost 5 km between 2000 and 2012 and by 4 km between 2012 and 2015. In epoch 2 the average SECR and VCR could not be estimated reliably because of significant gaps on the terminus in the 2015 DEM due to phase unwrapping errors (Fig. 1b); nevertheless, thinning rates higher than 40 m a⁻¹ could be observed at the 2012 front. Next to HPS 12 the highest loss rates at basin scale are observed on Jorge Montt Glacier (SECR = −4.01 and −4.95 m a⁻¹ during the two epochs) and on Upsala Glacier (SECR = −3.33 m a⁻¹ and −3.04 m a⁻¹).

The only glacier with positive mass balance in both epochs is Pio XI Glacier, showing a significant increase in VCR from 0.52 km³ a⁻¹ in epoch 1 to 1.26 km³ a⁻¹ in epoch 2 (Fig. 6). SEC rates in the elevation zones up to 1500 m a.s.l. were positive during both epochs. During epoch 1 the elevation zone between 100 and 400 m a.s.l. accounted for the main contribution to the total gain in ice mass. During epoch 2 an additional source of significant mass gain was the elevation zone between 1000 and 1500 m on the ice plateau (Fig. 10).

On the western sector south of HPS 12 (49.6^∘ S) and on the eastern sector south of Upsala Glacier (49.9^∘ S) the average loss rates are smaller than on the northern sector, but all glaciers covering areas >35 km² and the majority of smaller glaciers show negative SECR during epoch 1 (Table S9 and Fig. S14). On the main ice plateau the surface elevation was either stable or the SECR was slightly negative during epoch 1, becoming slightly positive during epoch 2. During epoch 2 the mass balance of several glaciers of the southern sector switched from negative to slightly positive values. However, the termini of the majority of glaciers were thinning during both epochs. The largest contributors to the SPI mass deficit in the southern sector during epoch 1 were Tyndall Glacier (VCR = −0.79 km³ a⁻¹) and Grey Glacier (VCR = −0.44 km³ a⁻¹). On both glaciers the volume loss rate decreased significantly during epoch 2, on Tyndall Glacier (VCR = −0.48 km³ a⁻¹) mainly due to decrease in losses above 700 m a.s.l. (Fig. 10) and on Grey Glacier (VCR = −0.07 km³ a⁻¹) at all elevations (Figs. 7, 10). Other glaciers with distinctly different hypsometric VCR between the two epochs are Penguin, Europa, Amalia and HPS 41 (Fig. S14).

Figure 8 shows a map of daily SECR on SPI during summer 2011/2012 based on DEM differencing spanning the periods 18 December 2011 to 26 March 2012 (99 d), 7 December 2011 to 15 March 2012 (99 d) and 29 December 2011 to 31 January 2012 (33 d). During summer the signal of surface lowering in the accumulation areas is mainly related to firn densification and melting of the top snow layers (Rivera, 2004). On the firn plateau, at elevations ≥1200 m, the average SECR in summer 2011/2012 was about −0.03 m d⁻¹ (blue curve in Fig. S3). In the ablation areas ice melt and dynamic downwasting (varying from glacier to glacier) are the main factors. High loss rates (SECR $<\sim -\mathrm{0.08}$ m d⁻¹) refer to areas that are subject to significant dynamic thinning, such as the lower terminus of Upsala and Viedma glaciers. Average summer melt rates for ice on the lower terminus of Perito Moreno Glacier (at 300 m altitude) are about 0.05 m d⁻¹ (Stuefer et al., 2007). On a glacier in a balanced state, surface lowering due to melt is in summer partly compensated for by uplift due to emergence.

The reported volume change and mass balance do not include subaqueous ice volume changes. Subaqueous losses are negligible with respect to the mass change of NPI since there are no large frontal retreats on water bodies. On SPI, the main glaciers and also many smaller ones terminate in proglacial lakes or in oceanic fjords. For the main retreating glaciers of SPI, Abdel Jaber (2016) estimated the subaqueous ice VCR at $-\mathrm{0.73}\pm \mathrm{0.22}$ km³ a⁻¹ for the period 2000 to 2011/2012. This number is obtained by measuring or estimating various parameters at the glacier front, including the glacier width, the water depth, the freeboard height on the two dates and the retreat distance. A bulk error of 30 % is assigned to the total subaqueous volume change rate, accounting also for unsurveyed glaciers. For the basal cross section at the calving front the shape of a semi-ellipse is assumed, except for four glaciers for which bathymetric data are available, enabling more precise estimates. For these glaciers a bulk error of 20 % is assumed for the subaqueous volume changes, amounting for the whole period to $-\mathrm{2.80}\pm \mathrm{0.56}$ km³ on the main front of Upsala Glacier, $-\mathrm{0.68}\pm \mathrm{0.14}$ km³ on Jorge Montt Glacier, $-\mathrm{0.59}\pm \mathrm{0.12}$ km³ on Tyndall Glacier and $-\mathrm{0.05}\pm \mathrm{0.01}$ km³ on Ameghino Glacier. The estimated subaqueous volume changes for the period 2000–2011/2012 are reported in Table S10, together with the frontal retreat distance.

5 Discussion

5.1 Spatial and temporal pattern of surface elevation and glacier volume change

Patagonian glaciers and ice fields experienced area retreat and shrinkage since the Little Ice Age which accelerated during recent decades associated with tropospheric warming (Davies and Glasser, 2012). Our estimate of mass loss for both ice fields during the period 2000 to 2016 is equivalent to 0.047±0.003 mm a⁻¹ eustatic sea level rise. This corresponds to 6 % of the ensemble mean contribution to sea level rise of glaciers and ice caps for the period 2005–2016 of 0.74±0.18 mm a⁻¹, based on global mass balance estimates from various sources (WCRP Global Sea Level Budget Group, 2018). Between epoch 1 (2000–2012) and epoch 2 (2012–2016) the rate of mass loss of SPI and NPI combined decreased by 9 % with a contrasting temporal trend between the two ice fields. The spatially detailed maps of SECR during the two subsequent epochs, derived from bistatic InSAR DEMs, provide a sound basis for studying the heterogeneous pattern of glacier response on NPI and SPI.

Regarding the ice bodies at large, on NPI the average loss rate increased by 31 % from epoch 1 to epoch 2 (VCR = $-\mathrm{4.26}\pm \mathrm{0.20}$ and $-\mathrm{5.60}\pm \mathrm{0.74}$ km³ a⁻¹, respectively). This was reversed on SPI, where the loss rate decreased by 20 % (VCR = $-\mathrm{14.87}\pm \mathrm{0.52}$ km³ a⁻¹ and $-\mathrm{11.86}\pm \mathrm{1.99}$ km³ a⁻¹, respectively). Reasons for the different behaviour are temporal changes of calving velocities, in particular on SPI, as well as a north–south gradient of air temperature increase in epoch 2 compared to epoch 1. Air temperatures, based on the European Centre for Medium-Range Weather Forecasts Interim Re-Analysis (ERA-Interim) (Dee et al., 2011; Berrisford et al., 2011) show for the ERA grid point 47.25^∘ S, 73.5^∘ W (NPI) at 850 hPa a mean annual temperature of +1.9 ^∘C during the period 2000–2011 and +2.3 ^∘C during 2012–2015. The corresponding values at the grid point 50.25^∘ S, 73.5^∘ W (southern SPI) are +0.7 and +0.8 ^∘C. The temperature difference between the two epochs was slightly larger during the main ablation period (1 November to 31 March): +4.1 ^∘C (summer 2000/2001 to 2011/2012) and +4.8 ^∘C (summer 2011/2012 to 2015/2016) on the NPI grid point and +2.4 and +2.7 ^∘C on southern SPI. The NCEP/NCAR Reanalysis 850 hPa mean temperature (Kalnay et al., 1996) is about 1 ^∘C lower but shows a similar temporal and spatial trend. Over an area extending from 48.00 to 51.75^∘ S, 72.75 to 74.25^∘ W, covering SPI, the mean annual precipitation, derived from ERA-Interim data, was slightly higher (8.4 %) in epoch 1 than in epoch 2 (Langhamer, 2017; Langhamer et al., 2018).

A main factor for the increased mass losses during epoch 2 on NPI is the higher air temperature compared to epoch 1, in particular during the main ablation period. Assuming a degree-day factor of 0.7 cm d⁻¹ on ice areas (Stuefer et al., 2007), the melt loss for an increase in surface temperature by 0.7 ^∘C during November to March corresponds to an additional loss of 0.74 m water equivalent. The hypsometric plot for the whole ice field shows changes of SECR by about −0.7 m a⁻¹ up to elevations of 1200 m a.s.l., indicating increased melt losses during epoch 2 not only on glacier termini but also on lower sections of the NPI plateau. At higher elevations the rate of surface lowering in both epochs, including the additional contribution during epoch 2, decreases gradually with elevation, reaching balanced state at about 1800 m in epoch 1 and about 2100 m in epoch 2 (Fig. 9). On NPI surface melt is the dominating process for mass depletion. During the period 2000 to 2009 the ice export due to calving amounted to about 20 % of the annual mass depletion by surface melt (Schaefer et al., 2013).

On lower sections of the main calving glaciers temporal variations in flow velocities are a main factor for the differences in SECR during the two epochs. In order to support the interpretation of differences in SECR between the two epochs, we derived maps of surface velocity gridded at 50 m for main glaciers from TerraSAR-X 11 d repeat pass data on various dates between 2010 and 2016, applying the offset tracking technique. The uncertainty of the velocity magnitude of these products is 0.05 m d⁻¹ (Wuite et al., 2015). Plots of velocities along central flow lines, extracted from the velocity maps, are shown in Fig. 11 for Jorge Montt, Pio XI, Upsala and Viedma glaciers.

Figure 11 Surface velocities along the central flow lines of Jorge Montt, Pio XI, Upsala and Viedma glaciers (SPI) on different dates, derived from TerraSAR-X repeat pass data. The distance origin refers to the ice front position on the first date: for Pio XI Glacier to the front of the southern branch and for NF to the front of the northern branch. Pio XI Glacier, dashed lines: velocities along the central flow line of the southern branch. Insets: TerraSAR-X images with location of central flow lines.

Flow velocities near the calving front of San Rafael Glacier reached magnitudes in excess of 18 m d⁻¹ in April 2007 (Willis et al., 2012a), dropping to 16 m d⁻¹ in May 2012 (Abdel Jaber et al., 2014). Mouginot and Rignot (2015) report a temporal peak in 2005 and a decrease by about 20 % until 2014 for the velocity at 10 km from the ice front. The drop in velocity between epoch 1 and 2 is reflected in the hypsometric curve of SECR, showing reduced loss rates below 800 m elevation during epoch 2 (Figs. 3 and 9). San Quintin Glacier, the largest glacier of NPI, reaches its maximum speed of about 3 m d⁻¹ 30 km from the front (Abdel Jaber et al., 2014; Mouginot and Rignot, 2015). Our analysis of TerraSAR-X data shows a mean increase in velocity by 10 % between May 2012 and June 2016 on the terminus. However, this caused only a minor additional increase in surface lowering (Fig. 9) because for this glacier the ice export due to calving accounts only for a very small part of total mass turnover (Schaefer et al., 2013).

On SPI, calving fluxes play a larger role for mass turnover than on NPI. This is reflected in the change in the average hypsometric curve of SECR of the ice field between the two epochs (Fig. 10). In spite of similar air temperatures during epoch 2 the average rate of surface lowering decreased at elevations below 400 m. Between 400 and 1000 m elevation the differences between the two epochs are very small. On the ice plateau, between 1000 and 2000 m, the loss rate decreased slightly, mainly brought about by minor changes on the southern sector of the ice field. Local increase in snow accumulation may play a role.

For six glacier basins, the VCR between the two epochs changed by more than +0.2 km³ a⁻¹, summing up to a combined decrease in volume losses by 2.20 km³ a⁻¹ (Table 3). The change of VCR from epoch 1 to epoch 2 amounted to +0.74 km³ a⁻¹ for Pio XI, to +0.37 km³ a⁻¹ for Grey and Dickson, to +0.33 km³ a⁻¹ for Upsala and Cono, to +0.30 km³ a⁻¹ for Tyndall, to +0.24 km³ a⁻¹ for Europa and to +0.22 km³ a⁻¹ for Penguin.

The behaviour of Pio XI Glacier, with frontal advance and positive mass balance has for many years been opposed to the general trend of SPI glaciers. The recent frontal advance trend started at the northern section of the terminus in 2006 and at the southern section in 2000 (Wilson et al., 2016). Between 2000 and 2014 a slowdown of velocity was observed on the southern section of the terminus (Mouginot and Rignot, 2015; Wilson et al., 2016). The slowdown went on until 2016, whereas the velocity of the northern section more than doubled between 2013 and 2016 (Fig. 11). Bathymetric data show shallow water with ridges running across the fjord at the present position of the ice front (Dowdeswell and Vásquez, 2013). This impedes calving at the southern ice front, causing during epoch 1 a main increase in surface elevation on the southern section, later on shifting towards the northern section that calves into Lago Greve (Fig. 6).

On Upsala Glacier the front retreated by 4 km between 2000 and 2014. The calving velocity reached a maximum in 2009/2010 (Abdel Jaber et al., 2012; Mouginot and Rignot, 2015) and decreased significantly afterwards, dropping from 8 m d⁻¹ in March 2011 to 5.9 m d⁻¹ in August 2014 and 4.8 m d⁻¹ in August 2016 (Fig. 11). This caused a major decrease in the thinning rate of the lower terminus during epoch 2 (Fig. 4).

The hypsometric curves of Grey and Tyndall glaciers show little change in SECR on the lower terminus close to the calving front and decreasing loss rates in the upper reaches of the terminus and in the accumulation area, an indication for surface mass balance as main cause for the change in SECR (Fig. 10). This is in line with TerraSAR-X surface velocity results between December 2011 and August 2016 showing only modest changes near the ice front and slowdown upstream. On Tyndall Glacier the velocity on the central flow line 0.5 km from the front was 0.96 m d⁻¹ in December 2011, 0.88 m d⁻¹ in October 2013 and 0.96 m d⁻¹ in August 2016. On Grey Glacier at the central flow line 3 km from the front, where the glacier splits into three branches, the velocity was 1.13 m d⁻¹ in December 2011, 1.11 m d⁻¹ in October 2013 and 1.02 m d⁻¹ in April 2016. Further upstream the velocity decreased by about 20 % between 2011 and 2016 on both glaciers. Weidemann et al. (2018) computed the surface mass balance of both glaciers and estimated the calving flux as residual of mean surface mass balance and geodetic mass balance over the period 2000 to 2014, pointing out that ice loss by surface ablation exceeds ice loss by calving. On Europa and Penguin, featuring steep narrow tongues, the SECR switched from slightly negative values to slightly positive values on the ice plateau above 1000 m elevation, also indicating a change in surface mass balance.

There are three glaciers with major increase in losses during epoch 2 (VCR becoming more negative by ≥0.2 km³ a⁻¹): the change of VCR for Jorge Montt is −0.36 km³ a⁻¹, for Viedma it is −0.28 km³ a⁻¹ and for Bernardo it is −0.20 km³ a⁻¹. Jorge Montt Glacier experienced a frontal retreat by 11 km between 1990 and 2011 (Rivera et al., 2012) and a further retreat by 2 km until 2016. The hypsometric profile shows high loss rates on the terminus at elevations up to 1000 m, with losses increasing during the second epoch (Table 3, Fig. 5) associated with major flow acceleration between 2007 and 2014 (Mouginot and Rignot, 2015). Figure 11 shows flow acceleration from 2010 to 2015/2016 extending 30 km up-glacier. On Viedma Glacier the increased mass loss during epoch 2 is caused by a major increase in the thinning rate on the glacier terminus below 1000 m elevation, an indication for changes in dynamic downwasting. This in accordance with a increase in ice velocity between 2010 and 2016 on the glacier terminus, notably in the lowest 5 km (Fig. 11). The calving velocity increased from 2 m d⁻¹ in August 2010 to 3.8 m d⁻¹ in August 2016.

The heterogeneous spatial pattern of elevation change on the two ice fields and its temporal evolution result from complex interdependencies between surface mass balance, responding directly to climate change signals, and effects of flow dynamics. Differences in surface elevation change and mass balance between individual glaciers and their temporal trends are particularly pronounced on SPI, where the calving fluxes represent a main component of mass turnover for most glaciers. The elevation dependence of the SEC reveals that ice dynamics exert a main control on topography not only on the glacier tongues but also on parts of the main ice plateau.

5.2 Comparison with previous estimates

A comparison of published results on volume change rates of SPI and NPI is reported in Table S11 for different epochs between 1968 and 2017, based on various methods including differencing of optical and/or interferometric DEMs and gravimetric time series of the GRACE mission. Similar comparisons are found in Malz et al. (2018) and Foresta et al. (2018). Our results are in line with geodetic mass balance results of NPI and SPI published by other authors, which suggest an overestimation of the mass losses retrieved from gravimetric time series such as those found in Chen et al. (2007), Ivins et al. (2011) and Jacob et al. (2012), the latter referring to Patagonia in general.

Our result for NPI during epoch 1 (VCR = $-\mathrm{4.26}\pm \mathrm{0.20}$ km³ a⁻¹) complies with the numbers reported by Abdel Jaber (2016) for the period 2000 to 2014 ($-\mathrm{4.40}\pm \mathrm{0.13}$ km³ a⁻¹) and Willis et al. (2012a) for 2000 to 2011 ($-\mathrm{4.06}\pm \mathrm{0.12}$ km³ a⁻¹), the latter based on SRTM and ASTER DEMs. Willis et al. (2012b) recomputed their previous estimate applying a 2 m offset to the SRTM DEM to account for signal penetration, which results in larger losses (VCR = $-\mathrm{4.9}\pm \mathrm{0.3}$ km³ a⁻¹). This correction is not comprehensible given the wet status of the snow surface during the summer acquisition of SRTM, as evident from the backscatter data (Sect. 3.2.2 and Abdel Jaber, 2016). Braun et al. (2019) report a VCR $-\mathrm{4.65}\pm \mathrm{0.17}$ km³ a⁻¹ for NPI over the period 2000 (SRTM data) to 2011/2015 (TDM data).

Our VCR for epoch 1 is slightly lower than the results of Dussaillant et al. (2018), who applied two methods: differencing of SPOT and SRTM DEMs (VCR = $-\mathrm{4.55}\pm \mathrm{0.41}$ km³ a⁻¹) and derivation of temporal elevation trends from ASTER DEM time series (VCR = $-\mathrm{4.72}\pm \mathrm{0.34}$ km³ a⁻¹). Our hypsometric curve of SEC shows a similar behaviour as their ASTER_trend results up to 2800 m of elevation, although with slightly lower losses at most elevations. Above 1000 m Dussaillant et al. (2018) report 35 % and 22 % of unsurveyed area for the SPOT-SRTM analysis and ASTER_trend, respectively, mostly due to the lack of contrast or the presence of clouds in the optical stereo images. For the same elevation band the unsurveyed area in our 2000–2012 SECR map of NPI is 6 %. On glaciers larger than 100 km² the SEC rates with both methods applied by Dussaillant et al. (2018) agree with our results within error bars. On two medium-sized glaciers, Exploradores (86 km²) and Grosse (67 km²), the average SECR of their two methods differs by more than 1.0 m a⁻¹, their ASTER_trend being ∼0.8 m a⁻¹ higher than our SECR and ∼0.6 m a⁻¹ higher than those of Willis et al. (2012a).

On SPI Willis et al. (2012b) estimate a VCR of $-\mathrm{21.2}\pm \mathrm{0.5}$ km³ a⁻¹ for the period 2000–2011, a much larger value compared to the value reported here for epoch 1 (VCR = $-\mathrm{14.87}\pm \mathrm{0.52}$ km³ a⁻¹) and to that of Abdel Jaber (2016) for 2000–2011/2012 (VCR = $-\mathrm{14.59}\pm \mathrm{0.37}$ km³ a⁻¹). The discrepancy largely exceeds the 10 % VCR contribution they attribute to the 2 m correction for signal penetration in the SRTM DEM.

Malz et al. (2018) present SECR maps and mass balance of SPI for the period 2000–2015 based on SRTM and several TDM DEMs of December 2015. We used the same raw data at the end of our epoch 2. They do not account for missing summer days and report a VCR of $-\mathrm{13.2}\pm \mathrm{3.6}$ km³ a⁻¹. Scaling our VCR results over the two epochs and accounting for the missing summer days in order to cover a period of 16 years, from mid-February 2000 to mid-February 2016, we obtain $-\mathrm{14.1}\pm \mathrm{0.9}$ km³ a⁻¹. The difference can probably be explained by the missing 48 to 76 summer days required for spanning a full period of 16 years. Applying the method of Sect. 3.1.3 for the missing summer days, we obtain an ice-field-wide average SECR value of −0.12 m a⁻¹, corresponding to a VCR of −1.5 km³ a⁻¹, which is not taken into account by Malz et al. (2018). For the southern sector of SPI, Malz et al. (2018) show SECR maps and hypsometric curves for the periods 2000–2012 and 2012–2015, based on the same TDM raw data used in this study (scenes 7 and 8 and scenes 13 and 14). The absence of a correction for 53 or 59 summer days at the end of the 4-year period leads to lower loss rates compared to our numbers for epoch 2.

Average SEC rates for single glaciers are reported by Willis et al. (2012b) and Malz et al. (2018). On several main glaciers, including Bernardo, Tempano, Occidental, Greve, Chico, Europa and Guilardi glaciers, a direct comparison is not possible because of different glacier outlines. Among main glaciers with similar areas our SECR estimates are in general lower than those of Willis et al. (2012b). Among glaciers >200 km², average SEC rates deviating by more than −1.0 m a⁻¹ from our results are reported for Tyndall, Pio XI and Perito Moreno glaciers.

Foresta et al. (2018) compute the geodetic mass balance of NPI and SPI for six glaciological years between 2011 and 2017 from SEC maps using swath-processed CryoSat-2 (CS2) interferometric data with sub-kilometre spatial resolution. The acquisitions dates vary spatially for different pixels. The authors explain that seasonality biases are avoided due to the regular flight path of CS2 ensuring data acquisition within each pixel at the same epochs in each glaciological year. The data coverage is relatively poor (46 % for NPI, 50 % for SPI), in particular on lower sections of glacier tongues. Termini of several main glaciers are not covered at all and the SECR data appear to be relatively noisy. To fill data gaps hypsometric average models are applied, using the values of polynomials (degree 1 to 3) fitted to the observed hypsometric SECR. This is performed for NPI at large and for SPI separately for six of the main glaciers and two large sub-regions, each of which is comprised of many glaciers with different sizes and physiographic features. Therefore, the resulting maps of SECR do not reflect the complexity and spatial variability of SECR derived from high-resolution geodetic data, particularly at lower elevations.

A comparison between VCR estimates of Foresta et al. (2018) and our results for 2012–2016 is provided in Table S12. The two datasets do not cover exactly the same period. The ERA-Interim 850 hPa data over NPI and SPI show an agreement within 0.1 ^∘C compared to 2012–2016 for the average air temperature of the main ablation period (November to March) 2011 to 2017, suggesting similar rates of surface melt during the two epochs. Foresta et al. (2018) report a mass balance of $-\mathrm{6.79}\pm \mathrm{1.16}$ Gt a⁻¹ for NPI and $-\mathrm{14.50}\pm \mathrm{1.60}$ Gt a⁻¹ for SPI using 900±125 kg m⁻³ for volume-to-mass conversion. This corresponds to a VCR of $-\mathrm{7.54}\pm \mathrm{0.75}$ and $-\mathrm{16.11}\pm \mathrm{1.43}$ km³ a⁻¹, respectively. These numbers for volume loss are significantly higher than our results for epoch 2 in all sub-regions, being overall 35 % higher for NPI and 36 % for SPI. Foresta et al. (2018) also show time series of cumulative mean observed elevation change for the nine sub-regions. For NPI and five sub-regions of SPI the plots show minima of the annual elevation during mid-winter in some years. This is not compatible with both the annual cycle of surface mass balance and the seasonal variation in flow velocities on glacier tongues showing a trend for higher velocities in summer compared to winter (Stuefer et al., 2007; Minowa et al., 2017).

6 Conclusions

We report on a detailed study focussing on the climate-sensitive northern and southern Patagonian ice fields, where high-resolution maps of surface elevation change were obtained for the epochs 2000–2012 and 2012–2016 from bistatic InSAR DEMs, allowing us to derive the total net mass balance of most of the glacier basins. We rely on a reprocessed version of the SRTM C-band DEM featuring improved absolute height calibration and on a series of TanDEM-X Raw DEMs, processed with a robust phase unwrapping method, leading to almost complete coverage, including narrow glaciers and high altitudes. Significant effort was dedicated to reduce systematic errors, especially critical for the mass balance of vast regions: a precise co-registration of the DEMs was performed, seasonal biases due to gaps in the full annual cycle were corrected based on a complementary TDM summer SEC map, the backscatter coefficient of all acquisitions (including SRTM) was analysed to assess signal penetration. A comprehensive uncertainty estimation including all main error sources of the SEC maps and of the mass balance was also performed.

A similar pattern of elevation change is found on NPI for the two epochs, with lowering on most of the termini and well into the main ice plateau with increased loss rates during epoch 2. Being mass depletion mainly driven by surface melt on NPI, this trend is likely due to higher average air temperatures during epoch 2. The estimated volume change rate increased by 31 % from $-\mathrm{4.26}\pm \mathrm{0.20}$ km³ a⁻¹ in epoch 1 to $-\mathrm{5.60}\pm \mathrm{0.74}$ km³ a⁻¹ in epoch 2.

On SPI the spatial pattern and the temporal trend of SECR are more complex. For 63 % of the glaciers larger than 35 km² in area the loss rate decreased or turned into positive values from epoch 1 to epoch 2, whereas for 37 % the loss rates increased. The overall volume change rate decreased by 20 % from $-\mathrm{14.87}\pm \mathrm{0.52}$ km³ a⁻¹ during epoch 1 to $-\mathrm{11.86}\pm \mathrm{1.99}$ km³ a⁻¹ during epoch 2. The only glacier with positive mass balance in both epochs is Pio XI Glacier, with an increase in the volume change rate from 0.52±0.04 to 1.26±0.25 km³ a⁻¹. The more complex behaviour of SPI glaciers is caused by the higher relevance of calving fluxes as a source of mass turnover in which the effect of ice dynamics on surface elevation changes extends to the main ice plateau. Different temporal trends between individual glaciers are commonly associated with opposing trends in calving velocities. On the accumulation areas south of ∼49.5^∘ S the SECR was either stable or slightly negative during epoch 1, turning to slightly positive in epoch 2. Air temperature remained relatively stable in the south of SPI, meaning a north–south gradient was present. This, coupled with a possible local increase in snow accumulation may be the cause of the decreased loss rates at elevations above 1000 m. Significant frontal retreat was observed on SPI during epoch 1, Abdel Jaber (2016) reported a coarse estimation of subaqueous volume change of $-\mathrm{0.73}\pm \mathrm{0.22}$ km³ a⁻¹ for the period 2000–2011/2012.

The eustatic sea level rise contribution of both ice fields, excluding subaqueous changes, was estimated to be 0.048±0.002 mm a⁻¹ in epoch 1 and 0.043±0.005 mm a⁻¹ in epoch 2. Behind these numbers lies a complex interplay between surface mass balance, responding directly to climate change and ice flow dynamics, mechanisms which regulate the heterogeneous spatial pattern and temporal evolution of the SEC on NPI and SPI.

This study confirms the potential of bistatic InSAR and particularly of the TanDEM-X mission for accurate, detailed and almost gapless mapping of surface elevation changes of large ice fields even for small basins and tongues. We recommend the use of TanDEM-X data – with an appropriate co-registration and care for radar signal penetration – to map SEC of all types of glaciers, as also recently shown in the northern Antarctic Peninsula (Rott et al., 2018). We hope that our results will encourage the development of remote-sensing missions capable of repeated bistatic InSAR observations allowing regular worldwide SEC mapping and mass balance estimations with improved temporal sampling.