High-temporal-resolution water level and storage change data sets for lakes on the Tibetan Plateau during 2000–2017 using multiple altimetric missions and Landsat-derived lake shoreline positions

The Tibetan Plateau (TP), known as Asia’s water tower, is quite sensitive to climate change, which is reflected by changes in hydrologic state variables such as lake water storage. Given the extremely limited ground observations on the TP due to the harsh environment and complex terrain, we exploited multiple altimetric missions and Landsat satellite data to create high-temporal-resolution lake water level and storage change time series at weekly to monthly timescales for 52 large lakes (50 lakes larger than 150 km2 and 2 lakes larger than 100 km2) on the TP during 2000–2017. The data sets are available online at https://doi.org/10.1594/PANGAEA.898411 (Li et al., 2019). With Landsat archives and altimetry data, we developed water levels from lake shoreline positions (i.e., Landsat-derived water levels) that cover the study period and serve as an ideal reference for merging multisource lake water levels with systematic biases being removed. To validate the Landsat-derived water levels, field experiments were carried out in two typical lakes, and theoretical uncertainty analysis was performed based on high-resolution optical images (0.8 m) as well. The RMSE of the Landsat-derived water levels is 0.11 m compared with the in situ measurements, consistent with the magnitude from theoretical analysis (0.1–0.2 m). The accuracy of the Landsat-derived water levels that can be derived in relatively small lakes is comparable with most altimetry data. The resulting merged Landsat-derived and altimetric lake water levels can provide accurate information on multiyear and short-term monitoring of lake water levels and storage changes on the TP, and critical information on lake overflow flood monitoring and prediction as the expansion of some TP lakes becomes a serious threat to surrounding residents and infrastructure.


Introduction
The Tibetan Plateau (TP), providing vital water resources for more than a billion population in Asia, is a sensitive region undergoing rapid climate change (Field et al., 2014). There are more than 1,200 alpine lakes larger than 1 km 2 on the TP where glaciers and permafrost are also widely distributed. With little disturbance by human activity in this area, lake storage changes may serve as an important indicator that reflects changes in regional hydrologic processes and responses to climate 30 change. Wang et al. (2018) showed that global endorheic basins are experiencing a decline in water storage whereas the endorheic basin on the TP is an exception. Given the fact that TP lakes have been expanding for more than 20 years (Pekel et al., 2016), quality data sets on lake water level and/or storage could be the basis for investigating its causes (e.g., climate change/variability) and interactions with the water/energy cycles and human society (e.g., increasing risks of inundation and overflow floods). 5 As an important component of the hydrosphere, terrestrial water cycle, and global water balance, millions of inland water bodies such as lakes, wetlands, and reservoirs have been investigated globally and their total storage was estimated to be 181.9 × 10 3 km 3 based on statistical models (Lehner and Döll, 2004;Messager et al., 2016;Pekel et al., 2016). Lake storage changes that play an important role in the regional water balance, can be derived from changes in lake water level and area (Frappart et al., 2005). Lake water levels and areas are mostly derived from satellite remote sensing due to the scarcity of in 10 situ data across the TP where the harsh environment and complex terrain make in situ measurements difficult to perform and costly (Cré taux et al., 2016;Song et al., 2013;Yao et al., 2018b;Zhang et al., 2017a). Lake water levels can be monitored using satellite altimeters initially designed for sea surface topography or ice sheet/sea ice freeboard height measurements.
Satellite altimeters determine the range between the nadir point and satellite by analysing the waveforms of reflected electromagnetic pulses . 15 There are two major categories of satellite altimeters, i.e., laser and radar. Laser altimeters, e.g., the Ice, Cloud, and land Elevation (ICESat) satellite, operating in the near infrared band has a smaller size of footprints and generally higher accuracy than radar altimeters, facilitating applications in glacier/ice mass balance studies (Neckel et al., 2014;Sandberg Sørensen et al., 2011). Rader altimeters, operating in the microwave band, have larger footprints and are more likely to be contaminated by signal from complex terrain when applied to inland water bodies. Nevertheless, it is possible to remove these impacts 20 with waveform retracking algorithms (Guo et al., 2009;Huang et al., 2018;Jiang et al., 2017). Zhang et al. (2011) mapped water level changes in 111 TP lakes for the 2003-2009 period using ICESat data that have a temporal resolution of 91 days.
ICESat data have relatively denser ground tracks but a lower temporal resolution than most of other altimetric missions. This means that ICESat covers more lakes but provide few water levels for each lake. After ICESat was decommissioned in 2010, CryoSat-2 data starting from 2010 were adopted in related studies (Jiang et al., 2017), due to its similar dense ground tracks 25 and competitive precision as ICESat. Other altimetric missions, such as Topex/Poseidon (T/P), Jason-1/2/3, the European Remote Sensing (ERS-1/2) satellite, and Envisat also have some but relatively limited applications in monitoring changes in lake water level on the TP due to sparse ground tracks. In this study, multisource altimetry data (i.e., Jason-1/2/3, Envisat, ICESat, and CryoSat-2) were combined if available for lakes in this study, with the Landsat-derived water levels developed in this study as a critical reference to increase the water level observations and merging data from multiple altimetric 30 missions.
Changes in lake area can be captured by optical or synthetic aperture radar (SAR) images from medium or high spatial resolution remote sensing data, such as Landsat and Sentinel series. Extraction of lake water bodies can be manually (Wan et al., 2016) or automatically (Zhang et al., 2017b) achieved. Automatic water extraction methods based on the water index and auto-thresholding are more efficient in dealing with a mass of remote sensing images. Even so, acquisition and preprocessing of such a large amount of historical data (~10TB) covering TP lakes are still intractable for researchers with limited computational resources. With the help of cloud-based platforms, such as the Google Earth Engine (GEE) that significantly reduces data downloading and preprocessing time, tens of thousands of images may be processed online in days instead of 5 months (Gorelick et al., 2017). In this study, more than twenty thousand Landsat images were processed online using GEE to extract lake water bodies based on the water index (McFeeters, 1996) and Otsu algorithm (Otsu, 1979).
There have been studies focusing on changes in lake water storage on the TP over recent decades, e.g., Zhang et al. (2017a) examined changes in water storage for ~70 lakes from the 1970s to 2015 with ICESat altimetry data and Landsat archives.
An individual lake area data set from the 1970s and annual area data after 1989 were used. Due to the short time span of 10 ICESat, they used the hypsometric method to convert lake areas into water levels. Yao et al. (2018b) used DEMs and optical images to develop hypsometric curves for lakes on the central TP, and estimated annual changes in water storage for 871 lakes from 2002 to 2015. These studies have a wide spatial coverage of lakes but relatively lower temporal resolution and little spaceborne altimetric information, which may limit the accuracy of trends in lake water level/storage in some cases and short-term monitoring of lake overflow floods. The Laboratoire d'Etudes en Gé ophysique et Océ anographie Spatiales 15 (LEGOS) Hydroweb provides a lake data set, including multisource altimetry-based changes in lake water level and storage as well as hypsometric curves for 22 TP lakes (Cré taux et al., 2016;Cré taux et al., 2011b). The data set incorporates more spaceborne altimetric information and has a higher temporal resolution. However, there may be a remaining bias when different sources of altimetric data are merged, due to the lack of some important reference that can be derived from optical remote sensing to be shown in this study. We term the reference data as the 'Landsat-derived water level' to be introduced in 20 section 3.2. Here, we list recent studies and data sets (Table 1) to provide a concise summary on remote sensing monitoring of water levels and storage changes over lakes on the TP. The overall objective of this study was to examine multiyear and short-term changes in water level and storage across 52 lakes with surface areas larger than 150 km 2 on the TP by merging multisource altimetry and optical remote sensing images to generate more coherent high temporal resolution lake water level and storage change data sets ranging from weekly to monthly timescales during 2000-2017 and the hypsometric curve (i.e., the lake area-water level relationship) for each lake.
To investigate changes in lake storage, lake water levels and areas need to be derived from multisource remote sensing. 5 First, water levels from various satellite altimeters ( Figure 1) for each lake as well as lake shoreline positions and lake areas from optical remote sensing images (i.e., Landsat) were derived. Second, systematic biases between different altimetry data were removed by either comparing the mean water levels during the overlap period ( Figure 1) or comparing the two water level time series with lake shoreline positions, depending on the length of the overlap period (details can be found in section 3.1). Lake shoreline position-derived water levels, termed as the Landsat-derived water levels in this study can serve as a 10 unique source of information reflecting water levels as well as a data merging reference. We will show that after deriving two or three regression parameters, lake shoreline positions can well reflect lake water levels with comparable accuracy as altimetry-derived water levels. Third, with information on lake water levels and areas derived from altimetry data and optical remote sensing images, the hypsometric curve that describes the relationship between the lake water level and storage changes can be derived. Fourth, the integral of the hypsometric curve was performed to convert lake water levels into storage 15 changes.
Results of this study provide a comprehensive and detailed assessment of changes in lake level and storage on the TP for the recent two decades, and short-term monitoring of lake overflow floods for some lakes. This study could largely benefit more detailed investigations into lakes, lake basins, and regional climate change, because the generated data sets have the highest temporal resolution during the study period with systematic biases well removed. To ensure the data quality, field experiments were carried out and in situ data were collected to examine the uncertainty in the Landsat-derived water levels. 5 Users are free to access the data set described in this paper at https://doi.org/10.1594/PANGAEA.898411 (Li et al., 2019). Figure 1. Spatial (the number of lakes covered) and temporal coverage and their overlap periods of multiple satellite altimetric missions used in this study, including Jason-1/2/3, Envisat, ICESat, and CryoSat-2.
2 Study area and data 10

