Limited Contribution of Glacier Mass Loss to the Recent Increase in Tibetan Plateau Lake Volume

The Tibetan plateau plays an essential role in the water supply to Asia’s large river systems and, as the largest and highest mountain plateau in the world, it drives the Asian monsoon and influences global atmospheric circulation patterns. The increase in the Tibetan plateau lake volume since the mid-1990s is well documented, however the drivers of lake growth remain largely unexplained. In this study we investigate changes in lake and glacier volumes, together with changes in precipitation and evapotranspiration at basin scale. We calculate the contribution of glacier mass loss to the lake volume increase for the period 1994–2015. We demonstrate that glacier mass loss does have a limited contribution to the lake volume increase (19 ± 21% for the whole Tibetan plateau). Glacier mass loss is thus insufficient to explain all of the lake volume gain, and despite large spread in various products that estimate precipitation and evaporation, we suggest that an increase in precipitation excess (precipitation - evapotranspiration) may be sufficient to explain the lake volume gain.


INTRODUCTION
Lakes are widespread features of the Tibetan plateau, with more than 1,000 lakes larger than 1 km 2 (Zhang et al., 2017a). As these lakes are located mainly in endorheic basins, and are minimally influenced by human interventions, they are expected to be very sensitive to climate fluctuations. Most lakes have undergone a striking increase in area and volume since the mid-1990s (Lei et al., 2013;Song et al., 2013;Lei et al., 2014;Zhang et al., 2017a;Zhang et al., 2017b;Yao et al., 2018;Qiao et al., 2019b). Many recent studies assess changes in lake volume using remote sensing techniques, typically with a combination of optical satellite images, such as Landsat, and topographic data (e.g., Crétaux et al., 2011). The lake delineation is not always consistent among the different studies, and neither are the study periods, however all of them conclude that the lakes have gained mass at a pace of 6-9 Gt/yr from the mid-1990s to the mid-2010s. This increase in lake storage is confirmed by an increase in the total land water storage, as observed in Gravity Recovery and Climate Experiment (GRACE) data Xiang et al., 2016;Zhang et al., 2017b;Yao et al., 2018;Loomis et al., 2019). During the same time period, the Tibetan glaciers have been losing mass (Gardner et al., 2013;Shean et al., 2020).
In order to explain these lake changes, two main hypotheses are presented in the literature: glacier mass loss and precipitation increase. Many studies estimated the glacier mass loss to be the primary contributor of lake expansion in Tibet (e.g., Yao et al., 2007). At the scale of one lake or catchment, glaciers were found to contribute between 10 and 75% to a given lake mass increase based on direct observations or modeling (Zhu et al., 2010;Song and Sheng, 2016;Tong et al., 2016;Zhou et al., 2019). This contribution could be 22-51% for the northwestern Tibetan plateau (Huang et al., 2011), and 60% for the whole plateau (Qiao et al., 2019a;Qiao et al., 2019b). The lake changes are also attributed to an increase in precipitation after approximately 1995, which is observed both in meteorological data (station data and reanalysis products) and indirectly through the greening of the Tibetan plateau (Zhang et al., 2017c;Yao et al., 2018;Treichler et al., 2019). Various observations favor the precipitation hypothesis. For instance, glacier-fed and non-glacier-fed lakes were found to follow the same relative area increase (Song et al., 2014), and the spatial pattern of lake volume increase does not match the pattern of glacier mass loss as revealed by the Ice, Cloud, and land Elevation Satellite (ICESat; Li et al., 2014;Treichler et al., 2019). However a quantitative estimate of the glacier contribution to lake mass increase at the scale of the Tibetan plateau is currently missing in the literature (Zhang et al., 2020b).
Due to harsh environment and low population density, direct measurements of precipitation and evaporation remain sparse over the plateau, and consequently their representativeness is questionable. Reanalysis and satellite-based precipitation datasets show a large spread in the amount and evolution of precipitation and evaporation over the Tibetan plateau (e.g., Yoon et al., 2019). As a consequence, it is more straightforward to test the glacier hypothesis instead of the precipitation hypothesis. Most previous studies quantified the glacier contribution at catchment scale, except two recent studies, which quantified it for the whole plateau (Zhang et al., 2017b;Qiao et al., 2019a). Qiao et al. (2019a) estimated the contribution based on the difference between glacier-fed and non glacier-fed lakes, but without considering glacier mass changes. Zhang et al. (2017b) relied on glacier mass change data from ICESat between 2003 and 2009, but these glacier mass changes might not be representative of the past few decades.
As a large part of the Tibetan plateau is endorheic, lake volume changes are a direct result of changes in the storage of snow and ice in respective watersheds, or a result of shifts in precipitation and evapotranspiration in the lake basins. Although previous work points to both glacier mass loss and precipitation shifts as possible explanations of lake growth, the evidence remains largely inconclusive. In this study, we use recently published comprehensive lake volume change data (Treichler et al., 2019) and glacier mass change data (Shean et al., 2020) in order to quantitatively assess the glacier contribution to endorheic Tibetan plateau lake mass increase for the early twenty-first century. Even though some ground-based measurements of lake and glacier changes are available (e.g., Yao et al., 2012;Lei et al., 2019), we decided to rely solely on the remote sensing estimates as they are the only measurements covering the whole Tibetan plateau in a comprehensive and consistent way. We then use multiple climate reanalysis datasets to assess whether the remaining lake volume increase can be explained by recent changes in precipitation excess (precipitation minus evapotranspiration (P − ET)). Treichler et al. (2019) estimated Tibetan plateau lake volume change for over 1,300 lakes between 1990 and 2015. Here, we limit the analysis to a subset of lakes over Tibet interior endorheic basins (1,009 lakes covering 37,700 km 2 ). The lake volume change estimate is based on time series of lake area derived from the Global Surface Water dataset (Pekel et al., 2016), which are converted to volume changes with empirical area-volume scaling relationships derived from the digital elevation model (DEM) shoreline extraction method for each lake (Treichler et al., 2019). We compute lake volume differences from the 1994 and 2015 lake volume data to maximize spatial coverage and time interval (Supplementary Figure S1). Most previous studies detected the onset of the lake volume increase in the mid-1990s (Zhang et al., 2017a;Qiao et al., 2019a). We calculate the lake volume change as the difference between the estimated volumes in year i and j, divided by (j − i). The lake volume change is converted into lake mass change (ΔM l in Gt/yr), assuming a density of 1,000 kg/m 3 When no data are available for a given lake in 1994, we use the lake volume from 1993 or 1995, if available. Consequently, the overall coverage is 86% of the total lake area (Supplementary Figure S1 and Table S2).

