Remote Sensing of Environment Heterogeneous and rapid ice loss over the Patagonian Ice Fields revealed by CryoSat-2 swath radar altimetry

The Northern and Southern Patagonian Ice Fields (NPI and SPI) in South America are the largest bodies of ice in the Southern hemisphere outside of Antarctica and the largest contributors to eustatic sea level rise (SLR) in the world, per unit area. Here we exploit swath processed CryoSat-2 interferometric data to produce maps of surface elevation change at sub-kilometer spatial resolution over the Ice Fields for six glaciological years between April 2011 and March 2017. Mass balance is calculated independently for nine sub-regions, including six individual glaciers larger than 300km 2 . Overall, between 2011 and 2017 the Patagonian Ice Fields have lost mass at a combined rate of 21.29 ± 1.98Gta − 1 , contributing 0.059 ± 0.005mma − 1 to SLR. We observe widespread thinning on the Ice Fields, particularly north of 49° S. However the pattern of surface elevation change is highly heterogeneous, partly re ﬂ ecting the importance of dynamic processes on the Ice Fields. The Jorge Montt glacier (SPI), whose tidewater terminus is approaching ﬂ oatation, retreated ~2.5km during our study period and lost mass at the rate of 2.20 ± 0.38Gta − 1 (4.64 ± 0.80m we a − 1 ). In contrast with the general pattern of retreat and mass loss, Pio XI, the largest glacier in South America, is advancing and gaining mass at 0.67 ± 0.29Gta − 1 rate.


