Lake storage variation on the endorheic Tibetan Plateau and its attribution to climate change since the new millennium

Alpine lakes in the interior of Tibet, the endorheic Changtang Plateau (CP), serve as ‘sentinels’ of regional climate change. Recent studies indicated that accelerated climate change has driven a widespread area expansion in lakes across the CP, but comprehensive and accurate quantifications of their storage changes are hitherto rare. This study integrated optical imagery and digital elevation models to uncover the fine spatial details of lake water storage (LWS) changes across the CP at an annual timescale after the new millennium (from 2002–2015). Validated by hypsometric information based on long-term altimetry measurements, our estimated LWS variations outperform some existing studies with reduced estimation biases and improved spatiotemporal coverages. The net LWS increased at an average rate of 7.34 ± 0.62 Gt yr−1 (cumulatively 95.42 ± 8.06 Gt), manifested as a dramatic monotonic increase of 9.05 ± 0.65 Gt yr−1 before 2012, a deceleration and pause in 2013–2014, and then an intriguing decline after 2014. Observations from the Gravity Recovery and Climate Experiment satellites reveal that the LWS pattern is in remarkable agreement with that of regional mass changes: a net effect of precipitation minus evapotranspiration (P-ET) in endorheic basins. Despite some regional variations, P-ET explains ~70% of the net LWS gain from 2002–2012 and the entire LWS loss after 2013. These findings clearly suggest that the water budget from net precipitation (i.e. P-ET) dominates those of glacier melt and permafrost degradation, and thus acts as the primary contributor to recent lake area/volume variations in endorheic Tibet. The produced lake areas and volume change dataset is freely available through PANAGEA (https://doi.pangaea.de/10.1594/PANGAEA.888706).


Introduction
Lakes are important hydrological components in alpine environments, where water budgets are highly vulnerable to climate change [1,2]. One of the world's largest groups of alpine lakes are located in the remote Changtang Plateau (CP) (figure 1). Here, surface water is landlocked due to arid climate and topographic barriers, forming a cluster of endorheic basins in the northwestern Tibetan Plateau (TP). Different from other arid endorheic basins, the CP possesses a high lake density of 4.6% [3], accounting for ∼60% of the total lake water storage (LWS) in the entire TP [4]. Under negligible human disturbance, these alpine lakes act as 'sentinels' of regional climate change [5]. Influenced by strengthened westerlies through advection of heat and moisture, the CP has undergone evident wetting [6] and warming [7,8] during the recent decades, which posed inevitable impacts on the water budget in its alpine lakes [9,10].
Using optical and altimetric satellites, previous studies have revealed that many large lakes on CP experienced dramatic increase in both area and volume during the past couple of decades, especially since the new millennium [11][12][13][14][15][16][17]. However, existing approaches for estimating volume variations across the CP are restricted to a limited number of lakes or multi-year to decadal temporal frequencies [13,15,16,[18][19]. For example, Zhang et al [16] applied ICESat altimetry and satellite imagery to quantify volume variations in 68 Tibetan lakes annually from 1989-2015. Yang et al [19] used satellite imagery and Digital Elevation Model (DEM) to estimate volume variations in 114 lakes at [5][6][7][8][9][10][11][12][13][14] year intervals during 1976-2013, but included only three discrete years (2000,2005 and 2013) after 2000. Lake dynamics varies across the CP and exhibits different changing rates during different periods. Particularly, in the past few years (since ∼2013), expansions of many Tibetan lakes have been decelerated or partially reversed [3,16,18,20], which suggests the necessity of a continuous, long-term monitoring of their volume variations.
Different from area change, volume variation has a unit dimension consistent with water fluxes such as precipitation (P) and evapotranspiration (ET). Thus, an accurate monitoring of CP lake volume is fundamental for a quantitative attribution of lake changing mechanism. Specifically, we aim at applying estimated volume changes to answer: was the observed lake dynamics across the CP (e.g. rapid expansion from ∼2000 and then deceleration since ∼2013) predominantly driven by the net change in precipitation and evapotranspiration (thereafter net precipitation) or by warming-induced glacier melt or permafrost degradation? Recently, Zhang et al [16] suggested that net precipitation is likely the key driver for lake volume variations on CP during 2003-2009 (when the ICE-Sat observations were available). Nevertheless, the 68 lakes they studied cover ∼53% of the total lake area across the CP [21], and the impact of net precipitation after 2009, including the period of decelerated lake expansion, remained unquantified.
Our recent effort [3] applied multi-source satellite imagery (Landsat and Huanjing) to produce a detailed mapping of all lakes greater than 1 km 2 (accounting for ∼97% of the total lake area [4]) across the CP from 2009-2014. Despite a limited period, our applied mapping methods and revealed lake area dynamics (see section 2.2.1) provide a critical basis that allows for an extended and more thorough monitoring of LWS variations across the CP. Thus in this study, we synergized our recent mapping with extended archival imagery and existing global DEMs (SRTM DEM [22] and ASTER GDEM [23]) to estimate water storage variations in the 871 lakes greater than 1 km 2 across the CP annually from 2002-2015. This estimate, to our best knowledge, achieved some of the highest spatial and temporal details in LWS monitoring on the climate-sensitive TP. Calculated LWS changes were next validated against the estimates derived from long-term radar altimetry, and then used along with the total water storage (TWS) variations observed by the Gravity Recovery and Climate Experiment (GRACE) satellites to test the hypothesis that 'net precipitation is the dominant driver for recent lake dynamics across the CP'. Our goal is to further eliminate the uncertainty in net LWS changes and improve our understanding of the recent and ongoing climate impacts on surface water budgets in the endorheic Tibet.

Monitoring LWS variations
We calculated the annual time series of water storage variation for each studied lake from 2002-2015 through a synergy of water area changes mapped from satellite imagery and bathymetric information retrieved from freely-available DEMs. Detailed methods are explained below.

Mapping lake area dynamics
Imagery from the Landsat archive (5, 7 and 8; 16 day repeat cycle) and China's Huanjing satellites (1A and 1 B; 2 day repeat cycle) were combined to improve the temporal coverage of cloud-free observations through data blending and fusion. All imagery have the same resolution (30 m) in visible and near-infrared bands (see [3] for details), and were collected from September to December during which lake extents on Tibetan Plateau are generally stable [2,3,13,[25][26]. Combined imagery from Landsat and Huanjing were used to extract annual extents in 871 lakes larger than 1 km 2 from 2009-2015. Among them, lake areas from 2009-2014 were directly obtained from our previous mapping [3]. For the period before 2009 (2002-2008) when Huanjing imagery were unavailable, images from Landsat 5 and 7 were fused using the planetary-scale platform Google Earth Engine (GEE) [27]. GEE provides parallel computation for large amounts of satellite data, making fusion of optical images efficient. For each lake, images acquired from September to December were filtered by the SimpleCloudScore algorithm provided in the GEE algorithms API, and then mosaicked using median-value composite. Despite an improved coverage through Landsat image fusion, cloud-free observations were not always feasible for each year and each lake. As a compromise, annual lake areas prior to 2009 were mapped only for 126 lakes larger than 50 km 2 .
A self-adaptive lake mapping scheme (similar to Li and Sheng [28] and Sheng et al [21]) was modified to extract lake extents from the optical images (see Yang et al [3] for details). In brief, the High Resolution Water Index (HRWI) [29] was applied to enhance the contrast between the water and non- Figure 1. Studied alpine lakes across the CP, a total area of ∼700 thousand km 2 in the northwestern TP, home to 396 endorheic basins [24]. Lake positions and sizes were derived from Yang et al [3].
water pixels in each image. Then, an initial threshold (T 0 ) of HRWI (T 0 = 0) was applied to flag possible water bodies. For each possible water body, the HRWI threshold was further fine-tuned by an iterative buffering method which updates the water extents and its local buffering region until the water extents converge. After that, any remnant errors in the automation results were removed manually with assistance of an interactive editing tool [30]. Our previous results have shown that lake extents extracted from both Landsat and Huanjing imagery using the proposed approach are highly consistent (slope = 1.00, R 2 > 0.99) (refer to Yang et al [3] for details). As our studied minimum lake area (1 km 2 ) far exceeds the spatial resolution of both imagery (30 m), the error in mapped lake areas was considered to be trivial [20,31,32] and thus not included in the uncertainty propagation for lake volume changes (see section 2.3).

Calculating lake volume changes
A DEM may provide useful bathymetric information for inferring lake volume changes if the lakebed topography between the minimum and maximum water levels during the study period was well exposed at the acquisition time of the DEM [19,33]. As most lakes on the CP have expanded since at least 2000 [34], the exposed topography surrounding each expanded lake on the SRTM DEM, which was acquired in February 2000, reveals the bathymetry that suffices the estimation of this lake's volume changes during our study period from 2002-2015. In other words, the hypsometry (i.e. level-area curve) derived from the SRTM DEM covers the range needed for volume change calculation without any downward extrapolation. For the remaining lakes where minimum water levels in the study period were higher than those in February 2000, the ASTER GDEM acquired during 2001-2008 was further considered if it could reduce the uncertainties of hypsometry extrapolation.
Specifically, we used the SRTM DEM version 3.0 [22] and the ASTER GDEM version 2 [23], both under the same resolution (1 arc second or ∼30 meters). The inundation area of a lake in the SRTM DEM was obtained from the SRTM Water Body Dataset (SWBD) [35] while that in the ASTER GDEM was detected based on surface flatness as in Fujisada et al [36]. Such a lake area (hereafter referred to as 'base area') corresponds to the water level (hereafter referred to as 'base level') at the acquisition time of the DEM (figure 2). Given higher vertical accuracy and quality consistency [37,38], the SRTM DEM was prioritized: it was chosen if (1) the base area is lower than the minimum lake area during the study period, or (2) the base area is lower than that in ASTER GDEM. Otherwise, the ASTER GDEM was applied. Then, base level was simulated to increase at a step of 1 m (the precision for these two DEMs) and the corresponding area below the simulated level in the chosen DEM was calculated until the area is greater than maximum lake area during the study period (figure 2). These levelarea pairs were used to construct the hypsometry. If the number of pairs is less than six, base level was further extended until the number reaches six for a better curve fitting purpose. Monotonic cubic spline fitting from the 'splines' package of the R software was implemented to fit the hypsometry considering the cubic spline generally outperforms other models (e.g. polynomial) in fitting statistics [4,39]. Finally, lake volume variations were calculated by the integral of the fitted hypsometry (f) as in equation 1: where Δ denotes the water storage change from time (t1) with lake area (A 1 ) and water level (L 1 ) to time (t2) with lake area (A 2 ) and water level (L 2 ). The proposed approach using DEM hypsometry was applied to estimate volume variations in all studied lakes except Zhuonai Lake (35.55 • N, 91.94 • E) because there was a drastic decline in its water level caused by a moraine dam failure in September of 2011 [40]. The DEM hypsometry for Zhuonai Lake was, thus, extended by an additional pair of water area and level after the dam failure observed by Landsat 5 imagery and Croysat-2 altimetry, respectively, in order to reduce the error from directly extrapolating the DEM hypsometry.
As previously described, annual time series of volume variation for all 871 lakes larger than 1 km 2 were estimated from the period 2009-2015, while those for the 126 lakes larger than 50 km 2 were further estimated for 2002-2008 during which smaller lakes were not estimated because of a poor temporal coverage caused by the absence of Huanjing imagery before 2009. These 126 larger lakes account for ∼83% of the total lake area and ∼96% of the total lake volume across the CP [3,4]. Given such dominance, their total volume changes may well represent those in all CP lakes. This is corroborated by a strong linear relationship between the two during the period of 2009-2015 ( figure 3(b)).
This relationship was used as a scaling function to infer annual volume variations in all CP lakes from 2002-2008 from those in the 126 large lakes (also see section 3.1). The best-fit linear regression was applied to calculate the trends in LWS time series (hereafter referred as 'DEM-derived trends'), and the uncertainty analysis is given in section 2.3.

Attributing LWS variations
We investigated the mechanism of monitored LWS variations across the CP using a simple mass balance model, where the change in TWS (ΔTWS), is the residual of precipitation (P), evapotranspiration (ET), and runoff (R).
On the CP where surface outflow is landlocked within each endorheic basin, the net runoff through each basin is considered to be zero (R ≈ 0), meaning that ΔTWS equals the difference between P and ET (i.e. net precipitation). If we partition ΔTWS on the CP into changes in LWS and other water components (integrating glaciers, permafrost, accumulated snow, soil moisture, and groundwater; hereafter non-lake water storage or NLWS), coupled with equation 2, LWS change (ΔLWS) equals net precipitation subtracted by NLWS change (ΔNLWS): Lacking in-situ meteorological measurements, direct acquisitions of P and ET are practically difficult in the remote CP [41]. For this reason, we substituted the changes in net precipitation by those of TWS observed by GRACE satellites [42][43][44]. Here we used monthly mass anomalies from the Jet Propulsion Laboratory mascon solution (JPL-RL05M version 2), which outperforms the conventional harmonic solutions with improved spatial localization and less signal leakage error [45][46][47]. The storage time series in each mascon can be considered to be fairly independent from those in the neighboring mascons [48]. Limited by the coarse mascon size (3 degrees, approximately 100 thousand km at the equator), TWS changes were retrieved only for the entire endorheic CP (∼700 thousand km 2 ) and four sub-regions (figure 1; sections 3 and 4) using a simple area scaling (see sections 2.3 and 4 for uncertainty analysis). The contribution of net precipitation to LWS changes was computed as the trend in GRACE-observed TWS time series (with removal of climatology) using best-fit linear regression, and the contribution of NLWS was calculated as the residual between net precipitation and LWS trends as in equation 3.

Accuracy assessments
DEM-derived LWS trends from 2002-2015 were validated against the trends during the same period but calculated using available hypsometric curves provided in the LEGOS Hydroweb database (http://hydroweb.theia-land.fr; hereafter Hydrowebderived trends) [49]. Many hypsometric curves provided by Hydroweb were constructed from multidecadal records of water levels observed by radar altimeters and water extents mapped from satellite imagery [18,49]. A total of 18 lakes on the CP with Hydroweb hypsometry were used for the validation. Volume of these lakes varies from ∼0.2 to ∼100 gigaton (see table S1); collectively, they account for approximately 70% of the total lake volume across the CP [4]. For a comparison purpose, we also validated the LWS trends from 2002-2015 using hypsometric curves constructed from ICESat altimetry and Landsat imagery by Zhang et al [16] (hereafter ICESat-derived trends). This validation is limited to 14 lakes where hypsometric curves are available in both Zhang et al [16] and Hydroweb.
Following standard error propagation rules, uncertainties of our estimated annual LWS changes were integrated from two major sources: (i) the bias of DEMderived trends, calculated as their discrepancy from Hydroweb-derived trends, and (ii) the uncertainty of estimated volume changes in small lakes during 2002-2008, calculated as the root mean square error (RMSE) of the rescaling relationship between volume variations in large lakes and all lakes ( figure 3(b)). Uncertainties of estimated annual P-ET were directly estimated from the storage errors provided in the GRACE mascon data, and the additional uncertainties induced by signal leakage in the fringe mascons are further explored in the Discussion. Propagated uncertainties were next applied to infer 95% confidence intervals of the calculated linear trends using a Monte Carlo approach similar to Wang et al [50].

Accuracy of estimated LWS trends
Our DEM-derived LWS trends from 2002-2015 are validated against the trends derived from Hydroweb hypsometry ( figure 3(a)). DEM-derived trends in the 18 assessed lakes appear to be highly consistent with the Hydroweb-derived results (slope: 1.008, R 2 : 0.968). Their aggregated storage change (the histogram in the upper left corner of figure 3(a)) shows a minor bias less than 5% of that derived from Hydroweb hypsometry. ICESat-derived trends, by contrast, are less consistent with the Hydroweb-derived trends (slope: 0.762, R 2 : 0.936). The aggregated storage change (the histogram in the lower-right corner of figure 3(a)) is ∼20% lower than that derived from Hydroweb hypsometry. This poorer consistency between ICESat-and Hydroweb-derived storage trends is likely caused by the extrapolation of hypsometric curves constructed from ICESat altimetry levels that are only available during 2003-2009 [18]. This contrast suggests that our proposed approach, using bathymetric information exposed in the DEMs, yields a comparable result with that produced by longer-term radar altimetry measurements, and is fairly reliable for estimating lake volume dynamics across the CP.
As shown in figure 3(b), calculated annual volume changes from 2009-2015 in all studied 871 lakes (>1 km 2 ) are closely correlated to those in the 126 large lakes (>50 km 2 ) with a slope of 1.153, R 2 of 0.993, and RMSE of 0.856 Gt. Therefore, net annual volume changes in the missing smaller lakes (<50 km 2 ) from 2002-2008 were inferred by this scaling relationship between volume changes in all lakes and large lakes, and the scaling errors (RMSE) were integrated for propagating LWS change uncertainties (refer to sections 2.3 and 4.3).

Spatial distribution of LWS changes
Alpine lakes across the CP experienced a widespread storage increase from 2002-2015 ( figure 4(a)). The magnitude of volume increase generally reduces along an east-to-west direction, which is consistent with the spatial gradient of decreasing lake abundance and size. As summarized in figure 4(b), LWS increase is proportional to lake size. For example, the largest lake Selingco (31.81 • N, 89.07 • E) in the southeastern CP experienced the fast volume increase of 1.22 ± 0.10 Gt yr −1 , accounting for 17% of the total LWS increase across the CP. About 78% of the total storage increase fed into lakes greater than 100 km 2 ( figure 4(b)), while only 9% occurred in lakes with size between 50-100 km 2 and 12% (estimated from the scaling relation as in figure 3(b)) occurred in smaller lakes (<50 km 2 ) ( figure 4(b)).  4(c)).

Trends in LWS across the CP and associations with net precipitation
The trajectory of net LWS across the CP exhibits three distinct phases (     precipitation in 2013. This time lag implies that NLWS (e.g. warming-induced glacier melting and permafrost thawing) may have compensated for the LWS reduction caused by immediate net precipitation decrease in the early period, but was unable to completely offset the long-term effect of net precipitation decline. As further summarized in table 1, net precipitation on the CP increased at an average rate of 4.11 ± 0.19 Gt yr −1 during our study period (2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015), which explains 56.0 ± 5.4% of the concurrent increase in LWS (7.34 ± 0.62 Gt yr −1 ). From 2002-2012, net precipitation increased at 6.31 ± 0.27 Gt yr −1 , accounting for 69.7 ± 5.8 % of the rapid LWS increase of 9.05 ± 0.65 Gt yr −1 , while net precipitation declined at a rate of −13.61 ± 2.03 Gt yr −1 from 2013 to 2015, which completely explains the concurrent LWS decline (−8.09 ± 3.37 Gt yr −1 ). Given these calculations, net precipitation appears to be the first-order contributor to the recent LWS changes across the CP.

Discussion
Our results indicate a dominant role of net precipitation on the recent decadal lake dynamics on the CP, which is consistent with several existing studies [2,13,16,17,[51][52]. However, this finding may contradict those of some other studies. For instance, Li et al [53] found limited impacts of glacier retreat and hypothesized that permafrost thawing was the primary cause for Tibetan lake expansion in recent decades, although Zhang et al [16] concluded that the magnitude of permafrost thawing is not substantially more than that of glacier retreat. Jiang et al [54] also suggested the important role of permafrost on rapid lake expansion in the northeastern CP. Thus, we here compare the relationships between LWS and net precipitation in four sub-regions of the CP ( figure  6). In northeastern, northwestern, and southwestern CP, rapid lake expansion accompanied with dramatic increase in P-ET during 2002-2015 (figures 6(a)-(c)), which confirms the dominant contribution of net precipitation. However, the covariation between LWS and net precipitation is not evident in the southeastern CP ( figure 6(d)). This may be caused by the substantial contribution from glacier retreat, as indicated in previous studies [55,56].
Several error sources limit our results and contribute to the uncertainty of this study. First, we ignored the additional uncertainty from hypsometry extrapolation for 6% of studied lakes (supplementary figure S1 available at stacks.iop.org/ERL/13/064011/ mmedia), where neither of the two DEM products could reveal the complete bathymetry required for storage change recovery during the study period. However, these lakes only account for less than 4% of the total lake area on CP [21], so their hypsometry extrapolation should not substantially affect our estimated net LWS variations. Second, except the errors provided by the RL05M mascon data, we did not   S2). It is worth noting that, despite a partial recovery of the signal variation within each mascon, the scale factors may not be suitable for deriving TWS trends at sub-mascon resolutions. This is because (1) the CLM does not include lake or glacier components so the simulated surface storage variations on the CP may be highly uncertain, and (2) the least-square correction involved in the factor calculation tends to be dominated by the annual cycles of water storage variations [45,57]. For these reasons, we did not apply the scale factors in the calculations as reported in section 3, but only used them for inferring possible uncertainty scales induced by signal leakage in the fringe mascons. Results of both solutions are consistent with our previous findings. The increase in net precipitation (TWS) accounts for most (∼70% or more) of the LWS increase from 2002-2012 and the net precipitation decrease fully explains the LWS loss from 2013-2015. Therefore, potential uncertainties due to signal leakage are not likely to alter our primary conclusion that LWS changes are predominantly attributed to the variations in net precipitation.

Summary and concluding remarks
This study provides a comprehensive estimate of LWS variations across the CP from 2002-2015, by synergizing satellite imagery (Landsat and Huanjing) and freely-available DEMs (SRTM DEM and ASTER GDEM). The sheer number of lakes (871) analyzed in our estimate, which account for 97% of total lake area across the CP [21], is substantially greater than any number recently studied for lake volume changes on the Tibetan Plateau (e.g. 68 lakes in Zhang et al [16] and 114 lakes in Yang et al [19]). Compared with existing estimates [15,16] using hypsometric curves extrapolated from short-term ICESat observations (available during 2003-2009), our LWS trends are more consistent with values derived from longer-term radar altimetry measurements (slope = 1.01, R 2 = 0.97). This implies that volume change estimates using extrapolated hypsometry may need to be interpreted with cautions. Given such improved spatiotemporal coverage and reduced estimation biases, this study advances our understanding of recent variations in lake water budget across the remote Tibetan Plateau. Our produced lake area and storage change data set is freely available at PANGAEA (https://doi.pangaea.de/ 10.1594/PANGAEA.888706). From 2002-2015, the net LWS increased at an average rate of 7.34 ± 0.62 Gt yr −1 (cumulatively 95.42 ± 8.06 Gt), manifested as a dramatic monotonic increase of 9.05 ± 0.65 Gt yr −1 before 2012, a deceleration and pause in 2013-2014, and then an intriguing decline after 2014. Using TWS anomalies from GRACE observations and a water balance model, we quantified that ∼70% of the monotonic LWS increase before 2012 was attributed to the increase in net precipitation (P-ET). Despite a smaller total influence, warming-induced NLWS changes, including glacier retreat and permafrost thawing, might compensate for the LWS reduction caused by the initial net precipitation decrease in 2013-2014, which is manifested by a 2 year time lag between net precipitation and LWS declines. However, the impact of NLWS was unable to offset the longer-term effect of net precipitation decline, leading to an evident LWS decrease in 2015: the first major reverse of a double-decadal lake expansion on the CP [16]. To this end, we conclude that net precipitation (i.e. P-ET) is a first-order contributor to alpine lake dynamics across the CP since at the least the new millennium.