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

The Northern and Southern Patagonian icefields (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 SAR interferometry data of SRTM and 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 icefields at large. We apply advanced TanDEM-X processing techniques allowing to cover 90 % and 95 % of the area 5 of NPI and 97 % and 98 % of SPI for the two epochs, respectively. Particular attention is paid to precisely coregistering the 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 icefields and reveal a heterogeneous spatial pattern of SEC and mass balance caused by different sensitivities in respect to direct climatic forcing and ice flow dynamics of individual glaciers. The estimated volume change rates for NPI 10 are −4.26± 0.20 km3 a−1 for epoch 1 and −5.60± 0.74 km3 a−1 for epoch 2, while for SPI these are −14.87± 0.52 km3 a−1 for epoch 1 and −11.86±1.99 km3 a−1 for epoch 2. This corresponds for both icefields to an eustatic sea level rise of 0.048±0.002 mma−1 for epoch 1 and 0.043± 0.005 mma−1 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 of flow velocities are a main factor for differences in SEC between the two epochs. Striking differences are observed 15 even on adjoining glaciers, such as Upsala Glacier with decreasing mass losses associated with slowdown of flow velocity, contrasting with acceleration and increase of mass losses on Viedma Glacier.

. Because the icefields 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 5 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 icefields are tidewater or freshwater calving glaciers, showing heterogeneous patterns 10 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 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 15 with considerable variability in magnitude and timing for individual glaciers. Only few glaciers advanced intermittently during recent decades. The most striking case is the Pio XI Glacier showing a large cumulative frontal advance since 1945, including recently a general advancing period of its southern and northern branches starting in , respectively (Wilson et al., 2016. Geodetic mass balance estimates of NPI and SPI have been derived from various sources. The first remote sensing based 20 estimates of NPI and SPI volume change and mass balance were reported by Rignot et al. (2003) March 2012, on the other hand by fitting pixel-based linear elevation trends over 118 DEMs calculated from ASTER stereo images acquired between 2000 and 2012. Icefield-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 icefields and estimated the mass balance for six years between April 2011 and March 2017. The maps 5 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 (shortly TDM) (Krieger et al., 2007), composed of the two formation-flying radar satellites TerraSAR-X and TanDEM-X, is operational since December 2010. The The paper presents the first spatially detailed analysis of surface elevation change and the derived total net mass balance over the Patagonian icefields 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 2 of NPI and > 9 km 2 of SPI. The results indicate a different temporal trend between the two icefields 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 10 interferometers.

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, 15 2018), whose performances are analyzed in Rizzoli et al. (2017) and Wessel et al. (2018), is the result of the combination of four 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 DLR-20 CAF, 2010) using ITP, the operational TanDEM-X processor (Breit et al., 2012;Fritz et al., 2011), in order to generate two elevation maps completely covering the icefields 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, 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 25 should be acquired at the same time of year in order to minimize seasonal changes, which can be significant on the Patagonian icefields. Since data availability restricted fulfilling the last criterium, the residual temporal gap had to be compensated.
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 the western termini of NPI where, due 30 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 5 coherence and flag mask indicating critical areas.