Lake and Glacier Mass Changes
We validate the lake volume and level rate of change estimates against higher resolution lake level change data, which rely on laser and radar altimetry (Supplementary Text S1; Crétaux et al., 2011;Li et al., 2019). From this comparison, we determine an empirical relation between lake level change uncertainty and lake area of individual and recombined lakes (Supplementary Text S1 and Figures S2 and S3). As the lake level uncertainty decreases with increased lake area, we group the data into regions containing more than 1,000 km 2 of lakes. We define 18 regions, hereafter named basins, based on the aggregation of neighboring endorheic basins from the HydroBASINS database (Lehner and Grill, 2013), until they contain more than 1,000 km 2 of lakes ( Table 1). Figure 1 shows the location of the 18 basins. The total lake mass change is calculated as the sum of individual lake mass change. To account for missing data in Treichler et al. (2019), we assume that unsurveyed lakes behaved like surveyed lakes in the same basin. For each basin, the total rate of lake mass change was calculated by attributing the basin area-weighted rate of level change to the lakes that were not surveyed, which represent 0-52.1% of the lake area per basin (Supplementary Figure S4). The uncertainty of the basin lake level change is calculated as a function of the lake area (Supplementary Equation S1), and multiplied by the total basin lake area to estimate lake mass change uncertainty in Gt/yr.
We use glacier mass change data derived from high-resolution DEM time series for the period from 2000 to 2018 (Shean et al., 2020). The glacier mass change (ΔM g in Gt/yr) is calculated as the sum of individual glacier mass balance estimates within a given basin. We combined these data with annual glacier change time series at the scale of the plateau (Supplementary Text S2; Zemp et al. 2019). Uncertainties in glacier mass change are calculated at the basin level as the area-weighted average of the individual glacier specific mass balance uncertainties (in m w. e./yr). This Frontiers in Earth Science | www.frontiersin.org November 2020 | Volume 8 | Article 582060 2 uncertainty is multiplied by the glacierized area in order to obtain a basin-wide uncertainty in Gt/yr. From ΔM l and ΔM g , we define two additional quantities: the ratio between −ΔM g and ΔM l is the glacier contribution to lake mass increase (C g , in %), and the sum of ΔM g and ΔM l (in Gt/yr) divided by the basin area is defined as the water excess (in mm/yr). C g is not defined when there is no lake change, which is never the case in this study. As ΔM g is calculated as the sum of the individual glacier mass balances, this implies that individual glaciers with positive mass balance have a negative contribution to lake change. At the basin scale, if ΔM g is positive (i.e., glaciers gaining mass), then C g is negative, but artificially set to zero for FIGURE 1 | Lake and glacier mass change for the 18 endorheic basins considered in this study. The blue bars represent the total mass changes and the black lines represent uncertainties, with all bars plotted using the same scale. The basins are shaded according to the ratio between the glacier area and the basin area. The basemap is a terrain layer from Stamen (http://maps.stamen.com). The inset shows the lake mass gain as a function of the glacier mass loss for each basin (black dots), with basin number labels in blue. For the time periods of lake, glacier and P-ET changes, refer to the text. Water excess corresponds to the changes in the lake and glacier water storage. Δ (P-ET) is calculated for ERA5 only.
Frontiers in Earth Science | www.frontiersin.org November 2020 | Volume 8 | Article 582060 3 basin 03 in Table 1. If the absolute value of ΔM g is larger than ΔM l , the glacier contribution is larger than 100%, which can happen given the uncertainties on both terms. This definition of C g is valid only within a closed basin budget approach and neglects a large number of processes, such as loss due to evaporation between glacier and lake reservoirs. It is not valid for single lakes because connectivity between glaciers and lakes is not considered in this study (Phan et al., 2013). It is also noteworthy that we investigate only the glacier contribution to lake mass increase, which is different from the glacier contribution to the lake water balance that has been defined and investigated in other studies, for example by the mean of hydrological modeling (e.g., Biskop et al., 2016).
At the scale of the plateau, we evaluate both the lake and glacier mass change uncertainties as the sum of the uncertainties at the basin level: where σ ΔM is the Tibet-wide uncertainty (in Gt/yr) and σ ΔMi the basin-level uncertainty. This presumes fully correlated uncertainties, corresponding to the upper bound of the uncertainty estimate. The uncertainties obtained are σ ΔM l 1.55 Gt/yr (22% of the total lake change) and σ ΔM g 1.48 Gt/yr (106% of the total glacier change). The lower bound of uncertainties is obtained by summing in quadrature (assuming independence), which gives σ ΔM l 0.37 Gt/ yr (5% of the total lake change) and σ ΔM g 0.50 Gt/yr (36% of the total glacier change). As the sub-basin level aggregation is performed assuming full correlation, we use the more conservative upper bound estimates. Assuming independence of the uncertainties for lake and glacier mass changes, the uncertainty for the glacier contribution (σ Cg ) can be expressed as the quadratic sum of the relative uncertainties on the lake and glacier mass changes:

Meteorological Data Used for Estimation of Precipitation and Evapotranspiration
Due to the poor accessibility and low population density there are a limited number of meteorological stations with accessible data for the inner Tibetan plateau (Supplementary Figure S5D). Consequently, we rely on spatially and temporally continuous reanalysis and/or satellite-based products to investigate recent changes in precipitation (P) and evapotranspiration (ET). There is no consensus about the best reanalysis data available for the entire plateau (Wang and Zeng, 2012;He et al., 2019;Yoon et al., 2019). To account for the uncertainty among different products, we consider several reanalysis products (Supplementary Table  S3). We selected reanalysis products which include P and ET estimates for the period 1979-2015, in order to capture the climatology before the lake volume increase. Our final selection includes (Supplementary Table S3): ERA5 (Copernicus Climate Change Service, 2017), ERA-Interim (Dee et al., 2011), MERRA-2 (Gelaro et al., 2017), GLDAS-CLSM (Rodell et al., 2004), GLDAS-NOAH (Rodell et al., 2004). In addition, we use output from the Weather Research and Forecast (WRF) model, which are available for the inner Tibetan plateau (de Kok et al., 2020), named WRF20k hereafter. Additionally, we extracted P time series from the Frequent Rainfall Observations on GridS (FROGS) database (Roca et al., 2019), Climatologies at High resolution for the Earth's Land Surface Areas (CHELSA) (Karger et al., 2017) and 108 meteorological stations from the China Meteorological Association (stations with less than 0.1% of missing data for the period 1979-2015, Supplementary Figure S5D). Within FROGS, we selected precipitation products that cover the period 1979 onwards, regardless of their origin (Supplementary Table S3). The FROGS and CHELSA monthly time series were averaged to create the FROGS-ensemble. We also use ET estimates from the GLEAM v3.3a product (Martens et al., 2017). All sources used to estimate P and ET are described in Supplementary Table S3.
All the monthly aggregates of P and ET from gridded products are processed in the same way. Raster masks of the 18 basins at native resolution are created using Geospatial Data Abstraction software Library (GDAL). For each basin and product, we calculate a precipitation anomaly (Δ(P) in mm/yr) as the difference between the annual average P for the period 1994-2015 (P 1994−2015 ), where lakes growth is reported, minus the annual average P for the reference period 1979-1993 (P 1979−1993 ). We calculate Δ(ET) and Δ(P-ET) in the same way.
Even for the endorheic basins considered here, Δ(P-ET) is not directly comparable to the water excess, because of the additional evaporation due to lake area increase is not taken into account by reanalysis products (Song et al., 2019). Consequently, we do a simple scaling with the lake area change. We calculate the cumulative water equivalent change (ΔM rea in Gt), assuming that the mean (P-ET) for the 1979-1993 period (P − ET 1979−1993 ) corresponds to balance conditions, and is completely balanced by lake evaporation, as suggested by the stable lake conditions during this period (Qiao and Zhu, 2019;Zhang et al., 2020b). We then scale the annual (P-ET) to the lake area [A l (i)] relative to the lake surface of 1994 where P − ET i is the mean P-ET for the year i. ΔM rea is defined for t ≥ 1995, and is set to zero for t 1994. We define ΔM rea for the whole plateau only (A TP is the inner Tibetan plateau area, 9.598×10 5 km 2 ), as annual lake area values have high uncertainties and many gaps at the basin scale. We also define ΔM rea,ref (t), which is defined similarly to ΔM rea and assuming a constant lake area [A l (1994) ΔM rea is directly comparable to the water excess, if we assume three hypotheses: 1) there are no changes in groundwater and soil moisture storage between the two periods, 2) there is no leakage to the deeper groundwater or outside the sub-basin, and 3) the overall water budget was in equilibrium during the period 1979-1993. All three hypotheses will be discussed in Section 4. We neglect the role Frontiers in Earth Science | www.frontiersin.org November 2020 | Volume 8 | Article 582060 of snow in this analysis, as the snow water equivalent is relatively small in Tibet (Smith and Bookhagen, 2018).