Study area
The TP can be generally divided into 12 major basins (Wan et al., 2016;Zhang et al., 2013), among which the Inner/Central TP (CP) is the only endorheic basin and home to most TP lakes including ~300 TP lakes larger than 10 km 2 . Therefore, it was chosen as the main study area. The endorheic basin covers an area of ~710,000 km 2 (~28% of total TP) with a mean elevation of ~4,900 m and has a semi-arid plateau climate with annual precipitation ranging from 96-295 mm (Li et al., 15 2017c). Most lakes in the endorheic basin were expanding under the influence of climate change/variability as opposed to other areas in the TP, e.g., Selin Co exceeded Nam Co in area and consequently became the largest lake in the endorheic basin between 2011-2012 and expanded by 26% over the past 40 years (Zhou et al., 2015), whereas Yamzhog Yumco (outside the endorheic basin, 350 km to the southeast of Selin Co) shrunk by ~11% during 2002-2014 according to Wan et al. (2016). Located in the southeast endorheic basin, the Nam Co basin covering about 10,800 km 2 , with 19% of the basin lake 20 water area and a mean lake elevation of ~4,730 m was chosen as a field experiment spot. The mean annual temperature and precipitation of Nam Co are 1.3 ℃, and 486 mm, respectively (Li et al., 2017a). The other experiment spot was Yamzhog Yumco, which has a mean lake elevation of ~4,440 m. Subject to steep mountainous terrain, the lake has a narrow-strip shape with complex shorelines. The basin of Yamzhog Yumco covers ~6,100 km 2 , with mean annual temperature and precipitation of 2.8 ℃ and ~360 mm, respectively (Yu et al., 2011). An overall map of experiment lakes is given in Figure 2. Yumco is located in the Yarlung Zangbo River basin (the Upper Brahmaputra River). Both lakes are close to Lhasa City.