Introduction
The Northern and Southern Patagonian Ice Fields are the two largest ice bodies in the Southern Hemisphere excluding Antarctica, with areas of about 4200 and 13,000 km 2 and volumes of about 1200 and 4300 km 3 (Carrivick et al., 2016), respectively, and elevation ranging from sea level to about 3900 m. They experience relatively warm and wet climatic conditions (Sagredo and Lowell, 2012) and lie on top of the narrow Andean mountain range, which forms an efficient barrier to the predominantly westerly winds and moisture rich air transported inland from the Pacific Ocean (e.g. Garreaud et al., 2013). The mountain and associated ice divide separates areas with contrasting climatic conditions. On the western side, orographic uplift of moist air produces extreme annual precipitation of up to 10 m a −1 (e.g. Lenaerts et al., 2014;Sakakibara and Sugiyama, 2014 and references within) as well as extensive and persistent cloud coverage with associated low shortwave and high longwave energy fluxes (Lenaerts et al., 2014). To the east of the divide, the Ice Fields are in the rain shadow and receive comparatively high amounts of shortwave energy (Lenaerts et al., 2014). Thus, the glaciers west and east of the ice divide are thought to be more sensitive to changes in precipitation and air temperature respectively (Rivera and Casassa, 1999;Warren and Sugden, 1993).
The two Ice Fields have been experiencing long-term thinning and retreat. Between the end of the Little Ice Age (LIA,~1870) and 2011, their area shrank by~12.5% on average (Davies and Glasser, 2012), associated with a combined mass loss of about 1.70 ± 0.25 Gt a −1 (Glasser et al., 2011). However, geodetic studies based on gravimetry (Chen et al., 2007;Ivins et al., 2011;Jacob et al., 2012;Gardner et al., 2013) and comparisons of Digital Elevation Models (DEM; Rignot et al., 2003;Willis et al., 2012aWillis et al., , 2012bJaber, 2016) estimated that between 1975 and 2012 the rate of mass loss of the Ice Fields has been in the range of 15 to 35 Gt a −1 , one order of magnitude more compared to the long term trend.
During the last 50 years, the Patagonian Ice Fields contributed an estimated 10% to the total mass loss from glaciers and ice caps, excluding those at the periphery of the Greenland and Antarctic ice sheets (Glasser et al., 2011 and references within), increasing to 15.4% in the first decade of the 21st century (Jacob et al., 2012;Gardner et al., 2013), second only to glaciers in Alaska and the Canadian Arctic, and larger than high mountain Asia (Brun et al., 2017) which all extend over areas~5-8 times larger. Currently, the Patagonian Ice Fields are the largest contributor to sea level rise per unit area in the world (Gardner et al., 2013;Carrivick et al., 2016).
Velocities of glaciers draining the Ice Fields (up to 10 km a −1 ) are amongst the fastest in the world (Sakakibara and Sugiyama, 2014;Mouginot and Rignot, 2015) and substantial ice flow acceleration has been observed, coincident with rapid frontal retreat, for a number of tidewater and lacustrine glaciers (Sakakibara and Sugiyama, 2014). These observations implicate the importance of the role of dynamics and tidewater glacier calving in the rapid wastage of the Ice Fields. In fact, > 80% of them terminate in proglacial lakes (mostly across the NPI and east of the SPI) or fjords (western side of the SPI) (Sugiyama et al., 2016 and references within).
Since 2010, the European Space Agency (ESA) radar altimetry mission CryoSat-2 (CS2) (Drinkwater et al., 2005;Wingham et al., 2006) has been acquiring topography data over land ice. Radar instruments are particularly suited to land ice applications since they can penetrate through clouds and do not depend on sunlight. Radar altimetry data have been previously exploited to map elevation change over ice caps (Rinne et al., 2011a(Rinne et al., , 2011b, but the technique has not been applied widely due to the limitation caused by the large radar footprint. CS2's state-of-the-art radar altimeter uses Synthetic Aperture Radar (SAR) along-track to reduce the footprint size as well as interferometry across-track to accurately locate the position of the surface reflection (SARIn mode; Wingham et al., 2006). Additionally, its orbit inclination of 92°and repeat cycle of 369 days provides an inter-track spacing of~5 km on average over the Patagonian Ice Fields. Finally, CS2's relatively short wavelength (2.2 cm; Ku band) restricts the penetration of the radar pulse in the snowpack, compared to, e.g., sensors working in C or X bands. These characteristics make CS2 better suited for monitoring changes in glacier areas with frequent cloud cover and considerable slopes. CS2 SARIn data have successfully mapped topographic changes over Arctic ice caps (McMillan et al., 2014a;Gray et al., 2015). Additionally, swath processing (Hawley et al., 2009) of CS2 SARIn data has been applied to generate high resolution DEMs of ice caps and selected areas of the Greenland and Antarctic ice sheets (Gray et al., 2013;Ignéczi et al., 2016;Gourmelen et al., 2017a) and to produce sub-kilometer maps of surface elevation change (Christie et al., 2016;Foresta et al., 2016;Gourmelen et al., 2017a;Gourmelen et al., 2017b), with a wide range of applications such as the identification of supraglacial lakes in NE Greenland (Ignéczi et al., 2016), subglacial lakes in West Antarctica and regions of subsidence in Iceland (Smith et al., 2017;Gourmelen et al., 2017a) as well as quantifying channelized basal melt under the Dotson ice shelf in West Antarctica (Gourmelen et al., 2017b) and volume and mass change of individual ice caps in Iceland (Foresta et al., 2016).
Despite their important contribution to ice mass loss and global SLR, studies quantifying mass changes of the Patagonian Ice Fields are limited in number and do not cover the most recent period. This paper focuses on quantifying the mass balance of the NPI and SPI during six glaciological years between April 2011 and March 2017. To this aim, we exploit swath processed CS2 SARin data to generate maps of surface elevation change rates at sub-kilometer spatial resolution, and convert them into estimates of glacier volume and mass change. For a number of large catchments on the SPI, such estimates are derived at the basin scale. Additionally, the dense L2swath elevation field enables the production of time series of elevation change for different sub-regions of the Ice Fields exhibiting contrasting patterns of change.

Data and methods
We exploit swath processed CS2 SARIn baseline C data (L2swath) as this technique (Hawley et al., 2009;Gray et al., 2013;Foresta et al., 2016;Gourmelen et al., 2017a) can generate up to two orders of magnitude more data than conventional Point-Of-Closest-Approach (POCA) processing and, importantly, provides more homogeneous spatial coverage over relatively small glaciated regions with considerable topography (Foresta et al., 2016;Gourmelen et al., 2017a). The L2swath processing scheme is similar to Foresta et al. (2016) and Gourmelen et al. (2017a), but we use a different procedure to discard noisy waveform samples before performing the phase unwrapping. This procedure, first developed for InSAR images (Weissgerber, 2016) and updated for CS2 SARIn data (Weissgerber et al., 2018), was shown to further increase the density of the L2swath elevation field and to improve the spatial coverage of the Jakobshavn glacier, Greenland (Weissgerber and Gourmelen, 2017) (Appendix A). The L2swath algorithm makes use of an external DEM to improve the precision of elevation measurements in the presence of slopes larger than 0.54°, where an entire waveform may be affected by a phase shift. Without this correction, observations may be wrong by several tens of meters in elevation and a few kilometres in geo-location (Gourmelen et al., 2017a). It is not straightforward to predict the accuracy needed for the DEM (Gourmelen et al., 2017a). However given the magnitude of the geolocation and elevation shifts, the DEM need not be extremely accurate. We used the SRTM C band DEM (Farr et al., 2007) acquired in 2000 as a reference for elevation, after including a correction to account for the elevation change occurred between 2000 and 2011 (Appendix B). L2swath data are then used to compute rates of surface elevation change for six glaciological years between April 2011 and April 2017 using a plane-fit algorithm (e.g. McMillan et al., 2014b).
One glaciological year is defined as the period between 1st April in year n and 31st March in year n + 1. CS2's acquisition dates vary spatially for different pixels due to the satellite's orbital path as well as to the local topography, so that the temporal resolution at the pixel scale is non-uniform and longer than monthly. However, seasonality biases are avoided given the regular flight path followed by CS2, which ensures that data within each pixel are acquired at the same epochs (within a few days) in each glaciological year. L2swath data are gridded at 500 m × 500 m spatial resolution and for each pixel, we model elevation z(x,y,t) using a linear relationship in space and time: where x, y and t are measured easting, northing, and acquisition time, respectively, and c 0 , c 1 , h, c 2 are the model coefficients. The time-dependent coefficient ḣretrieved from the model fit is the linear rate of surface elevation change for each given pixel. Each observation is assigned a weight according to the sample power as in Gourmelen et al. (2017a). We iteratively fit the model to the data using 3σ clipping until there are no more outliers. The formal uncertainty ε ḣo n each pixel's rate of elevation change ḣis extracted from the model covariance matrix M: where p is the vector of coefficients [c 0 c 1 ḣc 2 ] of the model parameters, z are the input elevations and G = [x y t 1] is the model matrix.
We simplify the data covariance matrix cov(z) to a variance matrix whose diagonal values are the squared elevation differences between the observed and modelled estimates (z-z′) 2 . The square root of the diagonal elements of P represents the standard deviations of the model parameters p.
Due to the complex topography (see Discussion), the ḣmaps do not have complete coverage. We use the relation between elevation and elevation change to model estimates for the gaps in the maps of surface elevation change rates (i.e. hypsometric averaging, see e.g. Moholdt et al., 2010;Nilsson et al., 2015a;Foresta et al., 2016), using the SRTM DEM for the elevation field. Polynomial models of order 1 to 3 are fitted to the data and used to generate elevation change rates for the individual pixels without an estimate. In order to avoid over-fitting the data, the F-test is used to evaluate if the improvement of the additional model parameters on the fit is statistically significant at the 99% confidence level. The median rate of elevation change is then computed in each 50 m elevation band k and multiplied to the area A k , extracted from the DEM, to produce elevation dependent volume change V̇k. The sum of these contributions represents the total rate of volume change. Uncertainties are calculated by error propagation using the same method as in Foresta et al. (2016), summarized in Appendix C.
This interpolation scheme is applied independently for the NPI and for different sub-regions of the SPI displaying highly contrasting patterns of change at similar elevations (Fig. 2). Finally, we assume that all changes relate to the gain or loss of ice of density ρ ice = 900 kg m −3 when converting the rate of volume change to mass balance. This simplification is based on the assumption that at least part of the observed changes are due to dynamics (see Discussion) and ignores possible differences in snow pack density below and above the firn line. To explore mass loss related to material with lower density, we calculate mass balance estimates using a dual density scenario. In this case the densities of glacial ice and firn are used when converting volumetric changes occurring, respectively, in the ablation and accumulation areas. We assign ρ firn = 600 kg m −3 (Malz et al., 2018 and references within). We used Equilibrium Line Altitude (ELA) values as reported in De Angelis (2014) and Barcaza et al. (2009), respectively, for the glaciers on the SPI and NPI. For each group of glaciers (SPI-G 1 , SPI-G 2 , NPI), we computed an average ELA from all glaciers with surface area larger than 100 km 2 in the given catchment.
Glacier outlines over the Ice Fields record their extent in -2001(RGI Consortium, 2017 and the Upsala and Jorge Montt glaciers (SPI) have retreated considerably since then and, for the latter, even during our study period. Their front positions in 2017 were manually digitized using Landsat8 scenes (Appendix D, Table D1) and their mass loss between 2011 and 2017 is calculated against their updated fronts. Area changes between 2011 and 2017 are not included in our estimates of mass loss. The only exception is Jorge Montt (SPI), which retreated considerably in this time period and for which we provide an additional estimate of mass loss due to area change. The front outline of the Jorge Montt and Upsala, as well as of Pio XI (SPI), was additionally digitized for a number of years between 2005 and 2017 (Appendix D, Table D1). This data is used for context in the Discussion (subsections Jorge Montt, Upsala and Pio XI) and is not employed for calculating mass balance.
Finally, we produce time series of mean observed glacier elevation change with the same methodology as Gray et al. (2015) and Foresta et al. (2016). The time series are generated at the catchment scale for each of the nine sub-regions with 90 (Pio XI, SPI-G 2 , SPI) or 120 days time step (Fig. 3).

Results
Swath processing of CryoSat-2 SARIn data provides 6.7 and 26.6 million valid observations of ice topography over the NPI and SPI, respectively, with the rate of elevation change for a single pixel being constrained by~1700 elevations (median) over a period of 5.6 years (median) between April 2011 and March 2017. For comparison conventional CS2 POCA delivers about 17,000 and 55,000 observations over the NPI and SPI respectively. Fig. 1 displays the maps of observed rates of elevation change over the Ice Fields. On the SPI, different catchments show distinct patterns of change over the study period ( Fig. 1). Given such heterogeneity, we apply the hypsometric averaging model independently for six large glaciers on the SPI, namely Jorge Montt, Viedma, Upsala, Pio XI, Grey and Tyndall. The spatial coverage of ḣestimates, at 500 m posting, ranges between 61 and 73% of total catchment areas (Table 1), with the exception of Grey and Tyndall (~52%). We combine data from the rest of the SPI in two groups of neighbouring glaciers, labelled SPI-G 1 and SPI-G 2 . The former includes all glaciers north of Pio XI and Viedma excluding Jorge Montt, while the latter is composed of all the glaciers west and south of Upsala excluding Grey and Tyndall (Fig. 1). Over the NPI, all glaciers are analysed together.
Combining data from different glaciers is needed if coverage is insufficient to provide a representative figure of elevation change in each and is justified providing that they show a similar pattern of change.
The observed median rates of elevation change for these nine subregions are shown as a function of cumulative surface aerial extent (10% steps; Fig.1, side panels) and as a function of elevation (50 m steps; Fig.2, side panels). Widespread thinning is occurring in the northern part of the SPI across all elevations, with average rates of 2 m a −1 on the plateau up to and above 1400 m elevation (SPI-G 1 ) and in excess of 10 m a −1 at the terminus margins of Jorge Montt (tidewater), Viedma and Upsala (both lacustrine) glaciers. Most glaciers in the south/southwest are close to balance (SPI-G 2 ), with the exception of Grey and Tyndall on the eastern side. Pio XI glacier is thickening at rates of~2 and 1 m a −1 below 1000 m and between 1000 and 1500 m elevation, respectively, and thinning by about 1 m a −1 above 1500 m altitude. Similar to the northern part of the SPI, the NPI is experiencing widespread thinning of up to 8 m a −1 with the exception of the ice divide close to the eastern margin of the ice field (Figs. 1 and 2). Hypsometric averaging is applied in each sub-region (Fig. 2, red lines) to generate maps of modelled elevation change rates for the Ice Fields ( Fig. 2), from which mass change is computed ( Table 1). The polynomial models (Fig. 2, red lines) compare well with the observed median rates of elevation change (Fig. 2, black lines with dots). A few exceptions are visible over the Tyndall and Upsala glaciers (at low and high elevation respectively) as well as over the Pio XI glacier below 1000 m and above 2500 m elevation. Model misfits have marginal impact on the glacier mass change when glacier area is negligible or data coverage is high (Fig. 2). For example the rate of mass loss of glacier Tyndall, assuming no elevation change below 350 m elevation, is reduced by about 9% (or 0.054 Gt a −1 ), which is well within its associated uncertainty (Table 1). Similarly for Pio XI glacier, the thinning predicted by the model above 2500 m elevation reduces the glacier net mass gain by only 6% (or 0.04 Gt a −1 ). However, at low elevation where the area of Pio XI glacier is not negligible and where there are no observations to constrain the elevation change (Figs. 1 and 2), the impact of the misfit on the glacier mass balance may be significant (see Discussion). Between April 2011 and April 2017, the NPI and SPI have been losing mass at rates of −6.79 ± 1.16 and − 14.50 ± 1.60 Gt a −1 , respectively, contributing 0.059 ± 0.005 mm a −1 to eustatic SLR. About 35% of the SPI mass loss is concentrated on glaciers in the SPI-G 1 group (−5.07 ± 0.79 Gt a −1 ), which represent~28% of the SPI surface. The Upsala glacier is the single largest contributor to the mass loss (−2.68 ± 0.40 Gt a −1 ) and is also the glacier with the second highest rate of loss per unit area (Table 1) after Jorge Montt. Both glaciers have retreated between 2011 and 2017, by about 0.5 and 2.5 km, respectively. Pio XI is the only glacier in the Patagonian Ice Fields with positive mass balance (0.67 ± 0.29 Gt a −1 ). Its southern tidewater and northern lacustrine termini have both advanced, respectively by about 500 m and 800 m during our study period. Using a dual density scenario, the rates of mass change in the nine sub-regions are lower by 11-19% compared to using the density of glacial ice at all elevations and the total rate of mass loss of the Ice Fields is 17.89 ± 2.03 Gt a −1 (Table 2). For most basins, dynamic processes are dominating the mass loss and hence the ice density scenario is the preferred option. However, in a few sectors the dual density scenario may be more accurate. For example, over the Pio XI glacier (SPI), where surface thickening is suspected (Malz et al., 2018), the mass change is 0.67 ± 0.29 Gt a −1 and 0.57 ± 0.25 Gt a −1 for the single and dual density scenarios, respectively.
Time series of mean observed glacier elevation change (Fig. 3) show negative trends for all sub-regions with the exception of Pio XI, which shows increasing elevation. Most sub-regions display a seasonal oscillation on the order of 1-3 m. The amplitude is highest (4 m) for the Grey Glacier, while it is less discernible for glaciers with the strongest mass losses per unit area (Jorge Montt and Upsala), possibly reflecting the importance of dynamic thinning also during the accumulation period.

Spatial coverage
The Patagonian Ice Fields are a challenging region for radar altimetry. The topography is similar to mountain glaciers, with elevation ranging from sea level to above 2000 m over distances of < 30 km. Furthermore, the flow of most glaciers on the Ice Fields is almost perpendicular to CS2's approximately north-south flight direction so that elevations change abruptly (> 1000 m) over short distances (400 m) along the flight track (Fig. 4), increasing the occurrence of loss-of-lock in the altimetric record and leading to gaps in the collected data (Dehecq et al., 2013). Over the Southern Patagonian ice field, conventional CS2 POCA altimetry provides about 30% spatial coverage at 500 m posting. Although swath altimetry is affected by loss-of-lock as much as POCA, enhanced spatial coverage is achieved because a swath of heights, rather than one single elevation, is acquired when the on board tracker correctly sets the range window. Swath altimetry thereby provides 61-73% surface coverage over the large (A > 400 km 2 ) glaciers in the northern part of the SPI and between 47% and 54% in other areas ( Table 1). The only exception is SPI-G 2 (39%), where the ice field is narrower and where there are no observations over a number of glaciers with relatively small surface area ( Fig. 1 and Fig. 2, inset). Despite the limited extent, their mass loss may be non-negligible. For example HPS12 (south of Pio XI) has an area of 165 km 2 and lost 0.63 Gt a −1 between 2000 and 2011/12 (Willis et al., 2012b). In comparison, work based on high resolution radar TanDEM-X DEMs have almost complete coverage at higher spatial resolution (Jaber et al., 2013;Jaber, 2016;Malz et al., 2018), although this data does not allow yet to generating time series of elevation change. Compared to TanDEM-X DEMs, optical ASTER DEMs achieve similarly high spatial coverage for decadal periods, decreasing to 57-73% for the entirety of the Ice Fields over shorter time periods comparable to that in this paper (Willis et al., 2012b). Despite CS2 L2swath's lower spatial coverage, we still capture in detail the various patterns of change. Furthermore, using a single sensor and frequent repeat measurements is advantageous as it limits penetration biases associated both with using multiple sensor types (e.g. optical vs radar or radar with varying wavelengths) and seasonal variations in surface mass density, thereby limiting the impact on surface elevation change estimates (Jaber et al., 2013;Jaber, 2016;Willis et al., 2012b;Malz et al., 2018). Given the similarity of the Patagonian Ice Fields to mountain glaciers, swath altimetry may also provide one additional tool for monitoring elevation change over these complex areas (Paul et al., 2015).     (Table 1), including the SPI as a whole (grey), in order of increasingly negative specific mass balance (top to bottom).

Rates of mass change
NPI and the northern part of the SPI (SPI-G 1 , Jorge Montt, Viedma and Upsala) are thinning very rapidly and account for 89% of the mass loss of the Patagonian Ice Fields between 2011 and 2017 (Table 1). The rate of mass loss of the Patagonian Ice Fields has increased in recent decades (Fig. 5a), with our estimated mass loss rate (21.29 ± 1.98 Gt a −1 ) being 46% higher than between 1944 and 1996 (Aniya, 1999), 42% higher than between 1975 and 2000 (Rignot et al., 2003) and 24% higher than between 2000 and 2012/14 (Jaber, 2016). However, for the period 2000-2011/12 Willis et al. (2012b) estimate a total rate loss of 24.39 ± 1.20 Gt a −1 , comparable to those based on gravimetry data (Chen et al., 2007;Ivins et al., 2011;Jacob et al., 2012;Gardner et al., 2013) but 30% more negative than that of Jaber (2016) (Fig. 5a). GRACE-based estimates rely on model predicted corrections for postglacial rebound and land water storage, which are a large source of uncertainty to the estimated rates of mass loss in Patagonia (e.g. Chen et al., 2007;Jacob et al., 2012). Additionally, mass loss from glaciers and ice fields in the vicinity (Möller and Schneider, 2008;Melkonian et al., 2013;Falaschi et al., 2017) may impact on the estimate since gravimetry methods are always sensitive to mass leakage effects from neighbouring areas (e.g. Sørensen et al., 2017). The disagreement between Willis et al. (2012b) and Jaber (2016) may be related to the 2 m elevation correction applied to the SRTM data by Willis et al. (2012b) in order to account for potential radar penetration through the glacier surface. However, analysis of the SRTM mean backscattering coefficient suggests wet surface conditions on the Ice Fields at the time of the SRTM acquisition (Jaber, 2016), which would be expected to prevent the radar signal from penetrating through the surface (Nilsson et al., 2015b). Additionally, Dussaillant et al. (2018) find that radar penetration over the NPI occurs only above 2900 m elevation, i.e. < 0.75% of the ice field's area. In absolute value, this correction has a larger impact on the SPI than on the NPI (Willis et al., 2012b), where the estimates from Willis et al. (2012b) and Jaber (2016) differ by only~10%. Separating the contributions of the Ice Fields (Fig. 5b) shows the difference in the progressive increase in mass loss between the NPI and SPI. Between 2011 and 2017 the NPI's rate of loss (6.79 ± 1.16 Gt a −1 ) is 70% more negative compared to the previous decade (about 4 Gt a −1 , Willis et al., 2012aWillis et al., , 2012bJaber, 2016;Dussaillant et al., 2018), which in turn was~37% higher than between 1975 and 2000 (Rignot et al., 2003) (Fig. 5b). Compared to the latter, Rivera et al. (2007) record higher rates of mass loss in a similar time period   (Fig. 5b), but their estimate is based on data mostly lying in the ablation zone of the NPI. In contrast, the mass loss over the SPI varies by just 8-10% between these three periods (Rignot et al., 2003;Jaber, 2016;Malz et al., 2018) (Fig. 5b) although its rate of loss is still more than twice that of the NPI. Estimates from Jaber (2016) and Malz et al. (2018), both based on comparing TanDEM-X data to the SRTM DEM, differ by 1.29 Gt a −1 although they agree within their uncertainties. The difference may be due to the 4 years longer time period analysed by Malz et al. (2018), who report positive elevation changes in the southernmost part of the SPI between 2011/12 and 2015/16. We observe only a slightly positive trend in elevation change in this time period, followed by a marked drop after 2015/16 (Fig. 3, SPI-G 2 ). However, our time series for SPI-G 2 is representative of an area roughly twice as large than that analysed by Malz et al. (2018) over multiple time periods. Finally, the estimated SPI's rate loss of 34.83 ± 3.96 Gt a −1 between 1995 and 2000 (Rignot et al., 2003), which is even more negative than any estimate for both Ice Fields combined (Fig. 5a), appears out of line with other values.

Glacier dynamics
We observe a sharp transition around 49°S (Figs. 1-2, dashed green line) from intense thinning in the north to a large area facing limited mass loss (SPI-G 2 ), which spans about 4800 km 2 or about 40% of the total surface of the SPI (Table 2). This pattern is in agreement with earlier observations between 2000 and 2012/16 (Willis et al., 2012b;Jaber et al., 2013;Jaber, 2016;Malz et al., 2018), therefore suggesting constant behaviour over decadal time scales. The topography of the northern and southern parts of the SPI has different characteristics, with a greater proportion of the northerly ice field lying at lower altitudes. In fact, about 72% of SPI-G 1 's surface lies below 1500 m elevation,~12% more than SPI-G 2 's at the same altitude (Fig. 1, insets SPI-G 1 and SPI-G 2 ). The Ice Field also narrows and steepens south of 49°S and even at low elevations the southern SPI shows only moderate thinning (Fig. 6). Glaciers Grey and Tyndall, at the southeastern tip, are the exception to this pattern. However the latter lies almost entirely below 1500 m altitude (Fig. 1, inset Tyndall) and both glaciers receive scarce precipitation due to their location east of the ice divide. The NPI has similar area-altitude distribution as SPI-G 1 (Fig. 1, inset SPI-G 1 and NPI) and the two areas show comparable mass losses per unit area (Table 1).
However, the northern part of the SPI contains some of the fastest flowing glaciers on the Ice Fields, including Jorge Montt and Upsala (Sakakibara and Sugiyama, 2014;Mouginot and Rignot, 2015) which accelerated significantly (> 500 m a −1 ) in the period 1984-2000 (Jorge Montt) and , coincident with rapid frontal retreat (Sakakibara and Sugiyama, 2014). These observations confirm the importance of dynamics in impacting the overall mass balance of the northern part of the SPI (e.g. Sakakibara and Sugiyama, 2014;Mouginot and Rignot, 2015), where three of the largest glaciers (Jorge Montt, Upsala and Viedma) are thinning very rapidly (Table 1 and Fig. 1-3).

Jorge Montt glacier
Jorge Montt, a tidewater glacier at the northernmost tip of the SPI, has been retreating since 1898 when it reached its LIA maximum extent . Its recession has been linked to fjord water depth, with periods of stable front positions corresponding to shallow depths and underwater pinning points . Additionally, water temperatures at depth (> 100 m) have been shown to be as high as 8°C in summer 2012 only 1 km from the glacier front (Moffat, 2014), which may further destabilise the glacier through submarine melting (Straneo and Heimbach, 2013). By 2011 Jorge Montt had retreated almost 20 km w.r.t 1898, with the highest rates of recession between 1990 and 1997 (993 m a −1 ) occurring when water depths beneath the glacier increased sharply to~300 m . The recent retreat history of Jorge Montt's reveals a slowdown to about 100-300 m a −1 in the early 2000s, followed by increased retreat after 2009 initiated at a location where bathymetry data reveals the deepest trough in the fjord Moffat, 2014; Fig. 7).
Between 2010 and 2011, Jorge Montt retreated almost 1 km Rivera et al., 2012b) and calved at a rate of 2.4 km 3 a −1 , when the terminus was floating (Rivera et al., 2012a). Manual delineation of the glacier front between 2011 and 2017 using Landsat optical data (Fig. 7) indicates that Jorge Montt retreated by an additional~2.5 km, likely through enhanced calving following retreat into deeper water (Rivera et al., 2012a). Given an average glacier freeboard height of 22 m above sea level at the terminus (Rivera et al., 2012a and width of 1.15 km, the glacier frontal retreat amounts to a mass loss rate of 0.07 Gt a −1 (~3% of the catchment's loss due to thinning, Table 1) between 2011 and 2017. This value is however likely underestimated since the slope of the glacier surface, and thus upglacier thickening, has not been considered. Due to the uncertainty associated with this calculation, and that at least part of the glacier terminus was floating in 2010/11 (Rivera et al., 2012a) and likely during our observational period, we report this loss separately in Table 1 and do not include it in our total estimate of glacier contribution to sea level rise. Between 2011 and 2017, Jorge Montt shows the highest mass loss per unit area, 4 times above the average for the SPI as a whole (Table 1). Its absolute rate of mass loss (2.20 ± 0.38 Gt a −1 ) is comparable to what reported by Jaber (2016) for the period 2011-2014 (2.59 Gt a −1 , uncertainty not reported), which increased by 50% compared to the 1.72 Gt a −1 (uncertainty not reported) rate of mass loss between 2000 and 2011 (Jaber, 2016).

Upsala glacier
Upsala, a freshwater calving glacier located on the north-eastern side of the SPI and draining into Argentino Lake, has also been retreating since the late 1970s (Naruse et al., 1997). Between 2008 and 2011, retreat rates quadrupled w.r.t the previous 8 years and the glacier retreated by almost 3 km (Sakakibara et al., 2013), while simultaneously speeding up by   20-50% (Sakakibara et al., 2013;Mouginot and Rignot, 2015) and thinning at a maximum rate of~40 m a −1 near the terminus (Sakakibara et al., 2013). The rapid retreat may have been caused by the glacier front reaching an area where the lake depth exceeds 560 m (Sugiyama et al., 2016). In early 2013 the glacier's velocity at the front was 2.9 m d −1 (Moragues et al., 2018),~33% lower compared to 2008 and more similar to values recorded in the early 2000s (Sakakibara et al., 2013). Moragues et al. (2018) report a doubling in maximum velocity between 2013 and 2014. However this increase is unlikely to have been sustained in time. In fact, between 2011 and 2017, the glacier front has been comparatively stable (Fig. 8), with a retreat rate of~85 m a −1 similar to the period 2000-2008 (Sakakibara et al., 2013). Coincident with a more stable front position, the average thinning rate within 16 km of the terminus decreased by a factor two from 13.4 m a −1 between 2006 and 2010 (Sakakibara et al., 2013) to 6.2 m a −1 between 2011 and 2017 (Fig. 8). Bertacchi Glacier, a tributary of Upsala, shows a similar pattern with current rates of elevation change decreasing to 8.5 m a −1 (Fig. 8) from~15 m a −1 between 2008 and 2011 (Sakakibara et al., 2013). We observe maximum thinning rates of~12 m a −1 5 km from the terminus of Upsala glacier, comparable to estimates at the front in the early 1990s (Naruse et al., 1997). Despite these reduced thinning rates, Upsala remains the glacier with the second most negative specific mass balance in the Patagonian Ice Fields after Jorge Montt, and is the largest single contributor to net mass loss amongst individual glaciers (Table 1 and Figs. 1-3).

Pio XI glacier
Pio XI, the largest glacier on the SPI and in South America, is the only glacier of the Patagonian Ice Fields to have experienced a net large advance since 1926 and the only known surge-type glacier on the SPI (Rivera et al., 1997a;Wilson et al., 2016). Published data of frontal changes, ice flow velocity at the termini, elevation change and mass balance, summarized in Fig. E1, reveal a complex history (Appendix E). Between 1951 and 1963, the glacier's westward and southward advance dammed a proglacial river originating from Greve glacier to the north, forming the current Lake Greve for at least the second time since 1926 (Rivera et al., 1997a). Since then, the glacier has been terminating in Ejre Fjord to the south and Lake Greve to the north. From 1945 to present, the tidewater terminus advanced~13 km and is currently at its Neoglacial maximum (Wilson et al., 2016;Fig. 9 and Fig. E1). Looped supraglacial moraines were used to identify up to six surge events since 1926 (Rivera et al., 1997a;Wilson et al., 2016), two of which were concurrent with front retreat, possibly due to enhanced calving at the tidewater terminus (Wilson et al., 2016;Fig. E1). The glacier has been thickening in the ablation area since the late 1970s, with the highest rates recorded at the termini, while the picture is more complicated in the accumulation area where data is sparse (Fig. E1). Between 2011 and 2017, we observe thickening at almost all elevations, by about 2.33 m a −1 and 0.57 m a −1 (median value) in the ablation and accumulation areas respectively ( Fig. 1; Pio XI inset) with thinning at the highest elevations above 1500 m altitude; findings which match those described in Jaber (2016). There is however no coverage in our data close to the two termini. Assuming a thickening trend at the two fronts, which was sustained for the last four decades (Rignot et al., 2003;Willis et al., 2012b;Jaber et al., 2013;Jaber, 2016;Wilson et al., 2016), we estimate that the Pio XI glacier is gaining mass at a rate of 0.67 ± 0.29 Gt a −1 between 2011 and 2017 (Table 1 and Fig. E1). However the rate is likely underestimated since the hypsometric averaging model for Pio XI glacier predicts less thickening compared to the observations below 1000 m elevation (Fig. 2, inset). The mass gain is likely a result of complex dynamics associated with both surge mechanisms and terminus calving processes, since Pio XI is the only advancing glacier within the Patagonian Ice Fields and air temperature has increased over the last 50 years (Rasmussen et al., 2007).

Conclusions
CryoSat-2 swath radar altimetry is employed successfully to map elevation change over the Patagonian Ice Fields at sub-kilometer spatial resolution. Despite the challenging topography, similar to that of mountain glaciers, the technique can be used to observe changes over individual glaciers or catchments with an area as small as 300 km 2 . The northern part of the SPI displays a high degree of complexity, although most of the area is thinning at all elevations, with Jorge Montt, Viedma and Upsala glaciers losing mass at rates higher than 2 Gt a −1 . Jorge Montt additionally retreated~2.5 km between 2011 and 2017, likely by enhanced calving in deeper fjord waters. The only exception to the overall pattern of thinning and retreat is the Pio XI glacier, which continues to advance at both its tidewater and lacustrine termini. The glacier, which is currently at its Neoglacial maximum, is estimated to have gained mass at a rate of 0.67 ± 0.29 Gt a −1 during our study period.
Between April 2011 and March 2017, the Ice Fields lost mass at a combined rate of 21.29 ± 1.98 Gt a −1 (equivalent to 0.059 ± 0.005 mm a −1 eustatic sea level rise), an increase of 24% and 42% when compared to the periods 2000-2012/14 and 1975-2000, respectively. We find that the NPI (−6.79 ± 1.16 Gt a −1 ), which is responsible for a third of the total loss, is losing mass 70% faster compared to the first decade of the 21st century. Given the ongoing and rapid wastage of the Patagonian Ice Fields, and their important contribution to the global budget of mass loss from glaciers and ice caps, continuous observations with excellent spatial and temporal resolution are essential. CS2 swath altimetry provides an important tool for monitoring these rapidly changing areas and quantifying their ongoing mass loss.

Acknowledgments
This work was performed under the European Space Agency's Support to Science Element CryoSat+ Mountain Glacier (contract 4000114224/ 15/I-SBo). L.F. acknowledges a Young Scientist Training grant from the European Space Agency's Dragon3 program. J.J.W. acknowledges a NERC PhD Studentship. The authors wish to thank ESA for providing free access to CryoSat-2 full bit rate data, NASA for providing free access to the SRTM DEM and USGS for providing free access to Landsat data. The authors are grateful to the Editor-in-Chief, Crystal Schaaf, and to two anonymous reviewers, whose comments have significantly improved the manuscript.

Appendix A. Filtering CS2 SARIn waveform samples
Selecting waveform samples based on fixed thresholds on coherence and power is an empirical approach which has been applied successfully to infer glacial topography and higher products based on it such as topography changes (Gray et al., 2013;Christie et al., 2016;Ignéczi et al., 2016;Foresta et al., 2016;Gourmelen et al., 2017a). However, in this paper we use a different procedure to discard noisy waveform samples before performing the phase unwrapping. This method relies entirely on the phase difference field and consists in identifying, within each waveform, groups of consecutive samples which can be modelled by a straight line. The original phase difference is divided into overlapping segments, their length being set to 64 samples and the overlap half of the length. The slope of the phase difference is then calculated independently in each segment. Instead of applying a linear regression, the algorithm applies a Fourier transform on the normalized complex coherence e (i ðɸ) , where ðɸ is the phase difference field. Compared to linear regression, this approach is both more efficient computationally as well as independent on phase wrapping. The Fourier transform enables to testing a large number of possible slopes and the one with the highest correlation with the input data is selected. The signal is oversampled to take into account that the slope of CS-2's phase difference can represent non-integer frequencies. Thus, each overlapping section has two possible slopes. A correlation is applied again to the data in each overlapping section, this time using only its two estimated slope values. Sections whose correlation is below a set threshold (for this work, 0.95) are considered noisy and discarded. Finally, the remaining segments are used to unwrap the phase difference. With this procedure, no smoothing is applied to the phase difference and no threshold is set on the power or coherence.

Appendix B. SRTM DEM correction
A number of freely available DEMs covering Southern Patagonia exist, namely: the SRTM (i) C and (ii) X band DEMs (Farr et al., 2007), the (iii) ASTER (Advanced Spaceborne Thermal Emission and Reflection Radiometer) GDEM2  and (iv) the ALOS (Advanced Land Observing Satellite) AW3D30 v1.1 Takaku et al., 2014;Tadono et al., 2016;Takaku et al., 2016). The latter is the most recent, but has large gaps over the SPI, which are filled using the SRTM C band DEM. Version 1 of the ASTER GDEM was known to be affected by large artefacts (Arefi and Reinartz, 2011) and, despite large overall improvements, version 2 still has high frequency noise, particularly over glacial terrain (Meyer et al., 2011). Visual comparison between the SRTM C band and ASTER GDEM2 DEMs shows evident noise in the latter, with differences at times on the order of tens of meters between neighbouring pixels. Finally, the SRTM X band DEM does not have complete coverage and gaps over the Ice Fields are significant. Therefore, we used the SRTM C band DEM as a reference for elevation, which fully covers the Patagonian Ice Fields, resampled at 300 m posting and referenced to the WGS84 vertical datum. The down-sampling of the DEM is mostly dictated by achieving a satisfactory performance in computing time while keeping the spatial resolution somewhat comparable to that of an individual elevation based on CryoSat-2 interferometric data (300 m in the along-track direction). We use linear interpolation when querying the DEM. The SRTM DEM is based on data acquired in February 2000, and a few areas at the margins of the SPI have thinned by at least 80 m since then (Willis et al., 2012b;Jaber et al., 2013;Jaber, 2016). This magnitude is comparable to the elevation offset caused by a 2π shift on CS2's phase (Data and methods, this paper; Fig. 3 in Gourmelen et al., 2017aGourmelen et al., , 2017b. Thus, over areas which experienced intense thinning rates, the swath algorithm may select an incorrect 2π multiple which best matches current observations to the glacier topography from 2000. In order to avoid that, the SRTM DEM needs to be registered to the beginning of our study period. To this purpose, we applied a first order correction of the SRTM DEM based on a visual inspection of results in Willis et al. (2012b) and assuming constant rates of elevation change between the SRTM and CS2 periods. This approach was sufficient to improve the phase unwrapping, leading to further pixels meeting the quality criteria for robust rates of surface elevation change. We refer to this corrected DEM simply as the SRTM DEM in this study. The correction of the SRTM DEM only affects the terminus areas of Jorge Montt and Upsala glaciers (SPI) since there are no CS2 swath altimetry observations over smaller glaciers on the Ice Fields which experienced similar thinning rates in the period 2000-2011/12 (e.g. HPS12, SPI; Willis et al., 2012b;Jaber et al., 2013;Jaber, 2016).

Appendix C. Error budget
The errors on the mass balance estimates are calculated as in Foresta et al. (2016). The uncertainties ε ḣo n the observed rates of elevation change for the individual pixels, extracted from the model covariance matrix (see Data and methods), are propagated when applying the hypsometric averaging method: where E ḣ( k) is the elevation change error in elevation band k and N(k) is the number of valid observations in the elevation band. A two-term decreasing exponential is used to interpolate values for elevation bands with no observations (e.g. at low elevation for bands with limited spatial extent). We multiply the area extent A(k) of the elevation band to the related E ḣ( k) and sum all contributions to estimate the total uncertainty on the rate of volume change: V k ḣẆ ith this method, the volume change uncertainty is only related to that area of the ice cap where there are valid rates of surface elevation change, but does not account for incomplete data coverage. The volume change uncertainty is therefore rescaled according to the data coverage (Table 1). This procedure generates a rather conservative (i.e. larger) error estimate since it assumes that the lack in data coverage has a direct impact on the total error estimate, which does not hold if the sampling is sufficiently uniform. Finally, we include an error on the density: when converting volume to mass change (e.g. Nilsson et al., 2015a).  E1 summarises published data of frontal changes, ice flow velocity at the fronts, elevation change and mass balance (Rivera et al., 1997a;Rivera et al., 1997b;Rivera and Casassa, 1999;Rignot et al., 2003;Lopez et al., 2010;Willis et al., 2012b;Jaber et al., 2013;Sakakibara and Sugiyama, 2014;Wilson et al., 2016).  (Rivera et al., 1997a) or retreat (Wilson et al., 2016) of the tidewater front are highlighted in light and dark grey, respectively. All relevant references are listed in the inset.