Lake and Glacier Changes
For the whole plateau, lakes and glaciers have contrasting temporal evolution (Figure 2). While glacier mass loss is relatively steady from 1990 to 2016 (Figure 2A), the lake mass increases from 1995 on-wards. Due to limited data availability in 1995, we assume that the beginning of the ongoing lake mass increase is 1994. The rate of lake mass increase slows around 2010 ( Figure 2A). The total lake mass increase between 1994 and 2015 is 7.19 ± 1.55 Gt/yr ( Table 1). All the basins have growing lakes (Table 1 and Figure 1). The spatial pattern of lake rate of level change is shown in Figure 3A. Basins 11 and 16 (Tu Co and Selin Co),  Frontiers in Earth Science | www.frontiersin.org November 2020 | Volume 8 | Article 582060 5 located in the south eastern plateau, have the largest mean rate of lake level change (0.34 ± 0.04 and 0.41 ± 0.04 m/yr, respectively; Table 2 and Figure 3), corresponding to a rate of mass increase of 0.81 ± 0.09 and 1.09 ± 0.10 Gt/yr, respectively. If we normalize the lake mass change by the basin area, the pattern of water gain is altered ( Figure 4A). For instance, basins 03 and 05 (Margai Caka and Jingyu Lake) have large rates of lake level increase, but they correspond to little mass increase at the scale of the basin, due to the small lake area relative to the basin area ( Figure 4A). The largest lake water gain per basin unit area is the southeastern plateau, and more specifically basin 16 (Selin Co), with an increase of 37.5 ± 3.4 mm/yr ( Figure 4A).
The total water contribution from glacier mass change for the period 2000-2018 is −1.39 ± 1.48 Gt/yr for the 18 basins together. With the exception on the northwestern plateau (basins 01, 02, 03 and 04), where glaciers are near balance (specific mass balances ranging from −0.10 to 0.05 m w. e./yr; Table 2 and Figure 3B), the Tibetan plateau is dominated by glacier mass loss. The largest specific glacier mass losses are located in the southeastern plateau with Bam Co basin reaching −0.66 ± 0.28 m w. e./yr ( Figure 3B). However, when converted to total glacier mass loss, no clear pattern emerges, as the highly glacierized basins are located where glaciers are close to equilibrium (Figure 4).
Water excess from lake and glacier changes together is largest on the southeastern inner Tibetan plateau ( Figure 4C), with a maximum value of 34.8 ± 3.4 mm/yr in the Selim Co basin (16). The pattern of water excess is dominated by the lake mass changes, as the water contributions from glacier mass changes are an order of magnitude smaller ( Figures 4A-C).
The overall contribution of glaciers to lake mass increase for the entire plateau is estimated to be 19 ± 21% (Table 1). In most basins, the glacier contribution lies between 0 and 25% ( Table 1 and Figure 1). Five basins have glacier contributions larger than 25%, either because of their relatively large glacierized areas (basins 01 and 06) or small lake mass increase (below 0.1 Gt/yr, basins 07, 12, 13). These large glacier contributions are thus associated with very large uncertainties. Basin 03 has a negative glacier contribution, because glaciers are gaining mass in the western Kunlun Shan mountain range (Kääb et al., 2015;Brun et al., 2017;Lin et al., 2017;Shean et al., 2020;Zhang et al., 2020a). For the southeastern basins, where the water excess is the largest, the glacier contribution is below 10%, with the notable exception of Nam Co (basin 18), where the glacier contribution is about 25%.

