# Lake and reservoir volume variations in South America from radar altimetry, ICESat laser altimetry, and GRACE time-variable gravity

Claudia Carabajal, Jean-Paul Boy

► To cite this version:

Claudia Carabajal, Jean-Paul Boy. Lake and reservoir volume variations in South America from radar altimetry, ICESat laser altimetry, and GRACE time-variable gravity. Advances in Space Research, 2020, 10.1016/j.asr.2020.04.022. hal-03037810

HAL Id: hal-03037810

https://hal.science/hal-03037810v1

Submitted on 5 Aug 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Lake and Reservoir Volume Variations in South America from Radar Altimetry, ICESat Laser Altimetry, and GRACE time-variable gravity

Claudia C. Carabajal (1) and Jean-Paul Boy (2)

(1) HEXAGON US Federal, Inc. @ NASA Goddard Space Flight Center, Geodesy and Geophysics Laboratory (Code 61A), Greenbelt, MD 20771, USA. (claudia.c.carabajal@nasa.gov)

(2) EOST/IPGS (UMR 7516 CNRS – Université de Strasbourg), 5 rue René Descartes, 67084 Strasbourg Cedex, France (jeanpaul.boy@unistra.fr)

## Abstract

Thanks to 25-years of radar altimetry measurements from TOPEX/POSEIDON, Jason 1, 2 and 3, ENVISAT and others, water level changes of major lakes and rivers can be measured regularly, remotely, with unprecedented precision, facilitating monitoring of continental water storage variations. In addition, the 18 laser altimeter measurement campaigns over 6 years by the Ice, Cloud and Land Elevation Satellite (ICESat) satellite complement these radar altimetry measurements and also validate them, expanding our monitoring capabilities. We extract derived water height estimates for various lakes in South America from radar and laser altimetry measurements. Using available imagery from the Moderate Resolution Imaging Spectroradiometer (MODIS), we estimate time-variable surface extent of these water bodies, and together with the altimetry, we derive corresponding estimates of relative volume variations for a number of lakes and reservoirs in South America. We see good agreement between our water extent and available products, with some differences when snow covered area are present. We find good agreement between the temporal variations of lake water heights obtained from ICESat, radar altimetry, and in situ measurements, when available. Lake-dependent biases still exist, likely due to differences in reference datum, choice of geoid model, and differences in data processing schemes. The is a good agreement in volume changes obtained from imagery and the various altimetry techniques. Furthermore, using the latest one-degree global iterated mascon solutions from NASA Goddard Space Flight Center, we successfully compare them to independent volume
changes assessments for these water bodies, derived from mass variations from time-
variable gravity measurements from 15 years of Gravity Recovery And Climate Experiment
(GRACE) observations. The agreement in water volume changes from GRACE with those
from altimetry and imagery is higher for lakes in arid and semi-arid regions than on
mountainous areas. GRACE volume estimates remain limited to the spatial resolution and
the quality of the hydrology models used to correct for soil-moisture and snow
contributions. Overall, studies like the one presented here further validate the power of
leveraging these data synergies to monitor continental water storage changes from space,
taking advantage of what each of the techniques has to offer. We show that water level
measurements from laser and radar altimeters have similar precision, of a few centimeters
over inland water bodies, and prove that GRACE global mascons allow to retrieve small lake
relative volume changes, even for small bodies, thanks to their higher sensitivity to small
wavelength mass variations.

Key words

Radar altimetry, Laser altimetry, time-variable gravity, inland water, ICESat, GRACE, MODIS.

1. Introduction

Surface freshwater, and in particular lakes and reservoirs, are an important component of the global water cycle (Huang, 2018). They can be monitored from space using imaging remote sensing (surface extent) and altimetry (water height). However, nadir altimeters on repeated orbits are missing a large number of targets due to their ground track spacing. Surface water extents are usually monitored using data in the visible and infrared bands (Tulbure et al., 2016; Huang et al., 2018). In the future, missions with this scientific goal, like the SWOT (Surface Water Ocean Topography) mission, optimized for these kind of applications, aim to observe almost all freshwater bodies thanks to its wide-swath (see, for example, Lee et al., 2010).

In this paper, we focus on lakes and reservoirs in South America. We chose this region
because the lakes are rather small compared to, for example, the Great Lakes in Northern
America or the East-African rift lakes (Moore and Williams, 2014), and their height variations are also small, posing more of a challenge for monitoring techniques. Their volume variations are also large enough to be derived from time-variable gravity measurements from the GRACE (Gravity Recovery And Climate Experiment) mission (Tapley et al., 2004). They are located in significantly different climatic and geographic areas, ranging from glacial lakes in the Andes to those in more temperate regions like Lake Chiquita. Our main goals are to increase the observations of water height measurements by adding laser altimetry to classical radar altimetry measurements with comparable accuracy, to show the feasibility of computing relative lake mass variations from time-variable gravity field measurements, and to compare relative volume variations derived from altimetry and imagery to time-variable gravity field estimates.

Thanks to nearly three decades of radar altimetry missions, radar altimeters have been very successful at monitoring inland water bodies, having demonstrated the ability to remotely detect water surface level changes in lakes and inland seas in sub-monthly, seasonal and inter-annual time scales with a precision up to a few centimeters (see, for example, Crétaux et al., 2011). Some of the advantages of this technique include the ability to observe under any weather conditions, night and day, and in various environments, including vegetated terrain and wetlands. Long term systematic observations along the repeat orbits follow closely a set of reference tracks, therefore facilitating the continuous monitoring of a particular water body. A fundamental difference between altimeter retrieved heights and traditional gauge measurements at a given location is that they represent an average footprint height, in the direction of the satellite motion, along a relatively narrow path (swath). These can represent height estimates for a relatively large area along track, from several meters to several kilometers. For all altimetry missions, the orbit configuration is a trade-off between the spatial and temporal coverage. The footprint size is determined by the instrument design, and the rate of the telemetry, in addition to the target's characteristics. Its geolocation accuracy is a function of the satellite orbit position and attitude knowledge, the knowledge of the altimetric range, including all geophysical corrections, and the bounce point (footprint) characteristics.
In addition to the classical radar altimeters, the available laser altimetry from ICESat (Ice, Cloud and Land Elevation Satellite) provides additional and complementary measurements of surface topography not only on the polar regions, but nearly globally, including water level measurements of lakes and reservoirs. The ICESat mission operated between 2003 and 2009 using the GLAS (Geoscience Laser Altimetry System). With the primary goal was to map the elevations of the polar regions to quantify ice sheet changes (Zwally et al., 2002 and 2011; Schutz et al., 2005), ICESat, operated in a primarily nadir-viewing orientation, day and night, along a series of reference ground tracks with coverage extending to ±86° latitude.

The Gravity Recovery And Climate Experiment (GRACE) mission, launched in April 2002 and followed by GRACE Follow-On in May 2018, has been very successful in tracking Earth's surface mass variations (Tapley et al., 2004) with spatial and temporal resolutions of typically a few hundreds of kilometers and one month, respectively. Although GRACE alone cannot separate contributing components to total water storage changes, surface water (lakes, reservoirs and even rivers) can be a significant contributor to the understanding of continental water storage variations (see, for example, Han et al., 2010; Longuevergne et al., 2013; Moore and Williams, 2014).

In the following sections, we describe the steps followed in this investigation. We present the study area in Section 2. In Section 3, we describe the processing of the MODIS NDVI (Normalized Difference Vegetation Index) data used for the retrieval of the surface extent of the various lakes and reservoirs in South America. In Section 4, we describe the analysis of the ICESat laser altimetry measurements, and compare them to the more classical radar altimetry products in Section 5. Section 6 is devoted to the analysis of the GRACE global iterated mascons (mass concentrations), and the retrieval of the lake mass variations. In Section 7, we compute the hypsometry for each individual body, i.e. the relation between height and surface. Comparisons between mass variations from remote sensing altimetry and imagery and the GRACE-derived estimates are provided in Section 8. Section 9 includes some concluding remarks and perspectives.
## 2. Selected study area

