Interannual variation in lake areas over 50 km² on the Tibetan Plateau from 1986 to 2020 based on remote sensing big data

ABSTRACT

Lake distribution on the Tibetan Plateau (TP) is extensive, and lake area changes are key indicators of the TP's climate change response. Many multisource remote sensing big data for the TP, particularly optical images, are unusable due to cloud cover. Therefore, an improved isobath interpolation-based lake area extraction method is proposed and applied to obtain annual average lake areas (≥ 50 km²) on the TP from 1986 to 2020 using remote sensing big data. The lake area result accuracy was verified using existing lake area and level datasets, yielding correlation coefficients of ∼0.9. The change points and segmented trends of each lake's interannual area sequence were obtained. The relationships between lake area and climatic variables were investigated. The positive accumulation of the total precipitation minus total evaporation explains the overall lake area expansion trend after 1995. The exorheic lake interannual area is related to precipitation more than that of endorheic lakes, but endorheic lake area changes are stronger. The shrinking of lakes on the southern TP may not be climate-driven but probably attributed to lake bottom leakage. We explore detailed interannual variation characteristics of lake areas on the TP and provide reference data for studying lake responses to climate change.

KEYWORDS:

• Tibetan Plateau
• lake area
• variation trend
• climatic variables
• remote sensing big data

1 . Introduction

The Tibetan Plateau (TP), with an average altitude of more than 4000 m, is colder than other regions except the polar regions and is known as the 'third pole of the Earth' (Qiu 2008). The TP has the largest ice and snow reserves outside the Antarctic and Arctic regions and is the headwaters of more than 10 major Asian rivers, known as the 'Asian Water Tower' (Yao et al. 2022). In addition, the TP is very sensitive to global climate change. Its warming rate (∼ 0.36 °C per decade) during 1970–2018 was approximately twice that of the global mean (∼ 0.19 °C per decade), and the TP also exhibited an obvious increase in precipitation from 1980 to 2018 (GISTEMP 2019; Zhang et al. 2021b; Zhu et al. 2019).

Lakes are an important natural synthesis that connects the atmosphere, biosphere, pedosphere and terrestrial hydrosphere by participating in the global/local water cycle and energy balance. Since 2018, there have been approximately 1400 lakes with an area larger than 1 km² on the TP, with a total area of approximately 5 × 10⁴ km² (Zhang et al. 2019a). The main factors affecting the changes in lakes on the TP include precipitation, evaporation, glacier melting, and permafrost degradation, which are closely related to climate change. Thus, lake area (Lei et al. 2014; Zhang et al. 2019a), lake level (Crétaux et al. 2011), lake volume (Yang et al. 2017), lake surface temperature (Zhang et al. 2014a) and lake ice phenology (Cai et al. 2019) are often used as important climate change indicators on the TP. The lake area can be effectively obtained through satellite remote sensing images due to the obvious spectral difference between water bodies and land surfaces. Commonly used satellite remote sensing data, such as Landsat, MODIS (Moderate Resolution Imaging Spectroradiometer), Sentinel and GaoFen (Chinese high-resolution optical/radar satellites) series images, are applied for lake mapping on the TP (Lu et al. 2017; Wangchuk and Bolch 2020; Zhang et al. 2020a; Zhang et al. 2017b; Zhao et al. 2022). Among them, Landsat images have the longest time span, which can be traced back to 1984 (excluding Landsat-1∼4), with a moderate temporal resolution (16 days) and high spatial resolution (30 m). Generally, in the field of surface water dynamics remote sensing, there exists a duality between high spatial resolution and low temporal resolution studies versus low spatial resolution and high temporal resolution studies (Sogno, Klein, and Kuenzer 2022). The use of Landsat images in this study facilitates the study of the response of long-term interannual changes of lake area to climate changes on the TP.

Commonly used water extraction methods based on optical remote sensing images include the single-band threshold method (Jupp et al. 1985), spectral relationship method (Zhou, Du, and Luo 1996), water index method (McFeeters 1996; Xu 2006), decision tree method (Du et al. 2001), and object-oriented analysis method (Shen et al. 2012). Among these methods, the water index is a fast, simple and efficient water identification method. A lake is a relatively closed water body with a certain area. To obtain the vector boundary of the lake, the main methods include edge detection (Liu and Jezek 2004), regional segmentation (Setiawan, Rusydi, and Pradityo 2017), active contour modeling (Kass, Witkin, and Terzopoulos 1988), and manual digitization. Canny edge detection is a fast and accurate method to obtain lake boundaries (Liu and Jezek 2004). It is a good choice to combine the water index and the edge detection methods to obtain the water surface and lake boundary.

However, when using optical remote sensing images to monitor the lake area, it is not feasible to extract the complete water body and boundary of lake due to contamination by clouds and cloud shadows, especially on the TP. It is necessary to select images during the melting period from May to October on the TP to reduce the influence of snow and ice cover on water body extraction. Taking 2020 as an example, a total of 1834 Landsat-8 surface reflectance Tier1 images were collected over the TP from May to October, among which 580 images with 0–10% cloud cover were collected, accounting for only 31.6% of the total number of images. There were more clouds over the TP from May to July, and the images with less than 10% cloud cover accounted for less than 21% of the total images of that month. The high cloud cover over the TP region further results in a limited number of usable images, which may lead to long time intervals between lake area results. In some cases, the images corresponding to certain lakes are contaminated by clouds for an entire year.

For lake images that are only partially contaminated by clouds, some researchers proposed an automatic correction method based on the Global Surface Water Data Set (GSWD) (Pekel et al. 2016; Zhao and Gao 2018). This method uses the probability that each pixel is a water body to restore contaminated images, and it has been applied to the dynamic evaluation of reservoir surface area. Moreover, there is a high-frequency time series method for constructing global lake and reservoir areas using Landsat images (Yao et al. 2019), which uses isobaths to perform interpolation estimations of complete lake areas in contaminated images. However, neither of these two methods provided a solution to area extraction of natural lakes with complex boundaries, and visual inspection and post-processing are required to ensure the quality of the automatically generated lake area time series. In Yao et al.'s algorithm (2019), the boundaries of natural lakes are automatically extracted from good images as isobaths. Nevertheless, there may be some topological errors due to the influence of small ponds around them, such as non-closure and hanging arcs. Complicated cases, such as a lake that may shrink from one whole to several parts or expand from many parts to one whole, have not been discussed.

The goal of this study is to improve the isobath-based interpolation method for extracting complete lake areas from cloudy and other contaminated optical satellite remote sensing images to make it suitable for complex lake situations with high extraction accuracy without visual inspection and post-processing requirements, so as to obtain the interannual changes in lake area on the TP during 1986–2020. This study (1) improved the interpolation method by extracting isobaths from the lake boundaries in the existing lake area datasets and conducting preprocessing, (2) extracted the area of lakes ≥ 50 km² on the TP from Landsat-5/7/8 images during 1986–2020 using the improved algorithm, (3) evaluated the extracted lake area results with the existing lake area and water level datasets, (4) analyzed the interannual lake area changes through the detection of multiple change points and segmented trend analysis, and (5) investigated the relationship between lake area changes and main climatic variables (such as air temperature, precipitation, and evaporation) based on ERA5-Land reanalysis data.

2. Study area and data

2.1. Study area