Precipitation and Evapotranspiration Changes
The reanalysis products show large discrepancies, both in terms of climatology ( Figure 5) and in terms of changes in P and ET ( Figure 6). The plateau-wide annual precipitation ranges from 200 to 600 mm/yr depending on the reanalysis product considered ( Figure 5). Five of the seven products show larger P and ET in the south (eastern), while MERRA-2 and GLDAS-CLSM shows the largest P and ET in basins 06-08 located in the northeastern part of the region ( Figure 5). The plateau P − ET 1979−1993 ranges from 20 to more than 312 mm/yr ( Figure 6), with very different patterns of aridity (defined as ET/P) ranging from 0.5 to almost 1.0 in the GLDAS-CLSM reanalysis product ( Figure 5). For instance, for some reanalysis products, the northwestern plateau is very arid (ERA5, GLDAS), whereas for others it is one of the wettest basins (WRF20k, FROGS-GLEAM). Similarly, changes in P, ET and P-ET exhibit a wide range of values. Most reanalysis show an increase in P, except the two GLDAS reanalysis. The eight members of the FROGS ensemble show a precipitation increase at the scale of the Tibetan plateau (Figure 7), with the largest increase located in the southeastern part ( Figure 6). Eleven out of the twelve China Meteorological Association stations, out of the twelve located within the study region measured an increase in P ranging from 0.1 to 59.9 mm/yr, with a mean increase in P of 15.9 mm/yr ( Figures 4D and 8). Despite very sparse spatial coverage, they also show a large increase in precipitation in the southeastern Tibetan plateau ( Figure 4D). All reanalysis show an increase in ET ranging from 1.0 to 23.5 mm/yr (Figure 6), which is roughly one third of the increase in P for the four reanalysis products with increasing P. Δ(P-ET) ranges −4.6-36.1 mm/yr, with positive Δ(P-ET) in four of the six reanalysis products (ERA5, ERA-Interim, WRF20k and MERRA-2; Figures 6 and 9). All the reanalysis and the combination of the FROGS ensemble and GLEAM product show very different patterns of Δ(P-ET), even though a wetting of the southern fringe of the plateau is a redundant pattern in ERA5, MERRA-2, WRF, ERA Interim and FROGS-GLEAM ( Figure 6).
In terms of cumulative water equivalent change, the reanalysis products show large discrepancies ( Figure 2B). The four reanalysis products with positive Δ(P-ET) have also positive ΔM rea,ref , but the magnitude of the latter is much larger than the cumulative lake and glacier changes. Accounting for the additional lake evaporation in a simple way leads to a better agreement, except for MERRA-2, which is the only dataset that shows a very strong increase in P around 2010 ( Figure 2B). For ERA5, ERA-Interim and WRF20k, ΔM rea remains much larger than the cumulative water change from the lakes and glaciers. The two GLDAS reanalysis products, which have negative Δ(P-ET), also have negative ΔM rea,ref and ΔM rea . For ERA5, the lake evaporation due to lake area increase leads to a decrease in the cumulative water equivalent change of about 30%, which would correspond to ΔM rea (t 2015)/Δt 14 mm/yr (to allow a comparison with Δ(P-ET) 20 mm/yr in Table 1).

