Monitoring Long-Term Lake Level Variations in Middle and Lower Yangtze Basin over 2002–2017 through Integration of Multiple Satellite Altimetry Datasets

: Satellite altimetry has been e ﬀ ectively used for monitoring lake level changes in recent years. This work focused on the integration of multiple satellite altimetry datasets from ICESat-1, Envisat and Cryosat-2 for the long-term (2002–2017) observation of lake level changes in the middle and lower Yangtze River Basin (MLYB). Inter-altimeter biases were estimated by using the gauged daily water level data. It showed that the average biases of ICESat-1 and Cryosat-2 with respect to Envisat were 6.7 cm and 3.1 cm, respectively. The satellite-derived water levels were evaluated against the gauged data. It indicated signiﬁcantly high correlations between the two datasets, and the combination of three altimetry data produced precise water level time series with high temporal and spatial resolutions. A liner regression model was used to estimate the rates of lake level changes over the study period after the inter-altimeter bias adjustment was performed. The results indicated that ~79% of observed lakes (41 / 52) showed increasing trends in water levels with rates up to 0.203 m / y during 2002–2017. The temporal analysis of lake level variations suggested that ~60% of measured lakes (32 / 53) showed decreasing trends during 2002–2009 while ~66% of measured lakes (79 / 119) exhibited increasing trends during 2010–2017. Most of measured reservoirs displayed rapidly rising trends during the study period. The driving force analysis indicated that the temporal heterogeneity of precipitation can be mainly used to explain the observed pattern of lake level changes. The operation of reservoirs and human water consumption were also responsible for the lake level variations. This work demonstrated the potential of integrating multiple satellite altimeters for the long-term monitoring of lake levels, which can help to evaluate the impact of climate change and anthropogenic activities on regional water resources.


Introduction
Lakes play a crucial role in the terrestrial hydrological cycle through complex processes at interfaces with the atmosphere and oceans and are valuable natural resources [1]. Surface water stored early radar altimeters, the quality and the number of lake height observations from different satellites can be inconsistent with the nominal accuracy and the nominal data volume, posing difficulties for the construction of multiple decades' lake level record. With careful selection of satellite altimeters, it is anticipated that a few decades of accurate lake level records can be achieved using the recent altimeter missions (i.e., ICESat-1, Envisat and Cryosat-2). The satellite-derived observations over lakes are noisy and usually affected by the hooking effect due to the large size of radar altimeter footprint and contaminations from the steep lakeshore or surrounding land [15,23]. In addition, there exist systematic errors among satellite altimeters due to inter-satellite biases, instrumental drifts, and different orbit realizations, which may also impede the long-term continuous observation of lake levels using multiple consecutive altimeter missions [24]. Therefore, sophisticated data processing methods like waveform retracking and inter-satellite calibration were suggested for the retrieval of lake levels using multiple satellite altimeters [21].
The middle and lower Yangtze Basin (MLYB) is home to the most densely distributed lake clusters in China. These lakes play important roles in supplying water resources and controlling flooding and in providing transit habitats for many rare or endangered migratory birds [25]. However, in recent decades, many lakes have undergone various transformations, largely as a result of natural factors such as climate, hydrologic, and ecosystem changes, as well as human activities [26,27]. In addition, the operation of numerous dams changed the natural water gradients at lake outlets/mouths meeting the rivers, and thus impacted inundation patterns of the downstream lake and wetland systems [28,29]. These changes have caused serious ecological and environmental problems, including eutrophication, reduction of lake water quantities, and destruction of ecosystems in the MLYB [27,28]. Therefore, it is important to monitor the long-term variations of lake levels and identify the driving forces for these changes. Previous studies on lake level monitoring in the MLYB only focused on a few major lakes (e.g., Tai Lake, Poyang Lake, and Dongting Lake) using a single satellite altimetry over several years span [1,30]. Hence, the spatio-temporal variation of lake water levels in the MLYB is still poorly known. Therefore, the primary aim of this study is to demonstrate the value of multiple satellite altimeters for generating long-term lake level time series in the MLYB and examine the driving forces for lake level changes. Specific key objectives are: (1) to evaluate the performance of multi-satellite altimetry; (2) to monitor the long-term lake level and identify the spatio-temporal pattern of lake level changes; and (3) to analyze relations between lake level variation and various driving forces.

Study area
The Yangtze River (YR) flows 6300 km through a wide variety of terrains from the Tibet Plateau to its mouth on the East China Sea. It can be divided into upper, middle and lower reaches according to the hydrological and geographical characteristics. The middle and lower reaches of YR extend from the Three Gorge Dam (TGD) to the river mouth, with a length of 1900 km. Formed by alluvial deposits from the YR and its tributaries, the MLYB is situated in middle and eastern China (106 • 7 -121 • 47 E, 24 • 30 -33 • 54 N), covering an area of 7.8 × 10 5 km 2 (Figure 1a). The MLYB can be further divided into six sub-basins (i.e., Han River, Yichang to Hukou, Dongting Lake, Poyang Lake, Lower Hukou, Tai Lake; Figure 1b) according to the spatial distribution of lakes and geographical differences [31]. Influenced by the subtropical monsoon, this region is characterized by a warm and humid climate, with an annual mean temperature of 14-18 • C and an annual mean precipitation of about 1300 mm.
The MLYB is a low-lying alluvial plain with elevations of less than 50 m, consisting of a large number of lakes and rivers. The statistics show that there are up to 100 lakes (area ≥10 km2) covering a total area of 10,593.4 km2 in this region, accounting for 38.17% of the total area of all freshwater lakes in China [32]. Among these lakes, a total of four of the five largest freshwater lakes of china are situated in this region, constituting the important components of the lake systems in the MLYB. Most of the lakes are distributed along the YR (Figure 1b), which can be categorized into two classes depending on their relations to the YR [28]. The class I lakes are those directly connected with the YR, including Poyang Lake, Dongting Lake, Shijiu Lake and Wu Lake. The class II lakes are those that are not directly Remote Sens. 2020, 12, 1448 4 of 19 connected to the YR. The MLYB is also spotted with numerous dammed reservoirs, with areas ranging from less than 1 km 2 to several thousand km 2 . The reservoirs are generally distributed far from the YR (Figure 1b). The operation of reservoirs influences the discharge feeding to lakes of downstream and YR.
Remote Sens. 2020, 12, x FOR PEER REVIEW 4 of 18 with areas ranging from less than 1 km 2 to several thousand km 2 . The reservoirs are generally distributed far from the YR (Figure 1b). The operation of reservoirs influences the discharge feeding to lakes of downstream and YR.