We select South American lakes in Bolivia and Peru (lakes Titicaca and Poopo), Chile (lakes General Carrera, Llanquihe and Todos Los Santos), Argentina (lakes Argentino, Chiquita, Colhué Huapi, El Chocón, Musters and Viedma) and Uruguay (Negro). Lake Titicaca is the largest lake in South America in terms of water volume (about 900 km³), located in the Altiplano at the border of Peru and Bolivia (Pillco et al., 2019). Its water inflow is mainly provided by 5 major rivers, and about 90% of the outflow is due to evaporation, the rest flows to lake Poopo through the Desaguadero River (Pillco and Bengtsson, 2006). Lake Poopo is a shallow lake with a surface varying greatly on seasonal basis, and almost disappeared in 2015 (Satgé et al., 2017). Lake Chiquita is located in the western part of the Parana-La Plata basin, in Central Argentina, in a subtropical and semi-arid region (Bucher and Stein, 2016). It is a closed, shallow, hypersaline-lake, with large changes in its water level over several decades, and it is considered a good regional and temporal water balance integrator of its catchment (Troin et al., 2010). Lakes Musters and Colhué Huapí are respectively deep and shallow lakes located in semi-arid central Patagonia, Argentina, and form the terminal stage of the Senguerr River endorheic basin (Rivas, 1978). The precipitation over the lakes reaches approximately 100 to 200 mm/year, mostly during winter months. Most of the outflow is due to evaporation (more than 1800 mm/year) (Scordo et al., 2018). Lakes General Carrera (also called Buenos Aires), Llanquihe, Todos Los Santos (Chile), Argentino and Viedma (Argentina) are deep lakes located in the Patagonia Andes, and are mainly fed by glacial melt water (Pasquini et al., 2008; Depetris and Pasquini, 2000). Their spatial extent varies from about 180 km² for Lake Todos Los Santos to almost 1900 km² for Lake General Carrera. Most of them experience a general growth in total water storage (Wilson et al., 2018) due to present-day ice melting of the Patagonian glaciers (Forestá et al., 2018). El Chocón and Negro are two reservoirs, and their corresponding surface areas are between 800 and 1000 km², respectively. El Chocón Dam is the largest hydroelectric plant in the Argentine Patagonia, built at the end of the 1970s along the Rio Negro, in the Neuquén province. The Rio Negro Reservoir (also called Rincón del Bonete or Gabriel Terra) is the largest fresh water body in Uruguay, built in the 1940s on a different river, also called Rio Negro, a tributary of the Uruguay River.
### 3. Water extent extraction from MODIS NDVI

MODIS (Moderate-Resolution Imaging Spectroradiometer) onboard Terra and Aqua satellites is imaging the Earth surface in 36 spectral bands, including infra-red and visible and has been widely used to derive inland water body spatial extent (see, for example, Buma et al., 2018; Omute et al., 2012). In this study, we use not only the NDVI (Normalized difference vegetation index) provided in the MOD13Q1 vegetation latest products (version 006), available with temporal and spatial resolutions of 16-days and 250-meters since March 2000 (Justice et al., 2002; Didan, 2015), but also the reflectance in band 1 (red), 2 (Near Infrared) and 3 (blue).

As a first approximation, negative values of NDVI identify water. NDVI values between 0 and 0.3 define dry-land, while values larger than 0.3 identify vegetation. However, we cannot use a simple threshold of NDVI to define water, as snow covered areas have also a low or negative NDVI. Only reflectances in the visible and NIR bands can separate snow (high reflectance) from water (lower reflectance). Therefore, we have defined as water cells, cells with NDVI values lower than 0.1, and red, blue and NIR reflectances lower than 0.2. An example for Lake Argentino (Argentina) is provided on Figure 1.

<figure><img src="image_2.png" /></figure>
Figure 1: NDVI (left) around lake Argentino for 2001/05/25 and corresponding middle infrared (MIR) (band 7), near infrared (NIR) (band 2), red (band 1) and blue (band 4) reflectances.

Figure 2 shows the full time series of surface extent from 2000 to 2018 for Lake Chiquita (Argentina), its comparison with Crétaux et al. (2011) estimates, and the spatial NDVI variations for the low (February 2014) and high (May 2003) extremes (shown as red dots on the upper curve).

<figure><img src="image_3.png" /></figure>

Figure 2: Surface area variations (km²) of lake Chiquita (Argentina): in black our estimates using MODIS NDVI products, in green Crétaux et al. (2011) results. Lowest (2014/02/02) and highest (2003/05/25) surface extents are shown with a red dot. These times correspond to the NDVI maps plotted at the bottom, which illustrate the extreme surface extents determined by our technique.
Our water extent estimates agree well with Crétaux et al. (2011, 2016) series for the common available period (2000-2010), with a bias of about 200 km². Khandelwal and Kumar (2019) also provide estimates of surface water area for various lakes and reservoirs worldwide, based on the Khandelwal et al. (2017) algorithm. However, their estimates for lake Chiquita have a significant bias of about 3000 km² compared to Crétaux et al. (2011, 2016) and our results.

The comparison between the different surface extent time series, e.g. Crétaux et al. (2011, 2016), Khandelwal and Kumar (2019) and our estimates is shown on Figure 3 for lake Poopo in Bolivia. Our time series agree well with Khandelwal et al. (2017), also based on MODIS data, but not with the Hydroweb time series. Our methodology, using a simple threshold in NDVI and red, blue and NIR reflectance values, is also validated for high altitude lakes, when land is partially covered with snow in winter.
<figure><img src="image_4.png" /><figcaption>Figure 3: Surface area variations (km²) of lake Poopo (Bolivia): in black our estimates using MODIS NDVI products, in green Crétaux et al. (2011, 2016) results, and in red Khandelwal et al. (2017) estimates (top). Surface water extent from Khandelwal et al. (2017) (middle) and corresponding NDVI maps (bottom) for two extremes (2001/03/06 and 2010/12/19); they correspond to the red dots in the time series.</figcaption></figure>
Figure 4 shows the comparison of water extent variations for 6 lakes in South America computed from MODIS/NDVI with Khandelwal and Kumar (2019) and Crétaux et al. (2011) products. In particular, it shows the underestimation of Lake Chiquita's water extent from Khandelwal and Kumar (2019) products compared to our time series and Hydroweb data (Crétaux et al., 2011). Except for Lake Titicaca, the linear fit between our estimate and the two other products is close to 1, with some small biases for Lake Viedma and Negro. For Lake Titicaca, we slightly overestimate the surface water extent during the winter season, due to the low NDVI of snow covered areas.

<figure><img src="image_5.png" /></figure>

Figure 4: Comparison of our surface water extent estimates using MODIS/NDVI (X-axis) with other estimates (Y-axis) for lakes Chiquita, Colhue Huapi, Viedma, Titicaca, Negro and Poopo. Relationship with Khandelwal and Kumar (2019) products is shown in green, and with Crétaux et al. (2011) products in red.
Among other sources, differences in extent values can be due to data sources and mask definitions (250 m MODIS/NDVI for our estimates, 500 m MODIS for Khandelwal and Kumar (2019), and Landsat for Crétaux et al. (2016)), and the ability to detect partially-filled pixels at the coastlines.

4. Water level estimates from laser altimetry using ICESat data

Laser altimetry (LIDAR, Light Detection and Ranging) relies on the measurement of the two-
way range distance from the satellite to the surface as derived from the timing of the
transmitted/returned optical laser energy distribution (laser transmit echo pulse/laser
return waveform), where the laser beam intercepted the surface (Harding and Carabajal,
2005). ICESat collected valuable altimetric measurements of surface topography in ice and
non-ice covered regions and vegetated terrain, water height measurements in land, and
ocean altimetry. The 3 lasers operated consecutively at infrared (1064 nm) and visible
(532nm) wavelengths (cloud lidar), pulsing at 40 Hz (~70 m footprints spaced ~170 m apart).
ICESat recorded the transmit laser pulse, and the infrared channel laser return waveforms.
During the first few weeks of operations, ICESat transitioned from an 8-day repeat check-out
period to its 91-day repeat mapping orbit (33-day near-repeat sub-cycle). Due to issues
experienced with the first laser, the mission operated mostly in campaign mode, collecting
data for nearly a month 3 times per year, and two times per year closer to the end of the
mission, during 18 campaigns (see Table 1 in Borsa et al., 2019).

Infrared laser measurements do not suffer from penetration into the surface of water, snow or ice. As other optical systems, ICESat was limited by cloud cover, which can lead to forward scattering effects that translate into delays in the range measurement. Corrections via the use of cloud flags, and derived from cloud products when such conditions where developed (Spinhirne et al., 2005) and were provided with the ICESat products. Slope broadening effects due to off-nadir viewing observations could also affect the accuracy of the range measurements. Saturation effects occur when the peak power of the returned signal from snow and ice surfaces is found to span a wider dynamic range than planned, often exceeding the linear dynamic range of the 1064nm detector assembly (Abshire et al., 2005).
In this work we have used an ICESat dataset developed as topographic Ground Control Points (GCPs) with sub-decimeter vertical accuracy and better than 10 m horizontal accuracy (Carabajal et al. 2010, 2011 and Carabajal and Boy, 2016), for which we apply a stringent editing criteria to yield a high quality global topographic control for land applications. GCPs do not include data contaminated by effects that result in degradation in the quality of the elevation products, such as return saturation, cloud contamination, or low energy returns. Laser returns that satisfy high accuracy elevation requirements have waveforms that show characteristics shared by waveform returns from water. The data processed was Release 33 (Zwally et al., 2011), current at the time of the last version of the GCPs processing. Later releases of the data did not include significant modifications to the elevations provided in the GLA14 products used for this analysis, where elevations were derived from ranges to the centroid of the laser waveforms. Inter-period laser biases corrections were later applied to the data for this work. We selected available returns from the GCPs for the water bodies studied as delineated by the water extent masks derived.