In this study, a 2500 m contour line was used to delineate the TP (approximately 25°−40°N, 67°−105°E) with an area of approximately 3 × 10⁶ km² and an average elevation of approximately 4000 m (Zhang et al. 2013). The westerlies, Indian summer monsoon, and East Asian summer monsoon bring abundant water vapor to the TP, affecting the precipitation in this area (Yao et al. 2013). Precipitation on the TP gradually decreases from southeast to northwest, concentrated from June to September, accounting for 60–90% of the total annual precipitation (Kang et al. 2010). More precipitation may also lead to more cloud cover during the same period. According to the influences of different large-scale atmospheric circulation systems, the locations of lakes can be divided into four regions: the northwestern TP (mainly affected by westerly circulation), southern TP (mainly affected by Indian monsoon circulation), central TP (roughly between 30°N and 35°N, with both westerly and Indian monsoon circulation as the main influences, also known as the transition region), and northeastern TP (mainly affected by East Asian monsoon circulation) (Figure 1) (Maussion et al. 2014; Zhang, Zhang, and Zhu 2020b).

Figure 1. The boundary of the Tibetan Plateau and the distribution of 180 lakes. (i) Northwestern TP (Westerlies domain), (ii) Central TP (Transition domain), (iii) Southern TP (Indian monsoon domain), (iv) Northeastern TP (East Asian monsoon domain).

From 1986 to 2020, there were 187 lakes on the TP with an area larger than 50 km² for at least one year (Zhang et al. 2019a). Here, 180 of these lakes were selected as study objects, excluding 7 salt lakes with significant changes in area and boundary due to human development (Figure 1). A total of 145 lakes are distributed in the Inner Plateau basin, accounting for 80.6% of all 180 lakes. Second, there are 8 lakes in the Qaidam River Basin, accounting for 4.4% of all 180 lakes, while other basins contain fewer lakes.

2.2. Data

The datasets used in this study mainly include Landsat-5/7/8 images, digital elevation model, existing datasets of lake area on the TP, a lake water level dataset, climate reanalysis data, and El Niño index data (Table 1).

Table 1. Multiple datasets used in this study.

Data Type	Data	Period	Temporal Resolution	Horizontal Resolution	Source
Optical remote sensing image	Landsat 5 Surface Reflectance Tier 1	1986–2011			United States Geological Survey
	Landsat 7 Surface Reflectance Tier 1	1999–2013	16 days	30 m
	Landsat 8 Surface Reflectance Tier 1	2013–2020
DEM	SRTM 30 m DEM	2000		30 m	United States Geological Survey
Existing datasets of lake area	China lake dataset (1960s–2020)	1960s–2020	20/10/5 years		National Tibetan Plateau Data Center
	Long-term sequence dataset of lake area on the Tibetan Plateau (1970–2013)	1970–2013	10/3 years
	The lakes larger than 1km² in Tibetan Plateau (v2.0) (1970s–2018)	1970s–2018	20/5/1 years
	The lakes larger than 1km² in Tibetan Plateau (v3.1) (1970s–2022)	1970s–2022	20/5/1 years
Lake water level dataset	High-temporal-resolution water level and storage change data sets for lakes on the Tibetan Plateau during 2000–2017	2000–2017	1 month/ 10 days		National Tibetan Plateau Data Center
Climate reanalysis data	The enhanced global dataset ERA5-Land Monthly Averaged	1986–2020	1 month	9 km	European Centre for Medium-Range Weather Forecasts
El Niño index data	Niño 3.4 index data	1986–2020	1 month		Climate Prediction Center of the National Oceanic and Atmospheric Administration

2.2.1. Landsat image and digital elevation model (DEM)

In this study, the Landsat-5, 7 and 8 Surface Reflectance Tier 1 datasets on the Google Earth Engine (GEE) platform were used as remote sensing data sources to extract the lake area, and a total of approximately 92,000 images were used. All the datasets were subjected to atmospheric correction and orthography correction (USGS 2019a; 2019b), and the C Function of Mask (CFMASK) algorithm (Foga et al. 2017) was used to generate masks of clouds, shadows and snow. The 30 m DEM in 2000 from the Shuttle Radar Topography Mission (SRTM) of the USGS (Farr et al. 2007) was used to exclude the impacts due to mountains and mountain shadows.

2.2.2. Existing datasets of lake area

To provide high-precision isobaths for the lake area extraction method, this study used the lake datasets published by the National Tibetan Plateau Data Center. For the China lake dataset (1960s–2020), both Landsat images and topographic maps were used to obtain lake areas (>1 km²) in the 1960s, 1970s, 1990, 1995, 2000, 2005, 2010, 2015, and 2020 over China by combining semiautomatic water extraction with manual visual interpretation (Zhang et al. 2019b). The long-term sequence dataset of lake area on the Tibetan Plateau (1970–2013) includes the areas of 364 lakes (>10 km²) from the 1970s to 2013 (Zhang et al. 2017a). This dataset mainly uses Landsat images from October, which have less cloud cover and more stable lake area, and they obtain lake area every three years (1976–1980, 1986–1990, 1991–1993, 1994–1997, 1998–2000, 2001–2003, 2004–2007, 2008–2010, 2011–2013) to maximize the available data while reducing seasonal variation. The lakes larger than 1 km² in Tibetan Plateau (v2.0) (1970s–2018) dataset covers the periods of the 1970s, 1990, 1995, 2000, 2005, 2010, 2013, 2014, 2015, 2016, 2017 and 2018 (Zhang et al. 2019a; 2014b). Because the abovementioned three existing lake area datasets have been visually inspected, manually edited, and fully verified, they are highly accurate. The three datasets for the TP provide the same lake area in the same year, but they complement each other in different years. Their union will be used for the extraction of isobaths. Recently, the v2.0 dataset was updated to v3.1, with new data for 2019, 2020, 2021 and 2022. The new version of the data is not used to extract the isobaths but is used to verify the accuracy of the lake area extraction results.

2.2.3. Lake water level dataset

Since the lake water level is highly correlated with the lake area, the high-temporal-resolution water level and storage change dataset for lakes on the Tibetan Plateau during 2000–2017 was used to verify the extraction accuracy of lake area (Li et al. 2019a). The dataset is based on multisource altimetry satellite data and Landsat images, and the systematic deviation between multisource altimetry water levels was removed using optical water levels with an accuracy of 0.1–0.2 m. The dataset provides water levels of 52 large and medium-sized lakes on the TP from 2000 to 2017.

2.2.4. Other data

The enhanced global dataset ERA5-Land (Muñoz Sabater 2019) was generated by replaying the land part of the fifth generation of European reanalysis (ERA5). ERA5-Land's data product 'ERA5-Land Monthly Averaged' provides monthly averages of 50 climate reanalysis parameters, including air temperature, total precipitation and total evaporation. These data can be used to calculate the annual average air temperature, annual total precipitation, annual total evaporation, and annual total precipitation minus total evaporation of each lake ROI (Region of Interest). Refer to Table S1 for the main abbreviations used in this article.

