The Rapid and Steady Mass Loss of the Patagonian Iceﬁelds throughout the GRACE Era: 2002–2017

: We use the complete gravity recovery and climate experiment (GRACE) Level-2 monthly time series to derive the ice mass changes of the Patagonian Iceﬁelds (Southern Andes). The glacial isostatic adjustment is accounted for by a regional model that is constrained by global navigation satellite systems (GNSS) uplift observations. Further corrections are applied concerning the e ﬀ ect of mass variations in the ocean, in the continental water storage, and of the Antarctic ice sheet. The 161 monthly GRACE gravity ﬁeld solutions are inverted in the spatial domain through the adjustment of scaling factors applied to a-priori ice mass change patterns based on published remote sensing results for the Southern and Northern Patagonian Iceﬁelds, respectively. We infer an ice mass change rate of − 24.4 ± 4.7 Gt / a for the Patagonian Iceﬁelds between April 2002 and June 2017, which corresponds to a contribution to the eustatic sea level rise of 0.067 ± 0.013 mm / a. Our time series of monthly ice mass changes reveals no indication for an acceleration in ice mass loss. We ﬁnd indications that the Northern Patagonian Iceﬁeld contributes more to the integral ice loss than previously assumed.


Introduction
The Southern Patagonian Icefield (SPI; 12,500 km 2 ) and Northern Patagonian Icefield (NPI; 4000 km 2 ) constitute, together, the largest temperate ice mass in the Southern Hemisphere. They are aligned in a north-south (N-S) direction along the crest of the Southern Andes (Figure 1). The icefields are strongly influenced by the Westerlies and the transport of wet air mass from the Pacific Ocean, receiving precipitation of up to 10 m/a water equivalent [1]. The quantification of present-day ice mass changes contributes to our understanding of ongoing climate change in the southern mid-latitudes and the southeastern Pacific. An observational assessment of the demise of the two Patagonian Icefields provides robust information for predicting the future of the regional cryosphere and the water cycle, as well as its contribution to the global sea level. Previous works suggest a disproportionally large contribution of the Patagonian Icefields to eustatic sea level rise when compared to their area, for example, 0.067 ± 0.004 mm/a [2], 0.105 ± 0.011 mm/a [3], and 0.059 ± 0.005 mm/a [4]. Bedrock global navigation satellite systems (GNSS) observations adjacent to and within the SPI reveal an uplift as high Southern Patagonia is characterized by a series of distinct Neoglacial advances and retreats during mid-to-late-Holocene times (0-5 ka before present) [10,11]. The most recent advance corresponds to the regional LIA. This history is inferred from datable moraine materials [11] and is Southern Patagonia is characterized by a series of distinct Neoglacial advances and retreats during mid-to-late-Holocene times (0-5 ka before present) [10,11]. The most recent advance corresponds to the regional LIA. This history is inferred from datable moraine materials [11] and is consistent with the regional paleo-climatological record [12,13]. While the climatic origins of these fluctuations are unknown, they are thought to be linked by global teleconnections [14,15]. More importantly, they establish a pattern of natural variability that is now perhaps in sync with anthropogenically driven warm late-summer conditions [16,17]. The latter exacerbates the negative net mass balance of the mountain glacier systems over the past 30-40 years, and that we know to be global in nature [18]. The future of the Patagonian Icefields may well be determined by the climatic conditions generated by the evolving state of the atmosphere and ocean on their Western Flank, as these are, in turn, linked to changes in both the El Niño and Southern Oscillation (ENSO) and the Pacific Decadal Oscillation (PDO), phenomena that we now suspect have been modulated by anthropogenic forcing since the 1990 s [19].
Our study of the 15-year mass evolution of the Patagonian Icefields provides new clues about the present-day land-ice trajectory of the higher latitude Southern Hemisphere. We discover that glacier mass loss is dominated by a strong, and relatively uninterrupted, linear trend, possibly lending new insights to some of the standing questions regarding the future of the Southern Ocean and its environs [20].

