Reconstructing Six Decades of Surface Temperatures at a Shallow Lake

: Lake surface water temperature (LSWT) plays a fundamental role in the lake energy budget. However, direct observations of LSWT require considerable e ﬀ ort for acquisition and hence are rare relative to a large number of lakes. In lakes where LSWT has not been covered su ﬃ ciently by in situ measurements, remote sensing and lake modeling can be used to produce a ﬁne spatio-temporal record of LSWTs. In our study, the Moderate-Resolution Imaging Spectroradiometer (MODIS) LSWT was used to compare with in situ data at the overpass times over the six sites in Lake Chaohu, a large shallow lake in China. MODIS-derived LSWT reﬂected the variation of lake surface temperature well, with a correlation coe ﬃ cient of 0.96 and a cool bias of 1.25 ◦ C. The bias was modiﬁed by an “Upper Envelop” smoothing method and then employed to evaluate the general lake model (GLM) performance, a one-dimensional hydrodynamic model. The GLM simulations showed good performance compared with MODIS LSWT data at an interannual time scale. A 57-year record of simulated LSWT was hindcast by the well-calibrated GLM for Lake Chaohu. The results showed that LSWT decreased by 0.08 ◦ C / year from 1960 to 1981 and then increased by 0.05 ◦ C / year. These trends were most likely caused by a cooling e ﬀ ect of decreased surface incident solar radiation and a warming e ﬀ ect of reduced wind speed. Our study promoted the use of MODIS-derived LSWT as an alternative data source, and then combined with a numerical model for inland water surface temperature, and also further provided an understanding of climate warming e ﬀ ect on such a shallow eutrophic lake. Key points: (1) Moderate-Resolution Imaging Spectroradiometer (MODIS) lake water surface temperature (LSWT) was validated with real-time in situ data collected at Lake Chaohu with high accuracy; (2) MODIS LSWT was modiﬁed by the bias correction and employed to evaluate a one-dimensional lake model at interannual and intraannual scale; The LSWT hindcast by a well-calibrated model at Lake Chaohu decreased by 0.08 ◦ C / year from 1960 to 1981 and increased by 0.05 ◦ C / year from 1982 to 2016.


Introduction
Lake surface water temperature (LSWT) is a fundamental driver of lake ecosystem structure and its energy budget [1]. LSWT affects the stability and the timing of stratification [2][3][4][5], the rate of lake evaporation [6], the gas exchange [7], mixing processes [8], and duration of ice coverage in some temperate regions [9]. Responding to climate changes, it has been reported that summer LSWTs have increased by 0.34 • C per decade [10,11], even faster than the rate of local air temperature [12][13][14].
In this study, we highlighted a direct overpass comparison between skin temperature measurements obtained from satellites and bulk temperatures of water column layers from observations. This comparison would enable the relationship between skin temperature and bulk temperature. We further proposed a novel strategy for adjusting remotely sensed skin temperature to the bulk temperature of surface layer, then used the data for calibration and validation in a 1-D hydrodynamic lake model. The scientific merit of this work is to provide an efficient bias correction technique for extrapolating remotely sensed skin temperature observation and produce an alternative data source for the lake model. Another scientific merit is that the well-calibrated lake model can forecast lake surface temperature with reasonable accuracy using key meteorological variables, such as wind speed and solar radiation. This is important because it could give a further understanding of how climate change affected a shallow and large lake, especially in an eutrophic condition. In our study, we selected a large shallow lake, Lake Chaohu, located at mid-latitudes in the east of China. The high-frequency in situ LSWT observation, MODIS temperature products, and a general lake model (GLM) were implemented to produce comparable long-term LSWT time series. We addressed the objectives in following steps: 1) evaluate the reliability of MODIS derived LSWT by high-frequency in situ water temperature observations; 2) assess the performance of the GLM simulation by in situ LSWT and MODIS LSWT at different time scales; 3) reproduce LSWT over a period of 57 years and quantify climatic variability of LSWT and its response to a changing climate.

Study Area
Lake Chaohu (31 • 43 -31 • 72 N, 117 • 29 -117 • 85 E) is a polymictic lake in Anhui Province in the lower reach of the Yangtze River ( Figure 1). A large portion of the inflow is from Nanfei River, Hangbu River and Yuxi River and one outflow connecting with Yangtze River. It is the fifth largest freshwater lake in China. The mean depth of the lake is 3 m (maximum depth is 7 m), and the surface area is about 780 km 2 . The nearest weather station is named Hefei (31 • 53 N, 117 • 15 E). The lake basin experiences a northern subtropical monsoon climate with an average annual rainfall of 1100 mm. The annual average air temperature is 15-16 • C.