The precipitation derived from the ERA5-Land reanalysis data may be overestimated relative to rain gauge data and satellite-based precipitation products (Xu et al. 2022; Liu et al. 2023a). However, the evaluation methods greatly depend on meteorological station data, but the distribution of meteorological observation stations on the TP is sparse (almost none in the central and western TP regions). The matching of stations and grids may result in loss of geographical detail, which may affect the evaluation results (Liu et al. 2023a). In addition, Wu et al. (2023) found the ERA5-Land reanalysis data can better reproduce the intra-annual variation and spatial–temporal pattern of precipitation on the TP. The multiyear (2003–2016) average evaporation amounts corresponding to the ROI of lakes from the ERA5-Land reanalysis data were compared with the published multiyear (2003–2016) average evaporation amounts of 75 lakes on the TP ($Eo$ and $Ea$ are the evaporation amount during the ice-free season and the annual evaporation amount, respectively) (Wang et al. 2020). Both datasets reveal proportional relationships ($r_{E,Eo}$≈ 0.51, $r_{E,Ea}$≈ 0.47) (Figure S1). The low evaporation within the lake ROI attributes to the influence of evaporation from the surrounding land. In general, ERA5-Land reanalysis data can provide long-term and spatially complete estimates of meteorological variables on the TP. Although climate variables from the ERA5-Land reanalysis data may exhibit overestimation or underestimation errors, their interannual variation characteristics are relatively reliable. These climate variables can be used to study the relationship between them and interannual lake area changes.

El Niño events are characterized by abnormal warming in the eastern equatorial Pacific Ocean and have global climate teleconnections. It is one of the most significant features of interannual climate variability (Yeh et al. 2009). Niño 3.4 index data (Trenberth 1997) were mainly used to study the relationship between El Niño or La Niña events and changes in lake area. El Niño or La Niña events are defined when the Niño 3.4 index exceeds +/−0.4 °C for a period of six months or more.

3. Methods

The flow chart of this study is shown in Figure 2, which is divided into two main parts: lake area extraction and accuracy verification, and interannual lake area trend analysis. The lake area extraction mainly includes (1) determining the three isobaths corresponding to the 180 lakes; (2) removing the contaminated parts and low-quality images; (3) calculating the modified normalized difference water index (MNDWI); (4) using the Otsu method to calculate the optimal water body segmentation thresholds; and (5) lake area interpolation. The verification of lake extraction accuracy mainly used the existing lake area datasets and the water level dataset to conduct correlation analysis with the results of this study. The Breaks for Additive Season and Trend (BFAST) change points detection, Mann-Kendall (M-K) trend test and Sen's slope methods were used to analyze the temporal variation trend in interannual lake area.

Figure 2. The overall technical flow chart of lake area extraction, accuracy verification, and interannual lake area trend analysis.

3.1. Improvement in the lake area extraction method

3.1.1. Determination of lake isobaths and ROI

The boundary of the lake can be considered the isobath of the lake, as the elevation of the whole lake surface is nearly the same. The lake boundaries corresponding to the maximum area, medium area and minimum area of the 180 lakes (≥ 50 km²) during 1986–2020 were selected from the abovementioned three existing lake area datasets, and the selected lake boundaries were recorded as the highest isobath, medium isobath and lowest isobath, respectively, as references for interpolation. Due to the high accuracy of the three existing lake area datasets mentioned above, the isobaths extracted in this study also have high precision.

Prior to using the derived isobaths to assist lake area extraction, there are several special situations that necessitate further processing and editing. (1) When three isobath lines partially intersect at the edge of some lakes with large slopes, to ensure the inclusion relationship between isobath lines of different heights, they need to be slightly edited. (2) For cases where multiple lakes expand and eventually connect into one lake or one lake shrinks and eventually separates into several smaller lakes, the isobaths of multiple smaller lakes need to be edited so that they are connected at the shortest distance and added very little area (Figure S2a and b). (3) For the presence of a large island in the center of the lake, the central island is excluded by establishing two lines at the shortest distance from the central island to the lake boundary (Figure S2c and d). These preprocessing methods for the derived isobaths can greatly improve the consistency of the research object and universality of the algorithm and ensure the suitability of the extracted isobaths for the subsequent fully automated process of lake area extraction, without visual correction and post-processing of the extraction results. Moreover, high-precision isobaths can result in the improvement in the final lake area extraction accuracy. The determination of lake ROIs is conducive to narrowing the research scope of lakes, reducing the amount of calculation, and facilitating the subsequent application of image segmentation algorithms (Figure 3). For each lake, the highest isobath region is selected to expand the buffer zone so that the area of the ROI is approximately twice the lake area.

Figure 3. Schematic diagram of lake ROI and three reference isobaths, taking Bamco Lake (31.26°N, 90.58°E, water level ≈ 4567.58 m) as an example.

3.1.2. Remove the contaminated parts and low-quality images

Based on the pixel quality assessment (pixel_qa) attribute contained in the Landsat-5/7/8 images and SRTM DEM 30 m data, all Landsat-5/7/8 images from May to October including each lake during 1986–2020, were masked by removing clouds, cloud shadows, snow and mountain areas. In addition, to avoid the reduction in the extraction accuracy of lake areas caused by excessive cloud cover and poor image quality, only the images of unmasked areas within the ROI accounting for at least 50% of the total ROI were retained for subsequent processing.

3.1.3. Calculation of the MNDWI

The MNDWI has a better water recognition effect than the normalized difference water index (NDWI) (Xu 2006), especially when water features are mixed with buildings or vegetation (Singh et al. 2015). The MNDWI can also reduce the influence of mountain shadows on water body extraction. By using the MNDWI, the distinction between water bodies and non-water bodies can be largely improved. (1) $\begin{array}{c}MNDWI={\displaystyle \frac{R_{Green}-R_{MIR}}{R_{Green}+R_{MIR}}}\end{array}$(1) where $R_{Green}$ is the surface reflectance of the green band, corresponding to band 2 of Landsat-5/7 and band 3 of Landsat-8; $R_{MIR}$ is the surface reflectance of the mid-infrared band, corresponding to band 5 of Landsat-5/7 and band 6 of Landsat-8.

3.1.4. Calculation of the optimal water segmentation threshold of the Otsu method

The Otsu algorithm, known as the maximum interclass variance method, is efficient in dividing gray images into binary images, as proposed by Otsu (1979). The Otsu method was used to segment the water and non-water areas in the MNDWI images by the optimized threshold determined by Fisher's linear discriminant analysis (Fisher 1936). In the MNDWI, pixels less than or equal to this threshold are identified as water bodies, while other pixels are identified as non-water bodies.

Then, Canny edge detection was used to obtain the lake boundaries from the water body distinction results, which were recorded as the known parts of the lake boundary from the part of the image that was not covered by clouds or other contaminated parts (Figure 3).

3.1.5. Lake area interpolation calculation

According to the known parts of the lake boundary extracted above, the two nearest isobaths are automatically selected from three reference isobaths, denoted as the lower isobath $l$ and the higher isobath $h$, respectively. Yao et al. (2019) assumed that between the two nearest isobaths, the change of water areas in the effective observation portion (unmasked) of a contaminated image is proportional to that of their full areas (Figure 4). Based on this assumption, $A_h$ and $A_l$ represent the full area corresponding to higher isobath h and lower isobath $l$, respectively, and $a_h$ and $a_l$ represent the partial area corresponding to higher isobath $h$ and lower isobath $l$ observed in the effective observation portion of the bad image, respectively. The variable $a$ is the partial lake area in the effective observation portion within the higher isobath range in the image, that is, the known part of the lake area identified in Section 3.1.4. $\hat{A}$ is the complete lake area corresponding to $a$, which can be estimated by the following formula. (2) $\begin{array}{c}\hat{A}\cong {\displaystyle \frac{A_h-A_l}{a_h-a_l}(a-a_l)+A_l}\end{array}$(2) By applying the abovementioned method to each Landsat image corresponding to each lake from 1986 to 2020, the lake area can be calculated from the available images. To analyze the interannual change in lake area results and reduce the seasonal difference within the year, it is necessary to calculate the annual average lake area.