Data Sources
The nadir-pointing radar altimeter (RA-2) onboard Envisat operated in Ku-band and S-band. The Ku radar is a pulse-limited altimeter which emitted pulses at 1800 Hz, but with a subsequently averaging of 100 return pulses onboard the satellite, resulting in an 18 Hz product being transmitted

Envisat/RA-2
The nadir-pointing radar altimeter (RA-2) onboard Envisat operated in Ku-band and S-band. The Ku radar is a pulse-limited altimeter which emitted pulses at 1800 Hz, but with a subsequently averaging of 100 return pulses onboard the satellite, resulting in an 18 Hz product being transmitted to the ground stations [17]. The ground-track has a spatial separation of about 85 km at the equator, with footprints of about 1.7 km in diameter and with an interval of about 390 m [33]. In this study, the Sensor Geophysical Data Record (SGDR) Full Mission Reprocessing (FMR) v3.0 product spanning from 2002 to 2012 was used to derive the lake levels. The SGDR data products were reprocessed for all cycles from 6 to 113 into a homogeneous standard in NetCDF format, which includes the data in the GDR product (RA-2 geophysical data, MWR data) and also RA-2 averaged waveforms (18Hz) and RA-2 individual waveforms (1800Hz).

CryoSat-2/SIRAL
CryoSat-2 carried a synthetic aperture interferometric radar altimeter (SIRAL). The SIRAL is a single Ku-band radar altimeter using the full deramp range compression, operating in three different modes: Synthetic Aperture Radar mode (SAR), SAR Interferometer mode (SARIn), as well as Low Resolution mode (LRM). The LRM is a conventional pulse-limited radar altimeter while the SAR and SARIn utilize a Delay/Doppler radar altimeter with finer along-track spatial resolution [34]. In the MLYB, CryoSat-2 operates in the LRM, which has a footprint diameter of approximately 1.65 km on the flat water surface, with an interval of about 300 m (sampling at 20 HZ). ESA provides different levels of data products. In this work, the Level-1b product from 2010 to 2017 was used. Level-1b data files contain geolocated and time-stamped waveforms with accompanying altitude, position, corrections, interferometric phase difference, coherence waveforms, etc. [34].

ICESat-1/GLAS
The GLAS onboard ICESat-1 is the first space-borne laser altimetry, which measures the two-way travel time of a pulse that is reflected by the ground at a frequency of 40 Hz with two channels, 532 nm and 1064 nm [35]. GLAS produces a series of approximately 70 m diameter footprints that are separated by nearly 170 m intervals along the track and 30 km across the track (at the equator) and 5 km at 80 • latitude. NASA provided 15 types of GLAS data product for scientific uses. In this paper, the GLA14 (the global land surface altimetry data) product release 34, acquired in L2 and L3 campaigns from 2003 to 2009, is used to monitor the lake levels. The data contain 95 different parameters in a record of different situations, which can be downloaded through the National Snow and Ice Data Center (NSIDC).

Auxiliary Data
The Tropical Rainfall Measurement Mission (TRMM) derived precipitation data were collected to examine the relationship between the rainfall and the lake water level over the study area. The TRMM is a collaborative effort between NASA and the Japanese Aerospace Exploration Agency (JAXA), designed to monitor and study tropical rainfall. Previous studies showed that there was a good correlation between TRMM precipitation data and in-situ data over Tai Lake Basin, with the determination coefficient (R 2 ) ranging from 0.89 to 0.94 [36]. In this work, the gridded TRMM 3B43 data from 2002 to 2017 was used, which had a calendar month temporal resolution and a 0.25 • by 0.25 • spatial resolution [37]. The dataset was downloaded from the NASA website.
The MODIS (Moderate Resolution Imaging Spectroradiometer) land-water mask (MOD44W) was used to mask out the lakes and reservoirs. The mask was based on MODIS 250 m imagery in combination with Shuttle Radar Topography Mission (SRTM) Water Body Data (SWBD) to create a global map of surface water at 250 m spatial resolution [38].
Daily water level observations were used to evaluate the performance of the three satellite altimetry measurements. The in situ gauged data for three lakes (i.e., Tai Lake, Chao Lake, and South Dongting Lake) were available, which were collected from the Hydrological Yearbook published by the Chinese Ministry of Water Resources. It is noted that the in-situ data were referred to a different height system from the satellite altimetry measurements.
The year-end reservoir storage data and human water consumption data of MLYB were collected from Changjiang and Southwest Rivers Water Resources Bulletin to investigate the influences on variations of lake water levels.

Data Processing
In order to obtain accurate radar satellite measurements over lakes, the advanced retracking method was employed to reprocess the waveforms [39,40]. The Narrow Primary Peak Threshold (NPPT) [41] was used in this paper, which was an extension to the widely used OCOG (Off Center of Gravity) method. The NPPT was used because it can give the most stable results over inland water [40,42]. The leading and trailing edge of the primary peak was identified by the retracker, and the OCOG method [43] was then used around the start and end gates to compute the amplitude of extracted sub-waveform [42]. In order to determine the optimum threshold of retracker, three frequently-used values (i.e., 30%, 50% and 70%) were examined for Tai Lake water level estimation using Cryosat-2 data. The results were shown in Figure 2, which suggested that the threshold of 30% outperformed the other two values in both the Standard Deviation (SD) and number of available observations. Hence, it was used to determine the retracking points of Envisat and CryoSat-2.
Remote Sens. 2020, 12, x FOR PEER REVIEW 6 of 18 The year-end reservoir storage data and human water consumption data of MLYB were collected from Changjiang and Southwest Rivers Water Resources Bulletin to investigate the influences on variations of lake water levels.