In Situ Lake Temperature Data
In situ temperatures in Lake Chaohu were measured at intervals of 15 min with loggers in 6 locations distributed over the lake at a depth of 0.3 m below the surface ( Figure 1, HOBO Onset, TidbiT; accuracy ±0.2 • C) [50]. The sites were distributed at the lake center, the bay center and bay mouth ( Figure 1). These in situ temperatures were acquired from 1 January 2016 to 31 December 2016 with a data gap from June to August due to malfunctioning. The daily average temperature was calculated at each site and then averaged over all the sites for the entire lake. Also, monthly observations were obtained from 2008 to 2016 from Kong, et al. [51], at a depth of 0.2~0.5 m near the lake center usually at 10:00.

Pre-Processing
Daily LSWTs were obtained from MODIS land surface temperature products (MOD11A1, MYD11A1, version 6), available on National Aeronautics and Space Administration (NASA) Earth Observation System Data and Information System [52]. The spatial resolution was 1 km and the split-window algorithm was used to retrieve surface temperature from band 31 and band 32 in products and cloud-contaminated temperatures have been removed from these products. The "land" surface temperature products estimated by satellite instruments include both land and inland water [29]. Because of diurnal variation of lake surface temperature [53], LSWTs were obtained from four overpasses per day: Terra overpasses at local times around 10:30 and 22:30, and Aqua overpasses at local times around 1:30 and 13:30. Overpass times could vary slightly and could be extracted from MODIS data files. MOD11A1 data (Terra) were available since 2001 and MYD11A1 data (Aqua) since 2002. In this study, we used data from MOD11A1 and MYD11A1 from 2002 to 2016.
The temperature map for the whole lake from MODIS temperature products was extracted by a lake shapefile mask from high-resolution data of Landsat Thematic Mappez (TM) image at 30 m with a 3-pixels buffer zone along the shore to avoid errors due to fluctuations of land-water interfaces or the variance of year-to-year lake boundary [23]. The temperature was valid only if the ratio of acceptable quality pixels to water pixels was larger than 50% in the whole mask. The pixels were used through the strict quality control according to quality control flag and quality assessment information provided by MOD11A1 and MYD11A1. The reliable quality pixels were accepted flagged with "good quality", "average LST error <= 1K", "average emissivity error <= 0.01", and "average emissivity error <= 0.02". The pixels with the other flags were rejected.

Comparison with In Situ Data
In our study, we compared MODIS-derived LSWT and in situ observation at the site scale and over the whole lake, respectively. The real-time comparisons were evaluated between MODIS-derived LSWT and in situ water surface temperature at six sites across Lake Chaohu. The site J was not included because it was located too close to the lake boundary. The 3 × 3 pixels matrix centered on sample sites was extracted from MODIS LSWT images. The area of this pixel matrix was large enough to represent the temperature estimated from MODIS images [26]. The exact overpass time information of MODIS LSWT was extracted and then the corresponding time of observed LSWT was selected from 15-min high-frequency in situ data. The mean temperature of the matrix was calculated. We further computed MODIS LSWT values over the whole area of the lake and then averaged all four records per day, which was identified as the daily temperature of LSWT from MODIS. The outliers of MODIS-derived LSWT were examined and rejected when the values deviated from climatological temperature by more than 3 times the variation. The climatological temperature was calculated by multi-year average temperature from MODIS.

Post-Processing
Although the cloud-contaminated pixels have been removed by the cloud mask [29], abnormally low values of LSWT were observed owing to the undetected cloud pixels which were much colder than lake surface temperature [29]. Also, systematic biases were found in MODIS-derived LSWT, such as cool bias from skin-bulk temperature difference, instrument noise, drift, sun glint, and split window coefficients for inland water [19,31,54,55]. To reduce the difference above and obtain a reliable daily time series LSWT, we adapted a novel method named "Upper Envelop" smoothing to minimize the difference. The upper envelop smoothing method adopted a simple linear interpolation method to primarily fill data gaps, caused by bad-quality pixels, and then adopted the three-point upper envelop smoothing method by Equation (1), according to Gu,et al. [56], to more effectively obtain the "best estimated" LWST from MODIS products. It was implemented to the solve negative noise reduction problem [56,57].
where LST u−e (t) represented MODIS temperature from Upper envelop method, LST(t), LST(t − 1) and LST(t + 1) represented the data on the t th day, the previous time-step and the following time-step.