Data and Methods
Natural conditions at the Patagonian Icefields limit the applicability of observational techniques that can be used to determine ice mass changes. Table 1 summarizes the previously published present-day ice mass change rates for the Patagonian Icefields. Conventional satellite radar altimetry is not appropriate because of the large footprint size compared to the narrow, fragmented, and often steep ice surface. Satellite laser altimetry provides insufficient data because of the frequent cloud cover and the unfavorable N-S orientation of the icefields. Recent advances using Cryosat-2 swath radar altimetry [4] are promising, but are limited in time to 2010 onward. The surface mass balance estimates derived from atmospheric modelling [17,21,22], or the changes derived from field glaciological measurements, suffer from a sparse distribution of meteorological stations and sampling sites that do not adequately capture the large variability over very short distances and rugged topography. Remote sensing techniques based on both optical [2,3,[23][24][25] and radar imagery [26][27][28] provide valuable data for quantifying the ice volume and ice mass changes of the SPI and NPI, in spite of having to deal with cloud cover or uncertain radar signal penetration. All of the remote sensing results (Table 1) have a partially reconcilable agreement on the net ice mass loss to both icefields. However, the published rates can differ by as much as a factor of 2, a discrepancy that cannot be explained by the differing imagery acquisition time spans. Mouginot and Rignot [29] present a comprehensive determination of the ice flow velocities over SPI and NPI based on satellite data covering 30 years (1984-2014). They reveal a heterogeneous behavior over time among the individual glaciers-some of the major glaciers accelerate significantly, others decelerate or show a reversal in acceleration/deceleration over a decade or so.
Time-variable global gravity field solutions provided by the gravity recovery and climate experiment (GRACE; [30]) satellite mission have proven an efficient basis for the quantification of ice mass changes (e.g., [31][32][33][34][35], as well as the references therein). Previous inversions of GRACE data targeted explicitly on SPI and NPI [36,37] had to rely on relatively short time spans compared to the complete mission record through 2017, along with larger uncertainties in the earlier GRACE Level-2 data releases and less refined correction models available at that time (e.g., for correcting the effects of GIA and the surface mass changes, and for improving the low-degree Stokes coefficients, i.e.: C 10 , C 11 , S 11 , C 20 , C 20 , and S 21 ). More recent ice mass change rates for southern Patagonia are presented as part of the global inversions that account for all mountain glacier and ice cap contributions to sea level rise [32,[38][39][40].
Here, we derive the ice mass changes of the Patagonian Icefields between April 2002 and June 2017 from the monthly ITSG-Grace2016 [41] gravity field solutions. The monthly SPI and NPI ice mass changes are derived from a GRACE Level-2 data analysis approach similar to the regional point mass inversion presented in Forsberg and Reeh [42] for Greenland. This approach has been refined and modified (e.g., [43][44][45]), and it has also been used to derive ice-mass change time series for Antarctica and Greenland for a reconciled estimate of ice sheet mass balance [34,35]. In our case, we used a-priori information on the location of the point masses and on their changes with time. This allows for the set-up of a very dense pattern of masses (SPI: 98, NPI: 76, in total 174), and we determined a common scaling parameter for the a-priori mass change rates.
There is no sharp spatial resolution limit of GRACE. The ability to separate the ice mass changes of the SPI and the NPI (with a 320 km distance between their aerial barycenters and a 80 km minimum distance between their areas) depends on a number of aspects, namely: (i) the noise level of the used GRACE products with respect to the specific geographic situation, (ii) the strength of the signal, and (iii) the exploitation of the a-priori information. The a-priori information on the geographic patterns of the mass change have been widely used to separate nearby signals. This has been shown to have great advantages for dealing with coastlines [46]. We explore the separability of the NPI and SPI signals, because each of these three factors may be readily exploited. (i) ITSG-Grace2016 has a reduced noise level w.r.t. to the previous GRACE Level-2 releases [47]. The higher frequency and shorter wavelength noise of the GRACE gravity field solutions is inherently lower for the along-track (N-S) accelerations detected by the K-Band Ranging (KBR) system, relative to the noise in an east-west (E-W) (track-normal) direction, wherein the KBR system does not directly provide gravity information [48]. (ii) The sustained mass change of the Patagonian Icefields is among the largest spatially concentrated long-term signals worldwide. (iii) We can exploit the a-priori information on the distribution of the glacier mass loss of the two icefields along with the extensive a-priori information on the disturbing signals outside the icefields, as explained below and in Appendix A. Table 1. Summary of ice mass change rates. The "comment" indicates the spatial extent for analyses using gravity recovery and climate experiment (GRACE) and Cryosat-2, and the data sources in the case of the remote sensing studies. SPI-Southern Patagonian Icefield; NPI-Northern Patagonian Icefield; SRTM-Shuttle Radar Topography Mission. In the following section, our inversion approach is described, and further details of the analysis procedure are given in the Appendix A. We used unconstrained monthly gravity field solutions of release ITSG-Grace2016 provided by TU Graz [41]. These are provided as Stokes coefficients to a spherical harmonic degree and order of 120. This time series comprises 161 monthly solutions that cover the period April 2002-June 2017. We added degree one coefficients (C 10 , C 11 , and S 11 ) derived from the ITSG-Grace2016 solutions [49]. The C 20 -term was replaced by the values obtained by satellite laser ranging [50]. The C 21 and S 21 coefficients were corrected for the pole tide, as recommended by Wahr et al. [51]. The monthly Stokes coefficients were corrected for GIA by applying the regional model A according to Lange et al. [6] (henceforth referred to as La-A), after removing the elastic contribution to the uplift. The Stokes coefficients were then further corrected for the effect of elastic loading deformation [44] (Equation (2)).

Reference
We retrieved GRACE pseudo-observations on a regional grid at the approximate GRACE satellite orbit altitude (500 km) from each monthly set of corrected Stokes coefficients [42][43][44][45]. Figure 1a shows the extension of this grid. The spacing is approximately 20 km, with a constant longitude increment and latitude increments that provide for an equal area of all of the grid cells. We chose to generate our pseudo-observations at the orbit altitude in order to avoid the amplification of noise in the downward continuation. No explicit filtering was applied to the GRACE data. For each grid node and month, 3D Cartesian gravity vector components were computed (instead of the scalar gravity disturbances used by Forsberg et al. [44] as pseudo-observations). These monthly sets of 3D gravity vectors were interpreted with respect to the mass redistributions at the Earth's surface. For comparison, an additional mass corresponding to a uniform layer of 1cm water equivalent of global extent at the surface (sea level) implies a change in the length of our gravity vectors at the orbit altitude of +7.22 nm s −2 .
The monthly gravity vector grids were corrected explicitly for five surface mass variation processes-the Antarctic ice sheet (AIS) [52], ocean mass change, changes in continental water storage [53,54], ice mass change rates in South America outside Patagonia [26], and water mass changes in the great Patagonian lakes [55]. Our ocean mass correction consisted of a gravitationally self-consistent distribution of the mass discharge from five major ice sheets/ice caps (AIS, Greenland ice sheet, Patagonian Icefields, Alaska, and the Canadian Arctic) and continental water storage over the ocean [56] (using ice mass change estimates [52,57] and continental water storage [53] as input), complemented by an empirical regional ocean mass change rate increment. The mass signals from outside our data grid domain may affect our pseudo-observations through leakage [44]. Our approach, with the global extent of the correction models (ocean mass and continental water storage) plus the additional solved-for ocean mass change rate increment, ensures that possible leakage effects are minimized. The time series of the pseudo-observations prior and after the correction, as well as of the individual mass-variation corrections, are shown for one data grid node in Figure 2. The individual impact of the applied corrections on the SPI and NPI ice mass change rate is given in Table 2. Figure 3 shows maps of the gravity rate components prior (top) and after (center) the application of all corrections. The relative change in gravity at the approximate GRACE orbit altitude prior (uncorrected, grey) and after (corrected, black) the application of the corrections described in the text is shown in the radial (top), north (center), and east (bottom) vector components. The applied corrections for the five mass-variation processes are shown in different colors. The location of this grid node is indicated in Figure 1 (yellow circle).

Figure 2.
Time series of the gravity signal and the applied corrections for a selected data grid node. The relative change in gravity at the approximate GRACE orbit altitude prior (uncorrected, grey) and after (corrected, black) the application of the corrections described in the text is shown in the radial (top), north (center), and east (bottom) vector components. The applied corrections for the five mass-variation processes are shown in different colors. The location of this grid node is indicated in Figure 1 (yellow circle).   Table 2. Impact of individual corrections applied in the GRACE data analysis scheme on the derived SPI and NPI ice mass change rate. Column 2 is the difference between our solution, including all of the corrections minus a modified solution, including all but the correction (or employing the alternative data set) specified in column 1. Column 3 indicates the percentage of column 2 with respect to the rate including all of the corrections (−24.4 Gt/a). Column 4 shows, for each modified solution, the ratio of the gravity vector RMS over the regional grid between the residuals after the adjustment and the corrected observations (for comparison, this ratio amounts to 0.152 for the solution including all corrections). We synthesized an a-priori pattern of the ice mass change at SPI and NPI. A scaling factor for this mass change pattern was derived from each monthly field of the corrected gravity vectors, simultaneously with the parameters of the empirical regional ocean mass change rate increment, by a least-squares adjustment without any regularization. In this way, a time series of the integral ice mass change was obtained ( Figure 4). First, an SPI and NPI a-priori pattern consisting of 174 point masses was derived from the average ice surface elevation rates for the accumulation and ablation zones of 87 glacier basins [2] (see Appendix A). However, these ice surface elevation rates were based on slightly different imagery acquisition time spans and processing algorithms for NPI [23] and SPI [2]. Therefore, a systematic bias between the icefields cannot be ruled out. We used the GRACE data to optimize the a-priori ice mass change pattern by deriving an individual scaling factor for each of the two subsets of point masses that represent SPI and NPI, respectively. Thus, two separate scaling factors were determined simultaneously, for which our adjustment yields a minimum residual gravity rate RMS. The monthly RMS values of the residual misfit between the gravity effect of the scaled optimum ice mass change pattern and the corrected gravity vector fields are included in Figure 4 (bottom).

Correction
A weighted fit of a linear plus periodic (annual and semi-annual sinusoid) model (red in Figure 4) to the time series yields the ice mass change rate for the SPI and NPI region. The GRACE Level-2 data quality degrades after March 2016, especially since the decommissioning of the accelerometer on the GRACE-B satellite in November 2016 [58] (see also Appendix A). We took this heterogeneity into account in the ice mass change rate estimation by an empirical downweightTing of the monthly ice mass change values after March 2016. For this purpose, the uncertainty of the monthly ice mass change estimates was calculated separately until and after that month. We obtained uncertainties (1σ) of ±22.7 Gt and ±53.5 Gt, respectively. This documents a noise increase by a factor of 2.4. The reciprocal of this value was introduced as weights for the monthly ice mass change values after March 2016 in the rate estimation. In addition, the Level-2 data for February 2015 was impaired by a two-day repeat orbital resonance [59]. Therefore, rate and accuracy estimates do not include this month.

Uncertainty Estimation
The uncertainty of the derived ice mass change rate is composed of the uncertainty of the individual monthly ice mass estimates and the uncertainties in the applied corrections. The residuals from the linear plus periodic model are used to estimate the precision of the ice mass change values [60]. The individual precision estimates for the two sub-periods, prior and after March 2016, lead to an uncertainty of the derived ice mass change rate of ±2.3 Gt/a.
The effect of all of the individual corrections applied on the ice mass change rate is summarized in Table 2. As a conservative estimate, we assume a 100% relative uncertainty for all of the corrections, that is, the corrections of the low-degree Stokes coefficients (lines 1-3 in Table 2) and the corrections for the surface mass change processes (lines 7-14 in Table 2). The total uncertainty of these corrections and the uncertainty of the monthly ice mass estimates (±2.3 Gt/a) is calculated as the root sum square of the individual components.
A crucial source of the uncertainty for the derived ice mass change rate is the applied GIA model. The GIA model La-A adopted in our correction is constrained by the GNSS observed uplift rates [6,7], and implies a (rock) mass gain equivalent to 9.1 Gt/a ( Table 2). Employing the alternative GIA model La-B [6] (which also satisfies the uplift data; 10.5 Gt/a) in our GIA correction results in a decrease in the derived mass change rate by −1.4 Gt/a ( Table 2). We adopted the standard deviation between both of the ice mass change rates as an estimate of the systematic impact of the uncertainty in the GIA correction. Adding this uncertainty contribution of ±1.0 Gt/a leads to a total uncertainty (1σ) of the ice mass change rate for the SPI and NPI region of ±4.7 Gt/a.

Results and Discussion
We obtained an ice mass change rate in the SPI and NPI region of −24.4 ± 4.7 Gt/a. The RMS of the residual gravity rate vectors over our grid domain (Figure 3, bottom) amounted to 0.38 nm s⁻² a⁻¹

Uncertainty Estimation
The uncertainty of the derived ice mass change rate is composed of the uncertainty of the individual monthly ice mass estimates and the uncertainties in the applied corrections. The residuals from the linear plus periodic model are used to estimate the precision of the ice mass change values [60]. The individual precision estimates for the two sub-periods, prior and after March 2016, lead to an uncertainty of the derived ice mass change rate of ±2.3 Gt/a.
The effect of all of the individual corrections applied on the ice mass change rate is summarized in Table 2. As a conservative estimate, we assume a 100% relative uncertainty for all of the corrections, that is, the corrections of the low-degree Stokes coefficients (lines 1-3 in Table 2) and the corrections for the surface mass change processes (lines 7-14 in Table 2). The total uncertainty of these corrections and the uncertainty of the monthly ice mass estimates (±2.3 Gt/a) is calculated as the root sum square of the individual components.
A crucial source of the uncertainty for the derived ice mass change rate is the applied GIA model. The GIA model La-A adopted in our correction is constrained by the GNSS observed uplift rates [6,7], and implies a (rock) mass gain equivalent to 9.1 Gt/a ( Table 2). Employing the alternative GIA model La-B [6] (which also satisfies the uplift data; 10.5 Gt/a) in our GIA correction results in a decrease in the derived mass change rate by −1.4 Gt/a ( Table 2). We adopted the standard deviation between both of the ice mass change rates as an estimate of the systematic impact of the uncertainty in the GIA correction. Adding this uncertainty contribution of ±1.0 Gt/a leads to a total uncertainty (1σ) of the ice mass change rate for the SPI and NPI region of ±4.7 Gt/a.

Results and Discussion
We obtained an ice mass change rate in the SPI and NPI region of −24.4 ± 4.7 Gt/a. The RMS of the residual gravity rate vectors over our grid domain (Figure 3, bottom) amounted to 0.38 nm s −2 a −1 compared with the 2.52 nm s −2 a −1 RMS of the observed gravity vectors with all of the corrections applied (Figure 3, center). This corresponds to a signal variance reduction of 98% and demonstrates the efficiency of the scaled ice mass change pattern to explain the observed gravity changes. The determination of a common factor for the initial SPI and NPI a-priori pattern as described in Section 2 yields a residual RMS value of 0.44 nm s −2 a −1 . A significant reduction of the RMS (down to the above-mentioned value of 0.38 nm s −2 a −1 ) is achieved when scaling the a-priori patterns of SPI and NPI separately. It is well-known that GRACE's ability to resolve masses close to each other is limited, and leakage effects may occur. In order to explore the sensitivity of the GRACE data concerning the separation of the SPI and NPI ice mass loss we carried out a forward computation of gravity rate vector components in our data grid according to the a-priori ice mass change pattern derived from Willis et al. [2]. This forward computation is repeated with systematically altered scaling factors for both SPI and NPI. For each pair of scaling factors (corresponding to a pair of mass change rates of the NPI and the SPI), the RMS of the residuals (gravity rate vector observed by GRACE after application of all of the corrections minus the computed gravity rate vector) is calculated. The results are shown in Figure 5. The RMS levels exhibit an elliptical shape with an optimum fit in the center. A movement from the RMS minimum along the direction of the semiminor axis of the ellipses would change the sum of both ice mass change rates and lead to a fast degradation of the fit. It is therefore an indication of the stability of our estimation of the sum of both ice mass change rates. A movement from the minimum along the semimajor axis of the ellipses would keep the sum of both of the mass rates constant, but would change the relation between the rates of the SPI and NPI. Here, a larger range of pairs of mass rates shows a reasonable fit. This demonstrates the degree of separability of the SPI and NPI ice mass change rates.
We interpret the tendency of our optimization, along with computational experiments applying different alternative a-priori ice mass change rate patterns, presented in Table A1 in the Appendix A, as an indication that the ice mass loss rates derived from the remote sensing results [2,26,28] (Figure 5) underestimate the ice loss contribution from NPI. Indeed, for the Benito glacier at the Western Flank of the NPI, a combination of field and satellite-derived elevation data suggests a dramatic surface lowering, attributed to a negative surface mass balance, with a mean rate of −3.0 ± 0.2 m/a between 1973 and 2017 [61].
The time series of the monthly ice mass change in the SPI and NPI region (Figure 4, top) is dominated by a steady decrease and an annual cycle. We estimated the annual and semi-annual amplitudes at 54.6 Gt and 0.8 Gt, respectively. The ice mass time series smoothed by the running annual mean shows small inter-annual variations and is dominated by a linear decrease. The long-term ice mass change is very well approximated by the linear model, in particular, we found no indication for an acceleration of the ice mass loss in the SPI and NPI region throughout the GRACE operation period (Figure 4). In fact, when including, in addition, a quadratic term in the model fitted to the ice mass change time series, this acceleration term was statistically not different from zero, according to a t-test with a 95% significance level. On the other hand, minimum observation periods of 10 years and 20 years are necessary for a statistically significant determination of ice mass loss acceleration of the AIS and Greenland ice sheet, respectively [62]. This makes a statistically sound inference of an acceleration of ice loss of the smaller Patagonian Icefields over the 15 year long GRACE observation period rather unlikely.
Prominent features in the time series of the residual gravity misfit (Figure 4, bottom) can be primarily explained by the orbital resonance appearing in the Level-2 solution for the month of February 2015 (RMS of almost 100 nm s −2 ). A slight increase in misfit after March 2016 (on average 23.1 nm s −2 ) compared to previous (on average 16.2 nm s −2 ) months also becomes evident. This, together with the increase in the deviation from the linear plus periodic fit, gives a quantitative measure of the inevitable degradation of the Level-2 GRACE product emerging from the Level-1a and 1b data obtained during the final year of the GRACE mission. In our study, the data quality was interpreted to mean the ability to fit our mass-change pattern, and depends only on those combinations of spherical harmonic coefficients that dominate shaping the monthly gravity changes within our grid domain. Such mainly short-wavelength signals correspond to a great extent to coefficients of a higher degree and order. A degradation concentrated on low-degree coefficients leads to a long-wavelength bias, which is largely mitigated by our joint adjustment of the empirical ocean mass increment and the scaling factor for the ice mass change pattern with the full exploitation of the directional information of the 3D gravity vectors. A more detailed discussion of the GRACE Level-2 data degradation and its impact on the SPI and NPI ice mass estimation is given in Appendix A.
Remote Sens. 2019, 11, x FOR PEER REVIEW 11 of 20 coefficients of a higher degree and order. A degradation concentrated on low-degree coefficients leads to a long-wavelength bias, which is largely mitigated by our joint adjustment of the empirical ocean mass increment and the scaling factor for the ice mass change pattern with the full exploitation of the directional information of the 3D gravity vectors. A more detailed discussion of the GRACE Level-2 data degradation and its impact on the SPI and NPI ice mass estimation is given in Appendix A. Figure 5. Separation of the ice mass loss at SPI and NPI. After the application of all corrections to the pseudo-observations derived from the GRACE Level-2 data, the residual gravity root mean square (RMS) is computed for pairs of ice mass change rates for SPI and NPI that are varied over plausible ranges. The residual RMS is shown in colors. White circle: ice mass change rates according to our pattern optimization; black symbols: ice mass change rates according to previous works; triangle: Braun et al. [26]; inverted triangle: Jaber [28]; diamond: Foresta et al. [4]; square: Willis et al. [2]. Table 2 shows the impact of the individual corrections applied in the course of our GRACE data analysis on the ice mass change rate for the SPI and NPI region. The GIA correction has by far the largest impact, contributing about one third to the derived ice mass change rate. Further corrections significant for the ice mass change rate are those for the mass changes in the ocean and the AIS. The ice mass loss and elastic solid Earth response at SPI and NPI cause a local lowering of the sea level adjacent to southern Patagonia through a loss of gravitational attraction coincident with the ice mass losses on land (Figure 1a). If this effect is not accounted for, the negative gravity pattern over the nearby ocean is absorbed into the ice mass estimate, while the long-wavelength sea level change is accommodated by the empirical regional ocean mass change rate correction. This effect has also been Figure 5. Separation of the ice mass loss at SPI and NPI. After the application of all corrections to the pseudo-observations derived from the GRACE Level-2 data, the residual gravity root mean square (RMS) is computed for pairs of ice mass change rates for SPI and NPI that are varied over plausible ranges. The residual RMS is shown in colors. White circle: ice mass change rates according to our pattern optimization; black symbols: ice mass change rates according to previous works; triangle: Braun et al. [26]; inverted triangle: Jaber [28]; diamond: Foresta et al. [4]; square: Willis et al. [2]. Table 2 shows the impact of the individual corrections applied in the course of our GRACE data analysis on the ice mass change rate for the SPI and NPI region. The GIA correction has by far the largest impact, contributing about one third to the derived ice mass change rate. Further corrections significant for the ice mass change rate are those for the mass changes in the ocean and the AIS. The ice mass loss and elastic solid Earth response at SPI and NPI cause a local lowering of the sea level adjacent to southern Patagonia through a loss of gravitational attraction coincident with the ice mass losses on land (Figure 1a). If this effect is not accounted for, the negative gravity pattern over the nearby ocean is absorbed into the ice mass estimate, while the long-wavelength sea level change is accommodated by the empirical regional ocean mass change rate correction. This effect has also been discussed in the context of the Antarctic Peninsula ice loss [63]. The far-field effect of the ice mass loss of the AIS manifests itself in a decrease in the N-S directed gravity component, with increasing intensity towards the south (Figure 2). Our correction model for the continental hydrology produces a signature of opposite orientation (intensity increasing towards N), yet smaller magnitude. Glaciers in the vicinity of the SPI and NPI are included in the derived ice mass change rate. The ice mass change of the extra-Patagonian glaciers in South America [26] results in a slight decrease in the SPI and NPI ice mass change rate, dominated by the glacial systems further south (Cordillera Darwin, Gran Campo Nevado). Mass variations in the adjacent lakes have a minor influence on the ice mass change rate; but are included because of their impact on the seasonal ice mass signal. The small magnitude of the empirical regional ocean mass change rate increment testifies to a high consistency of our explicit ocean mass correction with GRACE gravity. The last column in Table 2 demonstrates that the omission of any of the significant surface mass corrections degrades the fit of the adjusted model to the observations. The application of an identical data correction and analysis to alternative GRACE Level-2 CSR-RL06 harmonics (Center for Space Research, University of Texas [64]) over January 2003 to August 2016 (maximum degree and order 96) yields an SPI and NPI ice mass change rate within 3% of that derived from the ITSG-Grace2016 harmonics for the same period.
Other GRACE derived ice mass change rates as well as those based on Cryosat-2 swath altimetry, are essentially in agreement with our results (Table 1). A more substantive disagreement beyond the stated error bars is found with loss rates derived from remote sensing, especially for the SPI. The rates determined over the acquisition time spans not overlapping with the GRACE data period [3,25] shall not be discussed here. Furthermore, while remote sensing techniques are spatially limited, with separate results published for the two icefields, the GRACE derived rates represent integrals over a broader region. Therefore, the GRACE derived ice mass change rates can be compared only with those remote sensing results, which consistently sum the contributions of both icefields [2,26,28]. Synthetic-aperture radar interferometry (TanDEM-X) [26][27][28], in fact, seems to yield systematically smaller ice mass losses than optical imagery (ASTER, SPOT) [2,23,24]. Our results support the larger trend in ice loss as derived by Willis et al. [2] from the ASTER imagery.

Conclusions
Our analysis of the complete GRACE satellite gravimetry record between April 2002 and June 2017 yields an SPI and NPI ice mass change rate of −24.4 ± 4.7 Gt/a. The time series of the monthly ice mass change does not indicate an acceleration in ice mass loss throughout the GRACE period. This ice mass change rate corresponds to an average contribution to the eustatic sea level rise of 0.067 ± 0.013 mm/a, or about 3% of the mass contribution to the sea level implied by our ocean mass correction model. This confirms the disproportionally large impact of the Patagonian Icefields on sea level rise. The almost constant ice mass loss derived from GRACE over 15 years is also interesting in the context of the changes of ice flow velocities in the same time. In light of the very heterogeneous evolution in the glacier flow velocities [29], our result of a stable mass loss rate suggests that the differences among the individual glaciers basins in the non-linear mass balance variations cancel out to a great extent. We found strong indications that the NPI has a larger relative contribution to the total ice mass loss than previous remote sensing results have suggested [2,23,26,28].
We also show that the Level-2 data of the latest months of the GRACE mission provide valuable information on Patagonian ice mass changes. The derived ice mass change rate depends sensitively on the corrections for GIA and non-uniform ocean mass changes. Our new estimate of the SPI and NPI ice mass change rate is a step towards an improvement of the recent ice-unloading history, providing a strong link to the GIA and crustal deformation in this tectonically complex region, and, possibly, to a better determination of the effective viscoelasticity beneath southernmost South America.
Sustained and unabated ice mass loss over the GRACE observing period is remarkable, for there seem to be climate change processes (surface mass balance and ice discharge) that play out over a long time scale, especially when we consider the losses that can be traced into the middle of the 20th Century [65]. This fact forebodes ill for the environment of southernmost South America in nearly every way. The loss of mountain glacier water resources will lead, inevitably, to irreversible changes in the habitability of the regional flora and fauna. This is the most important societal ramification that may be derived from the 15 year-long GRACE record in Patagonia. Acknowledgments: We thank Petra Döll for providing updated grids of the continental water storage of the WaterGAP 2.2c model. The lake-level time series of the Patagonian lakes was retrieved from the BDHI data base of the Argentine Subsecretaría de Recursos Hídricos. We thank the German Space Operations Center (GSOC) of the German Aerospace Center (DLR) for providing, continuously, nearly 100% of the raw telemetry data of the twin GRACE satellites.

Conflicts of Interest:
The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript; or in the decision to publish the results.

GIA Correction
The monthly GRACE gravity fields are corrected for GIA according to the regional model La-A [6]. This model has been shown to fit best the surface deformation observed by GNSS [7]. The GIA uplift models are given as regional grids and include the elastic response to present-day ice mass changes [6]. This elastic contribution has to be removed before deriving the GIA correction models for GRACE. We applied surface densities according to the a-priori ice mass change rate pattern [2], as a load model and Green's functions of an elastic Earth model [66]. The elastic contribution to the vertical deformation is calculated for the grid nodes of the GIA model, and is subtracted from the model uplift rates. This reduced GIA model is expanded in spherical harmonics (degree/order 360). The deformation coefficients are converted to Stokes coefficients by an approximation [37] (Equation 9), enforcing mass conservation through a global layer of uniform mass outside Patagonia. These GIA Stokes rate coefficients are truncated to degree/order 2 to 120, multiplied by the time since the initial GRACE epoch (April 2002), and subtracted from the monthly GRACE Stokes coefficients sets.
Lange et al. [6] present, in addition, a second regional GIA model, termed model La-B, which fits, similarly well, the GNSS uplift observations, and has been deemed more plausible by those authors. An alternative GIA correction was derived from this GIA model La-B by the same procedure in order to quantify the impact of the uncertainty of the GIA models on the derived ice mass change (Table 2).

Surface Mass Variation Corrections
The monthly gravity vector grids are corrected for surface mass-variations of the AIS, the global ocean, continental water storage, and the great Patagonian lakes. In order to remove the far-field effect of ice mass change in Antarctica, the GRACE-derived time series of gridded Antarctic ice mass changes of the AIS_cci Gravimetric Mass Balance Product [52], generated in the scope of the European Space Agency's Climate Change Initiative project for the Antarctic Ice Sheet (AIS_cci), are applied. The latest months of the GRACE data period (November 2016 to June 2017) are not included in this product. These time series are therefore extrapolated for each AIS_cci grid cell, applying a linear plus periodic (periods: annual, semi-annual, and 161 days) model.
Our ocean mass correction consists in an updated solution of the gravitationally self-consistent distribution over the ocean of the water mass exchange with five major ice sheets/ice caps and the global continental hydrosphere [67] (following the pseudo-spectral approach [56]). The monthly ice mass change patterns of the AIS [52], Greenland ice sheet [52], and the Patagonian Icefields (this study, iterative solution), complemented by the ice mass rate patterns for Alaska and the Canadian Arctic [57], as well as the monthly patterns of the continental water storage according to the WaterGAP 2.2c model [53], are used as an input for the solution of the sea level equation. Because of the limited extent of the forcing models [52,53], the monthly ocean water mass distributions were computed through November 2016, and then extrapolated over the last GRACE mission months applying a linear plus periodic (annual, semi-annual) model. On a global average, this ocean mass model implies a relative sea level rise of 2.13 mm/a.
The effects of the changes in the continental water storage are removed by using total water storage time series on a global grid of the WaterGAP 2.2c model [53] with climate forcing WFDEI-CRU [54]. As these time series are available only until December 2015, they were extrapolated through June 2017, applying a linear plus periodic (annual and semi-annual) model. We took particular care of the effects of the hydrological mass variations associated with lake-level changes in Lago Buenos Aires/General Carrera (BAGC), Lago San Martín/O'Higgins, Lago Viedma, and Lago Argentino (Figure 1b). These large lakes receive the meltwater runoff from the adjacent NPI (BAGC) and SPI [55]. Brazo Rico and Brazo Sur of Lago Argentino were treated as an additional individual lake, because of the repeated damming by the Perito Moreno glacier during the GRACE data period. On short and seasonal time scales, we expect these lakes to produce a mass signal opposite to that of the icefields, but too close to be resolved separately by GRACE. We used tide gauge data in the lakes provided by the BDHI data base to derive the time series of the monthly mean lake levels. Before 2008, regular tide-gauge data are available only for Lago Argentino and Brazo Rico/Sur, the time series of the other three lakes were extrapolated backwards in time until April 2002, by stacking the annual lake-level signal. The mass signal of the lake-level changes was subtracted in the corresponding nodes of the global hydrological model grid prior to the application of the corrections to the GRACE data.
All of these surface mass-variation processes are represented by monthly time series of the surface density on global (ocean, hydrology) or regional (AIS, lakes) grids. For each of these processes and for each month, the corresponding 3D gravity vector components are derived for our regional grid and subtracted from the grids of the pseudo-observations.

Empirical Regional Ocean Mass Rate Correction
Our explicit ocean mass correction accounts for the mass exchange with five major glaciated regions and continental water storage. Additional sources may contribute to ocean mass variations in the region under investigation. Therefore, we apply an empirical, additional correction in order to accommodate the residual imperfections of our explicit ocean mass model. A set of point masses homogeneously distributed over the map extension in Figure 1a, spaced every 0.1 • in longitude and in latitude increments that provide for an equal area between the point masses, is used for this purpose. Over all of the point masses falling within the ocean, the residual linear mass changes are modeled by three parameters describing a mean mass change rate and mass-rate gradients in a north and east direction, respectively. These three parameters are solved, simultaneously with the monthly scaling factors for the ice mass change pattern, in a joint least-squares adjustment to the monthly gravity vector grids.
In fact, this empiric correction would accommodate long-wavelength effects not only originating from regional deviations from our explicit ocean mass correction, but also due to residual imperfections in other correction models or low-degree GRACE data. The empirical regional ocean mass rate increment produces only a small change in the ice mass change rate (1%, see Table 2), indicating a high consistency between the applied correction models and GRACE data.

A-Priori Mass Change Pattern
An a-priori SPI and NPI ice mass change pattern (Wi12) is derived based on published remote sensing results [2]. We use ice surface elevation rates dh/dt, determined by differencing ASTER digital elevation models (DEM) from 2012 (SPI) and 2011 (NPI), relative to a DEM derived from the SRTM in February 2000 [2]. The dh/dt are given as the average values for accumulation and ablation zones of 87 (SPI: 49; NPI: 38) glacier basins. We applied an ice density of 0.9 g cm −3 to convert the dh/dt to mass change rates. Glacier boundary polygons from the GLIMS data base [68] were used together with equilibrium line altitudes [2,23] and a SRTM derived DEM [69] in order to compute the barycenter coordinates for the 174 accumulation/ablation zones. The ice mass change rates are treated as point masses located in the barycenter of the corresponding accumulation/ablation zone (Figure 1b). This preliminary ice mass change rate pattern is then optimized by searching separate scaling factors for SPI respectively NPI that minimize the residual gravity RMS in our GRACE data adjustment.
Our approach is largely insensitive to the choice of the a-priori mass pattern, because the upward continuation to orbit altitude acts as a low-pass filter in the space domain, damping the resolution and detail of the surface mass configuration. This demonstrated computations employing the alternative a-priori patterns Ja16, Br19, and Fo18 (Table A1). All three of these alternative a-priori patterns consist of only two point masses each, in which the ice mass change rates derived from the remote sensing results [26,28] or swath altimetry [4] are concentrated. Although these models represent only 73%, 75%, and 91%, respectively, of the Wi12 total ice mass change rate, its impact on the adjusted SPI and NPI ice-mass rate does not reach 1.5% (Table A1). However, the application of either of the alternative patterns reduces the residuals, despite their coarser resolution. The systematic decrease in the residual gravity RMS among the different a-priori patterns with a decreasing SPI/NPI ice mass rate ratio is interpreted as an indication that the Wi12 model underestimates the ice-mass loss at the NPI relative to that at the SPI, and motivates our optimization of the Wi12 pattern by two separate scaling factors for SPI and NPI.

Impact of GRACE Level-2 Data Quality Degradation on SPI and NPI Ice Mass Determination
In order to assess the global temporal behavior of the GRACE data accuracy, we computed the residual equivalent water height anomalies for each spherical harmonic degree and order with respect to a model that includes linear, quadratic, and periodic terms. The standard deviation of these temporal residuals is interpreted as a measure of GRACE observation uncertainty. Figure A1 shows the degree amplitudes of these equivalent water height anomalies for the individual monthly Level-2 solutions. As expected, the maximum uncertainty is found at the highest degrees, while minimum standard deviations are found between degrees 10-30.
Among these degree amplitude curves, two outliers are observed, corresponding to the months of February 2015 (two-day repeat orbit; red) and March 2017. In addition, the deterioration of the data accuracy after March 2016 becomes evident, with degree amplitudes above the median (among all months of the entire GRACE data period), and substantially increased energy on the low degrees (<10).
In order to investigate the sensitivity of the SPI and NPI ice mass time series to this GRACE uncertainty evolution, the SPI and NPI region function (as applied in the classical regional integration approach) is expanded in the spherical harmonic coefficients. Unlike AIS and the Greenland ice sheet, which concentrate much of their signal on low degrees (<10), the small SPI and NPI region requiresa substantially broader band of degrees. Therefore, the uncertainty of those higher degrees is also more important. We derive degree amplitudes of the region functions for SPI and NPI, AIS, and the Greenland ice sheet, which are included in Figure A1. In general, only the portion of the error spectrum retained after filtering is integrated with the region function. In the present case, the filtering corresponds to the effect of using pseudo-observations at the orbit altitude. As the SPI and NPI region function distributes the energy over higher degrees, these are also less damped by the filtering (compared to AIS and Greenland). A comparison of the degree amplitudes of the observation uncertainty and the regional function reveals that the SPI and NPI region function has its strongest energy on the degrees for which the errors are smallest-degrees 10 to 30. Because of the small energy on the low degrees (<10), the accuracy degradation of the solutions after March 2016 has not had as efficient an effect on the SPI and NPI ice mass estimation, as, for example, for AIS or Greenland. Comparing the two outlier months reveals that February 2015 (repeat orbit) results in much larger residuals (green time series in Figure 4) than March 2017. This demonstrates that for SPI and NPI, the uncertainty on the higher degrees is much more important than the long-wavelength errors (unlike AIS or Greenland), but these higher degrees are not as severely affected by the accuracy degradation.
We conclude that the GRACE Level-2 data accuracy degradation affects mainly low degrees (<50), but the small SPI and NPI target is less sensitive to this low-degree degradation. Thus, also the latest months of the GRACE mission hold useful information about the Patagonian ice mass changes. In its exploitation, we took the data quality degradation into account by a cautious empirical downweighTing. An exclusion of the solutions affected by the data quality degradation has little impact on our final result-our time series over the period of April 2002 through March 2016 yields an ice mass change rate for the SPI and NPI region of −23.9 Gt/a.
A strategy is now in place at the Jet Propulsion Laboratory to model each of the thrusting events on the satellite lacking accelerometer data [70]. We therefore anticipate improvements in the Level-2 harmonics upon the next release that contains such modeling, and it will be interesting to further evaluate our scheme used here for the analysis of the late months of the GRACE mission. Greenland ice sheet, which concentrate much of their signal on low degrees (<10), the small SPI and NPI region requiresa substantially broader band of degrees. Therefore, the uncertainty of those higher degrees is also more important. We derive degree amplitudes of the region functions for SPI and NPI, AIS, and the Greenland ice sheet, which are included in Figure A1. In general, only the portion of the error spectrum retained after filtering is integrated with the region function. In the present case, the filtering corresponds to the effect of using pseudo-observations at the orbit altitude.
As the SPI and NPI region function distributes the energy over higher degrees, these are also less damped by the filtering (compared to AIS and Greenland). A comparison of the degree amplitudes of the observation uncertainty and the regional function reveals that the SPI and NPI region function has its strongest energy on the degrees for which the errors are smallest-degrees 10 to 30. Because of the small energy on the low degrees (<10), the accuracy degradation of the solutions after March 2016 has not had as efficient an effect on the SPI and NPI ice mass estimation, as, for example, for AIS or Greenland. Comparing the two outlier months reveals that February 2015 (repeat orbit) results in much larger residuals (green time series in Figure 4) than March 2017. This demonstrates that for SPI and NPI, the uncertainty on the higher degrees is much more important than the long-wavelength errors (unlike AIS or Greenland), but these higher degrees are not as severely affected by the accuracy degradation. We conclude that the GRACE Level-2 data accuracy degradation affects mainly low degrees (<50), but the small SPI and NPI target is less sensitive to this low-degree degradation. Thus, also the latest months of the GRACE mission hold useful information about the Patagonian ice mass changes. In its exploitation, we took the data quality degradation into account by a cautious empirical downweighTing. An exclusion of the solutions affected by the data quality degradation has little impact on our final result-our time series over the period of April 2002 through March 2016 yields an ice mass change rate for the SPI and NPI region of −23.9 Gt/a.
A strategy is now in place at the Jet Propulsion Laboratory to model each of the thrusting events on the satellite lacking accelerometer data [70]. We therefore anticipate improvements in the Level-2 harmonics upon the next release that contains such modeling, and it will be interesting to further evaluate our scheme used here for the analysis of the late months of the GRACE mission.