Temporal and spatial heterogeneity of recent lake surface water temperature trends in the Qinghai-Tibet Plateau

Abstract Lake surface water temperature (LSWT) monitoring is significant as it provides valuable information on climate changes and resultant changes in lakes. Here, monthly moderate resolution imaging spectroradiometer LSWT data of 364 lakes during 2001–2015 were analysed to determine spatiotemporal patterns of daytime (LSWTd), nighttime (LSWTn), daily LSWT (LSWTm) and their drivers in the Qinghai-Tibet Plateau at monthly, seasonal and annual timescales. We focused on spatiotemporal heterogeneity of LSWT trends. Results showed climatological LSWT presented a saddle-shaped distribution along the northeast-southwest direction except for LSWTd in March, while highly heterogeneous LSWT trends without spatial dependence were observed at all timescales. Diurnal asymmetry of LSWT trends was evident with 56.04–88.19% lakes showing smaller LSWTd trends than LSWTn trends except for August and September. Further, we found lake surface albedo trend could effectively explain spatial patterns of LSWT trends at different timescales, and its interactions with other variable trends generally were the largest or second-largest.


Introduction
Qinghai-Tibet Plateau (QTP), known as 'Third Pole of World' after the Arctic and Antarctic and 'Asian Water Tower', has become a focus area for climate change research, especially for the temperature change (Cai et al. 2017;Tian et al. 2020;Yao et al. 2021). During the last several decades, the change characteristics of near-surface atmosphere (T) and land surface temperature (LST) at different temporal and spatial scales as well as their drivers in the QTP have been intensively investigated, and an abundant of meaningful information has been obtained (Wang et al. 2002;Duan and Wu 2006;Yang et al. 2014;Fang et al. 2019;Tian et al. 2019;Xue et al. 2021). Lake is a fundamental component of the landscape and plays an essential role in land-air interaction (Williamson et al. 2008;Adrian et al. 2009;Williamson et al. 2009). As an important indicator of lake state and functioning, lake surface water temperature (LSWT) can provide detailed information on climate-induced changes in lakes (Adrian et al. 2009;Schneider et al. 2009;Schneider and Hook 2010) and can better characterize the long-term climate changes compared to T (Arvola et al. 2009;Torbick et al. 2016;Ye et al. 2019). QTP has the largest number of plateau lakes at the highest altitude worldwide because of the unique geological structure and topography (Qi et al. 2019). There were 1424 lakes larger than 1 km 2 in 2018, with a total area of 5.0 Â 10 4 ± 791.4 km 2 . However, the recorded impacts of climatic change on LSWT were relatively scarce in the QTP.
In situ monitoring is the most direct and effective way for ascertaining LSWT, while it is expensive and time-consuming, and generally has a high probability of undersampling (Alcantara et al. 2010). Consequently, in situ measurements were only available at suitable places and times such as the developed areas (Xiao et al. 2013;O'Reilly et al. 2015;Yang et al. 2019). This hampered the understandings of LSWT patterns at large temporal and spatial scales (Schwab et al. 1999). In contrast, the infrared images taken by remote sensing have relatively fine temporal resolution and large spatial coverage, which provide an unprecedented opportunity for dynamic LSWT monitoring around the world, especially in remote regions with harsh environmental conditions such as the QTP (Schwab et al. 1999). There were five commonly used satellite sensors for LSWT measurements, i.e. Terra/Aqua-moderate resolution imaging spectroradiometer (MODIS), Terra advanced spaceborne thermal emission and reflection radiometer (ASTER), NOAA advanced very high resolution radiometer (AVHRR), along-track scanning radiometer and Landsat series (Steissberg et al. 2005;Hodges et al. 2016;Virdis et al. 2020). Among them, MODIS images were better applicable to large-scale LSWT investigations due to their high temporal resolution and comparatively low spatial resolution (Yang et al. 2019).
Based on in situ and satellite-derived LSWT data, recent studies have documented that major lakes worldwide have experienced a warming LSWT trend with a larger rate in summer than that of local T in some regions in recent decades (Austin and Colman 2007;Hampton et al. 2008;Schneider et al. 2009;Schneider and Hook 2010;Fink et al. 2014;O'Reilly et al. 2015;Woolway et al. 2021). The drivers for LSWT variations are complex and various in both space and time. In addition to lake characteristics such as lake bathymetry, area and salinity, climatic factors also play a significant role in LSWT variations. For instance, Yang et al. (2019) used MODIS LST product to show that T was one of the determinant factors of lake warming in Yunnan-Guizhou Plateau during 2001-2017. In terms of the QTP, several studies have also reported the warming LSWT. Xiao et al. (2013) noted that Qinghai Lake was warming at a rate of 0.01 C/a during 2001-2010 by analysing data from the MODIS. Ke and Song (2014) took the Siling Co as a study area and found the lake-averaged LSWT has increased at a rate of 0.26 C/10a during 2001-2013 based on MODIS LST images, with a warmer rate of 0.34 C/10a at night than that at the daytime (0.18 C/10a). Zhang et al. (2014) investigated the trends in nighttime LSWT (LSWT n ) of 52 lakes using MODIS LST data and concluded that both warming and cooling LSWT n in the QTP can be the result of regional warming.
The above research provided an overview of LSWT changes in the QTP and proved the reliability of MODIS LST data in LSWT investigations. However, most studies primarily focused on long-term trends in summer and annual LSWT n at a specific lake or a limited region and comparative analysis between trends in LSWT and T. The diurnal, seasonal and monthly variations in daytime (LSWT d ), nighttime and daily LSWT (LSWT m ) were poorly understood. The effects of climate change on LSWT might be not equal among months and seasons, as well as among spaces. Additionally, little works have been devoted to the interpretation of spatial patterns of LSWT trends and their relationships to lake characteristics and climatic factors, such as lake surface albedo (A) regulating directly surface energy budget (Wang and Davidson 2007). What's more, previous studies preferred to use regression analysis to explore the relationships between LSWT and individual factor. These ignored the influences of interactions between factors on LSWT. In response to these concerns, the present study aims to present a geospatial analysis for changing LSWT in the QTP at different timescales (diurnal, monthly, seasonal and annual) and their drivers, focusing on potential spatiotemporal heterogeneity of LSWT trends. The outputs are expected to deepen understandings of climate changes and climatic responses of lakes in high mountainous areas.

Study area
QTP is the highest plateau around the world with an average elevation of over 4000 m a.s.l. It has a vast territory and generally can be divided into twelve river basins ( Figure 1) (Wan et al. 2017). QTP is chiefly characterized by plateau mountain climate, with a cool, rainy and humid summer, and a cold, snowy and dry winter. On average, from the southeast to the north and northwest, QTP becomes gradually higher, colder and drier . Until the remote Inner QTP, the average elevation reaches $5000 m a.s.l. and the resulting winter T can drop to 40 C below zero . Due to huge topography and unique geological structure, tens of thousands of glaciers and nearly half of the total number and area of lakes in China have developed in the QTP, with a concentrated distribution in the central and western part (Zhang et al. 2014). The lakes here are mainly fed on the meltwater of snow and glacier . Affected by a cold climate, most lakes are ice-covered periodically for the months generally from late November or December to the following April or May (Yao et al. 2012;Kropacek et al. 2013;Qiao and Zhu 2019).

Data
Monthly LSWT d (10:30 am) and LSWT n (10:30 pm), as well as geolocation data of lakes such longitude (Lon), latitude (Lat) and altitude (Alt) were taken from the LSWT dataset published publicly in international journal 'Scientific Data' (Wan et al. 2017) (Table 1). The dataset was generated on the basis of MOD11A2 version 6 LST product at resolution of 1 km with strict quality control and verification. It is the first comprehensive dataset including 15-year monthly LSWT d and LSWT n of 374 lakes (!10 km 2 ) in the QTP (Wan et al. 2017) and has been used successfully in subsequent LSWT studies (Wan et al. 2018). Among 374 lakes, 6 lakes fall outside the QTP delineated by Zhang et al. (2002) and 4 lakes have a smaller area than one LST pixel size after subtracting the area of 0.5 km buffer area offshore of lake from the lake area. Consequently, 364 lakes here were selected for data collection and analysis. Monthly LSWT m for a given lake then was obtained by averaging arithmetically the corresponding monthly LSWT d and LSWT n . The area-weighted mean monthly, seasonal and annual LSWT were prepared for the climatological analysis during the study period. According to Wan et al. (2017Wan et al. ( ), 2002Wan et al. ( , 2005Wan et al. ( , 2009 and 2014 lake areas were used, respectively, in the calculation of the area-weighted means during 2001-2003, 2004-2007, 2008-2011 and 2012-2015. For the characterization of climate conditions around the studied lakes, it is difficult to use the observations at nearby national meteorological stations in the QTP due to spatial mismatch between lakes and meteorological stations. Accordingly, we used China meteorological forcing dataset (CMFD) at 0.1 resolution instead ( Table 1). The dataset provides monthly T, precipitation (P), wind speed (WS), downward longwave radiation (DLR), downward shortwave radiation (DSR) and specific humidity (SH) (He et al. 2020). It was produced by merging multi-source data, including routine observations, satellite retrievals and reanalysis data (He et al. 2020). CMFD showed a better precision than the existing grid datasets (He et al. 2015;Long et al. 2018;) and thereby has been widely used in studies related to LSWT in the QTP Guo et al. 2018;Zhang et al. 2018;. Although CMFD has an inevitable error compared to the observations, it is sufficient for the present study whose main concerns are the trends in LSWT and climatic variables. In the analysis, SH was converted into relative humidity (RH) using standard equation (Wallace and Hobbs 2006).  (Table 1), which was produced daily based on 16 days of MODIS Terra and Aqua observations at 0.05 grid. The dataset provides both white-sky albedo and black-sky albedo of shortwave at local solar noon. MCD43C3 version 6 albedo has been considered reasonable compared to in-situ observations and widely used in climate and hydrology (Argaman et al. 2012;Lang et al. 2018). Given that black sky albedo of shortwave has similar value with white sky albedo (Zhu et al. 2011), here we only used black-sky albedo for the analysis. Additionally, the elevation data at 1 km resolution in the QTP were also used in this study to map Figure 1, which were sourced from the National Tibetan Plateau Data Center (Table 1).
Both CMFD and albedo data were resampled to a 1 km pixel size using bilinear interpolation in ArcGIS 10.2 to match with MOD11A2 version 6 LST product. Considering that LSWT data were provided monthly at the whole-lake scale, monthly lake-averaged meteorological data and MODIS albedo for each lake were extracted by lake boundary in ArcGIS 10.2. In data extraction, 2002, 2005, 2009 and 2014 lake boundaries were used for the periods of 2001-2003, 2004-2007, 2008-2011 and 2012-2015, respectively. The data processing was performed as the procedure of LSWT extraction described in Wan et al. (2017). Finally, the arithmetic sum (for precipitation) and mean (for the other variables) of parameters were calculated to obtain the corresponding seasonal and annual values.

Global Moran's index
Global Moran's index (GMI) is a commonly used statistic index in spatial autocorrelation analysis (Ren et al. 2016). It can be shown as follows: where n is the number of the studied lakes; X is the arithmetic mean of LSWT; X i and X j are LSWT values of the ith and jth lakes, respectively; W ij is the weight of the ith lake relative to the jth lake. The range of GMI is from À1 to 1. GMI >0 (<0) indicates a positive (negative) spatial autocorrelation and an obvious spatial cluster (spatial difference). If GMI is close to 1, this suggests greater spatial agglomeration; when GMI nears 0, a random spatial pattern is presented . In this study, GMI s of LSWT at different timescales were calculated using a spatial statistic tool in ArcGIS 10.2.

Geographical detector
Geographical Detector (GD, available freely at http://www.geodetector.cn/) is a novel spatial statistics method, which is developed to detect spatial heterogeneity in variables and their drivers (Wang et al. 2010). Unlike the conventional correlation and regression analysis methods, it is good for analysing the categorical data and identifying the interaction relationship between two independent variables without an assumption of linearity association (Wang et al. 2010). Therefore, GD has been widely used in various fields Hu et al. 2020;Liu et al. 2020;Wang et al. 2021). GD consists of four modules, i.e. risk detector, factor detector, ecological detector and interaction detector (Wang et al. 2010). The primary ideas for factor detector are spatial stratified heterogeneity and that if an independent variable X (i.e. climatic variable trend) impacts the dependent variable Y (i.e. LSWT trend), then Y would display a similar spatial or temporal pattern to X.
The higher degree of similarity between the two, the higher level of sensitivity of Y to X. The impact of X on Y is measured by Power of Determinant (PD) as follows： where h ¼ 1, … , L are the categories of X or Y; N h and N are the numbers of lakes in category h and the whole study area, respectively, and r h 2 and r 2 are the variances of Y, correspondingly. PD value ranges from 0 to 1, indicating that X explains PD Â 100% of the variations in Y. Moreover, the larger PD value represents a higher level of spatial heterogeneity of Y and the stronger influence of X on Y. When PD value is equal to 1, Y distribution has a perfect spatial heterogeneity and can be fully determined by X.
Interaction detector mainly determines whether two X variables (e.g. X 1 , X 2 ) have an interactive influence on Y and their interactive types, i.e. weakening, independent and enhanced relationships. Interaction results can be obtained by comparing PD(X 1 \X 2 ) with PD(X 1 ) and PD(X 2 ) as shown in Table 2. Before running GD, all variables were reclassified into eight categories using natural breaks classification, this ensures that the influences of climatic variable trends on LSWT trends were analysed and compared under the same category frame.
Additionally, the linear regression model with the least-squares method was employed to estimate monthly, seasonal and annual LSWT trends during 2001-2015 (Zhang et al. 2014;Wan et al. 2018;Wang et al. 2021); the Pearson correlation coefficients were calculated to qualitatively characterize the relationships between trends in LSWT and climatic variables (Shen et al. 2015). The significances of trends and correlation coefficients were evaluated by Student's t test.

Spatial patterns of LSWT at different timescales
Mean monthly, seasonal and annual LSWT in the QTP during 2001-2015 presented obvious spatial variations. Because of the similar climatological distributions of LSWT in different months within a season and the space limitation, March, June, September and December were selected, respectively, to represent the periods of March to May (i.e. Spring), June to August (i.e. Summer), September to November (i.e. Autumn), and December to February (i.e. Winter) in Figure 2. During the past 15 years, mean annual LSWT d , LSWT n, and LSWT m ranged from À0.94 to 25.98 C, from À12.20 to 7.28 C and from À4.84 to 12.90 C, respectively, the mean values were 6.19, À3.55 and 1.32 C, correspondingly. Except for the mosaic distribution pattern of high and low LSWT d values in March, mean monthly and annual LSWT d , LSWT n and LSWT m featured an approximate saddle-shaped pattern along the northeast-southwest direction. LSWT for most lakes in the central QTP (i.e. Inner E and F) was relatively low all the year-round compared to the northeast with high latitude and the southwest with a long ice-free duration (Wang Table 2. Interactive types between factors identified by geographical detector.

Interactive types Description
Enhance PD(X 1 \X 2 ) > PD(X 1 ) or PD(X 2 ), PD(X 1 \X 2 ) > PD(X 1 ) and PD(X 2 ) PD(X 1 \X 2 ) > PD(X 1 ) þ PD(X 2 ) Weaken PD(X 1 \X 2 ) < PD(X 1 ) or PD(X 2 ), PD(X 1 \X 2 ) < PD(X 1 ) and PD(X 2 ) PD(X 1 \X 2 ) < PD(X 1 ) þ PD(X 2 ) Independent PD(X 1 \X 2 ) ¼ PD(X 1 ) þ PD(X 2 ) et al. 2020). The finding accorded with ground investigations ). It was not only the direct result of locally high elevation but also was closely related to the relatively high A in the whole year: lakes in Inner E and F usually experienced a longer ice-covered duration than the other areas, the resulting mean annual A therein was generally higher (Lin et al. 2020). Further, the GMI s of mean monthly, seasonal and annual LSWT in the QTP were calculated, respectively (Table 3). Results showed that the GMI s at all timescales were positive (p < 0.01), indicating climatological LSWT in study area displayed obvious spatial LSWT n ¼ 73:17 À 1:20Lat À 0:08Lon À 0:01Alt R 2 ¼ 0:74 (4) LSWT m ¼ 74:01 À 1:10Lat À 0:11Lon À 0:01Alt R 2 ¼ 0:88 All regression coefficients were significantly negative (p < 0.01) and the determinable coefficients of three models were greater than 0.48, which meant that high LSWT values tended to distribute in the low latitude, longitude and altitude regions. These further showed that there were evident spatial dependence and spatial aggregation feature of mean annual LSWT in study area during 2001-2015. In comparison, Lat had maximum effects on mean annual LSWT distributions, indicative of dominant role of thermal conditions. The findings here were generally consistent with those obtained by Du et al. (2020) in the Eurasian continent.

LSWT trends at different timescales
Interannual trends in monthly, seasonal and annual LSWT presented obvious spatial heterogeneity. Linear regression analysis was performed to get the interannual LSWT trend for each lake (Figure 3). At monthly scale, over 68% of lakes displayed an increasing LSWT d during August to November, while an intensive decreasing LSWT d was detected during January to May. Especially in January to April, more than 87% of lakes showed a cooling trend and greater 29% of lakes were significant at p < 0.1. In contrast, although a considerable number of lakes presented a cooling LSWT n during the period when monthly LSWT d decreased widely, LSWT n primarily showed a warming trend in all months. The most intensive warming LSWT n occurred in October, which accounted for 82.7% of the studied lakes with 34.6% of lakes being significant at p < 0.1. The occurrences of interannual trends in monthly LSWT m were quite consistent with those of monthly LSWT d . Partly this might be because interannual fluctuations of LSWT d were larger than those of LSWT n . The frequency distributions of seasonal LSWT trends were consistent with those of monthly LSWT trends, with the prevailing cooling LSWT d and LSWT m in Spring and Winter, dominant warming LSWT d and LSWT m in Autumn and LSWT n in four seasons. The major differences between the two time scales were observed in spatial patterns of LSWT d and LSWT m trends. For instance, the number of lakes with cooling LSWT d in winter was noticeably smaller than that in December but higher than that in February, this discrepancy can be attributed to the arithmetical average method in which the changes in monthly LSWT were smoothed. At annual scale, LSWT changes were characterized by a general cooling LSWT d (71.2% lakes) and warming LSWT n (68.4% lakes). However, the number of lakes for warming LSWT m was comparable to that for cooling LSWT m .
In parallel, the GMI s of trends in LSWT d , LSWT n and LSWT m were calculated at different timescales (Table 3). On average, the GMI s of LSWT m trends were maximum (GMI 2 [0.09, 0.25]) relative to those of LSWT d and LSWT n at all timescales, and they for LSWT n trends were minimum with a range from 0.03 to 0.13. Although positive GMI s registered at a significant level of p < 0.01 at all timescales, most of them were comparatively small as a maximum value of 0.25, indicating that spatial aggregation tendencies in LSWT trends in the QTP at different timescales were substantially weak compared to those in LSWT. Use spatial distributions in annual LSWT trends as examples, the multiple linear regression analyses between trends in annual LSWT d , LSWT n , LSWT m and the geographical positions of lakes were carried out. The determinable coefficients of obtained regression models for trends in LSWT d , LSWT n and LSWT m were, respectively, only 0.13, 0.07 and 0.04, suggesting that annual LSWT trends in the QTP showed significant spatial heterogeneity but almost no spatial dependence. LSWT trends in lakes with very close distances might be also different markedly. To better understand spatial patterns in LSWT trends at different timescales, the following analyses will mainly focus on determinants of spatiotemporally heterogenous LSWT trends.

Relationships between the trends in LSWT and climatic variables
For a specific lake, LSWT variation is mainly affected by complex heat exchange processes at air-water interface (Alcantara et al. 2010). Lake water absorbs and reflects solar radiation (shortwave) and atmospheric radiation (longwave), and radiates heat to the atmosphere through evaporation, convection and outgoing radiation. When net radiation flux absorbed by lake water is greater (smaller) than total radiating heat flux emitted from lake water, LSWT will increase (decrease). As a consequence, trends in variables related to water and heat balance should explain the spatial patterns of LSWT trends. Given this, spatial correlations between trends in LSWT and climatic variables including P, T, DLR, DSR, A, RH and WS were examined to elucidate causes of spatial patterns in LSWT trends at diurnal, monthly, seasonal and annual scales, respectively (Figure 4).
The climatic variable trend whose spatial correlation with LSWT trend was the strongest was considered the most dominant factor (MDF) contributing to spatial pattern in LSWT trend in the analysis. It can be seen in Figure 4 that the trends in P, A and RH had negative effects on LSWT trend distributions and those of other factors exerted a positive role. Generally, the MDF s for spatial patterns in LSWT trends were various in different periods. For monthly LSWT d trend, it was T trend from April to May (p < 0.01), RH trend in June (p < 0.01) and A trend in other months (p < 0.01). For monthly LSWT n trend, A trend was identified as MDF at all months (p < 0.01), indicating that lakes with larger A trend tended to have smaller monthly LSWT n trend. For monthly LSWT m trend, it was RH trend in June and September, DLR trend in August and A trend in other months, and all the corresponding spatial correlation coefficients passed the significance test of p ¼ 0.01.
At seasonal scale, the MDF s presented a distinct temporal distribution, especially for spatial patterns in LSWT d trends. Specifically, they were respectively trends in T, A, P and A in Spring, Summer, Autumn and winter in terms of spatial distributions in LSWT d trends, with the corresponding spatial correlation coefficients (p < 0.01) of 0.23, 0.47, À0.19 and À0.39. Nevertheless, the MDF s were A trends for LSWT n trend distributions in four seasons, with a significant spatial correlation coefficient (p < 0.01) of À0.63, À0.31, À0.44 and À0.57 from Spring to Winter, respectively. The dominant roles of A trends were also found in spatial patterns in LSWT m trends with an exception in Summer. At annual scale, the MDF was RH trend for LSWT d trend, and A trend for both trends in LSWT n and LSWT m . A weak spatial correlation between trends in LSWT and solar radiation (i.e. DLR and DSR trends) registered at all timescales, which might be attributed to minor annual changes in solar radiation in study area during 2001-2015 ).

Geographical detector-based identification of major determinants on spatial patterns in LSWT trends
The factor detector was employed to further explore the relationships between spatial patterns in LSWT trends and those of trends in climatic factors. Figure 5 illustrated the power of determinant of trends in climatic variables on trends in LSWT d , LSWT n and LSWT m at different timescales. Similarly, the variable trend that showed the largest PD value was recognized as MDF contributing to spatial patterns in LSWT trends. As depicted in Figure 5, LSWT trends presented significant spatial heterogeneity (PD > 0.1) at different timescales. For monthly LSWT d trend, the spatial patterns primarily reflected the influences of A trends with PD values ranging from 0.14 to 0.41 in most months. For monthly LSWT n trend, the MDFs at all timescales were A trends with PD values varying from 0.10 to 0.42. For monthly LSWT m trend, although the MDF was changeable among months, the spatial patterns mainly depended on A trends in most months.
There were differences in MDF s contributing to spatial patterns of trends in seasonal LSWT d , LSWT n and LSWT m . The MDF s were, respectively, trends in T, A, P and A in Spring, Summer, Autumn and winter for LSWT d , and A trends for LSWT n and LSWT m in four seasons except for LSWT m in Summer. At annual scale, the MDF s leading to spatial distributions of annual LSWT d , LSWT n and LSWT m trends were trends in DLR, A and RH, respectively. Comparison of the MDF s identified by spatial correlation analysis and by factor detector at different timescales showed that they were correspondingly compatible with one another, especially for LSWT n trends. Moreover, it was found the differences between PD values of the MDF s determined by two methods were very small, suggesting that the identified MDF s were reliable and valid at different timescales.
The interactions between factors on spatial patterns in LSWT trends were analysed using the interaction detector. They were positive and significant, and the explanatory powers of some variables for LSWT trend distributions were improved significantly after interacting with another one (Table S1-S3 in Supplementary material). For example, the PD value of P trend increased significantly from 0.06 to 0.48 through its interaction with A trend. It can be concluded that LSWT variations were involved with the synergistic interactions of multiple variables. Generally, the interactions between any factors and A trends were the largest or second-largest, especially for LSWT n trends at monthly scale, suggesting that A trend played a dominant role in spatial distributions in LSWT trends in the QTP. This was in agreement with the conclusion obtained in Section 3.3.1. On average, the explanatory powers of the interactions between climatic variables were largest for spatial variations in LSWT m trends at different timescales. The largest PD value of the interactions registered between RH trend and A trend in July (0.60) for LSWT d trend, between DSR trend and A trend in April (0.58) for LSWT n trend, between RH trend and A trend in February (0.66) for LSWT m trend, indicating that the corresponding explanatory powers for spatial patterns in LSWT trends can reach over 58% via these interactions.

Discussion
In this study, we mainly investigated spatial patterns of LSWT trends and their drivers at different timescales from the perspective of spatiotemporal heterogeneity. Results confirmed an overall warming summer LSWT n in the QTP during 2001-2015 reported by Woolway and Merchant (2017). Compared with previous researches, however, one finding obtained herein was not fully reported, namely the significant spatiotemporal heterogeneity of LSWT trends in the QTP at different timescales. First, there was asymmetry between LSWT d and LSWT n variations. Except for August and September at monthly scale, LSWT d trends in most lakes (56.04-88.19%) were smaller than the corresponding LSWT n trends at all timescales ( Figure 6). These were in line with the asymmetric variations of the maximum and minimum T, which could be partially explained by increased cloud amount (Solomon et al. 2007). Additionally, the increasing A observed by MODIS in the frozen period from 2001 to 2015 and glacier and permafrost meltwater input to lakes can also contribute to cooling LSWT d (Zhang et al. 2014). At night, however, glacier and permafrost meltwater were negligible in the frozen period, LSWT n variations were closely related to energy storage in the daytime (Alcantara et al. 2010) and in the preceding periods. It can be assumed that water heat storage tends to increase under warming climate and lake expansion (Stefan et al. 1998). For the remaining months, especially from June to August, both increased evaporation and recharge by glaciers and permafrost meltwater mitigated LSWT d variations. The above suggested that multiple mechanisms acted together to produce diurnal asymmetry of LSWT trend. Secondly, we found that spatial patterns in LSWT trends were complex and dependent on time. There were no statistical correlations of LSWT trends with the geographical positions of lakes including longitude, latitude and altitude, which were different from Zhang et al. (2014) reporting the elevation dependence of LSWT trends. This discrepancy might be attributed to the different studied periods and lakes, illustrating the necessity of a considerable number of lakes and data for understanding LSWT variations in the QTP. By combining spatial correlation analysis and Geographical Detector, we found A trends at lake surface generally played a decisive role in spatial patterns of LSWT trends in the QTP at different timescales, especially for LSWT n . On average, the smaller A trend, the more rapid warming in LSWT, with an exception for LSWT d trend from June to September and Summer when a positive correlation between them was recorded. The positive correlation contradicts with physical principle, partly because of the inference by other factors (Vukovich et al. 1987). Burakowski et al. (2018) found that A had a negligible influence on LST at annual scale and high A induced cool LST in winter, which further supported our findings.
Previous research has documented that T was the most significant meteorological variable influencing on temporal variation of LSWT (Zhang et al. 2014;O'Reilly et al. 2015;Yang et al. 2019). To test whether it exists in the QTP during the study period, we further compute the Pearson correlation coefficients between LSWT d , LSWT n , LSWT m and climatic variables at 364 lakes at different timescales, respectively ( Figure 7). As seen in Figure 7, LSWT has respectively the most significant positive correlation with T but the most significant negative correlation with A at different timescales, and the absolute values of the correlation coefficients between LSWT and A were close to those between LSWT and T. This suggested that dynamic LSWT was primarily influenced by T and A at different timescales, and T was mainly responsible for warming LSWT. From the above analysis, it can be concluded that A played a dominant role in both temporal variations of LSWT and spatial patterns of LSWT trends. QTP has intensive solar radiation as its low latitude, high elevation and thin atmosphere. Therefore, minor A change can result in a significant effect on solar radiation flux into lake water, then leading to LSWT variation and resultant LSWT trend spatially. Because of a close link between A and lake characteristics such as water quality, the effects of A trends on LSWT trends represent those of lake characteristics to some extent.
Given that frequent heterogeneity in LSWT trends observed in various global regions, such as Scotland (Carvalho and Kirika 2003), eastern Finland (Voutilainen et al. 2014) and Wisconsin in USA (Winslow et al. 2017), it should be cautious about extrapolating observed LSWT trend for one period to other parts of the year and for a specific lake to other lakes. Especially, our study demonstrated that the rate of LSWT warming would be overestimated if only LSWT n data were used in the research. Additionally, in most existing models to simulate weather and climate (Thackeray and Fletcher 2016;Shoko et al. 2017), A has been recognized as an important input parameter. Because the available observations were scarce globally, especially in the QTP, most models either considered A as a specific constant or estimated it by the empirical relationship between it and LSWT. However, the direct correlation between A and LSWT is still a debatable point, sometimes it is rather weak . Therefore, climate variability simulated by the model based on the above schemes has many uncertainties . The stable and significant spatial relationships of LSWT trends with A trends at different timescales found in this study can provide an important constraint on A and LSWT model.

Conclusions
This study unravelled the spatiotemporal patterns of LSWT and their climatic drivers in the QTP at different timescales, focusing on spatiotemporal heterogeneity of LSWT trends. Climatological LSWT generally had obvious spatial clusters at all timescales with an approximate saddle-shaped distribution along the northeast-southwest direction (except for LSWT d in March), which can be explained well by the geographical locations of lakes. However, LSWT trends showed significant spatial heterogeneity but almost no spatial dependence. Affected by combined impacts of multi-factors, diurnal asymmetry of LSWT trends was evident in the QTP with generally smaller LSWT d trends for most lakes (56.04-88.19%) than LSWT n trends at different timescales except for August and September. Especially from winter to early summer, substantial lakes showed cooling LSWT d and warming LSWT n .
Further, it was found that spatial patterns in LSWT trends generally were associated with those of albedo trends rather than T trends and geographical locations of the lakes. What's more, we found that the interactions of A trends with trends in other variables overall were the largest or second-largest, suggesting A trend was mainly responsible for heterogeneous spatial patterns in LSWT trends at different timescales. Therefore, it should be cautious about extrapolating observed LSWT trend for one period to other parts of the year and for a specific lake to other lakes. These findings are expected to provide an insight into spatiotemporal patterns of LSWT in high mountainous areas and support the research of climate change and its effects. Future studies will focus on long-term LSWT monitoring practice and simulation to in deep investigate the underlying mechanism of heterogeneous LSWT trends in the QTP.

Disclosure statement
No potential conflict of interest was reported by the authors.