Data Processing
In order to obtain accurate radar satellite measurements over lakes, the advanced retracking method was employed to reprocess the waveforms [39,40]. The Narrow Primary Peak Threshold (NPPT) [41] was used in this paper, which was an extension to the widely used OCOG (Off Center of Gravity) method. The NPPT was used because it can give the most stable results over inland water [40,42]. The leading and trailing edge of the primary peak was identified by the retracker, and the OCOG method [43] was then used around the start and end gates to compute the amplitude of extracted sub-waveform [42]. In order to determine the optimum threshold of retracker, three frequently-used values (i.e., 30%, 50% and 70%) were examined for Tai Lake water level estimation using Cryosat-2 data. The results were shown in Figure 2, which suggested that the threshold of 30% outperformed the other two values in both the Standard Deviation (SD) and number of available observations. Hence, it was used to determine the retracking points of Envisat and CryoSat-2.
For determination of the water levels of lakes, it is necessary to preprocess the ICESat-1/GLA14 data. The potential elevation anomalies resulting from cloud reflections and atmospheric noise during the time of data acquisition were removed, using a quality control analysis based on four flags provided by NSDIC [44]. In addition, as elevation measurements from the three satellite altimetry datasets were referenced to different ellipsoids and vertical datum, the elevation data were converted to a common the reference system, horizontally in WGS84 and vertically in EMG96 (the Earth Gravitational Model 1996), using the methods provided in [13]. The 250 m MODIS land-water mask (MOD44W) was then applied to the three satellite altimetry datasets to obtain the footprints of lakes in the MLYB. The resolution of the mask was sufficient because only lakes and reservoirs with areas greater than 5 km 2 were examined in this study. The existence of the lakes identified from the MOD44W mask was further checked with Google Earth and appropriate Landsat images, which resulted in ~200 lakes/reservoirs with areas over 5 km 2 in the MLYB. The masked lake level observations were subject to outlier elimination in the next step.

Outlier Removal
Satellite altimetry observations are still quite noisy even after the retracking corrections and geophysical corrections. For example, Figure 3 showed the footprints of three satellite altimetry datasets over the Tai Lake, which illustrated that most campaigns were affected by anomalies. Hence, Outlier removal should be performed before computing the mean lake level. Standard outlier removal procedures, like RANSAC (Random Sample Consensus algorithm), were not effective, since For determination of the water levels of lakes, it is necessary to preprocess the ICESat-1/GLA14 data. The potential elevation anomalies resulting from cloud reflections and atmospheric noise during the time of data acquisition were removed, using a quality control analysis based on four flags provided by NSDIC [44]. In addition, as elevation measurements from the three satellite altimetry datasets were referenced to different ellipsoids and vertical datum, the elevation data were converted to a common the reference system, horizontally in WGS84 and vertically in EMG96 (the Earth Gravitational Model 1996), using the methods provided in [13].
The 250 m MODIS land-water mask (MOD44W) was then applied to the three satellite altimetry datasets to obtain the footprints of lakes in the MLYB. The resolution of the mask was sufficient because only lakes and reservoirs with areas greater than 5 km 2 were examined in this study. The existence of the lakes identified from the MOD44W mask was further checked with Google Earth and appropriate Landsat images, which resulted in~200 lakes/reservoirs with areas over 5 km 2 in the MLYB. The masked lake level observations were subject to outlier elimination in the next step.

Outlier Removal
Satellite altimetry observations are still quite noisy even after the retracking corrections and geophysical corrections. For example, Figure 3 showed the footprints of three satellite altimetry datasets over the Tai Lake, which illustrated that most campaigns were affected by anomalies. Hence, Outlier removal should be performed before computing the mean lake level. Standard outlier removal procedures, like RANSAC (Random Sample Consensus algorithm), were not effective, since we were dealing with large percentages of outliers. Therefore, a tailored outlier removal procedure based on the methods of Kleinherenbrink et al. [15] and Phan et al. [13] was proposed as follows: Remote Sens. 2020, 12, x FOR PEER REVIEW 7 of 18 we were dealing with large percentages of outliers. Therefore, a tailored outlier removal procedure based on the methods of Kleinherenbrink et al. [15] and Phan et al. [13] was proposed as follows: The first step is the estimation of the mode, which is the elevation bin where most measurements belong to. The step size is 1 m, large enough to cope with slopes in the water level due to various causes. Measurements outside the one meter range from the mode are considered as outliers and removed. In the second step, the mean lake level is calculated using the remaining elevations, and the Root Mean Square (RMS) of residuals between the elevations and the mean lake level is also computed. Those elevations will be regarded as outliers and removed if the absolute differences between the elevations and the mean lake level are greater than the RMS. It is noted that if the RMS of ICESat-1 altimetry data is less than 15 cm, it is set to be 15 cm, and if the RMS of the two radar altimetry data is less than 20 cm, it is set to be 20 cm (Figure 3b,c,d), considering the precisions of two types of altimetry data and high roughness of lakes' surface in winter [13]. In the third step, the outlier removal is continued by repeating the second step based on the remaining elevations. It will stop if the number of remaining elevations is less than 10 or 40% of all elevations in the lake track. The lakes' water level is then obtained by averaging all the remaining elevations.