Figure 4. Estimation of the full lake area from partial lake extents.

3.2. Interpolation of missing values in interannual lake area series

In the obtained interannual lake area series from 1986 to 2020, there are a small number of missing area values for some lakes in some years. To obtain the complete lake area series, the interannual lake area series are interpolated according to adjacent values. When the missing value is in the middle of the interannual area series, linear interpolation is performed by averaging forward and backward. When the missing value is at the beginning or end of the interannual area series, it is assigned the nearest area value.

3.3. Change points detection method

The BFAST method (Verbesselt et al. 2010) is commonly used for the decomposition and change identification of time series data, such as the normalized difference vegetation index (NDVI) (Geng et al. 2019) or air temperature (Li et al. 2019b). The BFAST method integrates the iterative decomposition of time series into trend, seasonal and noise components with methods for detecting change within time series, without the need to select a reference period, set a threshold, or define the trajectory of change (Verbesselt et al. 2010).

Multiple change points can be detected by applying the BFAST method without considering the seasonal term for the area series data of each lake from 1986 to 2020. The existence of abrupt changes is assessed by the ordinary least squares (OLS) residuals-based moving sum (MOSUM) test (Zeileis 2005). The trend components are estimated using a robust regression based on the M-estimation. The gradual change and mutation of variables are characterized by deducing the change time, amplitude and direction in the trend components of the time series. The positions of the detected change points usually correspond to the time points when the trend term of the lake area changes in direction and magnitude. Two parameters need to be set here. The parameter h represents the minimum segment size between potentially detected breaks in the trend model, given as a fraction relative to the sample size, and breaks represents the integer of the maximum number of breaks to be calculated. Many experiments that involve setting the values of h and breaks parameters were conducted. The parameters h = 0.2 and breaks = 3 were appropriate for 35-year sequences, and the detected change points were representative of shifts in long-term series data.

3.4. Trend analysis methods

The M-K trend test (Mann 1945), a non-parametric method that can be applied to time-series data with various distributions and tolerates outliers and missing values, was used to determine the trends and their significance for lake area and meteorological data. Sen's slope method (Sen 1968) was used to analyze the change rates of lake area and meteorological data. When the Sen's slope statistic β < 0, it shows a decreasing trend. When β = 0, there is no trend. When β > 0, there is an increasing trend. When the absolute value of the M-K test statistic Z is greater than or equal to 1.96 and 2.58, it means that it has passed the significance test at the 95% and 99% confidence levels, respectively. The significance values of the change trends are provided in Table 2.

Table 2. Significance of the change trends.

β	|Z|	Trends	Abbreviation
β < 0	|Z|≥2.58	Highly significant decrease	SD
β < 0	1.96≤|Z|<2.58	Significant decrease	SD
β < 0	|Z|<1.96	Not significant decrease	NS
β = 0	|Z|	No trend	NS
β > 0	|Z|<1.96	Not significant increase	NS
β > 0	1.96≤|Z|<2.58	Significant increase	SI
β > 0	|Z|≥2.58	Highly significant increase	SI

*‘SD’ refers to ‘Highly significant decrease’ or ‘Significant decrease’, ‘NS’ refers to ‘Not significant decrease’ or ‘No trend’ or ‘Not significant increase’, and ‘SI’ refers to ‘Significant increase’ or ‘Highly significant increase’. These abbreviations are used in the subsequent analysis of this study.

4. Results and discussion

4.1. Accuracy verification for the lake area series

A direct comparison of the annual total area of lakes over 50 km² obtained in this study with the annual total lake area from the existing dataset provided by Zhang et al. (2019a) is presented in Figure 5. The overall temporal variation trend in the total lake area between them is consistent in most years, and the Pearson's correlation coefficient (r) is approximately 0.997. The improved lake area extraction algorithm takes advantage of the long-term series of Landsat images. There is a slight difference between the lake area results obtained in this study and the existing lake area dataset in some years, mainly attributed to the seasonal fluctuations and interpolation errors in missing values that were induced by too few original images containing lake regions (Figure 5).

Figure 5. The variation trend of annual total lake area in existing datasets, total lake area in this study, the count of Landsat images containing lake and the count of lakes with void annual mean lake area.

Moreover, correlation analysis was conducted between the annual area sequence of each lake and the corresponding annual area sequence in the existing lake dataset, and the correlation coefficient of each lake is presented in Figure 6. The average value of r for 180 lakes is approximately 0.88, and 144 lakes (accounting for approximately 80%) yield an r greater than 0.8. This finding indicates the relatively high accuracy of the extracted lake area. To further explore the reason for the small r value of some lakes, the ratio of the maximum to minimum area of each lake during this period was calculated based on the existing dataset (Figure 6). The lakes with smaller r values mainly corresponded to those with smaller ratios of maximum to minimum area. This result suggests that when the lake area changes over a small range, some small extraction errors may lead to a reduction in the correlation coefficient.

Figure 6. The correlation coefficients between the annual lake area series extracted in this study and the existing lake dataset, as well as the ratio between the maximum and minimum lake areas in the existing lake dataset.

To further verify the accuracy of the lake area series extracted in this study, correlation analyses between the annual average area and annual average water level of each of 52 lakes from 2000 to 2017 (Li et al. 2019a) were performed. The average correlation coefficient of 52 lakes is approximately 0.9, among which 38 lakes (73.08%) have r ≥ 0.9, 45 lakes (86.54%) have r ≥ 0.8, and 50 lakes (96.15%) have r ≥ 0.6 (Figure 7a). The results revealed that the lakes with small r values are mainly those with water level changes by less than 3 m, and the ratio of the maximum to minimum area of these lakes is also small, which are mainly distributed in the southern part of the TP. Taking Qinghai Lake, the largest lake on the TP, as an example, the annual average area and annual average water level of Qinghai Lake from 2000 to 2018 and their correlation are shown in Figure 7b and c. The area results of Qinghai Lake in this study are highly consistent with the results of existing lake datasets, and the correlation coefficient with water level is as high as 0.987.

Figure 7. (a) The correlation coefficients between the annual lake area series extracted in this study and the lake water level dataset, and the range of variation in lake water level (maximum water level minus minimum water level). Taking Qinghai Lake as an example, the variation trends and correlation of annual average area and annual average water level are shown in (b) and (c), respectively.