Quantitative Evaluation
The quality of LSWT from different data sources was quantified by an evaluation matrix including bias, standard error (STD), and root mean square error (RMSE). Bias was the average difference between MODIS-derived LSWT model simulation and observation, which was expected to evaluate under-or over-estimation. STD was the standard deviation of the difference, which was expected to evaluate the heterogeneity or variance of two datasets. RMSE quantified the overall accuracy of estimated LSWTs, as it combined the notions of bias and standard error. The LSWT and the related atmospheric variables were also analyzed by linear trend with robust analysis, and Pearson correlation coefficient.
The change-point test of LSWT time series was detected by two methods which have been widely used in many studies [58][59][60][61]: the sequential t-test of regime shift [62] and a non-parametric and robust against outliers method, called CUSUM (the cumulative sum) [63], combined with the 1000-times bootstrap analysis to determine the confidence level.

Introduction of the General Lake Model
The general lake model (GLM, [64]), was a one-dimensional hydrodynamic model intrinsically assuming horizontal homogeneity. The model adopted a vertical Lagrangian layer structure and each layer had a unique density computed from the simulated salinity and temperature. The layer thickness was flexible, i.e., layers follow the vertical water movement to minimize numerical diffusion. Layers could split when they grow too voluminous or combine when sufficient energy became available to overcome density stratification between adjacent layers. GLM was particularly suited for long-term investigation, and required only little site-specific calibration [41]. The GLM model allowed users to simulate vertical profiles of temperature, salinity, and density. Mixing and surface layer dynamics were modeled at the confluence of adjacent layers and were dependent on a turbulent kinetic energy budget and potential energy required for mixing. Recent studies showed that GLM has been tested for water temperature in globally distributed lakes of 104 km 2 to 579,000 km 2 surface area [41] and more than 10,000 lakes in the United States [65].

Input Data
The input meteorological variables to force the GLM model (v2.4.0) were wind speed, air temperature, relative humidity, precipitation, shortwave radiation, and longwave radiation. We chose Water 2020, 12, 405 6 of 23 measurements of wind speed, air temperature, relative humidity, precipitation and sunshine duration from the weather station named Hefei, shown in Figure 1, provided from the China Meteorological Administration dataset during the study period of 1960-2016. The dataset was applied to the strict quality control and of good quality [66]. The climatological information of six variables could also be found in Figure 2.

Evaluation of MODIS Lake Surface Water Temperature (LSWT) with 15-Minute In Situ Data
MODIS LSWT was compared with in situ data at six sites across Lake Chaohu in 2016. The exact overpass time information of MODIS LSWT was extracted to make a comparison with high frequency observed LSWT in 2016. The comparison of MODIS LSWT and in situ data showed a good agreement (Figure 3) with the average correlation coefficient of 0.96 at six sites, and with the bias of 1.25 °C, and the STD of 1.45 °C, and the RMSE of 1.98 °C, respectively (Table 1). The incoming shortwave solar radiation (referred to as shortwave radiation in the following context) was calculated by observed sunshine duration according to a method by Yang, et al. [67], which was subsequently improved by Wang [68]. It has been proved that shortwave radiation derived from sunshine duration reflects the effects of clouds and aerosols on shortwave radiation [68,69], and performed well at a regional and global scale [70][71][72]. Calculated shortwave radiation showed good accuracy when compared with shortwave radiation observation. Due to longer availability, the derived shortwave solar radiation revealed a dimming before the 1980s and a brightening afterward [69,70,73]. The sunshine duration derived shortwave radiation agreed well with satellite retrievals, reanalysis, and climate model [69]. Furthermore, the sunshine duration derived shortwave radiation corrected the inhomogeneity of observed shortwave radiation caused by the sensitive drift and instrument replacement and could reflect variability from diurnal to decadal time scales [69].
The incoming longwave solar radiation (referred to as longwave radiation) was estimated following the method from Brutsaert [74] to determine the clear-sky radiation and then corrected to all-sky Water 2020, 12, 405 7 of 23 radiation by applying the cloud-sky equation from Crawford and Duchon [75]. This combination was adopted by the lake model internally and commonly used to estimate the longwave radiation [76][77][78][79].