Insufficient Glacier Mass Loss to Explain the Lake Mass Increase
Despite relatively large uncertainties in the glacier contribution to lake mass increase, the estimate of 19 ± 21% is much lower than some recent estimates of glacier contribution (Qiao and Zhu, 2019). These values would correspond to glacier mass balances of −0.89 and −0.48 m w. e. yr, respectively. These mass balances estimates are unreconcilable with recent geodetic observations covering the same period that are about −0.16 ± 0.16 or −0.14 ± 0.07 m w. e. yr (Brun et al., 2017;Shean et al., 2020). The discrepancy between the results is probably explained by the assumption by Qiao and Zhu (2019) of similar lake response to climate in neighboring glacier-fed and non-glacier fed basins. This study's estimate is comparable to the 13% estimate obtained when using ICESat data (2003)(2004)(2005)(2006)(2007)(2008)(2009) to constrain glacier mass change (Zhang et al., 2017b), despite different periods. Even though glaciers play a limited role in the lake expansion at the scale of the Tibetan plateau, their local contribution can be very large (Zhou et al., 2019). For basins 07 (Har Lake) and 12 (Ngangla Ringco), the absolute glacier mass loss is larger than the absolute lake mass change, though uncertainties are large. Our methodological framework does not consider any loss during transport between glaciers and lakes, as these are very difficult to quantify. We did not examine the glacier/lake FIGURE 4 | Lake (A), glacier (B) and lake + glacier (C) rate of mass change distributed over the basin area (in mm/yr). Δ(P-ET) from ERA5 reanalysis, and Δ(P) from the meteorological stations (dots) for the period 1994-2015 relative to 1979-1993 (D).
Frontiers in Earth Science | www.frontiersin.org November 2020 | Volume 8 | Article 582060 connectivity and look only at basin-scale aggregated values. The inconsistency between the glacier and the lake periods (2000-2018 vs. 1994-2015, respectively) should have little influence on the uncertainty, as the glacier mass balance is estimated at −1.4 Gt/yr for the period 1994-2015 (Supplementary Text S2), implying that there are limited changes in the glacier mass change rates between the two periods. To account for these unconstrained sources of error we estimate the contribution of glaciers to lake mass increase with large error bars, and chose the conservative assumption of fully correlated uncertainties when estimating the plateau-wide lake and glacier mass change.
Despite these large uncertainties, we can safely conclude that glacier mass loss is insufficient to explain the endorheic Tibetan plateau lake mass gain for the period 1994-2015. If the glaciers were fully responsible of the lake mass increase, it would imply a glacier specific mass balance of −0.89 m w. e./yr, which is much FIGURE 5 | Mean P and ET for the period 1979-2015 for each reanalysis product (with the exception of WRF20k, which ends in 2010, and the GLDAS reanalysis which ends in 2014). ET/P 1979−2015 is calculated as the annual mean ET for the period 1979-2015 divided by the annual mean P for the same period. The FROGS-GLEAM data correspond to P from FROGS and ET from GLEAM. The numbers in each subplot correspond to the area-weighted averages.
Frontiers in Earth Science | www.frontiersin.org November 2020 | Volume 8 | Article 582060 8 more negative than recent mass balance observations (Brun et al., 2017;Shean et al., 2020). However, the glacier contribution might increase or even dominate in the coming decades, as inner Tibet glaciers are expected to recede at increasing rates in the near future (Kraaijenbrink et al., 2017).