Inter-altimeter bias Adjustment and Accuracy Assessment
Bias adjustment is required for monitoring lake levels using multi-altimeter data, as there are systematic errors among different satellite altimeters. The most common method for bias adjustment was relative calibration among the altimeters [21]. It required different satellites to fly over the same lake at the same time, which was however not quite possible. Hence, the gauged daily water level data were used to facilitate the calibration. In the calibration, we assumed the lake surface is a flat plane and the satellite A and B fly over the lake at the same time. The bias (∆1 ) between the two satellites can be calculated as follows: where and are the observed lake levels by satellite A and B at j-th observation, and n is the number of observations. In fact, the satellites would not fly over the lake at the same time. Hence, ∆1 contains both the bias of inter-altimeters and the difference of lake levels (∆2 ) between the two observation time points. The latter can be obtained using the gauged lake levels.
where ℎ and ℎ are the gauged lake levels at the time points the satellite A and B fly over. Therefore, the inter-satellite bias can be determined by: The accuracy assessment of derived lake levels was carried out after the inter-altimeter calibration was performed. The assessment was conducted by comparing the satellite-derived lake levels against those of the gauged data. Due to the different reference systems between the gauged lake levels and altimetry data, a direct comparison of surface water elevation was not possible. Here, The first step is the estimation of the mode, which is the elevation bin where most measurements belong to. The step size is 1 m, large enough to cope with slopes in the water level due to various causes. Measurements outside the one meter range from the mode are considered as outliers and removed. In the second step, the mean lake level is calculated using the remaining elevations, and the Root Mean Square (RMS) of residuals between the elevations and the mean lake level is also computed. Those elevations will be regarded as outliers and removed if the absolute differences between the elevations and the mean lake level are greater than the RMS. It is noted that if the RMS of ICESat-1 altimetry data is less than 15 cm, it is set to be 15 cm, and if the RMS of the two radar altimetry data is less than 20 cm, it is set to be 20 cm (Figure 3b,c,d), considering the precisions of two types of altimetry data and high roughness of lakes' surface in winter [13]. In the third step, the outlier removal is continued by repeating the second step based on the remaining elevations. It will stop if the number of remaining elevations is less than 10 or 40% of all elevations in the lake track. The lakes' water level is then obtained by averaging all the remaining elevations.

Inter-altimeter bias Adjustment and Accuracy Assessment
Bias adjustment is required for monitoring lake levels using multi-altimeter data, as there are systematic errors among different satellite altimeters. The most common method for bias adjustment was relative calibration among the altimeters [21]. It required different satellites to fly over the same lake at the same time, which was however not quite possible. Hence, the gauged daily water level data were used to facilitate the calibration. In the calibration, we assumed the lake surface is a flat plane and the satellite A and B fly over the lake at the same time. The bias (∆1 j ) between the two satellites can be calculated as follows: where H A j and H B j are the observed lake levels by satellite A and B at j-th observation, and n is the number of observations. In fact, the satellites would not fly over the lake at the same time. Hence, ∆1 j contains both the bias of inter-altimeters and the difference of lake levels (∆2 j ) between the two observation time points. The latter can be obtained using the gauged lake levels.
where h A j and h B j are the gauged lake levels at the time points the satellite A and B fly over. Therefore, the inter-satellite bias can be determined by: The accuracy assessment of derived lake levels was carried out after the inter-altimeter calibration was performed. The assessment was conducted by comparing the satellite-derived lake levels against those of the gauged data. Due to the different reference systems between the gauged lake levels and altimetry data, a direct comparison of surface water elevation was not possible. Here, the correlation coefficient (R) between gauged data and altimetry-extracted data and the significant level (P) were calculated to indicate the accuracy. These parameters have also been used for accuracy assessment in similar studies due to inconsistent reference systems [9,22].

Lake Water Level Trend Estimation
The linear regression model was used to estimate the time series lakes' level change rates. The model was selected because it was less sensitive to random noises especially in case of the small number of footprints [22], and was also employed by many other previous studies [8,13,22,30]. Considering the estimates of a linear trend and annual amplitude, only lakes with at least 6 campaigns were included in the estimation. Hence, the rate represented the average change in lake level over the satellite acquisition period. The rate of lake level change was estimated by fitting a linear equation to the satellite observations using the least squares method [8], which was given as follows: where a is the slope which represents the annual rate of water level change, b is the intercept, x is the time (decimal year), and y is the time series water level of lakes.