Model Simulation
The model simulation required lake basic characteristics, input meteorological data, initial conditions, and a series of parameters. The model started on the first day of simulation period, such as 1 January 2016, and then ran forced by input meteorological data, producing daily water surface temperature. In our simulation, the heat effects of inflows and outflows were neglected, as inflow and outflow temperature differences were less than 1 • C according to previous observations (see Table 1 in Zhang, et al. [80] and Huang, et al. [81]), and residence times (~180 days [50]) were too long for a visible effect on water temperature. Also, the previous study showed the inflows and outflows in Lake Chaohu had a limited effect on the thermodynamics and currents in the lake, while energy exchange with atmosphere on lake surface had a more evident impact [82]. Hence, the effect of in-and outflows on lake surface temperature was not considered in our simulation. Table 1. Statistic information about the comparison between Moderate-Resolution Imaging Spectroradiometer (MODIS)-derived lake surface water temperature (LSWT) and observed LSWT from six sites in 2016. The pairs were count for the comparison, and the bias, standard error (STD), R square, and root mean square error (RMSE) were calculated. For most parameter settings, we used default values in our simulation, except for the light extinction coefficient and wind scaling factor. The Secchi depth (SD) was commonly around 0.3 m in spring and 0.48 in late autumn [83,84]. Hence a light extinction coefficient (named "K w " in GLM) in the range of 3-6 m −1 was appropriate, based on the relationship of K w and Secchi Depth from Liu, et al. [85]. The wind speed records were from the weather station at a considerable distance from the lake, and hence a wind scaling factor of 0.6-0.9 was applied based on our knowledge about the uncertainty of parameters and appropriate limits from previous applications of GLM ([65,86]. In our study, the best-fit parameters were calibrated and adjusted to minimize the RMSE based on the comparison of simulated data with the observed data in 2016. The observed data referred to both daily time series of MODIS-derived LSWT and in situ data. The validation and evaluation of the GLM model with best-fit parameters were performed from 2002 to 2016 by the statistical criteria RMSE, bias, STD.  (Table 1).   The correlation was at nighttime than at daytime better between MODIS-derived LSWT and observation (see Figure 3 and Table 1). The STD amounted to 1.16 • C to 1.53 • C at nighttime and 1.35 • C to 2.3 4 • C at daytime. This could be attributed to better spatial homogeneity of lakes surface temperature during nighttime, and the absence of solar heating during nighttime [29]. In addition, the cool bias and STD of LSWT were smallest in Aqua daytime, as the overpass time for Aqua daytime was around 13:30, and lay close to the peaks of the diurnal LSWT cycle.  It seemed that the LSWT from MODIS−UE data underestimated the surface summer during the summertime, especially when the temperature was above 27 °C. During summer, the monthly observation of LSWT was higher by 1.89 °C than MODIS-derived LSWT from 2008 to 2016. From June to September, the peak of shortwave radiation and lower wind speed was found in Lake Chaohu. The surface layers of lake would heat more and the skin effect through cooling would be As shown in Figure 4, a cool bias of MODIS was detected across all six sites when real-time match-ups were compared. The bias varied from −0.50 • C to −0.99 • C at daytime and from −1.22 • C to −2.12 • C at nighttime. The larger bias at nighttime was likely owing to temperature gradients in the shallow lake that the upper layer of water was stratified during the day and mixed well during the night [4]. Due to the high turbidity, Lake Chaohu absorbed most of incoming shortwave radiation used for heating water in the upper layers. During the daytime, the water temperature in the upper layers was quite mixed, but showed a larger gradient between a thin surface and a sub-surface layer. While at nighttime, the mixing was less complete and water temperature was stratified without heating from radiation, and hence the bias was larger between skin and bulk temperature. A similar situation has been documented in Lake Taihu, another shallow and turbid large lake further downstream in the Yangtze catchment [87].

Evaluation of MODIS
Based on the observation in six sample sites, the highest average temperature was 17.34 • C and the lowest value was 16.92 • C. Also, the bias varied from −0.50 • C to −0.99 • C at daytime and from −1.22 • C to −2.12 • C at nighttime between observation and MODIS data among six sample sites. Hence, the heterogeneous temperature was small on the surface layer. Then we retained 280 matching pairs to validate the MODIS-derived LSWT with in situ observations over the whole lake area in 2016. The results showed a significant correlation with the R 2 value of 0.98 ( Figure 5). There was also a noticeable cool bias of −1.76 • C when daily time series of MODIS-derived LSWT were compared with in situ measurements. curve of MODIS showed a good performance generally over the whole period, whereas an obvious cool bias remained in the summertime.

The GLM Model Evaluation
LSWT simulations compared well with MODIS-UE LSWT in 2016 (Figure 7). Bias and STD were 1.8 °C and 1.92 °C, and the correlation coefficient was 0.98. The variations were similar between MODIS-UE LSWT and GLM simulation. Then, LSWT simulation was also compared with MODIS-UE LSWT when excluding the data in the summertime, with the bias and STD value of 1.09 °C and 1.53 °C, respectively (as shown in Figure 7). The large bias of LSWT simulation was possibly contributed to the data in the summertime, as a cool bias was found clearly in the summertime. As the observation was missing during the summer of 2016, it was difficult to validate whether MODIS-UE LSWTs or simulations were corrected during the summer, especially the cool bias might have been larger than in other seasons due to the peak value of solar radiation and low wind speeds [54,88]. Hence, the monthly observations from 2008 to 2016 and the MODIS LSWT from 2002 to 2016 were also employed to validate the GLM simulations performance. When MODIS-derived LSWT was processed by the upper envelop method, and the cool bias reduced to 0.98 • C, as shown in Figure 5. The three-point "upper envelope" method worked well as all single measurements were either on the envelope or below it. The results revealed most of the temperature dynamics seen in in situ measurements qualitatively and quantitatively, provided enough support for the method.
The monthly LSWT observations from 2008 to 2016 were also used to validate the performance of MODIS and MODIS-UE (Upper Envelop line) LSWT, especially to bridge the data gap during summer (from June to August) due to malfunctioning of the high-frequency loggers. Monthly observations were obtained from Kong, et al. [51], sampled monthly at a depth of 0.2-0.5 m near the lake center usually at 10:00. As shown from Figure 6, the seasonal cycle of LSWT fitted well between MODIS derive LSWT and observation, with the correlation coefficient of 0.98 and the bias of −0.35 • C. The bias was −2.05 • C in summer, higher than other spring (−1.09 • C), autumn (−0.77 • C) and winter (−0.55 • C).  To validate the performance of GLM simulations during summer, the monthly observed LSWT from 2008 to 2016 were also used to assess the model results ( Figure 10). The seasonal pattern of simulations was consistent with the observed LSWT. The statistical result showed that the model could reflect the variation of temperatures well during the summer, with a bias of -0.03 °C and an STD of 1.82 °C ( Table 2). The performance of LSWT simulation looked more convincing in spring and winter than in autumn. The RMSE value was low (1.73 °C) in spring. The  It seemed that the LSWT from MODIS−UE data underestimated the surface summer during the summertime, especially when the temperature was above 27 • C. During summer, the monthly observation of LSWT was higher by 1.89 • C than MODIS-derived LSWT from 2008 to 2016. From June to September, the peak of shortwave radiation and lower wind speed was found in Lake Chaohu. The surface layers of lake would heat more and the skin effect through cooling would be more pronounced. In addition, the wind speed was low, and the skin layer could persist as it was not destroyed by wind mixing. Generally, a time series of upper envelope over the MODIS data (referred to as MODIS-UE LSWT), filled gaps properly and reduced uncertainty. The upper envelop curve of MODIS showed a good performance generally over the whole period, whereas an obvious cool bias remained in the summertime.

The GLM Model Evaluation
LSWT simulations compared well with MODIS-UE LSWT in 2016 (Figure 7). Bias and STD were 1.8 • C and 1.92 • C, and the correlation coefficient was 0.98. The variations were similar between MODIS-UE LSWT and GLM simulation. Then, LSWT simulation was also compared with MODIS-UE LSWT when excluding the data in the summertime, with the bias and STD value of 1.09 • C and 1.53 • C, respectively (as shown in Figure 7). The large bias of LSWT simulation was possibly contributed to the data in the summertime, as a cool bias was found clearly in the summertime. As the observation was missing during the summer of 2016, it was difficult to validate whether MODIS-UE LSWTs or simulations were corrected during the summer, especially the cool bias might have been larger than in other seasons due to the peak value of solar radiation and low wind speeds [54,88]. Hence, the monthly observations from 2008 to 2016 and the MODIS LSWT from 2002 to 2016 were also employed to validate the GLM simulations performance.
Water 2019, 11, x FOR PEER REVIEW 13 of 24 data and 0.34 °C compared to observation. In general, these comparisons showed that the GLM model could provide LSWT reliably and accurately. Thus, GLM could be used for long-term simulations to detect changes in LSWT.         To validate the performance of GLM simulations during summer, the monthly observed LSWT from 2008 to 2016 were also used to assess the model results ( Figure 10). The seasonal pattern of simulations was consistent with the observed LSWT. The statistical result showed that the model could reflect the variation of temperatures well during the summer, with a bias of -0.03 • C and an STD of 1.82 • C ( Table 2). The performance of LSWT simulation looked more convincing in spring and winter than in autumn. The RMSE value was low (1.73 • C) in spring. Table 2. Statistic information about the seasonal comparison between general lake model (GLM) simulation and observed LSWT from 2008 to 2016. The bias, standard error (STD), and root mean square error (RMSE) were calculated. The annual anomaly values of MODIS-UE LSWT data and GLM simulations from 2002 to 2016 showed a good agreement (Figure 11), with a correlation coefficient of 0.88, a bias of 0.75 • C, and a standard deviation of 0.13 • C. The annual anomaly observations from 2008 to 2015 were also compared with MODIS-UE LSWT and simulations. The STD value was 0.13 • C compared to MODIS data and 0.34 • C compared to observation. In general, these comparisons showed that the GLM model could provide LSWT reliably and accurately. Thus, GLM could be used for long-term simulations to detect changes in LSWT.

Long-Term Trends of LSWT and Atmospheric Forcing
With calibrated parameters, GLM was forced by meteorological variables from 1960 to 2016 to simulate LSWT in Lake Chaohu. As a result, annual LSWT showed a slightly increasing trend of 0.005 °C/year over the whole period. A change point occurred in 1980 (see Figure 12). There was a clear decreasing trend of 0.08 °C/year before 1980, and a significant increasing trend of about 0.05 °C/year after 1980.
The decadal variations of air temperature and wind speed showed a close relationship with the trend of simulated LSWT. The Pearson coefficient was high, with a value of 0.60 (air temperature), −0.73 (wind speed), respectively. Before 1980, solar radiation dropped sharply by 1.32 W/m 2 /year and the wind speed increased by 0.02 m/s/year (see Figure 2). Air temperature and longwave radiation showed a decreasing trend. These changes could have a cooling effect on the lake surface and would lead to a decreased LSWT in Lake Chaohu. Over the period of 1981-2016, a rapid increase of air temperature (0.046 °C/year), a decrease of wind speed (0.029 m/s/year), a sharp increase of longwave radiation (0.288 W/m 2 /year) and an insignificant decrease of shortwave radiation (−0.203 W/m 2 /year) became visible. The increasing energy from radiation contributed to heating the lake surface and the combination of decreasing wind speed and increasing air temperatures enhanced the stability of lake stratification.
To further quantify the effect of meteorological variables, which were mainly responsible for the simulated LSWT trends, we detrended the forcing variables by multi-annual mean value. Simulations were performed for every single detrended variable while all other input variables were retained. Three forcing variables were selected: air temperature, wind speed, and shortwave radiation. LSWT simulated by detrending shortwave radiation showed a significant positive trend of 0.14 °C/10year ( Figure 13). The LSWT driven by detrending air temperature and wind speed produced a negative trend of −0.12 °C/10year and −0.10 °C/10year, respectively. The simulations for any other detrended input, such as relative humidity, did not show any significant difference to the original simulation. Furthermore, one simulation was performed where wind speed, air temperature, and shortwave radiation were detrended together, which resulted in an insignificant decrease of −0.01 °C/10year of LSWT (Figure 13d).

Long-Term Trends of LSWT and Atmospheric Forcing
With calibrated parameters, GLM was forced by meteorological variables from 1960 to 2016 to simulate LSWT in Lake Chaohu. As a result, annual LSWT showed a slightly increasing trend of 0.005 • C/year over the whole period. A change point occurred in 1980 (see Figure 12). There was a clear decreasing trend of 0.08 • C/year before 1980, and a significant increasing trend of about 0.05 • C/year after 1980.  The decadal variations of air temperature and wind speed showed a close relationship with the trend of simulated LSWT. The Pearson coefficient was high, with a value of 0.60 (air temperature), −0.73 (wind speed), respectively. Before 1980, solar radiation dropped sharply by 1.32 W/m 2 /year and the wind speed increased by 0.02 m/s/year (see Figure 2). Air temperature and longwave radiation showed a decreasing trend. These changes could have a cooling effect on the lake surface and would lead to a decreased LSWT in Lake Chaohu. Over the period of 1981-2016, a rapid increase of air temperature (0.046 • C/year), a decrease of wind speed (0.029 m/s/year), a sharp increase of longwave radiation (0.288 W/m 2 /year) and an insignificant decrease of shortwave radiation (−0.203 W/m 2 /year) became visible. The increasing energy from radiation contributed to heating the lake surface and the combination of decreasing wind speed and increasing air temperatures enhanced the stability of lake stratification.
To further quantify the effect of meteorological variables, which were mainly responsible for the simulated LSWT trends, we detrended the forcing variables by multi-annual mean value. Simulations were performed for every single detrended variable while all other input variables were retained. Three forcing variables were selected: air temperature, wind speed, and shortwave radiation. LSWT simulated by detrending shortwave radiation showed a significant positive trend of 0.14 • C/10year ( Figure 13). The LSWT driven by detrending air temperature and wind speed produced a negative trend of −0.12 • C/10year and −0.10 • C/10year, respectively. The simulations for any other detrended input, such as relative humidity, did not show any significant difference to the original simulation. Furthermore, one simulation was performed where wind speed, air temperature, and shortwave radiation were detrended together, which resulted in an insignificant decrease of −0.01 • C/10year of LSWT (Figure 13d).

Discussion
Compared with our ground truthing of in situ LSWT, the MODIS data revealed a noticeable cool bias at the six sample sites, which also had appeared in many studies retrieving LSWT from MODIS products. Values ranged from 0.1 °C to several Celsius degrees, such as Lake Taihu (−0.63 °C~−1.53 °C) [87], Lake Urmia (−0.27 °C) [54], the Great Salt Lake (−1.5 °C) [88], Lake Tahoe (−0.1 °C) Figure 13. The plot of variation of simulation LSWT compared with detrending analysis from 1960 to 2016: the simulated LSWT forced by non-detrend meteorological input (black dot line) and forced by detrending air temperature (a), shortwave radiation (b), wind speed (c), and forced by detrending jointly variables (air temperature, wind speed, shortwave radiation) (d).

Discussion
Compared with our ground truthing of in situ LSWT, the MODIS data revealed a noticeable cool bias at the six sample sites, which also had appeared in many studies retrieving LSWT from MODIS products.
Values ranged from 0.1 • C to several Celsius degrees, such as Lake Taihu (−0.63 • C~−1.53 • C) [87], Lake Urmia (−0.27 • C) [54], the Great Salt Lake (−1.5 • C) [88], Lake Tahoe (−0.1 • C) [35], Lake Vänern and Lake Vättern (~−0.4 • C) [89]. One of the possible factors for the cool bias was skin-bulk temperature difference. Skin temperature was theoretically cooler than the bulk skin [13], because warming from solar irradiation penetrated deeper (in the order of a meter) into the water column and the heat, e.g., from evaporation, was instantly only taken from a very thin layer at the surface. Our data show a smaller cool skin effect during daytime than nighttime (Table 1). We argued the bulk layer became warmer than the skin layer as solar heating penetrated instantaneously deeper than surface cooling. Apparently, skin-water temperature effect on MODIS was stronger in summer, possibly because of the peak value of shortwave radiation and lower wind speed in summer [54,88]. It is possible that poor atmosphere and undetected cloud pixels could also contribute because cloud top temperature would be mistaken as LSWT value, but it was much colder than LSWT [29].
The cool bias from satellite temperatures had been analyzed over the oceans in previous studies, which were also caused by the heat effect [36,[90][91][92]. They were generally based on physical models and empirical parameterizations [36]. However, the characteristics of lakes were more heterogeneous than oceans when converting skin temperature to bulk temperature. Bias values varied with different conditions, such as the lake size, depth, and environmental conditions. To correct differences of bulk and skin temperature in lakes, the previous literature commonly adopted an estimated or an averaged bias value [87,88], or employed the smoothing filter (Loess or Gaussian filter [26,93], or Empirical orthogonal Function techniques [25]. However, the general smoothing methods were easily affected by abrupt drop values. In our study, the cool bias was corrected by an Upper Envelop smoothing method. This method was firstly used in LAI (leaf area index) research [56] for cool bias correction caused by poor atmospheric conditions (e.g., thick clouds) and undetected cloud pixels. Later it was accepted as an effective approach for reducing the effect of negatively biased noise [57,94]. The upper envelope method worked well as all single measurements were either on the envelope or below it. Hardly any data point was found significantly above the upper envelope. A comparison with in situ measurements revealed most of the temperature dynamics seen in in situ measurements qualitatively and quantitatively, as only some of the data points were affected by the cold bias and the others provided enough support for the envelope. Hence an upper envelope over the MODIS data provided a good solution to fill gaps and reduce uncertainty.
Calibrated by the upper envelop of MODIS-derived LSWT, the 1-D model GLM simulated LSWTs by well-fit parameters and exhibited a good agreement when validated by MODIS LSWT and observations at the monthly and annual scale. Especially, the simulations also showed a good performance in summer compared with monthly observations when MODIS under-estimated LSWT. The GLM model could reflect the variability of LSWT in Lake Chaohu accurately and could hindcast the long-term trend of LSWT with the long record of meteorological input data: the simulated LSWT initially decreased at a rate of 0.8 • C/10year from 1960 to 1980, and then increased 0.5 • C/10year after 1980 in Lake Chaohu. The rate of increasing LSWT in the latter phase was also consistent with but higher than the global lake warming trend (0.45 • C/10 year) from the observed and satellite records [10,12,95]. This strong increase of LSWT after 1980 had previously been found in shallow lakes in other regions, such as in northeastern United States [96,97], and northern Europa [61,98].
We detected a breakpoint in the LSWT simulation of Lake Chaohu in 1980, which has also been found in European lakes [61,99]. This shift corresponded to a climate regime shift in the 1980s [59], under the background of rapid warming [100], brightening and dimming [101], a decline in wind speed [102] at a global scale. In Lake Chaohu, air temperature increased by 0.46 • C/decade since 1980, at the rate twice of the global average (0.25 • C/decade) [100], and the wind speed decreased by 0.29m/s/decade, consistent with the atmospheric still phenomenon around the world [102]. Solar radiation decreased, under the background of a decrease in China (1.06 W/m 2 /decade) [70], contrary to the global brightening, especially in Europa [101]. Previous studies had discussed that higher air temperatures significantly contributed to recent warming of lakes [103,104]. In addition, the change in solar radiation [17,98], and wind speed [105,106] affected temperature trends in lakes. Our model results suggested that wind speed, air temperature, and solar radiation had an associated effect on the shift and the long-term trend of LSWT in Lake Chaohu, which helped us to understand the importance of meteorological drivers.
It had been reported that Lake Chaohu had suffered from eutrophication and serious cyanobacteria bloom after the 1980s [107,108], coinciding with the breakpoint of LSWT in 1980s. Apparently, the increased temperature played an important role [38,109,110]. The higher temperature (above 25 • C/) generally had led to increased cyanobacteria growth rates and reduced vertical mixing.

Conclusions
In this study, remote-sensing data (MODIS temperature products) and a well-calibrated lake model named GLM were jointly implemented to investigate the variability of LSWT in Lake Chaohu, a large and shallow lake in China. For ground truthing, in situ data from six locations in the lake were used to assess MODIS LSWT performance by real-time comparisons. According to the comparison, MODIS-derived temperature reflected the variation of lake surface temperature well, with a correlation coefficient of 0.96. However, we also detected a cool bias of 1.25 • C/ when validated by in situ records. As a consequence, MODIS-derived LSWTs were modified by a three-point "Upper Envelop" smoothing method, and then the time series of MODIS-derived LSWT with the bias correction was used as a data source to evaluate the lake model performance. A well-calibrated GLM was used to hindcast a long-term record of LSWT back to 1960 using meteorological data from the weather station. From 1960 to 2016, LSWTs decreased by 0.08 • C/year from 1960 to 1980 and then increased by 0.05 • C/year since 1980, mainly due to a combined effect of rapidly increased air temperature (0.05 • C/year), and decreased wind speed (0.03 m/s/year), and lower shortwave radiation (−0.22 W/m 2 /year).
In conclusion, one-year records of detailed LSWT measurements were sufficient to calibrate the lake model GLM with meteorological data from a nearby weather station. The simulation of LSWT could be extended back as far as the meteorological data were available. MODIS LSWT data was suited for validating the model back to 2002. In summer in particular, deviations could be attributed to MODIS cold bias. Numerical simulations could bridge data gaps in satellite observations and not be affected by atmospheric conditions which influenced both LWST and data collection in the case of satellite measurements. A similar approach may be implemented at other lakes of regional importance. Further work is needed to focus on exploring the links between a warming LSWT and algal concentration in such large and eutrophic lakes under climate change.