Lake area monitoring based on land surface temperature in the Tibetan Plateau from 2000 to 2018

Lake area change in the Tibetan plateau is an important indicator for climate change assessment. To overcome the temporal inconsistency of optical remote sensing-based lake area detections, a land surface temperature (LST)-based detection scheme was proposed by utilizing the big difference between land and water surface temperatures. A trend test conducted by the Mann–Kendall (MK) method was successfully applied to investigate lake area variation from 2000 to 2018 with the use of the annual mean temperature information derived from the Moderate Resolution Imaging Spectroradiometer (MODIS) LST daily product. A comparison with the monitoring results from Landsat images indicates the proposed method can provide spatial distributions of lake area change with high accuracy. More importantly, the temporal variation of annual mean LST provides a special way to detect the abrupt change year (ACY) in lake area.The ACYs of most lakes mainly occur from 2004 to 2012. For an individual lake, the ACY offers vital information about the lake area change process. In summary, this work demonstrates the good potential of the LST-based method for lake area monitoring and assessment.


Introduction
The Tibetan Plateau (TP) is the highest and most extensive highland in the world. It plays an important role in global climate and significantly affects the climate system over the surrounding areas. Meanwhile, it is the source of many major Asian rivers (e.g. the Indus, Brahmaputra, Mekong, Yangtze and Yellow rivers and others), sustaining more than 1 billion people living both in the region and downstream (Barnett et al 2005). The extreme high and cold environment of this special region develops many plateau lakes, which are widely distributed throughout the plateau (figure 1). These lakes play key roles in the hydrological cycle and greatly affect regional climate, water resource, and terrestrial and aquatic ecosystems.
Recently, the TP experienced evident climate changes (Yang et al 2014). The trend analysis of air temperature indicates an average increase of 0.28 • C per decade (Guo and Wang 2012). As a special landscape, the TP lake offers a very good proxy to investigate the impacts from climate change because of few disturbances from human activities. Obviously, lake water surface temperature is a direct indictor to present this effect (Wan et al 2018). Besides this variable, lake area is one of the most popularly investigated items because most lakes are endorheic lakes in this region. Up to now, many studies have conducted to investigate lake area changes in the TP (Song et al 2014b, Mao et al 2018, Zhang et al 2018, and generally extended to water volume or storage monitoring (Li et al 2019, Qiao et al 2019, Zhang et al 2019a. In these studies, remote sensing observations become key data sources to quantify these variations. The normalized difference water index (NDWI) (Mcfeeters 1996), developed based on surface reflectance at green band and near infrared band, was commonly used to separate water bodies from other land surfaces to analyze the temporal variation in the lake area.
However, optical remote sensing-based lake area extraction is seriously obscured by cloud cover, especially in the monsoon season. The cloud cover assessment based on the Himawari-8 satellite data indicated the averaged cloud-cover frequency in the TP has a typical annual cycle with high value in the summer monsoon season and low value in the winter and the mean value is above 35% (Shang et al 2018). Meanwhile, the lake water level also shows strong seasonal changing pattern in response to summer monsoon precipitation, snowmelt, and evaporation (Phan et al 2012, Song et al 2014a, and the seasonal pattern varies in different lakes (Lei et al 2017). Combined with the poor temporal resolution of the high-resolution satellite observations, the above issues jointly hamper the availability of optical remote sensing images with good temporal consistency. Consequently, it further influences accurate detections for the temporal changing pattern of lake area. How to resolve the temporal inconsistency on remote sensing data induced by cloud cover and get effective knowledge about lake area temporal change pattern are still challengeable. In this study, based on the significant differences in land surface temperature (LST) of water and land surfaces, an LST-based detection scheme based on the LST time series is proposed for monitoring lake area variation and its temporal characteristic in the TP.

Remote sensing data
Currently, there are many remote sensors providing LST observations via thermal infrared (TIR) bands, including Thermal Infrared Sensor (TIRS), Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER), Moderate Resolution Imaging Spectroradiometer (MODIS), and Spinning Enhanced Visible and Infrared Imager (SEVIRI). Compared with the poor temporal coverage of the high-resolution LST estimates from TIRS and ASTER, the MODIS LST daily product offers daytime and nighttime LSTs every day in most places with long-term records (since 2000). Although the spatial resolution of this product (1 km) is a little coarser than those of TIRS and ASTER, it could provide more detailed information than the geostationary observations. In addition, previous validation studies have indicated that the accuracy of the MODIS LST product is better than 1 K (Wan 2008, Li et al 2013, Duan et al 2019. Therefore, the MODIS LST daily product was a good data source for time series analysis, and the MODIS/Terra LST daily product (MOD11A1) from 2000 to 2018 was collected in this study.
In addition, Landsat images in 2000 and 2018 were used to evaluate the LST-based detection results. To alleviate the seasonal changes in lake area induced by the monsoon climate, the Landsat images acquired before the monsoon season were collected in this study (table 1). Similar as the previous studies, the threshold method with the use of NDWI was applied for lake area detection. The threshold value to delineate water surfaces was determined via the NDWI histogram of each image.

Annual mean land surface temperature estimation
Because of high temporal variability in LST, the LST values of a certain pixel at the same date but in different years cannot be directly compared to indicate whether there is a warming trend or not. Comparatively, the periodic information of LST, such as annual mean LST (LST m ), should be a good choice to perform such comparison because it covers general surface thermal property over one year.
It is well known that the annual temperature cycle (ATC) of LST can be depicted as a series of sinusoidal (harmonic) functions (Bechtel 2012, Fu and Weng 2016, Zhao et al 2019a. To ensure the simplicity and accuracy of the ATC simulation, the ATC model, expressed as two sinusoidal functions, was used in this study: T s (t) = T 0 + a 1 sin ωt + b 1 cos ωt + a 2 sin 2ωt +b 2 cos 2ωt (1) where T s (t) is the daily LSTs, t is the day of the year, a 1 , b 1 , a 2 , and b 2 are the coefficients for the primary component and the second component, and T 0 is the annual mean LST (LST m ). ω is a constant term which can be derived by 2π/d, where d represents the full cycle and it was set to 365 d (intercalary days were ignored). To reduce the impact from the LST estimation, the fitting process of each pixel was conducted with the LST values with good quality (QC = 0) according to the quality control file.
To present the fitting accuracy of the ATC model, four group sampling pixels as shown in figure 1, including both water pixel and its nearby land pixel, were selected as examples. These sampling pixels are widely distributed around the plateau with different altitudes and latitudes. The daily LSTs of these eight pixels in 2018 were collected and the cloud-cover days with filled value and the days with poor LST quality were excluded. As shown in figure 2, there are relatively few LST observations in summer due to frequent cloud cover at this monsoon season. The solid lines present the fitted annual cycle of LST. The fitted lines well capture the annual changing pattern of LST. The dashed lines indicate the annual mean LST values of these eight pixels. Obviously, the daily temperature difference between land and water pixels varies a lot at different seasons for each group. At annual scale, the temperature difference between water and land pixels at each group varies from 14.86 K to 24.47 K. For water pixels, the altitude and latitude differences present limited influences on the annual mean temperature. For example, water pixel E has the lowest LST m due to its highest altitudes (4904 m). However, the difference with other water pixels is relatively small (within 4 K). Comparatively, the LST m values of the land pixels are systematically higher than the water pixels but with higher variability. The variation is not only influenced by pixel geolocation, and land cover characteristic also plays an important role by varying surface thermal property. However, from all these sampling pixels, it can be concluded that the difference in annual mean surface temperature (more than 14 K) is large enough to discriminate water surface and land surface at most circumstances.

LST-based detection method 2.3.1. Lake change area identification.
From the significant differences in thermal characteristics of water and land surfaces, LST m will present strong responses in temperature increase or decrease under the circumstances of lake shrinkage or expansion. Therefore, the lake area change can be identified by performing trend analysis of the time series of LST m . The popularly Mann-Kendall (MK) trend test was selected as the primary tool to analyze this variation. As a non-parametric test, the MK test has been widely applied for trend detection in hydrometeorological variables ( The statistic S is calculated as: where n is the length of data, x i and x j are two generic sequential data values, and sgn(x j -x i ) is the signum function with following expression: The variance of the statistic is computed as:

18
(4) where m is the number of tied groups and t i denotes the number of ties of extent. In cases where the sample size n > 10, the standard normal test statistic Z S is computed as: Positive Z S value indicates an increasing trend and vice versa. No trend null hypothesis was rejected for an absolute value of Z S greater than Z 1 −α/2 and a significant trend exists in the time series. The α was set to 0.05 with Z 1−α/2 of 1.96.
In this study, a 10-km buffer zone was generated first for each lake to get a wide region covering the lake area and its nearby land areas. Then, the MK test was conducted for each pixel in the zone. For a pixel with a significant trend, its change rate (R) is derived via the Sen's slope (1968).As indicated by equation (6), the pixel with positive slope value was identified as lake shrinkage and the pixel with negative slope value was identified as lake expansion.
During the detection, the permanent water and land pixels during this period but with significant changing trends were manually masked out to exclude their influences.

Lake area abrupt change year detection.
The shifts in LST time series is closely connected with lake area change. To detection the start of the lake area changes in the period, the non-parametric Mann-Whitney (MW) test, developed by Pettitt (1979), was employed in this study. This method has been widely used in change point detection for hydrologic or climatic time series (Caloiero et al 2011, Salehi et al 2020). The key statistic of each year t in this method is expressed as: For LST time series with no change point, there will be a continually increasing trend for the absolute value of U t . When there is a change point at year K, there will be a local maximum value K n : The variation is significant when the probability (p) listed as follows lower than 0.05.

Performance of the LST-based change detection method
To reveal the performance of the proposed method, six lakes (Nam Co, Qinghai Lake, Selin Co, Ulan Ul Lake, Zharinam Co, and Ayakkum Lake) in the TP were selected as examples (table 1). These lakes have been widely analyzed in previous studies and they present different changing extents. Therefore, the evaluations in these lakes will be helpful to assess the accuracy of the proposed method at different changing situations. Table 1 lists the general information of these lakes, including the dates of the Landsat images used for the lake area extraction and the corresponding lake area in 2018. Figure 3 shows the annual mean land surface temperature of the six lakes in 2018. The temperature differences between water and land pixels are obvious in all lakes and the temperature boundary well follows the lake extent detected based on Landsat-8 images in 2018.
With the estimated LST m values of the lake areas from 2000 to 2018, the MK trend analysis was applied to detect the temporal variation in LST m and the pixels passing significant test (including increasing and decreasing) was figured out with the estimated change rate values. Figure 4 shows the lake area depicted via Landsat observations in 2000 and 2018 and the change rate of the detected change pixels of the six selected lakes. The pixels without significant variations are denoted as white background. The blue and red line represent the lake extents derived from Landsat observations in 2000 and 2018, respectively. For Nam Co as shown in figure 4(a), there is a very small variation in lake area. Therefore, the detection result indicates that there are almost no MODIS pixels with significant change during this period. Comparatively, the variation in lake area in other five lakes is significant, especially for Selin Co, Ulan Ul Lake, and Ayakkum Lake. All these five lakes generally show expansion pattern with negative slope values because of the decrease in LST from land surface to water surface. Moreover, no matter for places with relatively small variation such as Qinghai Lake and Zharinam Lake nor the region with big area variation, the LST-based detections are generally in line with the changes based on Landsat observations from spatial distribution. Only a few pixels present inconsistency with the variation in lake area and the reason should be partly attributed to the variation in thermal environment of the land pixels close to the lake boundary. Meanwhile, the temporal inconsistency between the annual temperature-based detections and the investigations before the monsoon with Landsat images is another reason. Additionally, the coarse resolution of the MODIS LST data partly reduces the sensitivity of surface temperature to lake area variation and then influences the spatial consistency. However, the overall distribution pattern generally reveals that LST variation is a good proxy to present lake area change and in turn good accuracy for lake variation identification. Figure 5 shows the detection result of the lake area variation in the TP with the proposed method. There are 238 lakes exhibiting significant changes in lake area during this period. Most of the changing lakes have an expanding status with red color, while the lakes with shrinkage status (noted with blue color) are quite limited. The location of the six selected lakes are also shown in this figure. The shrinking lakes are majorly located in the south part of the TP, such as the Yamzhog Yum Co, Pumo Yum Co, Como Chamling Lake, and Duoqing Co shown in figure  Meanwhile, the Zhuonai Lake in the north part is also detected with significant decreasing trend in lake area, which is in line with the finding in (Liu et al 2019). It partly demonstrates the good accuracy of the proposed method in detecting lake shrinkage.

Lake area variation from 2000 to 2018
From the detection results, the number of pixels located in the lake region with significant decreasing trend in LST m is 15 425, with the area more than 15 000 km 2 . Although the area cannot reflect the exact area of water surface change due to the mixed pixel effect in the MODIS coarse pixel, the results can precisely indicate the location of the variations. Similar condition should be appropriate to the lake area shrinkage. The number of pixels is only 1298 with significant increasing trend, largely smaller than that  of the opposite direction. Therefore, the TP generally experience a significant lake area expansion period, and the reason should be connected with the variation in precipitation, surface evaporation, and glacial meltwater (Qiao et al 2019).
Furthermore, the altitudinal distribution of the change area was analyzed based on surface elevation data. As shown in figure 6, the expansion lakes were majorly located at high elevation area from 4200 m to 5200 m, while the shrinkage lakes were mainly located from 4400 m to 4700 m. The distribution generally obeys the lake distribution in the TP as shown with the green line in figure 6. However, at different elevations, the magnitude of lake changing area varies a lot. The variation pattern of the shrinking lakes is not clear due to its relatively small amount. However, there is a general increasing trend with altitude in the lake expanding area as indicated by the arrow in figure 6. It implies that the lakes located from 4750 m to 4900 m present most significant variation. In turn, the impacts from climate change to lake area variation is the strongest at this elevation range.

Lake area abrupt change year from 2000 to 2018
In addition to lake area monitoring, the abrupt change point of lake area variation during this period was detected with the MW test. Figure 7 shows the abrupt change year (ACY) detection results of two typical pixels. The ACYs are denoted for each time series which exactly present the changing characteristics. With the detection method, the ACYs of the six selected lakes are shown in figure 8. Nam Co is not considered in this part because of its limited variation during this period. The color bar represents the variation of the change point from 2004 to 2012 because of the majority distribution in this range. The lake extent in 2018 is also added in this figure.
In figure 8, the ACY varies a lot at different lakes because of their different climatic environments, land cover patterns, and glacier or snow coverages. However, a common phenomenon is easy to be observed that ACY has an increasing trend from inner lake region to lake border region. These findings reveal the expansion characteristics that the inner lake  region firstly experiences the variation from land to water surface, resulting a significant reduction in LST. Then, the change area gradually spreads to the outside pixels of the original lake extent. This variation represents a homogenous changing pattern, which can be clearly observed in these lakes shown as the arrows. In Selin Co ( figure 8(b)) and Ayakkum Lake ( figure 8(d)), the variation is most significant because of their large change areas. The lake area investigation has found the Ayakkum Lake has a significant increase in 2008 (Jing et al 2018), and this result is consistent with the finding in this study. Most the ACYs in this lake is 2008. For Selin Co, 2004 and 2007 occupy a large amount, reflecting big variation around these two years. It is confirmed by the change detection conducted by Yi and Zhang (2015) and Zhang et al (2013) To recognize the overall temporal change characteristics, figure 9 shows the histogram plots of the ACYs for lakes with expansion and shrinkage, respectively. The ACY mainly changes from 2004 to 2012, with the majority happens from 2007 to 2010. The expansion and shrinkage have similar yearly distribution. The number implies that during these years, the changes in local climatic and environmental conditions greatly impact water surface area. The results are in line with the investigation with the ICESat and ICESat-2 observations that the lake level changes majorly happen in 2004 to 2018 (Zhang et al 2019a). Meanwhile, the periodic research also indicates that there is an increase rate (298.53 km 2 yr −1 ) in water surface area from 2005 to 2013 than that from 2000 to 2005 (Qiao et al 2019). The finding also confirms the detection in this study.

Discussion
The proposed method presents a novel way for monitoring lake area in the TP. The case study and the general spatial distribution pattern all indicate the good performance of this method in investigating lake area variation.
Unlike previous methods based on optical remote sensing data, this study introduces the use of land surface temperature, which presents significant differences for land and water pixels in the TP. To avoid the impact from cloud cover for daily LST observations, we applied the ATC model to derive periodic information from the time series of daily LST data. The annual mean LST was extracted from the fitting results to represent the annual temperature information and it well kept the obvious differences between land and water pixels. Therefore, the use of this term for lake area detection partly discards the uncertainty arose by the temporal inconsistency when conducting such detection with optical remote sensing data. The detected distribution pattern of lake area change indicates that most lakes are experiencing a significant expansion but few lakes in the southern TP are exhibiting shrinking trends. The results well agree with the findings in previous studies (Liu et al 2019, Qiao et al 2019. More importantly, the study provides key temporal information about the abrupt change information in lake area variation, which is rarely discussed in former analysis. This information will be more helpful to analyze the mechanism behind the changes. Besides the advantages, there are still some issues should be emphasized. The most important issue comes from the MODIS LST product. It has been widely used in many fields. However, the high vulnerability of thermal infrared remote sensing to cloud cover results many gaps in the daily LST data. The inclusion of available LST data from all days will assist more accurate derivation of the annual temperature parameters than the estimation with limited cloud free observations. The data fusion method may be a powerful tool to overcome this issue, in combination with land surface modelling-based temperature output (Martins et al 2019, Zhang et al 2019b, Long et al 2020. Additionally, the coarse resolution of the MODIS LST product brings mixed-pixel effect. It will reduce the sensitivity of LST to water surface changes for pixel with partial coverage resulting its relatively poor effect for monitoring small lake or lake border area with small variation in current study. Meanwhile, the coarse resolution brings big uncertainty in the change area statistics. The change area is tremendously larger than that conducted with 30-m Landsat images (Mao et al 2018, Liu et al 2019, Qiao et al 2019. Therefore, the application of the proposed method with higher spatial resolution LST product should be hopeful to solve this issue and the spatial downscaling method will be a good choice in further research (Hutengs and Vohland 2016, Bai et al 2019, Long et al 2019). The use of the downscaled LST data with similar spatial resolution as the popularly used optical remote sensing data will assist the direct evaluation of the proposed method in the accuracy of changing area detection.
The second issue is attributed to the uncertainty induced by the seasonal variation in lake area. Firstly, the detection was conducted based on the annual mean temperature values, the results cannot reflect the seasonal variation in lake area which have been identified by previous studies (Lei et al 2017, Zhang et al 2020. For example, the individual lake analysis against Selin Co indicated that the amplitude of the yearly lake area variation may reach more than 20 km 2 (Zhang et al 2020). The seasonal cycle is highly connected with the monsoon climate and glacier melting. On the other hand, the dramatic seasonal variation will influence the annual mean LST values for the lake boundary pixels and bring uncertainties to the trend analysis. Therefore, besides the analysis based on annual mean values, the seasonal values should be helpful to lighten the influences from the seasonal variation.
For the ACY detection, although its result has been partly confirmed by case analysis, ACY will present big variations even in nearby pixels when the lake area changes with big fluctuation but not a single expansion or shrinkage process. Meanwhile, the altitudinal distribution of the expanding lakes and the shrinking lakes at each year shown in figure 10 indicates quite different pattern for lake expansion and shrinkage. There is no obvious variation in average altitude of the expanding area at each year. Comparatively, the average altitude and its standard derivation of the shrinking area at each year present much more fluctuant. Additional, according to equation (6), the pixels with significant trend were identified as expanding or shrinking pixels. However, the changing rate was not yet considered in current analysis. The changing rate is also highly connected with the temporal change information of lake area. Therefore, combined with the detected ACYs and the change rates, a detailed analysis is required for investigating the driving factors of the temporal variation of each lake, which will be useful for simulating and assessing lake area variation in future (Liu et al 2017, Wu et al 2019a).
Finally, as an important attribute of lakes, lake area also provides important information for lake water storage assessment. The use of satellite altimetry (such as Jason, SARAL, Sentinel-3, and incoming SWOT mission) have demonstrates their advantages in deriving water level and storage changes (Munyaneza et al 2009, Maillard et al 2015, Huang et al 2018, Grippa et al 2019. The combined use of lake area detected from Landsat images and altimetry water level for water level investigation and water storage monitoring conducted by Li et al. (2019) provides a good reference for this study to further analyze the temporal variation of lake water level. It should be an important extension for current research.

Conclusions
Lake area monitoring in the TP is a popular topic in recent years due to its strong connection with climate change. Unlike the commonly used methods based on optical remote sensing data, this study firstly presented the monitoring of TP lake variations by using a new detection scheme based on the periodic information of LST time series data.
Compared with previous monitoring results, the findings revealed that both the change area (including expansion and shrinkage) and the abrupt change year of most lakes in the TP can be accurately monitored. The later one is rarely discussed in previous studies. The lake surface change pattern is consistent with previous investigation results that the TP lake mostly experienced an expanding status. Moreover, for each lake, the ACY result can also depict the lake surface change pattern during the period. The abrupt change year in lake area mainly occurs from 2004 to 2012.
Through integrating time series information of LST, the proposed method avoids the data absence due to cloud cover for optical remote sensing-based methods, and it shows good potential for lake area monitoring in the TP.