Inter-altimeter Calibration and Accuracy Assessment
In this study, the Tai lake was chosen for inter-altimeter calibration as the daily water level of the lake was available and the surface of the lake is flat with no gradient. As a result, 19 paired observations for Envisat and ICESat-1, and 17 paired observations for Envisat and Cryosat-2 were obtained, which were used to calculate the inter-satellite bias. Figure 4 showed the results of inter-altimeter calibration. It can be seen that the average bias of ICESat-1 and Cryosat-2 with respect to Envisat were 6.7 cm and 3.1 cm, respectively. The SD of ICESat-1 and Cryosat-2 with respect to Envisat were 12.1 cm and 3.7 cm, respectively. The obtained biases were used to correct the water level measurements by the two satellite altimeters to obtain a consistent water level measurement.
Three lakes were used for the accuracy assessment due to the availability of gauged lake level data. Figure 5 showed the correlation between the satellite altimeter measurements and the gauged water levels for the Tai Lake, Chao Lake and Dongting Lake. It was found that the satellite altimetry measurements were significantly correlated with gauged data (P < 0.001) with the correlation coefficients (R) ranging from 0.908 to 0.989, suggesting a high consistency between the two datasets.
observations for Envisat and ICESat-1, and 17 paired observations for Envisat and Cryosat-2 were obtained, which were used to calculate the inter-satellite bias. Figure 4 showed the results of interaltimeter calibration. It can be seen that the average bias of ICESat-1 and Cryosat-2 with respect to Envisat were 6.7 cm and 3.1 cm, respectively. The SD of ICESat-1 and Cryosat-2 with respect to Envisat were 12.1 cm and 3.7 cm, respectively. The obtained biases were used to correct the water level measurements by the two satellite altimeters to obtain a consistent water level measurement. Three lakes were used for the accuracy assessment due to the availability of gauged lake level data. Figure 5 showed the correlation between the satellite altimeter measurements and the gauged water levels for the Tai Lake, Chao Lake and Dongting Lake. It was found that the satellite altimetry measurements were significantly correlated with gauged data (P < 0.001) with the correlation coefficients (R) ranging from 0.908 to 0.989, suggesting a high consistency between the two datasets. The three satellite-derived water level time series were compared with those of the gauged. Figure 6 showed the two lake level time series for the Tai Lake. It can be seen that the satellite-derived lake level series followed the gauged data quite well (R = 0.910, P < 0.001), suggesting the good performance of the integration of multiple altimeters. While the obvious offset was ascribed to the different elevation reference systems between the satellite altimetry and gauged data. With the integration of multiple altimeters, sufficient measurements were obtained throughout the study period, allowing to capture the seasonal and annul cycles of lake level variations. The constructed long-term lake level time series provided a valuable tool for the hydrological and climatic studies. Correlation of Envisat, ICESat-1, Cryosat-2 measurements and the three combined measurements (x-axis) with the in-situ lake level data (y-axis) for (a-d) Tai Lake, (e-h) Chao Lake, and (i-l) South Dongting Lake, respectively. The correlation coefficient (R) and significant level (P) are also given. N is the number of observations. The three satellite-derived water level time series were compared with those of the gauged. Figure 6 showed the two lake level time series for the Tai Lake. It can be seen that the satellite-derived lake level series followed the gauged data quite well (R = 0.910, P < 0.001), suggesting the good performance of the integration of multiple altimeters. While the obvious offset was ascribed to the different elevation reference systems between the satellite altimetry and gauged data. With the integration of multiple altimeters, sufficient measurements were obtained throughout the study period, allowing to capture the seasonal and annul cycles of lake level variations. The constructed long-term lake level time series provided a valuable tool for the hydrological and climatic studies. also given. N is the number of observations. The three satellite-derived water level time series were compared with those of the gauged. Figure 6 showed the two lake level time series for the Tai Lake. It can be seen that the satellite-derived lake level series followed the gauged data quite well (R = 0.910, P < 0.001), suggesting the good performance of the integration of multiple altimeters. While the obvious offset was ascribed to the different elevation reference systems between the satellite altimetry and gauged data. With the integration of multiple altimeters, sufficient measurements were obtained throughout the study period, allowing to capture the seasonal and annul cycles of lake level variations. The constructed long-term lake level time series provided a valuable tool for the hydrological and climatic studies.

Lake Level Variation from 2002 to 2017
We initially tried to examine the water level changes of lakes with area greater than 5 km 2 in the MLYB from 2002 to 2017. After data processing, there were only 52 lakes and 18 reservoirs with more than six campaigns during each of the two sub-periods (2002-2009 and 2010-2017). Water-level change rates were determined for these lakes and reservoirs, as shown in Figure 7 and Table S1. Due to the combination of datasets from the three altimeters, sufficient campaigns were obtained for lakes, especially those large lakes. For example, the Poyang Lake, with an area of 2835.484 km 2 has 319 campaigns during 2002-2017 (Table S1), allowing an accurate estimation of lake-level variation.

Lake Level Variation from 2002 to 2017
We initially tried to examine the water level changes of lakes with area greater than 5 km 2 in the MLYB from 2002 to 2017. After data processing, there were only 52 lakes and 18 reservoirs with more than six campaigns during each of the two sub-periods (2002-2009 and 2010-2017). Water-level change rates were determined for these lakes and reservoirs, as shown in Figure 7 and Table S1. Due to the combination of datasets from the three altimeters, sufficient campaigns were obtained for lakes, especially those large lakes. For example, the Poyang Lake, with an area of 2835.484 km 2 has 319 campaigns during 2002-2017 (Table S1), allowing an accurate estimation of lake-level variation.
From Figure 7 and Table S1, It was clear that most observed lakes showed an increasing trend. To be specific, ~79% of lakes (41/52) revealed rising trend (18 are significant at 90% Confidence Level (CL)) with rates up to 0.203 m/y, while the remaining ~21% of lakes (11/52) exhibited falling trend (six are significant at 90% CL) with rates ranging from −0.102 to −0.007 m/y. However, the annual change rate of water level in most lakes was less than 10 cm/y. Their water levels seemed to be almost stable. In contrast, most of the reservoirs have the greater annual water level variations as they are controlled by humans. For example, the Dongjiang Reservoir took the first place in the annual water level change, with an increasing rate of 0.849 m/y. In specific, ~61% of observed reservoirs (11/18) exhibited rising rates (0.046 to 0.849 m/y) (7 are significant at 90% CL) while the remaining ~39% of reservoirs (7/18) showed falling rates (−0.513 to −0.104 m/y) (six are significant at 90% CL) in water levels. Spatially, the lakes close to the YR showed slight changes while the reservoirs far from the YR exhibited greater dynamics. Both the lakes and reservoirs showed varying annual changes with relatively small R 2 , less than 0.62.

Lake Level Variation in Different Periods
In order to examine the temporal characteristics of lake level variations, we divided the study period into two sub-periods:   Figure 7 and Table S1, It was clear that most observed lakes showed an increasing trend. To be specific,~79% of lakes (41/52) revealed rising trend (18 are significant at 90% Confidence Level (CL)) with rates up to 0.203 m/y, while the remaining~21% of lakes (11/52) exhibited falling trend (six are significant at 90% CL) with rates ranging from −0.102 to −0.007 m/y. However, the annual change rate of water level in most lakes was less than 10 cm/y. Their water levels seemed to be almost stable. In contrast, most of the reservoirs have the greater annual water level variations as they are controlled by humans. For example, the Dongjiang Reservoir took the first place in the annual water level change, with an increasing rate of 0.849 m/y. In specific,~61% of observed reservoirs (11/18) exhibited rising rates (0.046 to 0.849 m/y) (7 are significant at 90% CL) while the remaining~39% of reservoirs (7/18) showed falling rates (−0.513 to −0.104 m/y) (six are significant at 90% CL) in water levels. Spatially, the lakes close to the YR showed slight changes while the reservoirs far from the YR exhibited greater dynamics. Both the lakes and reservoirs showed varying annual changes with relatively small R 2 , less than 0.62.