Water returns, when saturation is not present, are generally narrow and closely resemble the transmit laser pulse. For ICESat, the transmit pulse was approximately 7 ns wide, about 1 m wide at full width at half the peak amplitude. Significant waveform width implies there is a significant height range encountered within the laser footprint due to relief (slope and/or roughness of the ground) and/or vegetation cover, or atmospheric forward scattering due to clouds effects. As one of the first discriminators, waveforms with widths from signal start to end of less than 5 m are selected, indicating that within-footprint relief is low, and data is from a relatively clear atmosphere. We also exclude waveforms with pulses narrower than the transmit pulse (widths > 0.5 m), and apply editing based on laser beam off-pointing and instrumental parameters, working with mostly nadir-looking data, with negligible saturation. Saturation broadens, distorts and truncates the received waveforms when the return energy exceeds the receiver's or the digitizer's dynamic range. We only keep data with saturation index (the number of waveform bins with an amplitude greater than the saturation threshold) of two or less (slightly saturated), and apply empirically derived corrections when available (Sun et al., 2016). To eliminate clouds, in addition to eliminating complex and severely attenuated waveform returns, we applied a cloud-clearing procedure that identifies
the outliers looking at the mean and standard deviations with Digital Elevation Models (DEMs) like SRTM (Shuttle Radar Topography Mission, Farr et al., 2007). We also apply editing based on instrumental parameters, not including low received energy returns (waveform amplitudes greater less than 0.15 Volts). We minimize data that may be affected by waveform broadening due to incidence angle between the laser beam vector and the normal to the surface slope by excluding data from off-nadir angles greater than 1 degree. Rigorous analysis has shown that for low relief locations the ICESat data meet the accuracy requirements of 6 m horizontal and 10 cm vertical (see Table 1 of Carabajal et al., 2011, and Borsa et al., 2019), corroborated by the mission using long-arc accuracy estimates from rigorous analysis of instrument calibration and validation schemes using ocean scan maneuvers and cross-overs.

The characteristics of each operating laser introduced small inter-period biases, correlated with laser transmit energies as they decayed throughout the life of the mission. We identified any mean biases per ICESat observation period separately, and applied the correction to the observations (Carabajal and Boy, 2016).

The altimetric heights are with respect to the TOPEX/Poseidon reference ellipsoid. For this work, the ICESat ellipsoidal heights are then converted into orthometric heights using EGM2008 geoid (Pavlis et al., 2012). Starting with the data provided in this database, we then extract valid shots for each track crossing all of the selected lakes and reservoirs (see Figure 5 for the example of Lake Chiquita for the L3G period). ICESat's longitudinal coverage at these latitudes during an approximately 30-day campaign is shown. Gaps in the data show the result of editing for clouds and saturation effects.
<figure><img src="image_6.png" /></figure>

Figure 5: ICESat tracks 55, 174, 345 and 1327 for the L3G operational period (Oct.-Nov. 2006) over Lake Chiquita, in Argentina (left), overlaid on the MODIS NDVI image for 2006/11/01. ICESat profiles along those tracks are shown on the right with different colors.

<figure><img src="image_7.png" /></figure>

Figure 6: ICESat laser tracks 137 and 323 crossing El Chocón reservoir, overlaid on the NDVI map for 2008/12/02 (left) and corresponding vertical profiles (right top). Individual ICESat periods are represented by a different color (L2B: 2004/02-03, L2C: 2004/05-06, L3A: 2004/10-11, L3B: 2005/02-03, L3D: 2005/10-11, L3G: 2006/10-11, L3I: 2007/10-11, L3J: 2008/02-03 and L2D: 2008/11-12). Mean and standard deviations for water measurements from all ICESat crossings are shown on the bottom right.
Figure 6 shows a selection of the available ICESat observations for El Chocón reservoir before editing for water returns. Elevation profiles for every observation period are also shown. The vertical axis was chosen to highlight the elevation measurements on the water (a few outliers remain). The consistency of the measurements along track, and the clear temporal variability of El Chocón water levels along each track for each observation period is evident, shown in the corresponding time series.

5. Comparison of ICESat water levels with radar altimetry available products

As illustrated in the previous examples, laser altimetry coverage after editing varies geographically and seasonally. For a particular water body, GCP data for each track is then re-edited based on the distribution of elevations over the largest water extent defined using the MODIS NDVI. We exclude any outliers using a 2-sigma rejection criteria. Figure 7 shows, as an example, the water height variations deduced from ICESat data and their comparison with radar altimetry products from Hydroweb (Crétaux et al., 2011), G-REALM (Birkett et al., 2011), and DAHITI (Schwatke et al., 2015).
<figure><img src="image_8.png" /><figcaption>Figure 7: Variation of Lake Chiquita water height means and standard deviations derived from ICESat (black) and radar altimetry publicly available series from Hydroweb (red), G-REALM (green), DAHITI (blue).</figcaption></figure>

For the common period (2003-2010), the 3 radar altimetry products agree well with our ICESat estimates, with a bias of about 1 meter. Prior to 2002, we observe a unique discrepancy between radar products, as DAHITI time series significantly differs from G-REALM and Hydroweb data. According to Bucher and Stein (2016), the decrease of water level seen in G-REALM and Hydroweb is also observed with in situ measurements.