Link Between Meteorology and Increased Water Storage
As mentioned in the method section, this work relies on three assumptions: 1) no changes in groundwater and soil moisture storage, 2) no leakage to the deeper groundwater and 3) equilibrium of the water budget during the period 1979-1993. The hypothesis of constant groundwater and soil moisture storage is supported by published trends in terrestrial water storage measured by GRACE. These include both surface and sub-surface water storage changes, thus a difference between numbers calculated in this study and terrestrial water storage trends would be indicative of changes in groundwater or soil moisture storage. Terrestrial water storage trends for the inner plateau for the period 2003-2016 range from 5.1 ± 1.0 to 7.3 ± 0.6 Gt/yr Yao et al., 2018), which matches Frontiers in Earth Science | www.frontiersin.org November 2020 | Volume 8 | Article 582060 9 ΔM l + ΔM g (as glacier mass is redistributed, not removed) of 5.80 ± 1.55 Gt/yr well, even though the large uncertainties on both terms could include some groundwater or soil storage changes.
The second hypothesis of no leakage to the deeper groundwater outside of the sub-basin is harder to justify. Very few estimates of the actual lake seepage are available, even though it could be a substantial contribution to lake water balance (Zhou et al., 2013). However, changes in seepage should be linked to changes in fault activities which, to our best knowledge, have not been reported for this study area. The third hypothesis of a balanced water budget during 1979-1993 is also challenging to validate due to the limited availability of remote sensing data in the past. Previous studies suggest that surface water storage was close to balance prior to 2000, with a lake mass change of −1.7 Gt/yr for  Frontiers in Earth Science | www.frontiersin.org November 2020 | Volume 8 | Article 582060 10 the period 1976-1990(Qiao et al., 2019b and moderate glacier mass loss for the period 1970-2000 (Zhou et al., 2018;Zemp et al., 2019). ΔM rea estimates are not very sensitive to the relaxation of this hypothesis (Supplementary Text S3 and Fig. S5).
However, we find that none of the reanalysis products used here can fully explain the lake mass change, and therefore close the water budget of these endorheic basins. Reanalysis products are notoriously challenging to validate in mountainous areas, with few in situ observations (Palazzi et al., 2013;He et al., 2019), and large biases (Orsolini et al., 2019). On top of the discrepancies between reanalysis mean values for P and ET (Figure 6), we can question the ability of the reanalysis products to close the water budget in regions where evaporation from lakes is a large component of the water budget and is not taken into account in the reanalysis. To illustrate this, we calculate the lake evaporation needed to meet the mean (P-ET) value at the basin scale. We find lake-equivalent evaporation ranging from 440 to 8,350 mm/yr, with none of the reanalysis products matching the range of observed values from 600 to 1,300 mm/yr (Supplementary Figure S6; Ma et al., 2016;Wang et al., 2020). This is a very serious limitation for the water budget closure, in particular for MERRA-2 and WRF20k products that show unrealistically high lake-equivalent evaporation (Supplementary Figure S6). Overall, the precipitation minus evapotranspiration balance is extremely sensitive in arid environments. Many additional processes could alter this balance. For example, the recent greening of the plateau could have dramatically increased the transpiration, something that is not captured by the reanalysis products (Zhu et al., 2016;Zhang et al., 2017c). Also an increase in the permafrost active layer depth and thawing permafrost could lead to enhanced evapotranspiration (Wu and Zhang, 2010).

Significance of the Lake Level Changes in Light of Climate Variability
Lakes act as low-pass filters of climate variability (e.g., Hasselmann, 1976;Mason et al., 1994;Huybers et al., 2016). Large changes in closed lake level and area can occur over decadal timescales under a constant climate in response to year-to-year stochastic climate fluctuations (Huybers et al., 2016). Our study focuses on the water budget at basin scale, and consequently does not investigate individual lake responses. However, we find that at the scale of the plateau, the increased lake mass corresponds to a sustained water excess of 6.0 ± 1.6 mm/yr, reduced to approximately 4 mm/yr if we account for the additional evaporation losses due to lake expansion (Song et al., 2019). In a qualitative approach, we can compare this value to the annual variability of the series of P and ET ( Table 3). The inter-annual variability in precipitation is about 41-46 mm/yr, and the interannual variability in evapo-transpiration is about 16-22 mm/yr, which are both higher than the water excess, even though the water excess is a 20-years average.
In a slightly more quantitative analysis, we calculate the response time (τ e ) of 34 of the largest lakes of the inner Tibetan plateau, which account for 38% of the total lake area (Supplementary Text S4 and Table S4): where dA l dl (in km 2 /m) is the lake area change in response to a level change, E l (in m/yr) is the lake evaporation and P l (in m/ yr). dA l dl is calculated by fitting linear fit of the relationship between A l and l (Treichler et al., 2019). We find response times FIGURE 9 | Annual precipitation (upper panel) and mean annual precipitation (lower panel) for the twelve China Meteorological Association stations located within the study region. Additionally, figure S7 shows the inter-annual variability of precipitation for each station.
Frontiers in Earth Science | www.frontiersin.org November 2020 | Volume 8 | Article 582060 ranging from 6 to 126 years (mean value 49 years). Given the calculated response times and observed climatic variability ( Table 3 and Supplementary Table S4), and considering the results of Huybers et al. (2016), we cannot rule out the possibility that the observed recent lake changes are a response to climate variability. Detailed modeling of water balance at lake scale and longer annual lake level records would be needed to further address this question.

Missing Processes
From this analysis, it is clear that glacier mass loss is not the main driver of the recent Tibetan lake volume increase. Most data sources (station data, the FROGS ensemble and four out of the six reanalysis products) point toward an increase in precipitation. However, the synoptic drivers for such changes are not clear, and could be linked to changes in the moisture sources (Gao et al., 2014) and/or to the influence of Westerlies (Mölg et al., 2017). There is a need to better identify these mechanisms. A better representation of evapotranspiration processes is also needed in reanalysis products. We still lack a comprehensive understanding of the Tibetan plateau water cycle in dry conditions, especially the role of snow and sublimation. In addition to total precipitation amount, the lake level changes might also respond to the intensity and frequency of the rainfall and snow fall and/or to changes in the precipitation seasonality. For instance, the increase in May precipitation observed in the south eastern edge of the plateau coincides with the location of the highest lake level increase (Zhang et al., 2017c). Consequently, the earlier onset on the South Asian monsoon might play an important role in lake volume increase .

CONCLUSION
By combining recent observations of lake and glacier mass change for the endorheic Tibetan plateau, we demonstrate that glacier mass loss has a relatively limited contribution of 19 ± 21% to the recent lake mass increase. The increase in surface water storage (lakes + glaciers) is largest on the southeastern edge of the plateau, which is also where the precipitation increase is largest. Overall, we do not manage to close the water budget, but we show that the increase in precipitation minus evapotranspiration after 1994 is three to five times larger than the increase in water stored in lake and glacier volume. We suggest that the remaining difference could be attributed to enhanced evaporation from the lakes due to warming and increased surface area, to enhanced evaporation from the ground due to greening of the land surface, and/or to changes in the seepage to the deep groundwater. Tibetan lakes are highly sensitive to changes in the water balance because they are sinks of large endorheic basins. Globalscale reanalysis and precipitation products show changes that are qualitatively, but not quantitatively, in agreement with the lake changes. Our results demonstrate the need for catchment-scale studies with local observations, to better understand the water budgets of the Tibetan plateau lakes.

DATA AVAILABILITY STATEMENT
The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

AUTHOR CONTRIBUTIONS
FB and WI designed the study. DT and DS provided data. FB analyzed the data. FB led the writing of the manuscript and all authors contributed to the writing.

ACKNOWLEDGMENTS
We thank Guoqing Zhang for valuable comments on the manuscript and four anonymous reviewers for critical comments that improved the original manuscript. We thank the Scientific Editor Christoph Schneider and the Chief Editor Regine Hock for handling the manuscript and for their careful reading.

SUPPLEMENTARY MATERIAL
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/feart.2020.582060/ full#supplementary-material P and ET are calculated as the mean value of the temporal average of the time series for the whole plateau for the whole study period  for each product. P′ and ET′ are calculated as the average of the standard deviation of each annual time series for the whole study period.