Data
Multisource altimetry data were used in this study as shown in Table 2. The earliest record dates back to 2002 (i.e., Envisat and Jason-1) and the latest record ends in 2017 (i.e., Jason-3 and CryoSat-2, Figure 1). Most of the 52 lakes examined in this study were covered by ICESat, Envisat, and CryoSat-2 data. ICESat data provided by the National Aeronautics and Space 10 Administration (NASA) were available on 42 lakes in this study. Envisat and CryoSat-2 data provided by the European Space Agency (ESA) were available on 35 and 51 lakes in this study. Jason-1/2/3 data provided by the Centre National d'Etudes Spatiales (CNES) were available only on 12 lakes in this study due to the relatively sparse ground tracks or data quality issues. Note that Jason-2 inherited the orbit of Jason-1 after its launch in 2008, whereas Jason-1 was shifted into an interleaved orbit and continued functioning until 2013, thereby increasing the spatial coverage of Jason altimetry series to some degree, e.g., Jason-1 data in Lake Qinghai, the largest lake on the TP, were only available after 2008 due to the orbit shift. ICESat and CryoSat-2 data have the largest spatial coverage but relatively long repeat cycles of 91 days and 369 days, respectively (Bouzinac, 2012;Zhang et al., 2011). The Envisat mission has a lower orbit than Jason-1/2/3 but higher than ICESat, resulting in a moderate spatial coverage and a temporal resolution of 35 days (Benveniste et al., 2002). To determine 5 if the altimetry data fall into the lakes, a lake shape data set generated by Wan et al. (2016) was used. An example of using the lake shape data set to determine altimetry data falling into the lake boundaries is given in Figure 3 (a), showing that data from all altimeters are available in Zhari Namco. *LSH stands for lake surface height. It should be noted that different altimeters vary with wavelengths of electromagnetic radiation and mechanisms. For instance, Jason-1/2/3 using the Ku and C bands and Envisat/RA-2 using the Ku and S bands work in the low resolution mode (LRM).
These dual-frequency radar altimeters can provide more accurate range corrections due to the ionospheric effect (Tournadre, 2004). The LRM is typical for the early version of satellite altimeters such as TOPEX/Poseidon. There are more advanced 15 modes, such as SAR and interferometric synthetic aperture radar (InSAR) for recent radar altimeters, which generally have smaller footprints than the LRM mode. CryoSat-2/SIRAL working at a single Ku band has three modes including LRM, SAR, and InSAR, which were designed to have an increasing resolution in turn and work in different zones. The InSAR mode uses interference phenomena so that shift of the nadir point across the track can be detected, improving the altimeter's performance on ice sheets with slopes (Bouzinac, 2012). The Geoscience Laser Altimeter System (GLAS) is the laser 20 altimeter carried by ICESat working in the near infrared band.
We used Landsat 5 TM (2000-2011), Landsat 7 ETM+(2000, and Landsat 8 OLI (2013-2017) surface reflectance data sets provided by GEE to generate information on lake shoreline positions and lake areas. Landsat 7 ETM+ was subject to sensor failure and all the Landsat 7 ETM+ images contain gaps after 2003 (Markham et al., 2004). There were more than 20,000 images processed and a half of them were excluded from the final results due to cloud contamination or gaps. We collected daily in situ water level measurements in Yamzhog Yumco for validation purposes with a pressure-type water level sensor. The in situ water level measurements spanned a half year from May to October 2018. We also performed unmanned aerial vehicle (UAV or drone) imaging over a 1 km-lake bank in Yamzhog Yumco and Nam Co for obtaining better knowledge on the experimental environment. 5 In addition, GaoFen-2 (GF-2, the China High Resolution Earth Observation System mission) images were used to performed a rigorous statistical analysis of uncertainty in the Landsat-derived water levels by taking the GF-2 derived lake shoreline positions as the ground truth to analyse the sub-pixel water area ratios of Landsat image pixels (see Section 4). GF-2 images have a spatial resolution of 0.8 m for the panchromatic band and pre-processing including orthorectification and radiometric calibration was performed. Before analysis, we performed an image to image registration with manually selected tie points 10 between GF-2 and corresponding Landsat 8 OLI images until the co-registration error reduced to ~2 m.

Satellite altimetry water level
The first step of deriving satellite altimetry water levels is to select correct ground tracks and valid footprints falling on the lakes. Because there is a random ground track shift at ~1 km in different cycles for most altimetry missions, it is uncertain that valid lake footprints can be obtained for each cycle, even though the nominal ground track seems to cross the lake. This 5 problem can be addressed by comparing geographic coordinates of the footprints with a lake shape data set (Wan et al., 2016). After picking out the valid footprints, the lake surface height can be calculated for each footprint. All radar altimetry data share a relationship: where LSH represents the lake surface height with respect to the geoid; Alt represents the altimeter height with respect to the reference ellipsoid; Range represents the distance between the altimeter and lake surface; and cor represents several range 10 corrections due to atmospheric effects, sensor design defects, or geophysical effects. Radar altimeters and laser altimeters need different corrections, given that they are working in different wavelengths and have different designs. For instances, corrections for radar altimeters include waveform retracking correction, wet/dry troposphere correction, ionosphere correction, pole/solid tide correction, and geoid correction. Laser altimeters also need atmospheric delay correction, geoid/pole tide correction, and geoid correction. Unlike radar altimeters, saturation correction instead of waveform 15 retracking correction is more important to laser altimeters.
The retracking correction plays an important role in removing the contamination of land signal when radar altimetry data are applied to inland water bodies. In this study, Jason-1/2/3 data were retracked using a classical waveform retracking algorithm, i.e., the improved threshold method (ITR), whereas CryoSat-2 data were retracked using the narrow primary peak threshold (NTTP) method (Birkett and Beckley, 2010;Cheng et al., 2010;Jain et al., 2015). Retracking corrections were not 20 performed for Envisat and ICESat data, because the Envisat/RA-2 product already provided a retracked range using the ICE-1 method and the ICESat GLAH14 product already included several corrections (such as saturation correction) that are sufficient for most applications including studies on inland water bodies (Zhang et al., 2011). The original idea of the NTTP, ICE-1, and ITR is quite similar. All of them adopt a threshold defined as the percentage of the waveform peak to determine the retracking gate, and then to convert the difference between the retracking gate and the nominal gate into range correction 25 by multiplying the gate range (c∆t/2, where c is the speed of light and ∆t is the time duration of a gate). The differences lie in the choice of thresholds as well as the calculation of waveform peaks. For instance, ICE-1 uses a 30% threshold whereas ITR uses a 50% threshold.
For each cycle of an altimeter, it is common that more than one footprint fall on a lake, thereby providing several lake surface height (LSH) observations on the same day. After removing outliers with the three-sigma rule, frequency 30 distributions of the LSHs from the same cycle were calculated. The mean value of the histogram bin with the highest frequency was selected to represent the LSH for the cycle. Meanwhile, the frequency of the chosen histogram bin was reserved to evaluate the data quality for the cycle, e.g., a cycle was marked as high quality if the frequency was higher than 0.8, moderate quality if it was only higher than 0.5, and poor quality if the frequency was lower than 0.5. The LSH from each cycle constituted the water level time series for a lake. LSHs that were marked as poor quality and obviously deviated 5 from the moving average were removed from the altimetry-based lake water level time series.
It is not uncommon that systematic biases exist in different altimetry data sets due to variations in orbit, the discrepancy between correction models, errors associated with sensors, and even the choice of the reference datum. After deriving lake water level time series for each altimeter, we first merged the Envisat and ICESat water levels if both were available for a lake, because they have the longest overlap period (Figure 1). We chose Envisat-derived water levels as the baseline and 10 removed the difference of the mean values of the two products during the overlap period from the ICESat data, because Envisat data are generally denser and longer than ICESat data. A similar process was applied to Jason-1/2/3, as there are two overlap periods connecting the three altimeters together. Figure 3 (b) shows a result of merged altimetry data when all sensors are available. There are tradeoffs between CryoSat-2 and Jason-2/3 data in terms of spatial coverage and time span.
CryoSat-2 data are available for all lakes in this study but they only have an overlap period with Jason-2/3 data, whereas 15 Jason-2/3 data are only available for 12 lakes. For most lakes without Jason-2/3 data, we merged CryoSat-2 data with either ICESat or Envisat using Landsat-derived water levels spanning from 2000 to 2017, because there is no overlap period between these altimetry water levels ( Figure 1). Details on how Landsat-derived water levels aid in merging the altimetry water levels are shown in section 3.2.

Landsat-derived water level 20
For most lake basins, it is possible to find a relatively flat portion of lake banks with an average slope of 1/30 or even smaller, where obvious interannual or intra-annual changes in lake shoreline positions can be detected using Landsat images (30 m).
These locations can be found by comparing lake images from the first year and the last year of the study period if the lake shows a clear expanding/shrinking trend. Otherwise, we can compare images acquired in early summer when the LSH is at lowest with those acquired in late autumn when the lake expands to its limit. In this study, we assumed that the selected lake 25 bank was flat enough such that the relationship between the lake water level and shoreline position can be depicted in a linear or quasi-linear (parabolic) way. Thus, we can transform the lake shoreline positions into Landsat-derived water levels by fitting with altimetry water levels. The validity of this assumption can be evaluated with the coefficient of determination (R 2 ) for each lake as shown in Table 3. For most lakes, the goodness of fit is higher than 0.7, suggesting the generally good fitting relationship between the lake water levels and shoreline positions. 30 Though there were ~500 Landsat images obtained for the selected lake banks during the study period, many of them were largely affected by cloud or cloud shadow. All the images were processed online using GEE application programming interface (API). Preprocessing such as radiometric calibration, atmospheric correction, as well as geometric correction was already performed in the production of the data sets. In addition, the failure of Landsat 7 sensor SLC left all the Landsat 7 ETM+ images with gaps after 2003 (Markham et al., 2004), making the available images even fewer. We managed to make use of some images with gaps in generating lake shoreline positions. By choosing the region of interest (ROI) that is parallel to the image gaps, we made most of the Landsat 7 ETM+ images useable. However, the width of ROI must be reduced to 5 avoid shifting gaps as shown in Figure 4 (b). The gaps may vary with time but are more like vibration around the mid-point.
The ROI did not fill the interval of gaps, because the wider the ROI, the higher possibility of shifting gaps cross it. Lake shoreline positions were characterized by water area ratios detected in the ROI. To automatically extract water areas from a mass of Landsat archives on GEE, the water index and Otsu threshold method were jointly used. We calculated the normalized difference water index (NDWI) and the modified normalized difference water index (MNDWI) of the images and compared their performance in different seasons. It was found that the MNDWI tends to be more sensitive to shallow water in summer, but is less effective than the NDWI when the lake bank is covered by snow in the cold season as shown in Figure 5 (c). Therefore, the two water indices were jointly used by applying the MNDWI to images acquired during May to October and applying the NDWI to the rest months. The NDWI and MNDWI can be calculated as follows (McFeeters, 1996;5 Xu, 2005) , and refer to surface reflectance of bands 2, 4, 5 for Landsat 5/7 TM/ETM+ images and bands 3, 5, 6 for Landsat 8 OLI images. After calculating the water index, the grayscale image was binarized using the Otsu method. If the selected ROI comprises ~50% water and ~50% land, the performance of the method is good, as the distribution of digital numbers of the grayscale 5 image is close to the assumption of the bimodal histogram implicit in the Otsu algorithm (Kittler and Illingworth, 1985;Otsu, 1979). The binarized images were further processed to provide the water area ratio in the ROI, which represents the lake shoreline position. The lake shoreline position time series were then converted into Landsat-derived water levels using linear regression or second-order polynomial fit with altimetry-derived water levels (Figure 4 (c)-(d)). For most cases, we only used linear regression, and only for 2 lakes with Jason-1/2/3 data we performed the second-order polynomial fit, because a 10 higher-order regression requires more input information to ensure the reliability. However, cloud, cloud shadow, and shifting gaps may contaminate the ROI and cause errors in the Landsat-derived water levels. Therefore, the QA band of the Landsat surface reflectance product was used to filter the images. Data were excluded if the fraction of the cloud or cloud shadowcovered area in the ROI was higher than 5%. For every Landsat 7 ETM+ image acquired after 2003, the pixel number of the ROI was counted and compared with those acquired before 2003. If the loss of pixels exceeded 2%, the ROI was considered 15 to be affected by gap and the data were consequently excluded from the subsequent analysis.
A critical function of Landsat-derived water levels is to aid in merging altimetry water levels when there was no overlap period between altimeters or the overlap period was too short. For lakes without Jason-1/2/3 data, lake shoreline positions were firstly translated into Landsat-derived water levels by fitting with CryoSat-2 data functioning as extrapolation of CryoSat-2 to 1-2 years. Then, we applied the same method of merging Jason-1/2/3 to merge the extrapolated CryoSat-2 data 20 with either Envisat or ICESat data. In doing so, we were able to remove all systematic biases between multisource altimetry water levels. After merging the altimetry water levels, we performed regression analysis for the second time between the Landsat-derived water levels and merged altimetry water levels to check if the linear relationship is stable during the entire study period and at different elevations. If the linear relationship was stable, the Landsat-derived water levels would be merged with the altimetry water levels using the linear fitting parameters from the second regression analysis. Otherwise, 25 there may have a change in the lake bank slope and, therefore, the extrapolation of CryoSat-2 data with Landsat-derived water levels was not suited. In this case, we reselected the ROI to extract lake shoreline positions and redid altimetry data merging until the Landsat-derived water levels and merged altimetry water levels agreed well with one another in the second linear regression. Detailed analysis about the potential extrapolation issue can be found in the supplementary file.
In summary, the basic idea of removing systematic biases of different altimetry water levels is to calculate the means of two 30 altimetry water level time series during the overlap period. The difference between the means is removed from one altimetry water level time series to make both altimetry water level time series consistent and to form a longer time series. This process was subsequently applied to all water level time series with overlap periods to merge them into a single time series for each lake. However, the overlap period may not be long enough, such as Envisat and CryoSat-2 (e.g., there are limited data points (e.g., 1-2) during the overlap period), or does not exist at all, such as ICESat and CryoSat-2. On these cases, Landsat-derived water levels are used to extend or create an overlap period that links two altimetry water level time series.

Hypsometry
We derived the hypsometric curve for each lake by polynomial fitting of the lake area and level time series. The lake area 5 comprises two parts: the inner invariable part and the outer variable part. As the variable water area was of more concern in this study, ROIs for extracting changes in lake area only cover the lake shoreline and its neighbouring areas as shown in Figure 6. The inner part of the water body was calculated only once and considered as invariant, making the calculation more efficient on GEE. Meanwhile, more images are available as the area of ROI becomes smaller, because the possibility of clouds covering the ROI is reduced compared with an ROI covering the entire lake. Landsat 7 ETM+ images after 2003 were 10 not included in this part of calculation as gaps negatively affected the ROI for lake area extraction. Similar to section 3.2, we selected images with less than 5% cloud cover on an ROI to generate time series of changes in lake area, obtaining 20-30 data points on average for regression. R 2 values for each lake are listed in Table 3, indicating that most lake basins agree well with the parabolic hypsometric curve.

Field experiment
Most Tibetan lakes are located in remote and inaccessible regions, resulting in the scarcity of ground-based in situ measurements that are, however, vital for data quality assessment. We made some in situ measurements in two lakes to validate the data quality of Landsat-derived water levels. The data quality of satellite altimetry on lakes or rivers has been 5 widely investigated and thus it is beyond the scope of our study. Many studies used in situ water levels to calculate statistical metrics, e.g., root mean squired error (RMSE). However, results provided by different studies vary, which could be associated with the cross-section width of the study water body in the ground track panel . This means that these results may not be comparable due to their unique applications. In addition, it is not rigorous to use in situ data of only one lake to represent the overall performance in the uncertainty assessment for altimetry water levels. Instead, we used 10 the standard deviation of valid footprints acquired in the same cycle as an estimate of uncertainty in satellite altimetry water levels. In contrast, the applicable condition of Landsat-derived water levels is not so variable as that of satellite altimetry data. Derivation of Landsat-derived water levels requires relatively flat bank as well as some altimetric information, which were available in all lakes. Since these selected bank slopes were similarly small (~1/30), it was possible to use a few typical lakes to represent all lakes. Therefore, we carried out a field experiment in Yamzhog Yumco and Nam Co to validate the 15 derived Landsat-derived water levels.
There were two main goals in our experiments: (1) collecting daily in situ water level data in a TP lake to validate the corresponding Landsat-derived water levels statistically; (2) testing the performance of extracting lake shoreline positions from high-resolution optical images (GF-2) to provide a theoretical uncertainty analysis of Landsat-derived water levels. On Yamzhog Yumco, we installed a pressure type water level sensor (type H5110-DY, manufactured by Shenzhen Hongdian technologies Co., Ltd.), which measured water pressure and temperature at the installation depth that were converted into water depths with a relative accuracy of ~0.1%. The device was carried onto the lake and put ~20 m below the water surface and 0.5 m above the lake bottom, suggesting an absolute error of ~2 cm. As for GF-2 images, the spatial resolution of the 5 panchromatic band is 0.8 m, which is able to provide a very accurate reference of lake boundaries for assessing water classification results for Landsat images. We used three GF-2 images acquired at different seasons (two in Jul and Sep, 2015 and one in Feb 2016) and different places on the TP to better represent the local conditions when extracting Landsat-derived water levels or areas. Image co-registration was performed to make sure that there was no obvious spatial shift between the GF-2 images and corresponding Landsat images. The accuracy of the image co-registration was ~2 m. 10

Uncertainty analysis of Landsat-derived water levels
Based on the in situ water level measurements made by the pressure-type water level sensor, we evaluated the accuracy of 15 Landsat-derived water levels statistically. We first calculated anomalies of in situ water levels and Landsat-derived water levels, and then water levels from the two sources acquired on the same day were used for analysis. There were 16 Landsatderived water level records available for the comparison against the in situ measurements, indicating an RMSE of the water level anomaly of 0.11 m. The linear fit shows a slope close to 1 and an R 2 of 0.89, suggesting the consistency between the in situ water level measurements and the derived Landsat-derived water levels (Figure 8(b)). It should be noted that the Landsat-derived water levels used for validation here were translated from lake shoreline positions using parameters derived from fitting with CryoSat-2 data, i.e., there is no in situ information involved in generating the Landsat-derived water levels shown in Figure 8. Furthermore, we performed a theoretical uncertainty analysis of the Landsat-derived water levels by looking at the original optical data and the generation process with the help of high-resolution GF-2 images. First, we took GF-2 images (after coregistration with the Landsat image for the same period and the co-registration errors were ~2 m) as the ground truth to 10 determine the accurate position and shape of the lake shoreline. Second, we performed water classification from the Landsat 8 OLI image for the same period jointly using the water index method and Otsu algorithm to derive the binarized image.
Landsat image pixels where the lake shorelines from the GF-2 images cross were delineated and marked as shoreline pixels as shown in Figure 9(a). Then the water area in each shoreline pixel was calculated.
Given that these shoreline pixels were classified as either water or land, a relationship between the water area ratio of the 15 shoreline pixel and the probability of the pixel being classified as water can be derived. This relationship generally describes the function of the water classification method by telling how likely a pixel will be determined as water, given the water area ratio of the pixel. Based on the observations of totalling 4128 Landsat shoreline pixels, a power function was chosen to represent the water classifier as Eq. (4) shows: where x represents the water area ratio in the shoreline pixel; f(x) represents the probability of the shoreline pixel being 20 classified as the water pixel; and n is the parameter that determines the shape of the curve. Parameter n was determined using the maximum likelihood method (Figure 9 (e)). As expected, the probability of the pixel classified as water increases with the water area ratio in the pixel (Figure 9(e)). The enclosed area of the fitting curve (y = x 1.43 ) is smaller than that of y = x on [0, 1], suggesting that there may be a lower probability of the occurrence of water pixels that is associated with a systematic bias of the lake shoreline detection. Note that the systematic bias can be removed when linearly fitting the lake shoreline positions and altimetry water levels as long 10 as the bias is stable. Therefore, uncertainty in Landsat-derived water levels developed in this study arises mainly from the variation in this systematic bias.
To describe the variation in the systematic bias, a new random variable X was introduced to represent the bias between the classified water area and the real water area in a shoreline pixel. Given the shape and position of the lake shoreline, the real water area in each shoreline pixel is a complex function of the relative position between the pixel and the shoreline. To 15 simplify the derivation, we assumed that the water area ratio in a shoreline pixel is uniformly distributed on [0,1], meaning that the probability of any value between 0 and 1 is equal. If we use X0 to represent the true water area ratio in the shoreline pixel and X1 to represent the classified results based on the water area ratio, the random variable X can be expressed as: where X1 can take on 1 or 0 (i.e., the classified results only tell us whether a pixel is water pixel or not), so X can only take on either 1-X0 or -X0. Because the range of X0 is [0,1], it is obvious that the range of X is [-1,1]. A derivation of F(X), i.e., the probability density function (PDF) of X can be found in the supplementary file (part 2).
Overall, F(X) describes how the bias between the classified water ratio and real water ratio in shoreline pixels is distributed as shown in Figure 10. If there are N shoreline pixels in an ROI, we can take them as N independent observations of X and 5 calculate the mean value ̅ . This value ̅ can represent an average shift of the detected lake shoreline from the real lake shoreline in the unit of one-pixel width (30 m). As we mentioned above, the systematic bias can be removed in the regression between the lake shoreline positions and altimetry water levels. As such, it is the variation of the bias that determines the accuracy of the Landsat-derived water levels. We can calculate the standard variation of ̅ to represent the uncertainty in lake shoreline position. Note that there is a simple relationship between σ ̅ and σ : 10 One only needs to calculate σx: If the slope of the shoreline is known, e.g., tanθ, the uncertainty of the Landsat-derived water level can be expressed as: where σho is the uncertainty of Landsat-derived water levels and d is the spatial resolution of the satellite image (30 m). In 5 this study, a typical width of ROI for deriving Landsat-derived water levels is ~10-pixel width, meaning that N is ~10. In addition, lake shores used for generating Landsat-derived water levels here generally have a relatively mild slope of ~1/30 or even smaller, which can be rounghly estimated from the maximun shoreline change and altimetry water level change within a year. Here if we use 1/30 as the slope tanθ, the uncertainty of the Landsat-derived water levels can untimately be estimated to be ~0.12 m, which is very close to the RMSE of 0.11 m based on the comparision between the opitcal water levels and in 10 situ water level measurements mentioned earlier.
However, for most cases we do not know the exact lake bank slope tanθ, which is the reason why we performed the regression analysis between the lake shoreline positions and altimetry-derived water levels. Information on the real lake bank slope is implicitly expressed in the linear fitting slope β (if the fitting line is y = βx + α). Uncertainty in altimetric information could evolve into the fitting parameters and impact the accuracy of the generated Landsat-derived water levels. Given that 15 the observed lake shoreline position is X1 (e.g., X1 = 5.6 meaning that the observed lake shoreline position is 5.6 Landsat pixels away from the initial postion corresponding to the lowest water level, different from Eq. (5), X1 here can be a rational number because it is determined by averaging all shoreline pixels in the ROI, whereas in Eq. (5) we focused on only one shoreline pixel), combining Eq. (5), the Landsat-derived water level (y) can be expressed as: where ( 1 − 1 ̅̅̅ ) denotes the observed lake shoreline change (in the unit of a Landsat pixel); 1 ̅̅̅ denotes the mean of observed lake shoreline positions used for linear regression; ̅ denotes the mean of altimetry water levels used for linear regression; ( 0 − 0 ̅̅̅ ) denotes the real lake shoreline change; ( − ̅ ) denotes the variation of the Landsat-derived shoreline position caused by the water extraction method; and β is the linear fitting slope. It is obvious that the expected value ( − ̅ ) is 0. As we discussed earlier, a systematic bias does not affect the accuracy of the Landsat-derived water level, but the 5 variation of the systematic bias does.
Based on Eq. (11), the overall uncertainty of Landsat-derived water level σy can be given as: where β and can be derived from the linear regression analysis; ̅ is given in Eq. (9); and ̅ is the uncertainty of the mean altimetry water level which can be estimated from the altimetry data. For a typical lake like Yamzhog Yumco, β = 0.35 m, = 0.02 m, Max(| 1 − 1 | ̅̅̅̅ ) = 11, ̅ = 0.13, ̅ = 0.015 m, which gives a maximum of 0.22 m. Note that ( 0 − 0 ̅̅̅ ) 10 is assumed to be the ground truth so there is no error associated with this term. This relationship shows that the uncertainty in the Landsat-derived water level increases with the distance from the center point ( 1 ̅̅̅ , ̅ ) represented by ( 1 − 1 ̅̅̅ ) 2 .
Interpretation of this phenomenon is that extrapolation of Landsat-derived water levels (far from the center point) may cause some errors and should be used with caution. More detailed discussion on the extrapolation can be found in the supplementary file. 15 Overall, the uncertainty quantification of the Landsat-derived water levels developed in this study indicates clearly that the accuracy of Landsat-derived water levels depends on the width of an ROI, e.g., the number of pixels/observations, slope of the lake shore, the effectiveness of the water classification method, and the uncertainty in altimetry water level used for regression. One of the advantages of the Landsat-derived water level is that an ROI does not necessarily cover a large area of lake shores, which maximizes the potential of optical remote sensing images to increase the spatial coverage and temporal 20 resolution of lake water level estimates that may not be realized by using satellite altimetry alone. Optical remote sensing images provide important complementary information on altimetry water levels and can subsequently facilitate lake water storage estimation.

Cross validation with similar products
We compared our product with a widely used lake water level/storage data set provided by the LEOGS Hydroweb, 25 indicating that the two products are, on the whole, consistent with each other (shown in Figure 11) but our product may perform better in terms of the temporal continuity as well as the temporal resolution (shown in Section 6.2). Both advantages are important in improving our understanding of responses of lakes to climate change. There are 21 same lakes in both our study and LEOGS Hydroweb. Annual trends in water level and lake storage during 2003-2015 are compared in Figure 11, indicating the overall consistency of the two products in terms of R 2 of the linear fit.

Data set description and availability
The derived TP lake water levels, hypsometric curves, and water storage changes are archived and available at 10 https://doi.org/10.1594/PANGAEA.898411 (Li et al., 2019). The data sets cover large 52 large lakes (50 lakes with a surface area larger than 150 km 2 , and 2 lakes are 100-150 km 2 ) on the TP. The data sets consist of two parts: (1) a table containing hypsometric curves and corresponding regression statistics (R 2 and the number of data pairs) for each lake, with parameters of the hypsometric curves listed in separate columns for the convenience of batch processing; (2) time series for each lake archived as 52 entities with geographic information (i.e., latitude, longitude, and size of the lake) that can be checked in an 15 online map provided by PANGAEA, avoiding the confusion of lake names. The time series of each lake include lake water levels and lake storage changes.
For data points in the water level time series, satellite or sensor type is shown (i.e., from Jason-1/2/3, Envisat, ICESat, CryoSat, or optical images). Uncertainty was calculated using the standard deviation of valid footprints in the cycle (only for altimetry data). The lake water storage time series were transformed from water level time series using the hypsometric 20 relationship so that they have equal data size. The lake water storage time series represent changes in lake storage with respect to a reference water level, which is listed in the corresponding hypsometric curve table as a parameter. The overall uncertainty of Landsat-derived water levels within the regression range (the range of altimetry water levels) is 0.1-0.2 m based on the experiment and analysis in this paper. The extrapolation of Landsat-derived water levels may occur during the time gap between altimetric missions and before 2002. The average uncertainty of altimetry water levels is 0.11 m.

Spatiotemporal analysis of changes in lake water storage in the Tibetan Plateau
Based on the lake water storage changes we derived, spatial patterns of lake storage trends during 2000-2017 were shown in 5 Figure 12. In the endorheic basin of the TP, similar to some reported results (Yao et al., 2018b;Zhang et al., 2017a), most lakes have been expanding rapidly, e.g., Selin Co (89.00 E, 31.80 N) gained 19.7±2.0 km 3 of water during the study period, Lake Kusai (92.90 E,35.70 N) experienced an abrupt expansion due to flood and gained 2.2±0.2 km 3 of water in 2011, as reported in related work . In contrast, some lakes in the southern part of the TP experienced shrinkage, e.g.,  However, spatial proximity cannot fully explain the intricate trend distribution in the Selin Co basin, where large lakes such as Selin Co were expanding whereas smaller adjacent lakes showed an opposite decreasing trend, e.g., Urru Co (88.00 E,31.70 N),Lake Co Ngoin (88.77 E,31.60 N),and Goren Co (88.37 E,31.10 N). In fact, we found that the decreasing trends in some small lakes like Goren Co were not detected in Yao et al. (2018b), which is likely due to the lower temporal resolution as shown in Figure 13. The three shrinking lakes are located in the upstream region and feed Selin Co through two small rivers. One of the rivers links Lakes Goren Co, Urru Co, and Selin Co, whereas the other links Lakes Co Ngoin and Selin Co. 5 Figure 13. Discrepancy of lake storage trends in Goren Co between Yao et al. (2018a) and our study.
A possible explanation of the disparity of changes in lake water storage in the Selin Co basin could be the principle of minimum potential energy. If we simplify the basin with the tank model and take the upstream small lake as a tank with a leaking hole, the storage of the small lake is mainly controlled by the height of the leaking hole. Given that surface water of 10 the small lake increased, most of the increased water would flow into the large lake (a lower tank), and the outflow discharge of the small lake at higher elevations would increase accordingly. The height of the leaking hole would decline (erosion) so as to increase the overflow capacity, which eventually results in the decrease in small lake storage. Another possible situation is that the height of the leaking hole remains the same and the water surface height of the small lake increases, but this situation is not consistent with the minimum potential energy principle, as more water potential energy is stored in the 15 small lake. This phenomenon shows that river-lake interactions may cause complex patterns of the regional surface water distribution. Therefore, decreases in small lake water storage and increases in water storage of Selin Co in the basin detected by our study seem reasonable. Increases in small lake water storage in this basin reported in some published studies may be associated with the sparse sampling of lake water levels.
We averaged the total lake water storage change in each season to generate time series shown in Figure 14 (a). The overall 20 storage change in the 52 lakes is 98.3±2.1 km 3 . The total lake water storage was increasing rapidly during the first 12 years but became relatively stable since 2012. Intra-annual variation in the TP lakes can also be investigated using the densified lake water level time series generated by our study. We removed the linear trend (sometimes there were multiple linear trends for a lake in different periods, which were removed in a stepwise fashion) and calculated the mean monthly water level anomaly for each lake over the study period. Then the intra-annual water level change was represented by the difference between the maximum and minimum values of the monthly water level anomaly. The histogram of the intraannual water level change in Figure 14 (b) shows that most of the TP lakes have water level variations ranging from 0.3-0.75 m in a year on average. Similar work was performed by Lei et al. (2017) but only a small number of lakes were 5 investigated in their study.

Quality assessment of similar data products
Some obvious discrepancies between the two data sources can be noticed, e.g., water levels of Taro Co. Both Hydroweb data and our estimation used ICESat and CryoSat-2 data. The difference lies in the fact that our CryoSat-2 product was more updated with a longer time span but Hydroweb used an additional altimetry satellite SARAL. Because the systematic biases of both products were removed, it is possible that we chose different baselines that resulted in the overall shift as shown in 15 Figure 15 (a). For instance, we may use different sets of ellipsoid and geoid models. In addition to the overall shift, some time-dependent discrepancy can be found in Figure 15, e.g., periods highlighted by shading areas.
The black curve shows the Landsat-derived water level we derived, which is a critical reference for connecting two different altimetry data time series without an overlap period. The Landsat-derived water level shows that the last two samples of ICESat data should not be lower than the first few samples of the CryoSat-2/SARAL data (see the dashed boxes). However, 20 it is apparent that Hydroweb data display a reverse relationship, showing that the last two ICESat measurements are smaller than the first a few CryoSat/Saral measurements. It is likely due to an unremoved systematic bias between ICESat and CryoSat/Saral time series from Hydroweb data in Taro Co.
Even though the Landsat-derived water levels were generated by linearly fitting the lake shoreline positions with altimetry data, the relative magnitude of water levels during different periods should not be largely affected by the fitting parameters, e.g., if Landsat-derived water levels show that Ha >= Hb, where Ha (Hb) means water levels acquired in period A (B), the Ha >= Hb relationship would not change with the fitting parameters used to generate the Landsat-derived water levels. This is the main reason for us to use Landsat-derived water levels as reference. Therefore, Hydroweb data may overestimate the 5 increasing trends in the water levels of Taro Co as their ICESat data are ~0.3 m lower than the SARAL/CryoSat data. A similar issue can be observed in Zhari Namco and Ngoring lakes shown in Figure 15 (b)-(c), and the explanation is similar to that of Taro Co. This problem may also exist in some similar studies when multisource altimetry data without overlap periods were used. *LSH represents lake surface height Figure 15. Similarities and differences between water level time series from the LEGOS Hydroweb database (Cré taux et al., 2011a) and our study. (a) Taro Co (84.12 E, 31.14 N); (b) Zhari Namco (85.61 E, 30.93 N); and (c) Ngoring Lake (97.70 E, 34.90 N). Shading areas highlight the differences between the two data sets.
As shown in Figure 16, optical data can be less noisy than altimetry data in certain lakes and significantly improve the 5 continuity of lake level and storage change monitoring. In addition, a more apparent seasonality in lake level change can be seen from the generated lake level time series. These advantages would largely benefit a better understanding of responses of TP lakes to climate change and facilitate hydrologic modelling of lake basins, regional water balance analysis, and even hydrodynamic analysis of lake water bodies. 10 *LSH represents lake surface height Figure 16. Lake water level (left y axis) estimates from our approach for six TP lakes. Black lines represent optical data and red dots represent altimetry data.

Lake overflow flood monitoring
As mentioned earlier in Section 5.1, Lake Kusai experienced an abrupt expansion in 2011, resulting from dike-break of an 5 upstream lake (Hwang et al., 2019;Liu et al., 2016;Xiaojun et al., 2012), named Lake Zhuonai (91.93 E, 35.54 N). The outburst of Lake Zhuonai occurred on Sep 14 (Liu et al., 2016) and 2.47±0.06 km 3 of water leaking into the Kusai River (as shown in Figure 17(b)), the main inflow of Lake Kusai. The water level of Lake Kusai increased by up to 7.9±0.5 m within 20 days (from Sep 11 to Oct 1 in 2011) based on Jason-2 data, and then started to drop as water overflowed from the southeast corner into Lake Haidingnuoer (93.16 E, 35.55 N) and Lake Salt (93.40 E,35.52 N). Lake Salt, the lowest part of 10 the basin close to the basin boundary, has gained 3.0±0.1 km 3 of water since 2011 and became a critical threat to the surrounding residents and railway ~10 km southeast to the boundary. Note that there are few satellite altimetry data available for Lake Salt except several CryoSat-2 observations, where Landsat-derived water levels can provide a near real-time monitoring of changes in lake water level and storage that are crucial to flood early warning and risk management. Aided by the high temporal resolution lake water level series, it was possible to estimate the height of the outlet of Lake Kusai, an important parameter for overflow estimation. The overflow of Lake Kusai can help predict the water level rise in Lake Salt and even serve as an indicator of flood forecast, as Jason-3 data with a 10-day revisit cycle are now available on 20 Lake Kusai. Several pairs of Landsat 8 OLI images and lake water levels for the same period were compared to provide a range of possible outlet heights, which are likely to be 4483.9 m to 4484.1 m, as shown in Figure 18 (a). Then we measured the mean width of the outlet from high resolution optical images provided by Planet Explorer (Team, 2017), which is relatively stable in Dec 2011 at 31.5±2.3 m in recent years. Given lake water levels, the outlet height and width, estimation of overflow can be made using the broad crest weir formula: Eq 13 where C is a parameter mainly reflecting geometric characteristics of the weir that mainly varies from 0.3-0.4; b is the width of the weir; H is the water head with respect to the top of the weir; and g is the acceleration of gravity. 5 We determined C (~0.3) by using stage 1 shown in Figure 19 as a calibration period. Details can be found in the supplementary file. Then we applied this result to stage 2 shown in Figure 19 to estimate the total overflow from Lake Kusai and compared the overflow with total water gain in stage 2 in Lake Salt. Since Lake Salt mainly relied on the replenishment of Lake Kusai during that period, with little precipitation input and negligible glacier meltwater in winter, the outflow of Lake Kusai can be comparable with the water gain in Lake Salt derived from remote sensing, though there was a small 10 amount of evaporation loss. This relationship can provide a straightforward validation of our developed method. However, it was not available in stage 1, because the outflow of Lake Kusai first replenished Lake Haidingnuoer until the later began overflowing. Results based on Eq. (13) indicate that the total outflow from Lake Kusai in stage 2 ranged from 0.21 -0.22 km 3 , whereas the water gain in Lake Salt from remote sensing was 0.19±0.01 km 3 . This indicates that our high temporal resolution lake water level time series are valuable in monitoring and predicting lake outflow flooding that is crucial for the 15 safety of downstream residents and infrastructure. 5 Figure 19. Changes in the water level of Lake Kusai after receiving the outburst flood from Lake Zhuonai. Stage 1 was used to determine the range of parameter C in Eq. (13). Stage 2 was used to compared the simulated lake outflow from Kusai Lake based on Eq. (13) with the water gain estimate from remote sensing of Lake Salt downstream during the same period; and (b) changes in water storage of Lake Salt derived from remote sensing using our developed method. There was 0.19 km 3 of water gained in stage 2, which was comparable to the outflow estimate of Lake Kusai (0.22 km 3 ) based on Eq. (13).

7 Conclusion
In this study, we develop high temporal resolution (i.e., weekly to monthly timescales) water levels and storage change data sets for 52 large lakes on the TP during 2000-2017 by combining multiple altimetric missions and optical remote sensing images. Generated from lake shoreline positions and regression analysis with altimetry data, the Landsat-derived water level serves as a unique reference covering the entire study period, enabling a more consistent merging of multisource altimetry 10 time series. Multisource altimetry water levels are first extracted separately from spaceborne altimetry products and then combined into a longer and denser altimetry water level time series with systematic biases well removed using Landsatderived water levels as reference. The combined altimetry and Landsat-derived water levels increase the overall sampling frequency to sub-monthly regardless of the lake size.
By comparison with a widely used LEGOS Hydroweb data set, we show that without Landsat-derived water levels as a 15 reference, there may be a remaining bias in the combined altimetry water levels in certain lakes. Our study has considerably improved the temporal resolution of the monitoring of lake water level and storage changes in the TP. For most lakes examined in the published studies, to our best knowledge, the estimates from our study provide the observations of the highest temporal resolution that can better reveal the interannual and intra-annual variability and trends in lake water level and storage, even in some relatively small lakes whose annual trends may, however, be incorrectly estimated by sparse 20 sampling of lake water levels. The developed data sets can also facilitate the monitoring of some rapidly expanding lakes with overflow risks and provide important information on flood prediction and early warning.
We evaluate the uncertainty in the Landsat-derived water levels by field experiments and rigorous uncertainty analysis. Both methods are consistent that the magnitude of the uncertainty is ~0.1 m, which suggests that Landsat-derived water levels are often more efficient and less noisy than altimetry data when the altimeter footprints on the lake surface are insufficient, 25 especially for small lakes. Based on our estimates, 52 large TP lakes accounting for ~60% of the total TP lake area have gained 98.3±2.1 km 3 of water during the past 18 years. Lakes in the endorheic basin on the TP have mostly expanded. The complex spatial pattern of lake storage changes in the Selin Co basin was quantified and a possible explanation was proposed in this study. Note that the quality of the Landsat-derived water levels before 2002 may not be as good as those after 2002, because no altimetry data before 2002 are used in this study. Extrapolation of the relationship between lake 30 shoreline positions and water levels may not be stable if the water level during 2000-2001 was much lower or higher than those from 2002-2017. Discussions on how the extrapolation may affect the data quality can be found in the supplementary file.

Author contribution
LD and LX designed the research. LX, LD, HQ, and ZF developed the approaches and data sets. LX, HQ, HP, and LD carried out the field experiment. LX, LD, and YW contributed to the analysis of results and writing of the paper. 5

Competing interests
The authors declare that they have no conflict of interests.

Acknowledgements
This study was supported by the National Natural Science Foundation of China (Grant Nos. 91547210 and 51722903). Mr.
Mingda Du's assistance in the field experiments, discussion on the use of satellite altimetry data for lake monitoring with Dr. 10 Gang Qiao from Tongji University, and efforts in improving and archiving the data sets made by Dr. Daniela Ransby from PANGAEA data publisher are acknowledged here.