Crétaux et al. (2011) estimate the water height along all crossing tracks, correcting for any geoid or mean lake surface slope. The final time series is computed by combining all altimetric data on a monthly basis. Birkett et al. (2011) estimate water level series, from a segment of a single altimeter track over the target, and the time series are then composed of data from consecutive altimeter missions measured along the same ground track, without attempting to combine contemporaneous missions. Therefore, G-REALM data are produced
only using Topex/Poseidon, and its successors Jason-1, -2 and -3. The most recently available database is the Database for Hydrological Time Series over Inland Waters (DAHITI, http://dahiti.dgfi.tum.de) launched in 2013. It includes estimates of water level time series for rivers, lakes, reservoirs and wetlands. It applies an extended outlier rejection methodology, using a Kalman filter approach, combining a variety of altimeter missions, implementing an optional waveform retracking, with full error propagation estimates for each measurement (Schwatke, et al. 2015).

Among other sources of discrepancies, the observed biases between radar altimetry products are partly due to the choice of the vertical datum, including different choice of geoid, and the geophysical corrections applied. In the Hydroweb products, the slope of the geoid is adjusted for each track from an a priori model, which does not fit some of short wavelength variations over the lake surface (Crétaux et al., 2011). The G-REALM products are not provided in a geoid-based datum, but remain in the reference ellipsoid system, with a vertical datum specific to each lake based on a mean water level computed from the first nine years of Topex/Poseidon altimeter (Ricko et al., 2012). In contrast to Hydroweb products, the geoid model is assumed to be errorless in DAHITI products, so no along-track geoid slope is adjusted when correcting for the geoid (Schwatke et al., 2015).

Water level estimates from radar altimetry are classically validated with comparisons with in situ gauge measurements when available. One challenge in this validation process is often the lack of knowledge about the reference vertical datum of the gauges, leading to height offsets between altimetry and in situ measurements. Therefore, only relative comparisons are performed in most cases. Crétaux and Birkett (2006) estimated RMS differences between in situ and radar altimetry observations for lakes Chiquita and Titicaca of 13 cm and 17 cm, respectively. Schwatke et al. (2015) provided RMS differences between radar altimetry and in situ measurements for lake Argentino of 14.6 cm (DAHITI) and 21.9 cm (Hydroweb), and lake General Carrera of 19.0 cm (DAHITI) and 29.4 cm (Hydroweb). In addition to those presented in Schwatke et al. (2015), more detailed assessments of the accuracy of lake and reservoir water level measurements from radar altimetry products can be found in Crétaux and Birkett (2006) and Ricko et al. (2012), among others.
As for lake Chiquita, there is also a good agreement of the temporal variations of water height from ICESat, radar altimetry and in situ measurements, when available, for all other lakes considered in this study (see Figure 8). Lake-dependent biases between the time series can be observed, even between radar altimetry series.

<figure><img src="image_9.png" /><figcaption>Figure 8: Comparison between water level heights from available radar altimetry products, ICESat water level and relative in situ gauges (an arbitrary offset has been added) for lakes General Carrera, Musters, Viedma and Argentino (Argentina) from the “Base de Datos Hidrológica Integrada” (BIDH, http://bdhi.hidricosargentina.gob.ar/).</figcaption></figure>

We provide in Table 1 some statistics of the differences between lake water level computed from ICESat and from radar altimetry. Except for lake Poopo, there is a high correlation (larger than 0.8) between laser and radar altimetry water levels. Biases between time series can reach more than 1 meter, and can be the result of the different vertical datum, choice of the geoid model if used, and particulars of the processing scheme used. Standard deviations are lower than a few tens of centimeters.
<table><thead><tr><th>Lake</th><th>Radar product</th><th>Correlation with ICESat</th><th>Mean bias (m)</th><th>Standard Deviation (m)</th><th>Nb of samples</th></tr></thead><tbody><tr><td rowspan="3">Argentino</td><td>DAHITI</td><td>0.954</td><td>-0.992</td><td>0.374</td><td>16</td></tr><tr><td>G-REALM</td><td>0.960</td><td>-1.243</td><td>0.357</td><td>16</td></tr><tr><td>Hydroweb</td><td>0.964</td><td>-1.454</td><td>0.453</td><td>6</td></tr><tr><td rowspan="3">Chiquita</td><td>DAHITI</td><td>0.967</td><td>-0.701</td><td>0.176</td><td>39</td></tr><tr><td>G-REALM</td><td>0.958</td><td>-0.170</td><td>0.230</td><td>40</td></tr><tr><td>Hydroweb</td><td>0.972</td><td>-0.631</td><td>0.160</td><td>23</td></tr><tr><td rowspan="2">General Carrera</td><td>DAHITI</td><td>0.814</td><td>-0.895</td><td>0.331</td><td>17</td></tr><tr><td>Hydroweb</td><td>0.802</td><td>-0.373</td><td>0.333</td><td>20</td></tr><tr><td>Poopo</td><td>Hydroweb</td><td>0.200</td><td>-0.902</td><td>0.441</td><td>26</td></tr><tr><td>Titicaca</td><td>Hydroweb</td><td>0.821</td><td>0.100</td><td>0.329</td><td>25</td></tr></tbody></table>

Table 1: Correlation between ICESat water level and available radar altimetry products for 5 lakes in South America, mean bias and standard deviation (in meters) of the ICESat minus radar altimetry series.

## 6. Lakes mass changes from GRACE data

The Gravity Recovery and Climate Experiment (GRACE), launched in April 2002, collected time variable gravity measurements through October of 2017. GRACE consisted of two twin satellites on a same polar orbit at an altitude of about 450 kilometers, and separated by about 220 km along track. The main goal of the mission is the recovery of the temporal variations of the Earth's gravity field, through the precise measurements of the inter-satellite differential speed using a K-band microwave ranging system, of each satellite accelerations using electrostatic accelerometers and of each satellite orbits (see Tapley et al., 2004 for more details).

GRACE solutions are generally expressed as monthly sets of spherical harmonic Stokes coefficients up to degree 60 or 90, i.e. with corresponding spatial resolutions of about 300 or 225 km. The data require in addition some post processing filtering, which ultimately reduces their spatial resolution. In addition to the classical unconstrained spherical harmonic
solutions, several other processing schemes have been developed, including localized
constrained solutions, or mascon solutions (Rowlands et al., 2010). These mascon solutions
do not require any post-processing filtering and tend to have improved spatial resolution, as
shown, for example, by Luthcke et al. (2013) in Greenland and Alaska.

We use the latest version (v02.4) of the 1-degree equal area global iterated mascon solution
provided by NASA Goddard Space Flight Center from January 2003 to July 2016 (updates of
Luthcke et al., 2013; Loomis and Luthcke, 2017). A mascon (mass concentration) represents
an excess or a lack of surface mass compared to the a priori mean gravity field in a
predefined region, directly inverted from the GRACE measurements. The excess or lack of
mass is represented as a uniform water layer, expressed in centimeters of equivalent water
over this area (Rowlands et al., 2010).

Mascon solutions developed at NASA Goddard Space Flight Center differ from the classical
constrained or unconstrained spherical harmonic solutions, as they allow to introduce
spatial constraints in the inversion of the K-band range rate (KBRR) residuals (Rowlands et
al., 2010), and to increase the spatial resolution of the solution, as post-processing filtering is
no longer required. The iteration scheme introduced in Luthcke et al. (2013) and Loomis and
Luthcke (2017), where the solution determined at each iteration is used as a forward model
in the next iteration, significantly reduces the KBRR residuals globally. This solution has been
shown to be very successful in exhibiting small wavelength mass variations, for example the
ice-mass loss in Alaska (Luthcke et al., 2013), or oceanic mass variations in the Red Sea
(Loomis and Luthcke, 2017).

The forward modeling used in this GSFC solution differs from the official GRACE products
available from the mission. Atmosphere and ocean de-aliasing products are computed using
ECMWF (European Centre for Medium-range Weather Forecasts) operational model
available at 3-hour resolution and TUGO-m barotropic ocean models (update of Carrère and
Lyard, 2003). Ocean tides are modeled using GOT4.7 (Ray, 1999). The C20 term is replaced
with the GRACE TN-07 product (Cheng et al., 2013), and geocenter is introduced using
Swenson et al. (2008) methodology.
Over land, GRACE is sensing the total mass changes, including solid Earth and continental
hydrology contributions, which cannot be separated without the help of other
measurements or models. To extract lake mass variations from GRACE, all other sources of
mass changes have to be removed, including mass variation from the solid Earth, snow and
glaciers, soil-moisture and ground water. We remove these hydrological contributions using
global hydrology models, as they do not generally include surface waters, such as lakes and
reservoirs. Figure 9 shows the linear trend, annual amplitude and variability of this GRACE
solution, after correcting from the solid earth contribution from Glacial Isostatic Adjustment
(GIA) using A et al. (2013) model over South America, and the comparison with two state-of-
the-art global hydrology models, e.g. GLDAS2.1/Noah (Rodell et al., 2004) and MERRA2
(Gelaro et al., 2017).
<figure><img src="image_10.png" /></figure>

Figure 9: Linear trend, annual amplitude and variability of GRACE (v02.4) iterated global mascon solution (Luthcke et al., 2013; Loomis and Luthcke, 2017), after correcting for Glacial Isostatic Adjustment (GIA) (A et al., 2013), compared to the output from GLDAS2.1/Noah (Rodell et al., 2004) and MERRA2 (Gelaro et al., 2017) hydrology models, computed over the entire GRACE time series (2003-2016).

Except for the Amazon basin, there is a good agreement between GRACE seasonal and high frequency derived terrestrial water storage variations and GLDAS and MERRA2 hydrology models at seasonal timescales, but also at higher frequency. For the Amazon basin, the
differences are mainly due to the lack of the river storage component in most hydrology models, as it represents almost half of the full variations (Han et al., 2010).

Linear trends of hydrology models and GRACE significantly differ over the Southern part of the Andes mountains; the decrease in GRACE trends are mostly due to the present-day Patagonia ice fields melting (see, for example, Li et al., 2019 and Foresta et al., 2018). Permanently ice-covered regions are generally not well taken into account in global hydrology models.

The 1-degree mascon cells have not been optimized for the retrieval of lake mass variations, therefore some of the lakes may cover several cells. Even if the leakage errors are minimized with the mascon technique compared to the spherical harmonic solutions, we cannot neglect them. Leakage errors are due to the contamination of the signal from the surrounding areas because of the truncation of the high spatial frequency parts of the gravity spectrum. For this purpose, we add mass variations of each 1-degree mascon cells with a distance from its center to the lake less than 350 km for each lake. This 350-km threshold appears to be the best compromise, considered the effective GRACE spatial resolution. In order to assess the lake mass variations, we need to remove other geophysical signals remaining in GRACE solutions, mainly Glacial Isostatic Adjustment (GIA) and other hydrological contributions, such as snow, soil moisture, groundwater, etc. We choose the A et al. (2013) GIA model and we consider two global up-to-date hydrological models: GLDAS2.1 (Rodell et al., 2004) and MERRA2 (Gelaro et al., 2017). The differences between the GLDAS2.1 and MERRA2 are almost negligible for all lakes. For this reason, we only consider the lake volume estimates using GLDAS2.1 in the later part of the paper. The lake mass variation is then equal to GRACE observations averaged over 350-km around the lake minus the GIA and hydrology models in the same area.

## 7. Lakes water levels versus surface area

There are some challenges to estimate absolute lake water volume from a combination of these datasets. On one hand, the computation of the total volume of each lake or reservoir
requires the knowledge of the bathymetry and topography. On the other hand, GRACE does not allow to retrieve the full mass of the lake, but only its variations. For these reasons, we decide to compute and compare the relative lake volume variations. In order to derive the relative volume variations ΔV(t) of the various lakes and reservoirs, we need to integrate the following equation (see, for example, Abileah et al., 2011):

$$ \Delta V(t) = \int_{h_0}^{h(t)} S(h)dh $$

where h(t) is the temporal evolution of the lake level, h₀ a reference water level and S(h) is the surface water extent as a function of the water height. Because radar altimetry products offer a more complete time span and temporal resolution, we only use the radar altimeter (and not laser altimetry) to derive the hypsometry. Figure 10 shows the surface variations computed from MODIS NDVI as a function of water elevation from three different publicly available radar altimetry products.

<figure><img src="image_11.png" /></figure>
Figure 10: Lake surface variations (in km²) derived from MODIS NDVI as a function of the water level (in m) from three different radar altimetry products, e.g. Hydroweb in red (Crétaux et al., 2011), G-REALM in green (Birkett et al., 2011), and DAHITI in blue (Schwatke et al., 2015).

For Lakes Chiquita (Argentina) and Poopo (Bolivia), there is high correlation (greater than 0.95 and equal to 0.70, respectively) and a clear linear relation between water extent and water level, with estimates between 816 km²/m (DAHITI) and 1083 km²/m (G-REALM) for Lake Chiquita, and equal to 177 km²/m (Hydroweb) for Lake Poopo. The correlation for lake Colhué Huapi (Argentina) is also equal to 0.70, but there is a quadratic relation between height and surface extent. For the other lakes in this study, the relationship between water extent and water height is not clear, as evidenced in Figure 10. Therefore, we chose to consider the surface of those lakes as a constant surface.

## 8. Comparison of lake volume variations from the different techniques

Figure 11 shows the relative lake volume variations estimated from altimetry and MODIS, and from the GRACE mascons, for three Argentine glacial lakes in the Andes. As shown on Figure 9, the GRACE and hydrology model linear trends significantly differ in this area. For this reason, we had to remove this trend in the GRACE time series before the comparison. These different trends are likely caused by the mismodeling of the present-day ice mass loss in the Patagonian ice fields (see, for example, Foresta et al., 2018) in GLDAS2.1 or MERRA2 hydrology models. Even if the two estimates (altimetry + remote sensing, and GRACE - hydrology model) have the same order of magnitude, there are some differences that remain, especially for the phase of the annual variations. GRACE maxima is reached at the end of winter, and altimetry maxima at the early part of the winter. In addition to the significant mass decrease seen in GRACE, this may indicate that the hydrology models do not model properly the seasonal and long term snow and ice mass variations in the Andes.
<figure><img src="image_12.png" /></figure>

Figure 11: Comparison of relative volume variations for glacial lakes in Argentina and Chile from MODIS and radar altimetry (blue: DAHITI) or laser altimetry (black: ICESat). GRACE time series is shown in green, after removal of a linear trend computed over the 2002-2016 period.

Figures 12 and 13 show similar comparisons, but for other lakes located in different climatic areas, and different in size. For these water bodies, the volume variations derived from MODIS and altimetry and from GRACE agree well both at seasonal and longer timescales, even for the small Negro and El Chocón reservoirs, with a mean surface around 800 km². We prove here the ability of the NASA GSFC global mascons to retrieve small wavelength mass variations in continental areas.
Long term evolution of water volumes of Lakes Chiquita (Troin et al., 2010) and Titicaca (Pillco et al., 2019), shown on Figure 12, is linked to changes in precipitation and evaporation. Lake Poopo is partly filled from flows originating from Lake Titicaca and both lake levels are highly correlated. Water volume variations are linked to the inflow from Lake Titicaca, precipitation and evaporation changes (Pillco and Bengtsson, 2006).
<figure><img src="image_13.png" /><figcaption>Figure 12: Comparison of volume variations for Lakes Titicaca and Poopo (Bolivia), Negro and El Chocón (Argentina) from MODIS and radar altimetry (red: Hydroweb) or laser altimetry (black: ICESat). GRACE time series is shown in green.</figcaption></figure>
<figure><img src="image_14.png" /><figcaption>Figure 13: Comparison of relative volume variations for lake Chiquita (Argentina) from GRACE mascons (green), from radar altimetry (red) and ICESat (black) using our hypsometry model, and the Hydroweb estimates (blue).</figcaption></figure>

Figure 13 shows a more detailed comparison of the different relative volume variation estimates for lake Chiquita (Argentina). We especially show the differences between the conversion of the same water level height (radar altimetry measurements from Hydroweb) into volume using different hypsometry, ours and the ones used by Hydroweb. The volume increase in 2008-2009 in the Hydroweb database is clearly not observed in our volume estimates using water level measurements from ICESat (black) and radar (red) altimetry when we use our surface extent variations, and it is neither observed in the GRACE derived series (green). In general, our estimates are in better agreement with volume changes retrieved by GRACE (see, for example, the small volume increase in 2007 of a few km³). The large increase in the GRACE derived water volume estimates for Lake Chiquita in 2010 may be due to a mis-modeling of the intense precipitation in the GLDAS2.1 hydrology model that was taken out, which may not accurately capture extreme events like floods. The general increase on the total water storage in the basin is clearly seen by GRACE (see Abelen et al., 2015).
When converting GRACE lake mass changes into volume variations, we assume a constant water density of 1000 kg/m³. For large lakes, especially with strong seasonal temperature variations, this may introduce a small error in the GRACE estimates, and therefore explain some of the seasonal amplitude differences between GRACE and altimetry-derived volume changes at seasonal timescales.

Some statistics of the lakes surface extent variations and comparisons of the volume variations derived from DAHITI, Hydroweb, and ICESat altimetry and MODIS, using our hypsometry, and GRACE are provided in Table 2. In comparison to the time series provided by radar altimetry and imagery (2000-2018), the ICESat time series (2003-2009) and GRACE time series (2003-2016) are much shorter. ICESat campaign-based observations scheme to manage laser lifetime, operating during 3 or 2 months per year (February-March, May-June and October-November), optimizing monitoring ice sheet and sea ice mass changes in the polar regions, makes the comparison between long-term time series from radar altimetry and ICESat not trivial. Gaps in the table represent lack of available data for radar altimetry. The statistics show that, with the exception the volumes computed for the lakes at high elevations, there is a good agreement in amplitude and phase of the seasonal variations, as well as for the standard deviations of lake volume changes computed from radar and laser altimetry and MODIS/NDVI, and the GRACE estimates.
<table><thead><tr><th rowspan="2">Lake Volume</th><th colspan="3">MODIS</th><th colspan="3">DAHITI</th><th colspan="3">Hydroweb</th><th colspan="3">ICESat</th><th colspan="3">GRACE</th></tr><tr><th>Long (°)</th><th>Lat (°)</th><th>Min</th><th>Max</th><th>stdev</th><th>ampl</th><th>phase</th><th>stdev</th><th>Ampl</th><th>phase</th><th>stdev</th><th>ampl</th><th>phase</th><th>stdev</th><th>ampl</th><th>phase</th><th>stdev</th></tr></thead><tbody><tr><td>Lake Argentino*</td><td>-72.50</td><td>-50.25</td><td>1944.00</td><td>2794.06</td><td>199.61</td><td>2.79</td><td>-48.6</td><td>2.14</td><td>2.81</td><td>-48.1</td><td>2.09</td><td>3.65</td><td>-65.4</td><td>3.11</td><td>2.05</td><td>-162.2</td><td>3.75</td></tr><tr><td>Lake Chiquita</td><td>-62.60</td><td>-30.50</td><td>4430.60</td><td>6767.31</td><td>483.42</td><td>0.43</td><td>20.5</td><td>7.71</td><td>0.81</td><td>7.03</td><td>7.13</td><td>0.15</td><td>-84.5</td><td>7.12</td><td>0.76</td><td>-70.8</td><td>4.15</td></tr><tr><td>Lake Chocón</td><td>-69.00</td><td>-39.50</td><td>853.69</td><td>965.94</td><td>17.84</td><td></td><td></td><td></td><td></td><td></td><td></td><td>0.53</td><td>-142.6</td><td>2.79</td><td>1.56</td><td>147.5</td><td>1.72</td></tr><tr><td>Lake General Carrera*</td><td>-71.75</td><td>-46.50</td><td>1983.88</td><td>2806.56</td><td>179.78</td><td>0.96</td><td>-62.9</td><td>0.96</td><td>1.26</td><td>-58.4</td><td>1.78</td><td>1.49</td><td>-64.4</td><td>1.39</td><td>2.53</td><td>-178.9</td><td>2.55</td></tr><tr><td>Lake Llanquihe*</td><td>-72.75</td><td>-41.00</td><td>577.12</td><td>868.38</td><td>57.32</td><td>0.17</td><td>137.7</td><td>0.17</td><td>0.18</td><td>143.1</td><td>0.18</td><td>0.10</td><td>-146.8</td><td>0.15</td><td>1.49</td><td>-152.3</td><td>1.83</td></tr><tr><td>Lake Mirim</td><td>-53.10</td><td>-32.85</td><td>4566.88</td><td>4943.38</td><td>76.38</td><td></td><td></td><td></td><td></td><td></td><td></td><td>1.80</td><td>103.5</td><td>5.44</td><td>2.18</td><td>-139.9</td><td>4.65</td></tr><tr><td>Lake Musters</td><td>-69.20</td><td>-45.42</td><td>504.62</td><td>566.94</td><td>9.05</td><td>0.29</td><td>167.5</td><td>0.34</td><td></td><td></td><td></td><td>0.34</td><td>177.5</td><td>0.14</td><td>1.62</td><td>167.4</td><td>1.75</td></tr><tr><td>Lake Negro</td><td>-56.00</td><td>-32.75</td><td>785.81</td><td>971.06</td><td>38.97</td><td></td><td></td><td></td><td></td><td></td><td></td><td>0.99</td><td>155.1</td><td>1.83</td><td>0.73</td><td>134.4</td><td>1.78</td></tr><tr><td>Lake Poopo*</td><td>-67.00</td><td>-18.75</td><td>1725.25</td><td>2029.50</td><td>57.41</td><td></td><td></td><td></td><td>0.72</td><td>-11.2</td><td>1.95</td><td>0.44</td><td>-78.9</td><td>1.35</td><td>0.99</td><td>-114.1</td><td>1.78</td></tr><tr><td>Lake Titicaca*</td><td>-69.25</td><td>-16.00</td><td>8384.75</td><td>8972.19</td><td>109.49</td><td></td><td></td><td></td><td>2.72</td><td>-9.0</td><td>6.23</td><td>2.39</td><td>6.02</td><td>6.49</td><td>3.91</td><td>-40.2</td><td>3.89</td></tr><tr><td>Lake Todos Los Santos*</td><td>-72.25</td><td>-41.25</td><td>127.94</td><td>238.31</td><td>20.70</td><td></td><td></td><td></td><td>0.07</td><td>76.1</td><td>0.11</td><td>0.16</td><td>88.0</td><td>0.11</td><td>1.05</td><td>-169.3</td><td>1.27</td></tr><tr><td>Lake Viedma*</td><td>-72.50</td><td>-49.50</td><td>1480.81</td><td>1930.94</td><td>106.38</td><td>1.41</td><td>-56.3</td><td>1.13</td><td></td><td></td><td></td><td>1.98</td><td>-59.2</td><td>1.98</td><td>2.14</td><td>-156.1</td><td>3.48</td></tr></tbody></table>

Table 2: Statistics of the lake surface changes (extremes and standard deviation), and annual relative volume variations (amplitude in km³ and phase in degrees) and standard deviation of the full time series from radar and laser altimetry and MODIS (DAHITI, Hydroweb and ICESat) and from GRACE global mascons. High elevation lakes are marked with a star after their name.

As illustrated above in Table 2, there are more differences with the GRACE-derived volume changes, even if the order of magnitude of the extremes and standard deviations are quite similar. The differences are larger for the Patagonian glacial lakes (Argentino, General Carrera, Llanquihe, Todo Los Santos and Viedma) where the hydrology models, which are needed to remove other hydrological contributions to GRACE mass changes, are known to have deficiencies.

For the lakes or reservoirs in arid or semi-arid regions (Chiquita, Negro and El Chocón), where all other hydrological contributions (soil moisture, snow, etc.) are small, there is a very good agreement between altimetry and imagery and time-variable gravity field
estimates. Because of the different sampling rates and duration, this is better shown on
Figures 12 and 13 than on Table 2.

Even though validation of GRACE products is difficult, and volume changes for lakes and reservoirs estimates challenging, we believe that it is still useful to use storage change estimates from GRACE iterated global mascon solutions to complement other available estimates. Time variable gravity solutions provide a global picture of water storage variations independent of other in-situ data. It is true that one of the main limitations is the quality of the hydrology models used to remove other contributions aside from that of surface water. Another way to model lake volume variations is to model precipitation minus evaporation in addition to river inflows and outflows. All these data may also suffer from significant errors.

9. Conclusions and perspectives

In the present study we compute relative volume variations from several lakes and reservoirs in South America using two different methodologies: the classical combination of remote sensing altimetry and imagery as in Busker et al. (2019), but also derived from GRACE global iterated mascons, after removal of soil-moisture and snow from global hydrology models. With ICESat observations from ~70 m footprints along track every 170 m, sub-satellite tracks repeated nearly monthly, spaced about 15 km at the equator, 11 km at mid-latitudes, and 2.6 km at 80-degree latitude, we are able to retrieve water level variations of small water bodies such as Lakes Chocón and Negro, from the laser altimeter. When available, the temporal variations of water levels from radar and laser altimetry agree at a few centimeter level. Some unresolved biases remain, and are even observed between different radar altimetry products. These biases are likely related to disparity of reference systems, geophysical corrections, and/or the use of different geoid models to convert the water elevations to orthometric heights.

With the exception of lake Poopo, correlation between ICESat water level estimates and
those from radar altimetry products have coefficients larger than 0.81 (Argentino, Chiquita,
General Carrera, and Titicaca). Biases between time series vary from 0.1 m than 1.5 m, likely due to vertical datum, choice of the geoid model, processing scheme differences. Standard deviations are below 0.45 m. Our surface water extent estimates compare as well as other products compare to each other. Except for Lake Titicaca, where extents are less accurate, the linear fit between our estimates and the two other available products is close to 1, with some small biases for Lakes Viedma and Negro.

Larger differences in volume changes derived from GRACE and those derived from altimetry and imagery are observed, especially for the Patagonian glacial lakes (Argentino, General Carrera, Llanquihe, Todo Los Santos and Viedma). This is likely due to deficiencies in the hydrological models. In semi-arid regions, with small soil moisture, snow, and other hydrological contributions (lakes Chiquita, Negro and El Chocón), there is a very good agreement between altimetry and imagery and time-variable gravity field estimates, with differences far below 1 km³.

As an integrative measurement, GRACE is not able to discriminate mass changes in each individual component of the water cycle. To access lake and reservoir mass variations, GRACE measurements have to be corrected for all other hydrological contributions (soil-moisture, snow, etc.). After removal of soil moisture and snow contributions using global hydrology models such as GLDAS2.1 (Rodell et al., 2004) or MERRA2 (Gelaro et al., 2017), isolating surface water contributions, it is possible to retrieve lake mass variations from the GRACE global 1-degree iterated mascons (Luthcke et al., 2013), even for small water bodies. When the quality of the hydrology model decreases, the GRACE estimates clearly differ from those obtained by the classical approach, using remote sensing altimetry and imagery. This is particularly the case for the high altitude lakes in the Andes, likely due to the mis-modeling of the present-day mass loss in the Patagonia ice fields in most hydrology models. Lake volume variations of small water bodies (surface less than 50,000 km²) cannot be individually retrieved from classical unconstrained or even constrained spherical harmonic solutions (see, for example, Longuevergne et al., 2013), showing the advantage of high resolution localized GRACE solutions to retrieve small wavelength mass variations at the Earth's surface.
No other laser altimeter has measured surface height variations from space since the end of the ICESat operations in 2009, and until the launch of the ICESat-2 satellite in September, 2018 (Markus et al., 2017, Neumann et al., 2019). In contrast to the single beam of the GLAS (Geoscience Laser Altimeter System), the instrument onboard ICESat, the ATLAS (Advanced Topographic Laser Altimeter System) onboard ICESat-2 has six beams arranged in three pairs (one strong and one weak beams), allowing a better estimation of the surface slope on ice-sheets. It also provides more ground coverage, with a nearly polar orbit and a 91-day repeat (33-day near repeat cycle). Preliminary analysis of ICESat-2 data has already demonstrated very promising results for inland water applications.

Current and upcoming laser altimetry mission data, with increased coverage, improved spatial and temporal sampling, like those provided by ICESat-2 and GEDI (Global Ecosystem Dynamics Investigation), with more precise measurements and advanced processing techniques and algorithms, including inland water products specifically formulated for these kinds of applications, will provide valuable data to combine with SWOT, and GRACE-Follow On mission data. Synergy between all these data products will help address the need for continuous monitoring of continental water storage variations from space measurements.
Acknowledgments

We thank Hydroweb (data available at http://hydroweb.theia-land.fr/), G-REALM (data available at http://ipad.fas.usda.gov/cropexplorer/global_reservoir/) and DAHITI (data available at http://dahiti.dgfi.tum.de/en/) research groups for providing lake and reservoir products used in this research. We also thank the Argentine government, the Ministerio del Interior, Obras Públicas y Vivienda and the Secretaría de Infraestructura y Política Hídrica for providing the BDHI in situ water level measurements (data available at http://bdhi.hidricosargentina.gob.ar/). We used the open source GMT (Generic Mapping Tools) software (Wessel et al., 2013) to prepare all figures. We would like to particularly thank the anonymous reviewers whose insightful comments have greatly improved this manuscript.
References

A, G., J. Wahr and S. Zhong (2013). Computations of the viscoelastic response of a 3-D
compressible Earth to surface loading: an application to Glacial Isostatic Adjustment in
Antarctica and Canada, Geophys. J. Int., 192, 557-572, doi:10.1093/gji/ggs030.

Abelen, S., F. Seitz, R. Abarca-del-Rio and A. Güntner (2015). Droughts and Floods in the La Plata Basin in Soil Moisture Data and GRACE, Remote Sensing, 7 (6), 7324-7349, doi: 10.3390/rs70607324.

Abileah, R., S. Vignudelli and A. Scozzari (2011). A completely remote sensing approach to monitoring reservoirs water volume, International Water Technology Journal, 1, 63–77.

Birkett, C.M., Reynolds, C., Beckley, B. and Doorn, B. (2011). From Research to Operations: The USDA Global Reservoir and Lake Monitor, chapter 2 in ‘Coastal Altimetry’, Springer Publications, eds. S. Vignudelli, A.G. Kostianoy, P. Cipollini and J. Benveniste, pp. 19-50, doi: 10.1007/978-3-642-12796-0_2.

Borsa, A. A., H. A. Fricker and K. M. Brunt (2019). Terrestrial Validation of ICESat Elevation Measurements and Implications, IEEE Transactions On Geoscience And Remote Sensing, 57 (9), 6946-6959, doi: 10.1109/TGRS.2019.2909739.

Bucher, E. H. and A. F. Stein (2016). Large Salt Dust Storms Follow a 30-Year Rainfall Cycle in the Mar Chiquita Lake (Córdoba, Argentina), PLOS ONE, 11 (6), e0156672, doi: 10.1371/journal.pone.0156672.

Buma, W. G., S. I. Lee and J. Y. Seo (2018). Recent Surface Water Extent of Lake Chad from Multispectral Sensors and GRACE. Sensors, 18 (7), 2082, doi:10.3390/s18072082.

Busker, T., A. de Roo, E. Gelati, C. Schwatke, M. Adamovic, B. Bisselink, J.-F. Pekel and A. Cottam (2019). A global lake and reservoir volume analysis using a surface water dataset
and satellite altimetry, Hydrol. Earth Syst. Sci., 23, 669-690, doi: 10.5194/hess-23-669-2019.

Carabajal, C. C., D. J. Harding and V. P. Suchdeo (2010). ICESat Lidar and Global Digital Elevation Models: Applications to DESDynl, Geoscience and Remote Sensing Symposium (IGARSS), 2010 IEEE International, Honolulu, HI, 1907-1910, doi: 10.1109/IGARSS.2010.5650201.

Carabajal, C. C., D. J. Harding, J.-P. Boy, J. J. Danielson, D. B. Gesch and V. P. Suchdeo (2011). Evaluation of the Global Multi-Resolution Terrain Elevation Data 2010 (GMTED2010) Using ICESat Geodetic Control, SPIE Proceedings, International Symposium on LIDAR and Radar Mapping: Technologies and Applications (LIDAR & RADAR 2011), Nanjing, China, 8266, 82661Y, doi: 10.1117/12.912776.

Carabajal, C. C. and J.-P. Boy (2016). Evaluation of ASTER GDEM V3 using ICESat Laser Altimetry, Int. Arch. Photogramm. Remote Sens. Spatial Inf. Sci., Volume XLI-B4, 117-124, https://doi.org/10.5194/isprs-archives-XLI-B4-117-2016.

Carrère, L. and F. Lyard (2003). Modeling the barotropic response of the global ocean to atmospheric wind and pressure forcing Comparisons with observations, Geophysical Research Letters, 30, 1275, doi: 10.1029/2002GL016473.

Cheng, M., B. D. Tapley and J. C. Ries (2013). Deceleration in the Earth's oblateness, J. Geophys. Res. Solid Earth, 118, 740-747, doi: 10.1002/jgrb.50058.

Crétaux, J.-F. and C. Birkett (2006). Lake studies from satellite radar altimetry, Comptes Rendus Geoscience, 338 (14-15), 1098-1112, doi: 10.1016/j.crte.2006.08.002.

Crétaux, J.-F., W. Jelinski, S. Calmant et al., (2011). SOLS: a lake database to monitor in the Near Real Time water level and storage variations from remote sensing data, Adv. Space Res., 47 (9), 1497-1507, doi: 10.1016/j.asr.2011.01.004.
Crétaux, J.-F., R. Abarca-del-Río, M. Bergé-Nguyen, A. Arsen, V. Drolon, G. Clos and P. Maisongrande (2016). Lake Volume Monitoring from Space, Surveys in Geophysics, 37 (2), 269-305, doi: 10.1007/s10712-016-9362-6.

Depetris, P. J. and A. I. Pasquini (2000). The hydrological signal of the Perito Moreno Glacier damming of Lake Argentino (southern Andean Patagonia): the connection to climate anomalies, Global and Planetary Change, 26 (4), 367-374, doi: 10.1016/S0921-8181(00)00049-7.

Didan, K. (2015). MOD13Q1 MODIS/Terra Vegetation Indices 16-Day L3 Global 250m SIN Grid V006, NASA EOSDIS LP DAAC, doi: 10.5067/MODIS/MOD13Q1.006.

Farr, T.G., P. A. Rosen, E. Caro, R. Crippen, R. Duren, S. Hensley, M. Kobrick, M. Paller, E. Rodriguez, L. Roth, D. Seal, S. Shaffer, J. Shimada, J. Umland, M. Werner, M. Oskin, D. Burbank and D. Alsdorf (2007). The Shuttle Radar Topography Mission, Reviews of Geophysics, 45, RG2004, doi: 10.1029/2005RG000183.

Foresta, L., N. Gourmelen, F. Weissgerber, P. Nienow, J. J. Williams, A. Shepherd, M. R. Drinkwater and S. Plummer (2018). Heterogeneous and rapid ice loss over the Patagonian Ice Fields revealed by CryoSat-2 swath radar altimetry, Remote Sensing of Environment, 211, 441-455, doi: 10.1016/j.rse.2018.03.041.

Gelaro, R., W. McCarty, M. J. Suárez et al. (2017). The Modern-Era Retrospective Analysis for Research and Applications, Version 2 (MERRA-2), J. Climate, 30, 5419-5454, doi: 10.1175/JCLI-D-16-0758.1.

Han, S.-C., I.-Y. Yeo, D. Alsdorf, P. Bates, J.-P. Boy, H. Kim, T. Oki, and M. Rodell (2010). Movement of Amazon surface water from time-variable satellite gravity measurements and implications for water cycle parameters in land surface models, Geochem. Geophys. Geosyst., 11, Q09007, doi:10.1029/2010GC003214.
Harding, D. J., and C. C. Carabajal (2005). ICESat Waveform Measurements of Within-footprint Topographic Relief and Vegetation Vertical Structure, *Gephys. Res. Lett.*, L21S10, doi: 10.1029/2005GL023471.

Huang, C. (2018). Seeing surface water from space, Eos, 99, doi: 10.1029/2018EO103123.

Huang, C., Chen, Y., Zhang, S., and Wu, J. (2018). Detecting, extracting, and monitoring surface water from space using optical sensors: A review. Reviews of Geophysics, 56, 333–360, doi: 10.1029/2018RG000598.

Jasinski, M. F. and J. D. Stoll. (2019). ATLAS/ICESat-2 L3A Inland Water Surface Height, Version 1, Boulder, Colorado USA. NASA National Snow and Ice Data Center Distributed Active Archive Center, doi: https://doi.org/10.5067/ATLAS/ATL13.001 (last access, August 27, 2019).

Jasinski, M., J. Stoll, D. Hancock et al. (2019). Algorithm Theoretical Basis Document (ATBD) for Inland Water Data Products ATL13, version 1 (May 28, 2019), https://icesat-2.gsfc.nasa.gov/sites/default/files/page_files/ICESat2_ATL13_ATBD_r001.pdf (last access, August 27, 2019).

Justice, C. O., J. R. G. Townshend, E. F. Vermote, E. Masuoka, R. E. Wolfe, N. Saleous, D. P. Roy and J. T. Morisette (2002). An overview of MODIS Land data processing and product status, Remote Sens. Environ. 83, 3–15, https://doi.org/10.1016/S0034-4257(02)00084-6.

Khandelwal, A. and V. Kumar (2019). Pre SWOT Hydrology Global Lake/Reservoir Surface Inland Water Area Extent V2. Ver. 2. PO.DAAC, CA, USA. 10.5067/UCLRS-AREV2, https://podaac.jpl.nasa.gov/dataset/PRESWOT_HYDRO_L3_LAKE_RESEVOIR_AREA_V2.

Khandelwal, A., A. Karpatne, M. E. Marlier, J. Kim, D. P. Lettenmaier and V. Kumar (2017). An approach for global monitoring of surface water extent variations in reservoirs using MODIS data, Remote Sensing of Environment, 202, 113-128, doi: 10.1016/j.rse.2017.05.039.
Kokot, R. R., (2004). Erosión en la costa patagónica por cambio climático, Rev. Asoc. Geol. Argent., 59 (4), 715-726.

Lee, H., M. Durand, H. C. Jung, D. Alsdorf, C. K. Shum and Y. Sheng (2010). Characterization of surface water storage changes in Arctic lakes using simulated SWOT measurements, International Journal of Remote Sensing, 31 (14), 3931-3953, doi: 10.1080/01431161.2010.483494.

Li, Z., J. L. Chen, S. Ni, L. Tang and X. G. Hu (2019). Long-term and inter-annual mass changes of Patagonia Ice Field from GRACE, Geodesy and Geodynamics, 10, 100-109, doi: 10.1016/j.geog.2018.06.001.

Longuevergne, L., C. R. Wilson, B. R. Scanlon and J. F. Crétaux (2013). GRACE water storage estimates for the Middle East and other regions with significant reservoir and lake storage, Hydrol. Earth Syst. Sci., 17, 4817-4830, doi: 10.5194/hess-17-4817-2013.

Loomis, B. and S. B. Luthcke (2017). Mass evolution of Mediterranean, Black, Red, and Caspian Seas from GRACE and altimetry: accuracy assessment and solution calibration, Journal of Geodesy, 91 (2), 195-206, doi: 10.1007/s00190-016-0952-3.

Luthcke, S. B., T. J. Sabaka, B. D. Loomis, A. A. Arendt, J. J. McCarthy and J. Camp (2013). Antarctica, Greenland and Gulf of Alaska land ice evolution from an iterated GRACE global mascon solution, J. Glaciology, 59 (216), 613-631, doi: 10.3189/2013JoG12J147.

Markus, T., T. Neumann, A. Martino et al. (2017). The Ice, Cloud, and land Elevation Satellite-2 (ICESat-2): Science requirements, concept, and implementation, Remote Sensing of Environment, 190, 260-273, doi: 10.1016/j.rse.2016.12.029.

Moore, P., and S. D. P. Williams (2014). Integration of altimetric lake levels and GRACE gravimetry over Africa: Inferences for terrestrial water storage change 2003–2011, Water Resour. Res., 50, 9696–9720, doi:10.1002/2014WR015506.
Nandigam, V. K. Lin, M. Phan, A. Borsa, S. J. S. Khalsa and C. Crosby (2018). Rapid Access and Visualization of Spaceborne Altimetry Data from ICESat and ICESat-2, IGARSS 2018 - 2018 IEEE International Geoscience and Remote Sensing Symposium, Valencia, 2018, 5270-5273, doi: 10.1109/IGARSS.2018.8518105.

Neumann, T. A., A. J. Martino, T. Markus, S. Bae, M. R. Bock, A. C. Brenner, K. M. Brunt, J. Cavanaugh, S. T. Fernandes, D. W. Hancock, K. Harbeck, J. Lee, N. T. Kurtz, P. J. Luers, S. B. Luthcke, L. Magruder, T. A. Pennington, L. Ramos-Izquierdo, T. Rebold, J. Skoog, T. C. Thomas (2019). The Ice, Cloud, and Land Elevation Satellite-2 mission: A global geolocated photon product derived from the Advanced Topographic Laser Altimeter System, Remote Sensing of Environment 233, 111325, doi: 10.1016/j.rse.2019.111325.

O'Loughlin, F. E., J. Neal, D. Yamazaki, and P. D. Bates (2016). ICESat-derived inland water surface spot heights, Water Resour. Res., 52, 3276-3284, doi: 10.1002/2015WR018237.

Omute, P. R. Corner and J. L. Awange (2012). The use of NDVI and its Derivatives for Monitoring Lake Victoria's Water Level and Drought Conditions, Water Resources Management, 26 (6), 1591-1613, doi: 10.1007/s11269-011-9974-z.

Pasquini, A. I., K. L. Lecomte and P. J. Depetris (2008). Climate change and recent water level variability in Patagonian proglacial lakes, Argentina, Global and Planetary Change, 63 (4), 290-298, doi: 10.1016/j.gloplacha.2008.07.001.

Pavlis, N. K., S. A. Holmes, S. C., Kenyon and J. K. Factor (2012). The development and evaluation of the Earth Gravitational Model 2008 (EGM2008), J. Geophys. Res., 117, B04406, doi: 10.1029/2011JB008916.

Pillco, Z. R., and L. Bengtsson (2006). Long-term and extreme water level variations of the shallow Lake Poopó, Bolivia, Hydrological Sciences Journal, 51 (1), 98-114, doi: 10.1623/hysj.51.1.98.
Pillco, Z. R., L. Bengtsson, R. Berndtsson et al. (2019). Modelling Lake Titicaca's daily and monthly evaporation, Hydrol. Earth Syst. Sci., 23, 657-668, doi: 10.5194/hess-23-657-2019.

Ray, R. D. (1999). A global ocean tide model from Topex/Poseidon altimetry: GOT99.2, NASA Tech Memo 209478, 58 pp., Goddard Space Flight Center.

Ricko, M., C. M. Birkett, J. A. Carton and J.-F. Crétaux (2012). Intercomparison and validation of continental water level products derived from satellite radar altimetry, Journal of Applied Remote Sensing, 6 (1), 1-24, doi: 10.1117/1.JRS.6.061710.

Rivas, A. (1978). Batimetría del lago Colhué Huapí. Informe Técnico del Centro Nacional Patagónico (Bathymetry of the Colhué Huapí lake. Technical Report of the National Patagonian Center), 18, 1-8, http://www.repositorio.cenpat-conicet.gob.ar/bitstream/handle/123456789/181/0018B-BatimetriaDelLagoColhueHu.pdf (last access, March 20, 2020).

Rodell, M., P.R. Houser, U. Jambor, J. Gottschalck, K. Mitchell, C. Meng, K. Arsenault, B. Cosgrove, J. Radakovich, M. Bosilovich, J.K. Entin, J.P. Walker, D. Lohmann and D. Toll (2004). The Global Land Data Assimilation System, Bull. Amer. Meteor. Soc., 85, 381-394, doi: 10.1175/BAMS-85-3-381.

Rowlands, D. D., S. B. Luthcke, J. J. McCarthy, S. M. Klosko, D. S. Chinn, F. G. Lemoine, J.-P. Boy and T. J. Sabaka (2010). Global mass flux solutions from GRACE: A comparison of parameter estimation strategies - Mass concentrations versus Stokes coefficients, J. Geophys. Res., 115, B01403, doi: 10.1029/2009JB006546.

Satgé, F., R. Espinoza, R. Z. Pillco et al. (2017). Role of Climate Variability and Human Activity on Poopó Lake Droughts between 1990 and 2015 Assessed Using Remote Sensing Data, Remote Sensing, 9 (3), 218, doi: 10.3390/rs9030218.
Schutz, B. E., H. J. Zwally, C. A. Shuman, D. Hancock, and J. P. DiMarzio (2005), Overview of the ICESat Mission, Geophys. Res. Lett., 32, L21S01, doi: 10.1029/2005GL024009.

Schwatke, C., D. Dettmering, W. Bosch and F. and Seitz (2015). DAHITI - an innovative approach for estimating water level time series over inland waters using multi-mission satellite altimetry, Hydrol. Earth Syst. Sci., 19, 4345-4364, doi: 10.5194/hess-19-4345-2015.

Scordo, F., V. Y. Bohn, M. C. Piccolo and G. M. E. Perillo (2018). Mapping and Monitoring Lakes Intra-Annual Variability in Semi-Arid Regions: A Case of Study in Patagonian Plains (Argentina), Water, 10, 889, doi: 10.3390/w10070889.

Sun, X., J. B. Abshire, A. A. Borsa, H. A. Fricker, D., Yi, J. P. DiMarzio, F. S. Paolo, K. M. Brunt, D. J. Harding and G. A. Neumann (2017). ICESAT/GLAS Altimetry Measurements: Received Signal Dynamic Range and Saturation Correction, IEEE Transactions on Geoscience and Remote Sensing, 55 (10), 5440-5454, doi: 10.1109/tgrs.2017.2702126.

Swenson, S., D. Chambers and J. Wahr (2008). Estimating geocenter variations from a combination of GRACE and ocean model output, J. Geophys. Res., 113, B08410, doi:10.1029/2007JB005338.

Tapley, B. D., S. Bettadpur, J. C. Ries, P. Thompson and M. M. Watkins (2004). GRACE Measurements of Mass Variability in the Earth System, Science, 305 (5683), 503-505, doi: 10.1126/science.1099192.

Troin, M., C. Vallet-Coulomb, F. Sylvestre and E Piovano (2010). Hydrological modelling of a closed lake (Laguna Mar Chiquita, Argentina) in the context of 20th century climatic changes, Journal of Hydrology, 393 (3-4), 233-244, doi: 10.1016/j.jhydrol.2010.08.019.

Tulbure, M. G., M. Broich, S. V. Stehman and A. Kommareddy (2016). Surface water extent dynamics from three decades of seasonally continuous Landsat time series at
subcontinental scale in a semi-arid region, Remote Sensing of Environment, 178, 142-157,
doi: 10.1016/j.rse.2016.02.034.

Wang, X., P. Gong, Y. Zhao et al. (2013). Water-level changes in China's large lakes determined from ICESat/GLAS data, Remote Sens. Environ., 132, 131-144, doi: 10.1016/j.rse.2013.01.005.

Wessel, P., W. H. F. Smith, R. Scharroo, J. Luis, and F. Wobbe (2013). Generic Mapping Tools: Improved Version Released, EOS Trans. AGU, 94 (45), 409-410, doi: 10.1002/2013EO450001.

Wilson, R., N. F. Glasser, J. M. Reynolds, S. Harrison, P. Iribarren Anacona, M. Schaefer and S. Shannon (2018). Glacial lakes of the Central and Patagonian Andes, Global and Planetary Change, 162, 275-291, doi: 10.1016/j.gloplacha.2018.01.004.

Wu, H., N. Wang, X. Jiang and Z. Guo (2014). Variations in water level and glacier mass balance in Nam Co lake, Nyainqentanglha range, Tibetan Plateau, based on ICESat data for 2003-09, Annals of Glaciology, 55 (66), 239-247, doi: 10.3189/2014AoG66A100.

Ye, Z., H. Liu, Y. Chen, S. Shu, Q. Wu and S. Wang (2017). Analysis of water level variation of lakes and reservoirs in Xinjiang, China using ICESat laser altimetry data (2003–2009), PLoS ONE 12 (9), e0183800, doi: 10.1371/journal.pone.0183800.

Zhang, G., H. Xie, S. Kang, D. Yi, S. F. Ackley (2011). Monitoring lake level changes on the Tibetan Plateau using ICESat altimetry data (2003–2009), Remote Sensing of Environment, 115 (7), 1733-1742, doi: 10.1016/J.RSE.2011.03.005.

Zwally, H. J., R. Schutz, W. Abdalati, J. Abshire, C. Bentley, J. Bufton, D. Harding, T. Herring, B. Minster, J. Spinhirne and R. Thomas (2002). ICESat's laser measurements of polar ice, atmosphere, ocean, and land, Journal of Geodynamics, 34 (3-4), 405-445, doi: 10.1016/S0264-3707(02)00042-X.
Zwally, H. J., R. Schutz, C. Bentley, J. Bufton, T. Herring, J. Minster, J. Spinhirne and R. Thomas (2011). GLA14: GLAS/ICESat L2 Global Land Surface Altimetry Data, Release 33; NASA National Snow and Ice Data Center Distributed Active Archive Center: Boulder, CO, USA.