Overall, the lake area on the TP extracted in this study is reliable, with an average correlation coefficient of approximately 0.9 with both the existing lake area and water level datasets. The existing lake datasets used in this study particularly refer to Zhang's lake area datasets, because they have been visually inspected, manually edited and fully verified, and their accuracy can be guaranteed. In Figure S3, Zhang's lake area datasets and lake water level dataset (Li et al. 2019a) were used as reference values to compare the lake area results obtained by the improved isobath interpolation-based lake area extraction method in this study with those obtained by the original isobath interpolation-based method (Yao et al. 2019). The lake area results of the two methods are highly correlated with Zhang's lake area datasets and water level dataset (all $\overline{r}$ > 0.9), and the lake area results of this study are closer to Zhang's lake area datasets. This may occur because we use the existing lake area datasets (Zhang's lake area datasets) as prior data, which facilitates automatic lake area extraction, greatly improving the efficiency of the algorithm. These prior data also provide a reliable reference for the lake area, reduce abnormal lake area results extracted by automatic algorithms, and improve the lake area extraction accuracy. Moreover, this study overcomes the challenges due to cloud and atmospheric influences (Sogno, Klein, and Kuenzer 2022) in extracting lake areas with a high temporal resolution using Landsat images.

4.2. Multiple change points detection and segmented trend analysis

The variation trend of lake area in the TP during 1986–2020 was not monotonously increasing. To further explore the interannual variation pattern of lake areas, this study carried out multiple change points detection and segmented trend analysis for each lake area sequence.

4.2.1. Multiple change points detection

The BFAST method was applied to the area series data of each lake from 1986 to 2020. This study set a maximum of 3 change points that could be detected in each lake, and the obtained change points all passed the significance test of p ≤ 0.05. Figure 8 shows the distribution of the number of detected change points. The change points are mainly distributed between 1992 and 2013, and there are two peak years: 1999 and 2013. As shown in Figure 5, from 1986 to 2020, the total area of 180 lakes increased from approximately 33,612.89 km² to 41,480.36 km², with a large magnitude of approximately 7,867.47 km² and a high rate of 231.40 km²/a. The lake area changes from 1986 to 1999 showed both increases and decreases, but with a small magnitude of approximately 204.62 km² and a low rate of 15.74 km²/a overall. From 1999 to 2013, the lake area consistently showed an increasing trend, with a large magnitude of approximately 5,544.15 km² and a high rate of 396.01 km²/a. The lake area from 2013 to 2020 initially showed a slight decreasing trend, and then exhibited a rapid increasing trend, with an overall magnitude of approximately 2,118.70 km² and a rate of 302.67 km²/a.

Figure 8. Distribution of BFAST change points count.

Previous studies found that there was a significant increase in lake area both after 1995 and 2015, and these two time points are often used as turning points for studying changes in lake area on the TP (Lei et al. 2019; Zhang et al. 2021b; 2019a). This difference might be caused by the use of the BFAST method, which determines change points based not only on the direction of the trend change but also on the magnitude of the trend change. In addition, the parameter settings of the method may cause slight differences in the results of the change points. Considering both the direction and magnitude of changes in lake area, the two peak years of the change points detected in this study are also reasonable.

4.2.2. Segmented trend analysis based on the multiple change points

Based on the detected change points in each lake area sequence, the segmented M-K test and Sen's slope methods were conducted on the annual lake area series from 1986 to 2020 for each lake. The segmented change trend of the annual lake area series of 180 lakes, classified according to the different sizes, locations, and types (endorheic/exorheic, glacier-fed/non-glacier-fed) of lakes, is presented in Figure 9. According to the average lake area of each lake during 1986–2020 obtained in this study, the lakes were classified into four categories: (1) area ≥ 1000 km² (3 lakes, accounting for approximately 1.7%), (2) 500 ≤ area < 1000 km² (11 lakes, 6.1%), (3) 100 ≤ area < 500 km² (71 lakes, 39.4%), and (4) 50 ≤ area < 100 km² (95 lakes, 52.8%, including 28 lakes with an average area < 50 km², which had an area greater than 50 km² during the period of 1986–2020). The 180 lakes in this study are distributed in the northwest (westerlies domain), central (transition domain), south (Indian monsoon domain), and northeast (East Asian monsoon domain) TP, with 48, 114, 6, and 12 lakes, accounting for 26.7%, 63.3%, 3.3%, and 6.7%, respectively. In addition, based on the geomorphological features around the lakes and the sources of water supply, lakes can be divided into endorheic lakes or exorheic lakes, glacier-fed lakes or non-glacier-fed lakes (Zhang et al. 2021a; Zhang, Zhang, and Zhu 2020b). Among them, there are 132 endorheic lakes and 48 exorheic lakes, accounting for 73.3% and 26.7%, respectively. There are 114 glacier-fed lakes and 66 non-glacier-fed lakes, accounting for 63.3% and 36.7%, respectively.

Figure 9. The segmented trends of the annual lake area series of 180 lakes from 1986 to 2020. (a) the lakes classified by size, (b) the lakes classified by region, (c) the lakes classified by whether they are endorheic or exorheic, and (d) the lakes classified by whether they are glacier-fed or non-glacier-fed.

In Figure 9, the horizontal axis is the lake ID, and the vertical axis is the year. The color of each column from bottom to top shows the change trend of each lake from 1986 to 2020. In Figure 9a, most lakes showed an ‘SI’ trend after 2000, while lakes showed a ‘SI’, ‘SD', or ‘NS’ trend in a similar proportion before 2000. Figure 9b shows the trends of lake area changes in four main regions (as shown in Figure 1). The most obvious finding is that the lakes in the southern TP generally showed a decreasing trend from 1986 to 2020. Most of the lakes with continuous ‘SD' or ‘NS’ trends from 1986 to 2020 are mainly exorheic lakes and glacier-fed lakes (Figure 9c and d).

For the characteristics of different change points of each lake, to uniformly analyze the change in lakes, trend analyses of lake area in the whole period from 1986 to 2020 were carried out (Figure 10a), and segmentation trend analyses were conducted according to the two peak years in 1999 and 2013 (Figure 10b–d). The number of lakes with different classifications and change trends in each period are shown in Table 3.

Figure 10. The trends of lake area variation in (a) 1986–2020, (b) 1986–1999, (c) 1999–2013 and (d) 2013–2020, respectively.

Table 3. The number of lakes with different classifications and trends of changes in each period.

Trends Classifications	1986–2020			1986–1999			1999–2013			2013–2020			Total
	SD	NS	SI	SD	NS	SI	SD	NS	SI	SD	NS	SI
>1000 km²	0(0.0%)	1(33.3%)	2(66.7%)	1(33.3%)	0(0.0%)	2(66.7%)	0(0.0%)	0(0.0%)	3(100.0%)	0(0.0%)	2(66.7%)	1(33.3%)	3
500–1000 km²	1(9.1%)	1(9.1%)	9(81.8%)	3(27.3%)	6(54.5%)	2(18.2%)	1(9.1%)	0(0.0%)	10(90.9%)	0(0.0%)	5(45.5%)	6(54.5%)	11
100–500 km²	4(5.6%)	8(11.3%)	59(83.1%)	25(35.2%)	28(39.4%)	18(25.4%)	5(7.0%)	13(18.3%)	53(74.7%)	3(4.2%)	34(47.9%)	34(47.9%)	71
50–100 km²	5(5.3%)	8(8.4%)	82(86.3%)	32(33.7%)	42(44.2%)	21(22.1%)	7(7.4%)	17(17.9%)	71(74.7%)	2(2.1%)	48(50.5%)	45(47.4%)	95
Northwestern TP	1(2.1%)	3(6.2%)	44(91.7%)	17(35.4%)	23(47.9%)	8(16.7%)	2(4.2%)	10(20.8%)	36(75.0%)	0(0.0%)	16(33.3%)	32(66.7%)	48
Central TP	4(3.5%)	11(9.7%)	99(86.8%)	38(33.3%)	43(37.7%)	33(29.0%)	6(5.2%)	19(16.7%)	89(78.1%)	5(4.4%)	62(54.4%)	47(41.2%)	114
Southern TP	4(66.6%)	1(16.7%)	1(16.7%)	1(16.7%)	4(66.6%)	1(16.7%)	5(83.3%)	0(0.0%)	1(16.7%)	0(0.0%)	6(100.0%)	0(0.0%)	6
Northeastern TP	1(8.3%)	3(25.0%)	8(66.7%)	5(41.7%)	6(50.0%)	1(8.3%)	0(0.0%)	1(8.3%)	11(91.7%)	0(0.0%)	5(41.7%)	7(58.3%)	12
Endorheic	3(2.3%)	8(6.0%)	121(91.7%)	51(38.6%)	47(35.6%)	34(25.8%)	5(3.8%)	9(6.8%)	118(89.4%)	5(3.8%)	51(38.6%)	76(57.6%)	132
Exorheic	7(14.6%)	10(20.8%)	31(64.6%)	10(20.8%)	29(60.4%)	9(18.8%)	8(16.6%)	21(43.8%)	19(39.6%)	0(0.0%)	38(79.2%)	10(20.8%)	48
Glacier-fed	8(7.0%)	14(12.3%)	92(80.7%)	42(36.9%)	47(41.2%)	25(21.9%)	9(7.9%)	21(18.4%)	84(73.7%)	2(1.7%)	50(43.9%)	62(54.4%)	114
Non-glacier-fed	2(3.0%)	4(6.1%)	60(90.9%)	19(28.8%)	29(43.9%)	18(27.3%)	4(6.1%)	9(13.6%)	53(80.3%)	3(4.5%)	39(59.1%)	24(36.4%)	66
Total	10(5.6%)	18(10.0%)	152(84.4%)	61(33.9%)	76(42.2%)	43(23.9%)	13(7.2%)	30(16.7%)	137(76.1%)	5(2.8%)	89(49.4%)	86(47.8%)	180

*The number of lakes (the percentage of lake numbers with different trends under each classification).

The lakes of different sizes showed similar variation trends during the same period (Table 3). From 1986 to 2020, most lakes (accounting for approximately 84.4%) exhibited the ‘SI’ trend regardless of the size of the lakes. From 1986 to 1999, approximately 42.2% of all lakes showed the ‘NS’ trend, and 33.9% showed the ‘SD' trend. Approximately 76.1% of all lakes showed the ‘SI’ trend from 1999 to 2013. For lakes ≥ 1000 km² and lakes (500–1000 km²), higher proportions of the ‘SI’ trend (100% and 90.9%) were found. However, from 2013 to 2020, approximately 49.4% of all lakes showed the ‘NS’ trend, and 47.8% showed the ‘SI’ trend, with similar variation trends for different lake sizes. During 1986–2020, 1986–1999, 1999–2013 and 2013–2020, the main variation trends in lake area (≥ 50 km²) on the TP can be summarized as ‘SI’, ‘NS’ or ‘SD', ‘SI’, ‘NS’ or ‘SI’ (Figure 10 and Table 3).

For lakes in different regions, from 1986 to 2020, most lakes in the northwestern TP (westerlies domain) (approximately 91.7%), central TP (transition domain) (86.8%), and northeastern TP (East Asian monsoon domain) (66.7%) showed the ‘SI’ trend, but in the southern TP (Indian monsoon domain), they mostly (66.6%) showed the ‘SD' trend (Table 3). This finding is similar to the results of most previous studies, such as Zhang et al. (2019a) and Song et al. (2014). Lakes in any region mainly tended to be ‘NS’ from 1986 to 1999, with the second most common trend being ‘SD' and the least common being ‘SI’. From 1999 to 2013, most of the lakes in the northwestern TP (75.0%), central TP (78.1%), and northeastern TP (91.7%) exhibited the ‘SI’ trend, while most of the lakes in the southern TP (83.3%) showed the ‘SD' trend. From 2013 to 2020, the lakes in the northwestern TP and northeastern TP mainly showed the ‘SI’ trend (66.7% and 58.3%, respectively), followed by the ‘NS’ trend (33.3%, 41.7%). In the central TP, the lakes mainly showed the ‘NS’ trend (54.4%) and the ‘SI’ trend (41.2%). All the lakes in the southern TP presented the ‘NS’ trend. The mid-latitude westerlies, Indian monsoon, and East Asian monsoon transport water vapor to the TP (Bothe, Fraedrich, and Zhu 2009), and all these large-scale atmospheric circulation anomalies can influence water vapor convergence (Zhang, Zhang, and Zhu 2020b). Among the four regions classified in this study, the interannual variations in lake area in the northwestern and central TP exhibit the most similar characteristics. This similarity may be attributed to the strong influence of the westerlies in these two regions. The results revealed that large-scale atmospheric circulation systems, such as the westerlies, Indian monsoon, transition in between, and East Asian monsoon, might control the pattern of interannual lake area variations in different regions (Figure 10 and Table 3).

The endorheic lakes showed a larger proportion of significant changes than exorheic lakes during the same period (Table 3). During 1986–2020, approximately 91.7% of endorheic lakes exhibited the ‘SI’ trend, 64.6% of exorheic lakes showed the ‘SI’ trend and 20.8% showed the ‘NS’ trend. From 1986 to 1999, ‘SD' and ‘NS’ trends were mainly found in the endorheic lakes (38.6% and 35.6%, respectively). The exorheic lakes mainly exhibited the ‘NS’ trend (60.4%), with similar proportions of the ‘SD' (20.8%) and ‘SI’ (18.8%) trends. From 1999 to 2013, most of the endorheic lakes had the ‘SI’ trend (89.4%), while the exorheic lakes mainly presented ‘NS’ (43.8%) and ‘SI’ (39.6%) trends. From 2013 to 2020, most endorheic lakes showed ‘SI’ and ‘NS’ trends (57.6% and 38.6%, respectively); for the exorheic lakes, 79.2% exhibited the ‘NS’ trend, and 20.8% exhibited the ‘SI’ trend. These results showed that the interannual area variation patterns of the endorheic lakes are more pronounced than those of the exorheic lakes (Table 3). This finding is consistent with the results of other studies: Zhang et al. (2021a) and Zhang, Zhang, and Zhu (2020b) reported that the multiyear and seasonal changes in endorheic lakes are stronger than those in exorheic lakes and may be attributed to differences between the closed nature of endorheic lakes and the water inflows and outflows of exorheic lakes.

Most glacier-fed (80.7%) and non-glacier-fed (90.9%) lakes exhibited the ‘SI’ trend from 1986 to 2020 (Table 3). From 1986 to 1999, the glacier-fed lakes mainly showed ‘NS’ (41.2%) and ‘SD' (36.9%) trends. Most of the non-glacier-fed lakes had the ‘NS’ trend (43.9%), with a similar proportion of the ‘SD' (28.8%) or ‘SI’ (27.3%) trends. From 1999 to 2013, both the glacier-fed (73.7%) and non-glacier-fed lakes (80.3%) mainly presented the ‘SI’ trend. From 2013 to 2020, the glacier-fed and non-glacier-fed lakes mostly exhibited ‘SI’ (54.4% and 36.4%, respectively) or ‘NS’ (43.9%, 59.1%) trends. These results revealed that the climatic conditions in different periods may lead to differences in the melting rates of glaciers, which in turn leads to some differences in the variation trends of glacier-fed and non-glacier-fed lakes (Table 3).

4.3. Relationship between lake area variation and climatic variables

The climatic variables that affect changes in lake area mainly include air temperature, precipitation, evaporation, and the teleconnection effect of El Niño (Yeh et al. 2009; Zhang et al. 2020a). According to the classification of lakes based on size, location, and type (endorheic/exorheic, glacier-fed/non-glacier-fed), the annual average area of different classified lakes (mean lake area /km²) was calculated. The average air temperature (T/°C), average total precipitation (P /mm), average total evaporation (E /mm), and average values of total precipitation minus total evaporation (P-E /mm) within the ROIs of different classifications of lakes are presented in Figures 11–13, respectively. Table 4 lists the trends of the area variation in all lakes and different classified lakes and corresponding climatic variables during 1986–2020. The correlation coefficients between the annual lake area series and the annual series of T, P, E, and P-E of all lakes and different classifications of lakes from 1986 to 2020 are shown in Figure 14.

Figure 11. The average values of lake area (mean lake area /km²), the average values of mean air temperature (T/°C), the average values of total precipitation (P /mm), the average values of total evaporation (E /mm), and the average values of total precipitation minus total evaporation (P-E /mm) of the total and different size classified lakes for each year, during 1986–2020. The yellow horizontal line represents the P-E value of 0. Additionally, the Niño 3.4 index curve from January 1985 to December 2021, and the red stripes represent El Niño years and the blue stripes represent La Niña years.

Figure 12. The average values of lake area and multiple climatic variables of the total and different region classified lakes for each year, during 1986–2020. (Refer to Figure 11 for others).

Figure 13. The average values of lake area and multiple climatic variables of the total and different types classified lakes for each year, during 1986–2020. (Refer to Figure 11 for others).

Figure 14. Correlation coefficients between mean area series of the total and each classified lakes, and climatic variables (T, P, E, and P-E), during 1986–2020. (a) All lakes and different size classified lakes. (b) All lakes and different region classified lakes. (c) All lakes and different types classified lakes.

Table 4. The variation trends of the area of all lakes and different classified lakes and corresponding climatic variables during 1986–2020.

Variables Classifications	Area	T	P	E	P-E
All lakes	Highly significant increase	Significant increase	Highly significant increase	Highly significant increase	Significant increase
>1000 km²	Highly significant increase	Highly significant increase	Not significant increase	Highly significant increase	Not significant increase
500–1000 km²	Highly significant increase	Highly significant increase	Highly significant increase	Highly significant increase	Significant increase
100–500 km²	Highly significant increase	Significant increase	Highly significant increase	Highly significant increase	Significant increase
50–100 km²	Highly significant increase	Not significant increase	Highly significant increase	Highly significant increase	Significant increase
Northwestern TP	Highly significant increase	Not significant increase	Highly significant increase	Highly significant increase	Not significant increase
Central TP	Highly significant increase	Not significant increase	Significant increase	Highly significant increase	Not significant increase
Southern TP	Highly significant decrease	Significant increase	Significant increase	Not significant increase	Not significant increase
Northeastern TP	Significant increase	Significant increase	Significant increase	Significant increase	Not significant increase
Endorheic	Highly significant increase	Significant increase	Highly significant increase	Highly significant increase	Significant increase
Exorheic	Highly significant increase	Significant increase	Highly significant increase	Highly significant increase	Significant increase
Glacier-fed	Highly significant increase	Significant increase	Highly significant increase	Highly significant increase	Significant increase
Non-glacier-fed	Highly significant increase	Not significant increase	Highly significant increase	Highly significant increase	Not significant increase

From 1986 to 2020, the mean area of all lakes increased highly significantly (Figures 11–13 and Table 4). Moreover, the corresponding climatic variables (mainly T, P, E, and P-E) of all lakes also showed a significant increasing trend. The correlation between the lake area and E is the strongest (r ≈ 0.61), followed by P (r ≈ 0.39), T (r ≈ 0.38), and P-E (r ≈ 0.22), as shown in Figure 14. Evaporation usually has a negative effect on the lake area, but the overall lake area is positively correlated with evaporation. The possible reason is that the change in lake area may be dominated by other factors and less influenced by evaporation, while the change in evaporation is more influenced by the lake area. The increase or decrease in lake area has an effect on the increase or decrease in E in the lake ROI, so the two are mostly positively correlated. Thus, E is not considered the main factor affecting the change in lake area in the subsequent study. For all lakes, the main climatic variables affecting lake area change are precipitation and air temperature. The effect of precipitation on the lake area is mainly through direct precipitation in the lake region and inflows of surface runoff. Additionally, air temperature usually affects the change in lake area by influencing summer precipitation, glacier/snow melting, permafrost degradation, and then changing runoff into lakes (Zhang et al. 2020a). In El Niño years, precipitation tends to decrease, and the lake area tends to decrease or slowly increase. In La Niña years, precipitation tends to increase, and the lake area tends to increase (Figures 11–13). After approximately 1995, the values of P-E within the lake ROIs were nearly all positive, and the lake area showed an upward trend from 1995. This result indicates that the positive accumulation of P-E promotes the expansion of lakes.

The interannual variation in the lake area exhibited a similar increasing trend in lakes of different sizes from 1986 to 2020 (Figure 11 and Table 4). The larger the lakes are, the larger the increasing amplitude is. Among them, the area variation trend of the lakes (100–500 km²) is highly consistent with that of the overall lakes. The climatic variables corresponding to different size classes of lakes all showed significant or nonsignificant increasing trends from 1986 to 2020. As shown in Figure 14a, among the correlation coefficients between the area of lakes (≥ 1000 km²) and climatic variables, the correlation coefficient with T is the highest (0.52), followed by E (0.51), P (0.13) and P-E (−0.01). For the lakes with areas of 500–1000 km² and 100–500 km², the correlation coefficients between lake area and E were the highest (0.48 and 0.61, respectively), followed by T (0.45, 0.39), P (0.42, 0.36) and P-E (0.28, 0.19). The area of lakes (50–100 km²) had the highest correlation coefficient with E (0.62), followed by P (0.44), T (0.35) and P-E (0.27). This finding indicates that the main climatic variable that has a great effect on the area change in the large lake (>100 km²) is air temperature, followed by precipitation. Precipitation is the main climatic variable affecting the area change in small lakes (50–100 km²), followed by air temperature. A more notable feature is that the larger the lake is, the greater the correlation between its area and T is, possibly because air temperature changes have a more significant influence on lakes with larger areas. The underlying reason may be that larger lakes have larger lake basins and reach peaks later than small lakes (Zhang, Zhang, and Zhu 2020b). Therefore, when the temperature rises by a certain degree, the area growth of large lakes is greater and more significant than that of small lakes.

The interannual variations in lake area of the southern TP (Indian monsoon domain) and northeastern TP (East Asian monsoon domain) lakes are quite different from those of all lakes. From 1986 to 2020, the lake areas of the northeastern TP had large amplitudes of changes, and those of the southern TP even showed a highly significant decreasing trend (Figure 12 and Table 4). The climatic variables corresponding to different location-classified lakes all showed significant or nonsignificant increasing trends from 1986 to 2020. For the lakes in the southern TP, T, P, E, and P-E had relatively high values, but negative correlations with the annual mean lake area were found at −0.28, −0.08, −0.18 and −0.03, respectively (Figure 14b). Based on the lake water balance, the decrease in lake area in the southern TP may be attributed to the excessive outflows of surface and subsurface runoff. Chen et al. (2021) and Liu et al. (2023b) reported that the shrinkage of lakes in the southern TP is related to water leakage at the lake bottom and the existence of leakage channels in the rift zones. From the perspective of remote sensing big data analysis, this study provides indirect evidence that the decrease of lake area in the southern TP may not be driven by climate, but is probably related to the leakage of the lake bottom. In the northwestern TP (westerlies domain), the correlation coefficient between the area and E is the highest (0.54), followed by P (0.45), T (0.29) and P-E (0.28). Similarly, the correlation coefficient between the area and E is the highest (0.62) in the Central TP (Transition domain), followed by T (0.37), P (0.35) and P-E (0.18). In the northeastern TP (East Asian monsoon domain), among the correlation coefficients between the area and climatic variables, P is the highest (0.31), followed by E and P-E (both 0.24) and T (0.2). These results indicate that the lakes in the northwestern TP and northeastern TP, where the westerlies or East Asian summer monsoon bring abundant precipitation, may be influenced more by precipitation than by other climatic variables. Additionally, the lakes in the Central TP, where the westerlies and Indian monsoon co-occur, may be mainly affected by air temperature and precipitation at the same time.

During 1986–2020, the interannual variations in lake area of endorheic and exorheic lakes all exhibited increasing trends, while the amplitude of the area changes in endorheic lakes was larger than that in exorheic lakes. The corresponding climatic variables in the endorheic and exorheic lakes showed similar significant increasing trends (Figure 13 and Table 4). For the endorheic lakes, the correlation coefficient between the area and E is the highest (0.62), followed by T and P (both 0.37) and P-E (0.19). Similarly, the correlation coefficient between the area and E is the highest (0.61) of the exorheic lakes, followed by P (0.55), P-E (0.43) and T (0.31) (Figure 14c). The results reveal that in the long-term interannual variation in lake area, endorheic lakes are mainly affected by precipitation and air temperature, while exorheic lakes are mainly affected by precipitation. The area of exorheic lakes is more related to precipitation than that of endorheic lakes, but the changes in endorheic lakes are stronger than those in exorheic lakes; that is, the amplitudes of changes are larger. This result is attributed to the difference between the natural characteristics (e.g. whether there is water flowing out) of the endorheic lakes and the exorheic lakes (Zhang et al. 2021a; Zhang, Zhang, and Zhu 2020b). There are no outflows of surface runoff in the endorheic lakes, and the inflows mainly caused by rivers, precipitation and glacier/snow meltwaters can be conserved. However, exorheic lakes have outflows of surface runoff; even if the inflows are extremely large, they will be discharged, so the changes in lake area are not strong.

The glacier-fed lakes and non-glacier-fed lakes presented similar increasing trends in interannual areas from 1986 to 2020. The amplitude of the area changes was larger for glacier-fed lakes than non-glacier-fed lakes. The climatic variables corresponding to glacier-fed lakes and non-glacier-fed lakes all showed significant or nonsignificant increasing trends (Figure 13 and Table 4). The annual mean area of the glacier-fed lakes is generally higher than that of the non-glacier-fed lakes, while all the corresponding climatic variables of the non-glacier-fed lakes are generally higher than those of the glacier-fed lakes. Among the correlation coefficients between lake area and climatic variables in both glacier-fed and non-glacier-fed lakes, the correlation coefficients between lake area and E are the highest (0.57 and 0.67, respectively), followed by P (0.39, 0.40), T (0.38, 0.34) and P-E (0.24, 0.21), as shown in Figure 14c. This finding indicates that the areas of both glacier-fed and non-glacier-fed lakes are mainly affected by precipitation, followed by air temperature. The influence of air temperature on glacier-fed lakes is greater than that on non-glacier-fed lakes because of its effect on glacier melting.

5. Conclusions

This study proposed an improved lake area extraction method based on isobath interpolation, which can effectively extract the complete lake area from optical remote sensing images, particularly when the lake is partly contaminated by cloud cover. The lake area extraction method primarily improves the acquisition scheme for isobaths and their accuracy, especially for complex lake conditions. More complete and accurate isobaths were obtained as prior information for the subsequent automatic extraction procedure. The proposed method effectively solves the cloud contamination problem and improves data utilization, calculation efficiency and accuracy. The annual results of lake area obtained in this study are conducive to the highly refined analysis of their change trends.

The interannual area series from 1986 to 2020 for lakes ≥ 50 km² on the TP were obtained using the proposed extraction method based on Landsat big data mining, along with the change points and segmented trends for each lake. Overall, most of the lakes on the TP showed significant increasing trends (84.4%) from 1986 to 2020. During 1986–1999, most of the lakes exhibited no significant change (42.2%) or significant decreasing trends (33.9%), while during 1999–2013, most of the lakes showed significant increasing trends (76.1%). From 2013 to 2020, most of the lakes exhibited no significant change (49.4%) or significant increasing trends (47.8%).

Combined with meteorological data for preliminary qualitative analysis, it was found that the corresponding climatic variables (mainly T, P, E, and P-E) of all lakes and different classified lakes showed a significant or nonsignificant increasing trend. Precipitation and air temperature are the two main climatic variables affecting interannual variations in lake area. El Niño and La Niña events affect precipitation on the TP, thereby influencing lake area changes. The overall expansion trend of lake area after 1995 can be attributed to the positive accumulation of P-E. There is a slight difference in the influence of climatic variables on lake area changes in different sizes, regions and types (endorheic/exorheic, glacier-fed/non-glacier-fed) on the TP. The areas of larger lakes (>100 km²) are more significantly affected by air temperature than those of smaller lakes (50–100 km²). The large-scale atmospheric circulation systems acting on different regions, such as the westerlies, Indian monsoon, transition in between and East Asian monsoon, were found to affect the interannual variation pattern of lake area. Importantly, the decrease in lake area within the southern TP may not be climate-driven but may be related to leakage from the lake bottom. The interannual area of the exorheic lakes is more related to precipitation than that of the endorheic lakes, but the interannual variations in the area of the endorheic lakes are stronger than those of the exorheic lakes. Glacier-fed lakes are more affected by air temperature than non-glacier-fed lakes, along with stronger interannual area variations.

In summary, this study provides a new method, data, and characteristic results for research on the interannual variations in lake areas on the TP, which can be used by research communities and various hydrological and water resource management departments for practical applications. The limitations of the method proposed in this study are that the isobaths need to be preprocessed, and the isobaths need to be updated if the study time range is extended, especially the isobaths corresponding to the maximum and minimum lake areas. In the future, it is necessary to further improve the isobath extraction method and fully utilize the existing high-precision lake area datasets and multisource remote sensing big data. More high-precision isobaths need to be obtained to further study the seasonal variations in lake area on the TP. Additionally, the relative contributions of various climate variables to the lake area need to be more precisely assessed through in-depth water balance analysis.

Acknowledgements

The authors would like to acknowledge all the scientists and engineers who contributed data at the United States Geological Survey, the National Tibetan Plateau Data Center, the European Centre for Medium-Range Weather Forecasts, and the Climate Prediction Center of the National Oceanic and Atmospheric Administration. The authors also thank Google for providing free access to the Google Earth Engine platform and Esri for providing the service of World Hillshade basemap, and finally, thank anonymous reviewers and the editors for supporting this work.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Data availability statement

The key codes for extracting lake area in this study are available to interested users: https://code.earthengine.google.com/17ce545ea551f0933587100715e9921d, https://code.earthengine.google.com/b9276a3fc78bfb8783ef38c8aaa049c9. The lake area dataset produced in this study was publicly published on the National Tibetan Plateau Data Center (https://doi.org/10.11888/Terre.tpdc.300883).

Additional information

Funding

This work was supported by the Strategic Priority Research Program of the Chinese Academy of Sciences [grant number XDA19070104], and the National Science and Technology Major Project of China's High Resolution Earth Observation System [grant number 21-Y20B01-9001-19/22].