Lake Level Variation in Different Periods
In order to examine the temporal characteristics of lake level variations, we divided the study period into two sub-periods:  Figure 8 illustrated the spatial pattern of lake level changes in the two sub-periods. It can be seen that the majority of lakes showed decreasing rates in water levels during 2002-2009. The lakes exhibited slight changes were those close to the YR, while those with great fluctuations were situated far from the YR, and most of which were reservoirs (Figure 8a). Specifically,~60% of lakes (32/53) showed falling rates (−0.215 to −0.004 m/y) (15 are significant at 90% CL) in water levels and~40% of lakes (21/53) exhibited rising rates (0.001 to 0.158 m/y) (four are significant at 90% CL), which indicated a generally declining trend of lake levels in the MLYB during this sub-period (Figure 8a). Meanwhile, nine reservoirs showed rising rates (0.158 to 1.568 m/y) (three are significant at 90% CL) and nine reservoirs exhibited falling rates (−0.818 to −0.004 m/y) (one is significant at 90% CL). Figure 8b showed the pattern of lake level changes during 2010-2017. It can be seen that most of the lakes showed an increasing trend in water levels. In specific,~66% of lakes (79/119) showed an increasing rate (0.002 to 0.481 m/y) (17 are significant at 90% CL) while~34% of lakes (40/119) showed a decreasing rate (−0.281 to −0.003 m/y) (four are significant at 90% CL) in lake levels. Meanwhile,~69% of reservoirs (36/52) showed increasing rates (0.037 to 4.21 m/y) (seven are significant at 90% CL) and~31% of reservoirs (16/52) showed falling rates (−0.790 to −0.045 m/y) (four are significant at 90% CL). Similarly, the lakes adjacent to the YR generally exhibited a slight increasing trend in water levels, while the reservoirs far from the YR displayed greater water level variations. The water level time series derived from multi-altimeters provided more than one-decade record of lake dynamics, revealing different evolution characteristics of lakes that cannot be captured by a single altimeter.   Figure 8 illustrated the spatial pattern of lake level changes in the two sub-periods. It can be seen that the majority of lakes showed decreasing rates in water levels during 2002-2009. The lakes exhibited slight changes were those close to the YR, while those with great fluctuations were situated far from the YR, and most of which were reservoirs (Figure 8a). Specifically, ~60% of lakes (32/53) showed falling rates (−0.215 to −0.004 m/y) (15 are significant at 90% CL) in water levels and ~40% of lakes (21/53) exhibited rising rates (0.001 to 0.158 m/y) (four are significant at 90% CL), which indicated a generally declining trend of lake levels in the MLYB during this sub-period (Figure 8a). Meanwhile, nine reservoirs showed rising rates (0.158 to 1.568 m/y) (three are significant at 90% CL)

Performance of Multi-Satellite Altimetry
The retracking is of key importance for the accuracy of satellite-derived lake level estimation. Over the inland water bodies the waveforms are generally contaminated by noise resulting from multiple land returns such as vegetation, bare sands, or steep shorelines [45]. The shape of return waveforms does not correspond to the theoretical Brown model showing, e.g., peaks or fast decaying trailing edges. Hence, the onboard tracking algorithm is not always able to position the 50% point in the nominal tracking gate. The retracking method was used to reprocess waveforms to improve the accuracy of lake elevation estimates. Several retracking algorithms were available, where each one displays unique advantages and disadvantages. For example, Uebbing et al. [45] compared six retrackers as for the reprocessing of the Jason-1/2 waveforms over the Volta Lake, Africa, and showed that the improved threshold (IT) and the Institute for Theoretical Geodesy (ITG) retrackers had better performance, with SD values of lake footprint elevations of only~10 cm. Nielsen et al. [40] evaluated the performance of four empirical retrackers for Cryosat-2 waveform processing over five lakes, and reported that the NPPT outperformed the others, with the SD values of lake footprint elevations of less than 10 cm. There exists no single algorithm that meets the diverse needs of the whole application range. In this work, the employed NPPT retracker improved the accuracy of lake footprints elevations, with SD value lower than 10 cm, which seemed to work well over the lakes of the study area.
The inter-altimeter calibration is critical for the construction of a consistent long-term water level record. For the calibration, the high redundancy was preferred, which can be achieved by the double measurements of the lake surface in all crossing tracks [24]. In this work, the duration of overlapping operation for ICESat-1 and Envisat, Cryosat-2 and Envisat were six years and two years, leading to the double observations of 19 and 17, respectively. The number of measurements seemed to be adequate to calculate the inter-altimeter biases as a similar number was used in the previous studies. For example, in the inter-altimeter calibration over the lake Ngangzi Co, eight paired-sample observations for T/P and Jason-1, and 18 paired-sample observations for Jason-2 and Jason-3 were used [21]. The identified biases of 6.7 and 3.1 cm were also comparable to similar works. For example, Wang et al. [21] reported that Jason-1 had a mean lake level bias of 2.5 ± 5.4 cm with respect to T/P, and the bias of Jason-3 with respect to Jason-2 is −2.8 ± 3.5 cm after removing an outlier. Bosch et al. [24] found the Envisat mission exhibited a significant offset of about 3 cm between the range biases of Side A and Side B, and the estimated range biases of individual periods of ICESat-1 varied significantly, due to operating only in short episodic laser periods. This may also explain the calculated SD value of Envisat-ICESat-1 was larger than that of Envisat-Cryosat-2 ( Figure 4).
The high accurate satellite-derived lake level estimations were obtained after these processing were applied. In this work, the satellite-derived lake levels was significantly correlated to the gauged data with the R value ranging from 0.908 to 0.989 ( Figure 5). The similar level of accuracy was also reported by previous works. For example, Song et al. [22] showed that the ICESat-1 and CryoSat-2 altimetry measurements were significantly correlated with gauge-based data, with the correlation coefficients (R) of 0.97 and 0.91 for Yamzhog Yumco, Tibetan. Qi et al. [9] found a high determination coefficient (R 2 = 0.95) between river gauge network data and Envisat derived water levels in the Great Brahmaputra River. Moreover, the altimeter-derived lake level change trends can also be confirmed by other observations. For example, Wang et al. [30] investigated water level changes in Chinese large lakes during 2003-2009, which found that several large lakes (i.e., Poyang Lake, Longgan Lake) in the MLYB showed a decreasing trend in lake levels. Jiang et al. [1] examined the water level and freshwater resources of over 1000 lakes and reservoirs in China from 2010 to 2016. It was found that most of the lake in the MLYB exhibited a rising trend. Overall, our results were consistent with the previous studies in both the accuracy and lake level change trend, suggesting the good performance of multi-satellite altimetry in the long-term monitoring of lake levels.
The integration of multiple satellite altimeters allows for the reconstruction of long time series of lake levels. To derive a consistent long-term data record, altimeter systems should operate in a sequence of different, partly overlapping time spans [20]. In this work, the employed ICESat-1, Envisat and CryoSat-2 have operated during [2003][2004][2005][2006][2007][2008][2009][2002][2003][2004][2005][2006][2007][2008][2009][2010][2011][2012], and 2010-current, respectively, making them temporally optimal combination to derive a long-term record. The strength of ICESat-1 was the small footprint size and high ranging accuracy and the only weakness perhaps was the coarse temporal resolution [12]. While the Envisat had a repeat period of 35 days, and a small cross-track distances, allowing the satellite to observe more lakes in a relatively higher temporal resolution [14]. The combination of these two altimeters enabled more campaigns over the lakes during 2002-2012, comparable to those from the CryoSat-2 during 2010-2017. The CryoSat-2 altimetry can therefore be considered as a reliable extension of ICESat-1 and Envisat altimetry missions for monitoring lake level changes ( Figure 6). The integration of these three altimeters allowed the construction of 2002-2017 measurements of lake levels with a high temporal resolution, demonstrating the potential for the long-term observation of lakes ( Figure 6). It was especially useful for the lakes in remote areas where no gauging stations were available.