SRTM
The SRTM (Farr et al., 2007;Rabus et al., 2003) was launched 11 February 2000 and produced in 9 days 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 applying the ScanSAR technique to four sub-swaths featuring different polarization (HH, 10 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 9 ascending and 9 descending datatakes covering the Patagonian icefields (Seal and Rogez, 2000) are listed in Table S3. The performance of SRTM was assessed among others by Rodriguez et al. (2005); Brown et al. (2005); Carabajal and Harding (2006);Wendleder et al. (2016). The main issue is the presence of long-wavelength height errors with magnitude up to ∼ 20 m globally and spatial variation scales of hundreds to 15 thousands of kilometres, mainly caused by residual roll errors due to the attitude adjustment manoeuvres of the 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 20 NASADEM (NASA JPL, 2018) as the elevation map of year 2000 for both icefields. The choice was done after comparing on a vast region surrounding the Patagonian icefields the NASADEM and the SRTM ver. 3 (SRTMGL1) (NASA JPL, 2013) to the TDM global DEM rescaled to 1 arcsec. The SRTMGL1 data set, 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 geo-location biases (∆h images are shown in Fig. S2) and a higher RMS when compared to the NASADEM. On the icefields the differences 25 between the two SRTM data sets are larger on NPI and in the very south of SPI.
We furthermore retrieved the SRTM radar brightness images (SRTMIMGR) (NASA JPL, 2014) for the sub-swaths covering the icefields (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.

Glacier outlines 30
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). The use of ITP to process the single 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 coregistration (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 . We 10 hence relied on 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 (maximum half of the HoA). The absolute phase of the data is then reconstructed by 15 summing to the unwrapped differential phase the phase simulated from the reference DEM, this way 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 depends. ITP allows to manually update the APO value and perform a new geocoding, for instance to fine-tune the coregistration with a reference DEM, as described in Sect. 3.1.2.

DEM coregistration
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, 25 coregistered to each other and to the SRTM DEM we coregistered 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 30 geocoding with ITP. This method assures high precision by exploiting the geometrical parameters of the SAR acquisition and allows avoiding critical aspects of the generic coregistration 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 icefields 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 5 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: 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 ( Values of δh off ranged in magnitude between 0 m and 1.8 m for the TDM Raw DEMs and was used to fine-tune the APO. 15 For the NASADEM it was equal to 0.3 m and 0.1 m in absolute value on NPI and SPI, respectively, and was subtracted after checking the horizontal coregistration with respect to the TDM reference DEM at 0.4 arcsec in the proximity of the icefields. Furthermore range and azimuth tilts caused by baseline errors (Hueso González et al., 2010) were verified and found to be negligible for all Raw DEMs. Height consistency between overlapping Raw DEMs was also checked in order to ensure a seamless elevation map. 20 The coregistration procedure partly compensates the crustal uplift rates due to the glacial isostatic adjustment affecting the region, characterized by rates up to 40 mm a −1 on the plateau of SPI and decreasing with distance as reported by Dietrich et al. (2010).

Seasonal correction
Seasonal variations of SECR should be taken into account for deriving annual rates of surface elevation change if the time span 25 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), respectively 12 years (2000 to 2012), taking seasonal differences in SECR for missing days 30 into account versus mean annual SECR without accounting for such differences.  Fig. S1a). For the main sections of SPI the percentage of missing days ranges from 3.6 % to 5 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 applied seasonal corrections also to this data set. 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 on the seasonal correction for the epoch 2012-2016 because of their larger impact. Depending on 10 the availability of additional TDM DEM data, the following procedures were applied for different sections of the icefields: -Three additional TDM acquisitions of December 2011 (Table S2)   to 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 Moreno Glacier (Stuefer et al., 2007). The applicability of the 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 icefield (Bravo et al., 2019). Further details are given in the Supplement, Sect. S4.  (Stuefer et al., 2007).
Finally the correction rasters were obtained by scaling pixelwise the daily correction rate by the temporal gaps in days ( Fig.   S4). To avoid biases of the mass balance we masked-out artefacts due to phase unwrapping, layover, shadow, etc. The elevation of the water surface subject to frontal retreat, usually decorrelated, was manually edited in order to correctly capture the freeboard SEC.

Derivation of SECR maps and estimation of mass balance
By multiplying the average SECR with the corresponding glacier area over elevation intervals of 50 m the altitude-dependent 25 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 on the entire icefields 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 −3 for the conversion of the VCR to mass change rate. This value is commonly 30 used for geodetic mass balance measurements and provides traceability for comparisons with other studies (Cogley, 2009).
The main mass losses on the Patagonian icefields refer to ice areas, and for the accumulation areas assumptions on changes of the vertical profiles of snow/ice density would be speculative.

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 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 5 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 C-and X-band the penetration depth drops rapidly below 0.2 m already with a LWC of approximately 0.5 % vol . We used the backscattering coefficient σ 0 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. 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 6:00 local time (UTC−4h), which is generally the coldest time of the day, although the plateaus of NPI and SPI usually feature limited daily variations of air temperature due to 20 the dense clouds and strong precipitation occurring most of the year (Garreaud et al., 2013;Schaefer et al., 2013).
The backscattering of SPI is more heterogeneous compared to NPI. The 2015 coverage features σ 0 ≤ −19 dB revealing wet snow on large parts of the plateau (particularly on the western margin). The σ 0 of the 2012 coverage is in average higher (especially on scene 4/5 acquired at the end of March). Over the main parts of the plateau σ 0 is still lower than −16 dB, an 30 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 σ 0 < −18 dB imputable to wet snow on most of the plateau. Some isolated regions with higher σ 0 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 which we considered prone to signal penetration in each DEM mosaic (Fig. S5). The outlining was performed manually in order to identify areas on rough ice surfaces and on radar fore-slopes, where high σ 0 is not an indicator for signal penetration. We 5 assigned a potential penetration height offset to each of these polygons according to its average σ 0 and look angle. This offset is taken into account in the error budget. The offsets are based on empirical observations of the relationship between σ 0 and height offset performed on multiseasonal TDM Raw DEMs of NPI, showing a mean penetration bias of 4 m for an increase of σ 0 by 10 dB from wet to dry snow (Abdel Jaber, 2016).

Assessment of SRTM backscatter 10
The SRTM absolute and relative radiometric accuracy nominal values are 3 dB and 1 dB, respectively (Farr et al., 2007). The SRTMIMGR product provides the radar brightness β 0 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 midlook angle of each sub-swath, introducing this way an error up to ±0.6 dB and computed the backscattering coefficient using the provided local incidence angle (θ loc ) mask as σ 0 = β 0 · sin θ loc .
15 Figure S6 shows the arithmetic mean (σ 0 ) and the standard deviation computed pixelwise from the sub-swath σ 0 images covering the icefields (4 to 7 stacked pixels are usually found). The measure of spread supports the interpretation ofσ 0 . While  (Table S3).
Values ofσ 0 < −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. S6). Values ofσ 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 25 snow with a rough surface (Nagler and Rott, 2000). Above 2300 mσ 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σ 0 , it can be concluded that the SRTM elevations are not affected by a bias due to C-band radar  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.

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. The random error of each SECR sample σ SECR was computed pixelwise 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: For TDM elevations the random error is given in the HEM raster, which expresses the interferometric standard error for each 10 sample computed assuming a normally distributed error as (Rossi et al., 2010): 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 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:

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 coregistration to the reference DEM (Sect. 3.1.2) was defined for each Raw DEM and for the 2. To account for signal penetration we used the penetration height offsets assigned to critical regions on each DEM mosaic 15 as local systematic errors, ranging between 1 m and 6 m according to σ 0 , look angle and radar frequency (Sect. 3.2).
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 20 GIA effects, residual tilts, unmasked local errors due to PU or layover, etc. components (ε reg , ε pen and ε add ) were estimated separately for the summer 2011/2012 SECR. 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 pixelwise scaling by the number of corrected days and by the appropriate ∆t in years was applied, leading to a fourth systematic error The total systematic error ε (x, y) of each SECR sample was obtained pixelwise as (omitting (x, y) for compactness): ε (x, y) = ε 2 reg + ε 2 pen + ε 2 add + ε 2 seas .
The mean values of ε (x, y) and of its components for the four SECR maps are reported in Table S7.