Impact Factors of Water Level Variations
Precipitation and evaporation are the two direct climatic factors affecting the lake levels [8]. Figure 9 was the annual precipitation in the MLYB based on monthly TRMM data over the study period. It can be seen that the amount of precipitation was increasing (12. Although the annual precipitation fluctuated during the study period, an obvious increase in the average precipitation can be identified between the two sub-periods. This was consistent with the observed lake level trends of the two corresponding sub-periods ( Figure 7). Thus, it suggested that the observed lake level changes can be attributed to the precipitation variations during the study period.
sequence of different, partly overlapping time spans [20]. In this work, the employed ICESat-1, Envisat and CryoSat-2 have operated during [2003][2004][2005][2006][2007][2008][2009][2002][2003][2004][2005][2006][2007][2008][2009][2010][2011][2012], and 2010-current, respectively, making them temporally optimal combination to derive a long-term record. The strength of ICESat-1 was the small footprint size and high ranging accuracy and the only weakness perhaps was the coarse temporal resolution [12]. While the Envisat had a repeat period of 35 days, and a small crosstrack distances, allowing the satellite to observe more lakes in a relatively higher temporal resolution [14]. The combination of these two altimeters enabled more campaigns over the lakes during 2002-2012, comparable to those from the CryoSat-2 during 2010-2017. The CryoSat-2 altimetry can therefore be considered as a reliable extension of ICESat-1 and Envisat altimetry missions for monitoring lake level changes ( Figure 6). The integration of these three altimeters allowed the construction of 2002-2017 measurements of lake levels with a high temporal resolution, demonstrating the potential for the long-term observation of lakes ( Figure 6). It was especially useful for the lakes in remote areas where no gauging stations were available.

Impact Factors of Water Level Variations
Precipitation and evaporation are the two direct climatic factors affecting the lake levels [8]. Figure 9 was the annual precipitation in the MLYB based on monthly TRMM data over the study period. It can be seen that the amount of precipitation was increasing (12. Although the annual precipitation fluctuated during the study period, an obvious increase in the average precipitation can be identified between the two sub-periods. This was consistent with the observed lake level trends of the two corresponding sub-periods ( Figure 7). Thus, it suggested that the observed lake level changes can be attributed to the precipitation variations during the study period.  A quantitative analysis between lake level and sub-basin precipitation was carried out for 16 large lakes in the MLYB, as shown in Figure 10. It can be seen that the annual water levels of 14 of 16 large lakes showed a good correlation with annual precipitation of their sub-basins, with R > 0.5. It was noted that most of these lakes showed significant relationships (P < 0.01), suggesting the important role of precipitation in lake level variations. In addition to precipitation, the evaporation also impacted the lake levels, which was closely associated with air temperature. It was reported that there was no obvious climate warming occurred over the MLYB in recent years [30], and therefore, the evaporation did not change much over the study period. Hence, the precipitation variations primarily contributed to the lake level fluctuations in the study area.
large lakes showed a good correlation with annual precipitation of their sub-basins, with R > 0.5. It was noted that most of these lakes showed significant relationships (P < 0.01), suggesting the important role of precipitation in lake level variations. In addition to precipitation, the evaporation also impacted the lake levels, which was closely associated with air temperature. It was reported that there was no obvious climate warming occurred over the MLYB in recent years [30], and therefore, the evaporation did not change much over the study period. Hence, the precipitation variations primarily contributed to the lake level fluctuations in the study area. Anthropogenic activities, such as the operation of dams or reservoirs and human water consumption, had direct impacts on lake level variations. There are about 1300 dams in the Yangtze River Basin [46]. Figure 11a showed the annually year-end water storage of the Three Gorges Reservoir (TGR) and the other large and medium reservoirs in the MLYB. It can be seen that the storage of TGR was increased from ~15 km 3 in 2003 to ~40 km 3 in 2017, with a rate of 2.01 km 3 /y. And the other reservoirs in the MLYB experienced a steady rise in water storage, with a rate of 2.61 km 3 /y. Anthropogenic activities, such as the operation of dams or reservoirs and human water consumption, had direct impacts on lake level variations. There are about 1300 dams in the Yangtze River Basin [46]. Figure 11a showed the annually year-end water storage of the Three Gorges Reservoir (TGR) and the other large and medium reservoirs in the MLYB. It can be seen that the storage of TGR was increased from~15 km 3 in 2003 to~40 km 3 in 2017, with a rate of 2.01 km 3 /y. And the other reservoirs in the MLYB experienced a steady rise in water storage, with a rate of 2.61 km 3 /y. Those dams usually store water from August to October and release water from November to May [46,47]. Due to complex interactions of rivers and lakes in the MLYB, these dams or reservoirs seasonally altered the water level of downstream rivers through seasonal regulation of discharges, which may impact the water levels of lakes. For example, during the impoundment of TGR, water levels along the downstream YR decreased by an average of~0.2-1.0 m [48]. The drop of river water level would increase the outflow from lakes to the river, leading to the lakes level decrease. It was reported that the Dongting Lake and Poyang Lake have experienced low water level during the impoundment season of TGR as a result of the increased outflow [47,49]. Hence the operation of dams had a seasonal influence on water levels of those class I lakes.
rate of 1.67 km 3 /y. Unlike the TGD's influence, the impact of human water consumption stayed negative in all seasons, contributing to the decrease of lake levels. In short, anthropogenic activities can affect the lake levels in different directions, which may further complicate lake level variations in addition to climatic forcing, possibly leading to small R 2 in water level changes. A clear knowledge of the driving mechanism of lake dynamics still requires more and longer extensive observations of lakes, climate, and human activities.

Conclusions
Satellite altimetry has been effectively used for monitoring the lake level changes in recent years. This work has demonstrated the potential of integrating multiple altimetry datasets derived from ICESat-1, Envisat and Cryosat-2 for the long-term (2002-2017) observation of lake level variations in the MLYB. The inter-altimeter calibration indicated that there were a few centimeter biases among these three altimeters, which were corrected to obtain a consistent measurement. The performance assessment suggested that the satellite altimetry measurements were significantly correlated with the gauge-based data, and the combination of three altimetry datasets produced unprecedented temporal and spatial resolutions for the long-term observation of lakes. The results of lake level observation in the MLYB showed that ~79% of measured lakes displayed an increasing rate in water levels during 2002-2017. Lake levels varied regionally. Lakes along the YR exhibited slight rising rates. In contrast, those reservoirs far from the YR showed marked fluctuations. Lake levels varied temporally as well. Most of the lakes showed declining trends in water levels during 2002-2009, and rising trends during 2010-2017. Most reservoirs showed rapidly rising rates during the study period. The different change patterns of lake levels were closely linked to the spatial heterogeneity of precipitation variations. Anthropogenic activities, such as the operation of dams or reservoirs, human water consumption also impacted the lake level variations. This research work reported herein is an important step towards integrating the wealth of multiple satellite altimeters to provide continuous long-term observations on lake level dynamics in the MLYB. Based on our results, we recommend Human water consumption was another important factor affecting the water levels of lakes. Figure 11b showed the amount of annual human water consumption in the MLYB. It can be seen that the total human water consumption increased from 131.99 km 3 in 2003 to 160.01 km 3 in 2017, with a rate of 1.67 km 3 /y. Unlike the TGD's influence, the impact of human water consumption stayed negative in all seasons, contributing to the decrease of lake levels. In short, anthropogenic activities can affect the lake levels in different directions, which may further complicate lake level variations in addition to climatic forcing, possibly leading to small R 2 in water level changes. A clear knowledge of the driving mechanism of lake dynamics still requires more and longer extensive observations of lakes, climate, and human activities.

Conclusions
Satellite altimetry has been effectively used for monitoring the lake level changes in recent years. This work has demonstrated the potential of integrating multiple altimetry datasets derived from ICESat-1, Envisat and Cryosat-2 for the long-term (2002-2017) observation of lake level variations in the MLYB. The inter-altimeter calibration indicated that there were a few centimeter biases among these three altimeters, which were corrected to obtain a consistent measurement. The performance assessment suggested that the satellite altimetry measurements were significantly correlated with the gauge-based data, and the combination of three altimetry datasets produced unprecedented temporal and spatial resolutions for the long-term observation of lakes. The results of lake level observation in the MLYB showed that~79% of measured lakes displayed an increasing rate in water levels during 2002-2017. Lake levels varied regionally. Lakes along the YR exhibited slight rising rates. In contrast, those reservoirs far from the YR showed marked fluctuations. Lake levels varied temporally as well. Most of the lakes showed declining trends in water levels during 2002-2009, and rising trends during 2010-2017. Most reservoirs showed rapidly rising rates during the study period. The different change patterns of lake levels were closely linked to the spatial heterogeneity of precipitation variations. Anthropogenic activities, such as the operation of dams or reservoirs, human water consumption also impacted the lake level variations. This research work reported herein is an important step towards integrating the wealth of multiple satellite altimeters to provide continuous long-term observations on lake level dynamics in the MLYB. Based on our results, we recommend the integration of altimetry data from ICESat-1, Envisat, and CryoSat-2 can be used for the long-term observations of lakes in remote areas where no gauging stations were available.
Author Contributions: P.L. was responsible for the processing, analysis, and interpretation of data and wrote the manuscript; H.L. led the conception, design and supervision of the work, and revised the manuscript; F.C. was involved in discussion and editing; X.C. was involved in review and editing. All authors have read and agreed to the published version of the manuscript.
Funding: This research was funded in part by the National Natural Science Foundation of China, grant numbers 41201429 and 41391240191, and the Chinese Geological Survey (DD20190263).