Geodetic mass balance error
The geodetic method was applied to estimate the average SECR on separate elevation bins and the corresponding volume and 5 mass change rates. The total error of the mean SECR on elevation bin b was computed by summing in quadrature the mean systematic error ε b on bin b and the standard error SE b of the spatial average on bin b (which is generally negligible compared to ε b ), obtained as: where N b is computed according to Eq. (5). 10 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 icefield, 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.

Results
The SECR maps of NPI and SPI after seasonal correction are shown in Fig. 1 Table   25 1 the SECR, the volume change rate (VCR), the mass balance and the contribution to sea level rise are specified for the entire icefields. Table 2 provides SECR and VCR for NPI glaciers larger than 2 km 2 , Table 3 for SPI glaciers larger than 35 km 2 and Table S8 in the Supplement for SPI glaciers with area between 35 km 2 and 9 km 2 . The tables report also 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. S10 and S11 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. S12 and S13. 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. S9 together with the 2012 TDM DEM average surface elevation.

5
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 icefield 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 2 show volume losses during both epochs except Leones Glacier revealing a modest increase in ice volume ( Table 2). The volume loss rate and −0.87 km 3 a −1 ), but the loss pattern changed (Fig. 3, Fig. 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 The only glacier with positive mass balance in both epochs is Pio XI Glacier, showing a significant increase of VCR from 0.52 km 3 a −1 in epoch 1 to 1.26 km 3 a −1 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 m 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 m and 1500 m on the ice plateau (Fig. 10).

5
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 2 and the majority of smaller glaciers show negative SECR during epoch 1 (Table 3 and S13). On the main ice plateau the surface elevation was either stable or the 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 30 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 is 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 −2.80 ± 0.56 km 3 on the main front of Upsala Glacier, −0.68 ± 0.14 km 3 on Jorge Montt Glacier, −0.59 ± 0.12 km 3 on Tyndall Glacier and −0.05 ± 0.01 km 3 on Ameghino Glacier. The estimated subaqueous volume changes for the period 2000-2011/2012 are reported in Table   S9, together with the frontal retreat distance.

Spatial and temporal pattern of surface elevation and glacier volume change
Patagonian glaciers and icefields experienced area retreat and shrinkage since the Little Ice Age which accelerated during  . 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 of 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-day repeat pass data on various dates 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 of the average hypsometric curve of SECR of the icefield between the two epochs (Fig. 10). For six glacier basins the VCR between the two epochs changed by more than +0.2 km 3 a −1 , summing up to a combined decrease of volume losses by 2.20 km 3 a −1 ( 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 of 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(Abdel Jaber et al., 2012Mouginot and Rignot, 2015) and decreased significantly afterwards, dropping from 8 m d −1 in March 2011 to 5.9 m d −1 in August 2014 and 4.8 m d −1 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 5 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  dependence of the SEC reveals that ice dynamics exerts a main control on topography not only on the glacier tongues, but also on parts of the main ice plateau.

Comparison with previous estimates
A comparison of published results on volume change rates of SPI and NPI is reported in Table S10 for different epochs between   1968  days and report a VCR of −13.2 ± 3.6 km 3 a −1 . 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 −14.1 ± 0.9 km 3 a −1 . 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 icefield-wide average SECR value of −0.12 m a −1 , corresponding to a VCR of −1.5 km 3 a −1 , 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/8 and 13/14). The absence of a correction for 53/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, for higher velocities in summer compared to winter (Stuefer et al., 2007;Minowa et al., 2017).

Conclusions
We report on a detailed study focussing on the climate-sensitive Northern and Southern Patagonian icefields, where high resolution maps of surface elevation change were obtained for the epochs 2000-2012 and 2012-2016 from bistatic InSAR DEMs allowing to derive the total net mass balance of most of the glacier basins. We rely on a re-processed 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 coregistration of the DEMs was performed, seasonal biases due to gaps in the full annual cycle were corrected based on a 5 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 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 icefields even for small basins and tongues. We recommend the use of TanDEM-X data-with an appropriate coregistration and care for radar signal penetration-to map SEC of all types of glaciers, as recently shown also in the northern Antarctic Peninsula (Rott et al., 2018). We hope that our results will encourage 30 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.