A Slight Temperature Warming Trend Occurred over Lake Ontario from 2001 to 2018

Satellite-derived lake surface water temperature (LSWT) measurements can be used for monitoring purposes. However, analyses based on the LSWT of Lake Ontario and the surrounding land surface temperature (LST) are scarce in the current literature. First, we provide an evaluation of the commonly used Moderate Resolution Imaging Spectroradiometer (MODIS)-derived LSWT/LST (MOD11A1 and MYD11A1) using in situ measurements near the area of where Lake Ontario, the St. Lawrence River and the Rideau Canal meet. The MODIS datasets agreed well with ground sites measurements from 2015–2017, with an R2 consistently over 0.90. Among the different ground measurement sites, the best results were achieved for Hill Island, with a correlation of 0.99 and centered root mean square difference (RMSD) of 0.73 K for Aqua/MYD nighttime. The validated MODIS datasets were used to analyze the temperature trend over the study area from 2001 to 2018, through a linear regression method with a Mann–Kendall test. A slight warming trend was found, with 95% confidence over the ground sites from 2003 to 2012 for the MYD11A1-Night datasets. The warming trend for the whole region, including both the lake and the land, was about 0.17 K year−1 for the MYD11A1 datasets during 2003–2012, whereas it was about 0.06 K year−1 during 2003–2018. There was also a spatial pattern of warming, but the trend for the lake region was not obviously different from that of the land region. For the monthly trends, the warming trends for September and October from 2013 to 2018 are much more apparent than those of other months.


Introduction
Lakes are always considered as effective and valid sentinels in relation to climate change. They can provide accurate indications for climate change analysis and its regulation [1][2][3][4][5]. However, some of the indicators, such as the lake surface water temperature (LSWT), which could be acquired by various methods, ranging from ground observations to satellite observations and numerical measurements, have not been thoroughly analyzed [6] due to various reasons. For example, sites used in ground networks are not representative of the neighboring locations [7], and there are many lake and weather related uncertainties in estimates from remote sensing. Nevertheless, remote sensing-based observations have gained increasing popularity for analyzing LSWT and have been used as a good proxy indicators for studying the impact of climate change by providing spatial and temporal temperature observations on lake surfaces [6,8].
The 2013 Intergovernmental Panel on Climate Change (IPCC) and some researchers claimed that the world is in a warming hiatus [9][10][11]. Early studies claimed that the global warming hiatus occurred during 1998 to 2012 [11,12], whereas recent studies have claimed the warming rate remained slow into the early 2010s [12]. However, not many details have been revealed through in situ or satellite observations. Fyfe et al. [13] found the recently observed global warming trend is less than that predicted from climate models. The warming hiatus in areas over the world has attracted researchers' interest. Meanwhile, some researchers debated about whether the warming was on a "hiatus" or "slowed down" [14][15][16].
In 2019, Cheng et al. [17] claimed that the observation-based measurements show a rapid warming of the Earth's oceans over the past few decades. How fast is Earth warming? What happened in specific parts of the world and with the various land/water surface types? Is there any possibility that we can use satellite-derived surface datasets themselves for warming trend analysis when the air temperature is not available?
The IPCC has also emphasized that the urgent need for long term satellite-derived land surface temperature (LST) (also known as surface skin temperature) data in global warming studies to improve the limits of near surface (2 m) air temperature observations [18]. The relationships between air temperature and surface temperature have been explored on daily and annual bases [19,20]. The daily minimum and maximum air temperatures could be acquired both from ground observations at meteorological stations and from geostationary satellite data through an empirical multiple-linear regression model [7]. Some researchers [21,22] also concluded that the relationships between surface temperature and air temperature are strongest during summer and autumn. Additionally, some studies showed that there was a high correlation between air temperature and surface temperature [18].
Lake Ontario is one of the Great Lakes of North America [23], and there are many big cities around it. The spatial and seasonal trends of Great Lake LSWT have significant impacts on regional and global water management and climate analysis [24]. Nalley et al. [25] used the monthly temperature datasets and concluded that a warming trend occurred in the southern of Ontario due to winter and summer warming, caused by the increases in dew point, air moisture, humidity and the urban heat island effects. They also found that the northern hemisphere has been experiencing more warming since the 1950s than in the southern hemisphere. The importance of Lake Ontario to North America and Earth is obvious. Mohsin et al. [26] used the sequential Mann-Kendall (SQMK) test to conclude that winter is contributing the most to the increase in annual minimum temperature in the Greater Toronto Area (GTA) on the shore of Lake Ontario. In addition, a gradually increasing temperature has been observed in western Canada during last few decades, and it is projected to increase into the mid-century [27][28][29].
Nevertheless, LSWT validation with MODIS surface temperature products and time series trend analysis for Canada have not been thoroughly studied for recent years, especially for inland water and lakes. The Moderate Resolution Imaging Spectroradiometer (MODIS) LSWT product has been validated for many places on Earth through a number of studies [8,[30][31][32][33]. However, the quality and feasibility of the satellite retrievals in the different study areas along with systematic errors are also important in the LSWT trend analysis, evaporation and water surface energy flux analysis [34][35][36]. How the satellitederived products could support the global trend evaluation for the long-term and how the satellite datasets could be used to fill the spatial and temporal gaps of in situ measurements in the Lake Ontario region are challenges for remote sensing researchers [17].
Therefore, this research intended to investigate the feasibility and potential of the MODIS-derived surface temperature products for warming trend analysis over Lake Ontario and its surroundings in eastern Ontario region. This paper first gives an evaluation of the relationships between the current commonly used MODIS-derived LSWT and the in situ measurements from the area near where Lake Ontario, St. Lawrence River and the Rideau Canal meet during the flooding seasons over 2015-2017. Then, the MODIS LSWT datasets obtained from 2001-2018 were used for LSWT recent trend analysis using the linear regression method with Mann-Kendall (MK) testing. Based on the validation and trend analysis, the MODIS data's compatibility and capability are discussed.

Study Area
Lake Ontario, one of the Great Lakes, which encompasses between 42.94-45.25 • N and 80.46-75.44 • W, including Lake Ontario, most of the large Canadian cities and most of the population in Canada was our study area. The average depth and total basin area of Lake Ontario are 86 m and 82,990 km 2 , respectively [24].
The six validation sites are located in and near lakes around where Lake Ontario, Rideau Canal and St. Lawrence River meet. Rideau Canal is an inland waterway between the Canadian capital of Ottawa and Lake Ontario at Kingston, Ontario, where the St. Lawrence River begins. There are many inland lakes and water bodies within this region. The land cover map of Lake Ontario and the six enlarged and specified validation sites are presented in Figure 1 (2017). The land cover map was obtained from MODIS land cover product from the Application for Extracting and Exploring Analysis Ready Samples (AppEEARS). The enlarged picture details the six ground sites around, surrounding Kingston, around the Lake Ontario outlet, St. Lawrence River, in southeastern Ontario, Canada. The six specified validations images were downloaded from Google Earth. It can be clearly seen in the ground sites images that all the validation sites are all covered with vegetation or water. Although the cover is not homogeneous, it can largely meet the validation criterion.

MODIS Datasets
MODIS, acquiring data in 36 spectral bands (or groups of wavelengths) is on board the Terra (a.m.) and Aqua (p.m.) satellites. Terra passes south to north over the equator in the morning, and Aqua passes the equator in the afternoon [8].
The MODIS datasets are impeded by cloud coverage; therefore, the available MOD11A1/ MYD11A1 data have 1-4 day intervals. The 8-day MOD11A2/MYD11A2 datasets produced from the average values of clear-sky MOD11A1/MYD11A1 datasets were used for the temperature trend analysis. They were downloaded from AppEEARS. The acquisition times for the datasets were different. The daytime passing times for MOD11A1 and MYD11A1 were about 10:30 a.m. and 13:30 p.m. local time, respectively; the nighttime passing times for MOD11A1 (MOD for short) and MYD11A1 (MYD for short) were about 22:30 p.m. and 1:30 a.m. local time, respectively.

Ground Air Temperature Datasets
The ground temperature measurements from six sites, surrounding Kingston, managed by the Queen's University Biological Station (QUBS), include air temperatures acquired at 15 min intervals from March to August in 2015-2017. However, not all sites have complete records due to the impacts of frozen ice in the winter and flowing waters in the spring melting season. Table 1 shows the major characteristics of the ground sites and the periods of each site with ground measurements each year. The accuracy of the ground-measured temperatures is +/−0.1 K at 283-303 K. around the Lake Ontario outlet, St. Lawrence River, in southeastern Ontario, Canada. The six specified validations images were downloaded from Google Earth. It can be clearly seen in the ground sites images that all the validation sites are all covered with vegetation or water. Although the cover is not homogeneous, it can largely meet the validation criterion.

MODIS Datasets
MODIS, acquiring data in 36 spectral bands (or groups of wavelengths) is on board the Terra (a.m.) and Aqua (p.m.) satellites. Terra passes south to north over the equator in the morning, and Aqua passes the equator in the afternoon [8].
The MODIS datasets are impeded by cloud coverage; therefore, the available MOD11A1/MYD11A1 data have 1-4 day intervals. The 8-day MOD11A2/MYD11A2 datasets produced from the average values of clear-sky MOD11A1/MYD11A1 datasets were used for the temperature trend analysis. They were downloaded from AppEEARS. The acquisition times for the datasets were different. The daytime passing times for MOD11A1 and MYD11A1 were about 10:30 a.m. and 13:30 p.m. local time, respectively; the nighttime passing times for MOD11A1 (MOD for short) and MYD11A1 (MYD for short) were about 22:30 p.m. and 1:30 a.m. local time, respectively.

Ground Air Temperature Datasets
The ground temperature measurements from six sites, surrounding Kingston, managed by the Queen's University Biological Station (QUBS), include air temperatures acquired at 15 min intervals from March to August in 2015-2017. However, not all sites have complete records due to the impacts of frozen ice in the winter and flowing waters in the spring melting season. Table 1 shows the major characteristics of the ground sites and the periods of each site with ground measurements each year. The accuracy of the groundmeasured temperatures is +/−0.1 K at 283-303 K.   The land cover was obtained from MCD12Q1-same source as Figure 1. Note that the land cover might not be consistent with Google Earth observations. Figure 2 shows the scheme used in this paper. The evaluation part is important in remote sensing data applications, especially in climate analysis. Therefore, the daily LSWT products were first compared with ground air temperature measurements (Figure 2a). The daily LSWT products' restrictions and feasibility could be obtained in this process. Based on the LSWT products' performance analysis, the LSWT datasets were then used in the  Figure 2 shows the scheme used in this paper. The evaluation part is important in remote sensing data applications, especially in climate analysis. Therefore, the daily LSWT products were first compared with ground air temperature measurements (Figure 2a) The daily LSWT products' restrictions and feasibility could be obtained in this process Based on the LSWT products' performance analysis, the LSWT datasets were then used in the temperature trend analysis (Figure 2b). After that the site's linear regression with an MK test was conducted and the results were analyzed.

Evaluation of MODIS Products with an In Situ Dataset
The evaluation using the flooding season datasets of QUBS air temperature was o Terra/MOD and Aqua/MYD LSWT products. It was accomplished for 2015-2017 ( Figure  2a). The matchups were selected by the coincident and collocation areas between QUBS sites and MODIS pixels. The satellite datasets included Terra/MOD11A1 (Day/Night) and Aqua/MYD11A1 (Day/Night), at a spatial resolution of 1 km. The daytime and nighttime evaluations are accomplished separately to perform restriction and feasibility analysis o the LSWT products easily and clearly.
The total statistics and relationships between QUBS and MODIS were calculated, in cluding the QUBS site's mean, standard deviation (SD), centered root mean square differ ence (RMSD) and correlation. The statistics were all scaled according to the associated datasets. The mean and SD values are for MODIS and QUBS, separately. The RMSD and correlation were calculated between them.
For the LSWT datasets (Tsi, i = 1 to n) and ground air temperature datasets (Tgi, i = 1 to n), the calculation methods are shown in the following equations: (1

Evaluation of MODIS Products with an In Situ Dataset
The evaluation using the flooding season datasets of QUBS air temperature was of Terra/MOD and Aqua/MYD LSWT products. It was accomplished for 2015-2017 ( Figure 2a). The matchups were selected by the coincident and collocation areas between QUBS sites and MODIS pixels. The satellite datasets included Terra/MOD11A1 (Day/Night) and Aqua/MYD11A1 (Day/Night), at a spatial resolution of 1 km. The daytime and nighttime evaluations are accomplished separately to perform restriction and feasibility analysis of the LSWT products easily and clearly.
The total statistics and relationships between QUBS and MODIS were calculated, including the QUBS site's mean, standard deviation (SD), centered root mean square difference (RMSD) and correlation. The statistics were all scaled according to the associated datasets. The mean and SD values are for MODIS and QUBS, separately. The RMSD and correlation were calculated between them.
For the LSWT datasets (T si , i = 1 to n) and ground air temperature datasets (T gi , i = 1 to n), the calculation methods are shown in the following equations: Land 2021, 10, 1315 6 of 16

Monthly Temperature Trends Analysis from the MODIS Product
The commonly used linear regression was performed in the temperature trend analysis of the study area using the following equation: where y denotes the LSWT; x is the time in years; a is the intercept; b is the slope (K year -1 ), i.e., the warming or cooling rate.
In the annual temperature analysis, the 8-day monthly MODIS, including MOD11A2 and MYD11A2 daytime and nighttime datasets, were first averaged. Then, the monthly average temperature trends for the study area were calculated with the linear regression method. One should note the uncertainty created by ice throughout December-April in the high latitude areas. The monthly analysis was intended to be a good complement to the annual analysis.

Mann-Kendall (MK) Test
Long-term trends of satellite-derived LST have been studied through a parametric statistical method, such as linear regression, in many studies [16]. However, the warming or cooling trend beginnings or hiatuses might not pass the significance test within linear regression. The widely used Mann-Kendall (MK) test was used in this study to observe the change trend and its significance [37][38][39]. As opposed to the linear regression method, the assumption of a normal distribution of the satellite-derived datasets is not a concern in the MK method [26]. The statistical s is first defined as follows: where x i and x j are sequential values in a series. Then the trend Z could be defined as: When Z > 0, the trend is ascending; when Z < 0, the trend is declining. When the absolute values of Z are larger than 1.28, 1.64 and 2.32, the respective data fell within 90%, 95% and 99% confidence intervals, respectively.

Evaluation Results of MODIS Products with In Situ Datasets
Air temperatures from six sites were first compared with MODIS datasets from 2015 to 2017 (Figure 3). To simplify, only the daytime (Figure 3a) and nighttime (Figure 3b) comparisons from the MYD11A2 are displayed. It should be noted that only datasets from the flooding season were compared to avoid the mixed-up lake ice and melt pixels. For example, all sites' first year datasets were from March 12 to 23 August 2015, except the missing HI ( Table 1). The time interval in those datasets is about 1-4 days normally, as mentioned before, due to the cloud obscuring the satellite observations. For HI, we had all other datasets from 2016 to 2017. In late March and early April, due to the cold weather, some frozen surfaces had low temperatures below 273.15 K (0 degree). In these time periods, the differences between the MODIS LSWT and ground air temperature are larger than in the warmer periods. The differences between the two datasets are about 3-4 K. There are a lot of scatter points off the y = x diagonal line in Figure 3a  LSWT accuracy was relatively high relative to the ground air temperature. The detailed statistics can be found in Table 2. some frozen surfaces had low temperatures below 273.15 K (0 degree). In these time periods, the differences between the MODIS LSWT and ground air temperature are larger than in the warmer periods. The differences between the two datasets are about 3-4 K. There are a lot of scatter points off the y = x diagonal line in Figure 3a compared to Figure  3b, which indicates that the daytime MODIS LSWT had a lower correlation with ground air temperature than the nighttime MODIS LSWT. However, the overall MODIS LSWT accuracy was relatively high relative to the ground air temperature. The detailed statistics can be found in Table 2.     Since we used the flooding season datasets in the comparison, it is reasonable that the mean temperatures were all above 280 K ( Table 2). As shown from Figure 1, HI sits south of the other five sites; hence, it had the highest temperatures among all six sites during both daytime and nighttime.
As shown in Table 2 and Figure 3, the MODIS LSWTs agree well with in situ air temperature measurements. Aqua/MYD nighttime datasets are highly correlated with ground measurements with a correlation of over 0.90 for 2015-2017. Meanwhile, the Terra/MOD daytime results are better correlated with the ground measurements than the Terra/MOD nighttime results. For the different sites, the best results are from HI with a correlation of 0.98 for Aqua/MYD nighttime. The worst results are those for HI and EL with a correlation of 0.88 for Aqua/MYD daytime. The nighttime LSWTs matched better with the ground air temperatures than the daytime products. Hence, the nighttime LSWTs were used in the following trend analysis of Lake Ontario seasonal temperature. Nevertheless, the differences between Terra/MOD and Aqua/MYD should be paid more attention in future analyses.

Temperature Trends Found by Linear Regression at Six Ground Sites from MODIS Datasets
To examine the temperature trend and its confidence interval, the linear trends of annual average temperature at the six sites during different periods were examined from the 8-day MODIS-derived nighttime LST product. Table 3 Figure 1, the surfaces of these two sites are fulfilled with lake water in the MODIS pixels, which means that the LSWT trend is not obvious, even during the apparent warming period of other sites and the whole area (Table 3). The highest Z value of 2.50 is for the EL site during 2003-2012, which passed the 99% confidence test. This means that the warming trend of 0.25 K year −1 for EL according to the MYD datasets is highly significant. This warming trend was also the largest found among all six sites. For the period of 2003-2010, the data at EL also passed the 95% confidence test with a warming trend of 0.19 K year −1 .
For the other time periods of MOD and MYD, no significant trends could be concluded. For the slope lower than zero (declining), the absolute value of Z is lower than 1.28, which means the data did not pass the MK test.

Spatial Patterns of Annual Temperature Trend Found by Linear Regression from MODIS Datasets
Further analysis was performed on the trend of total average temperature, with the data passing the 95% confidence threshold over the whole study region from MOD11A2 and MYD11A2 during 2001-2018 ( Figure 4). The lake and land were differentiated using the HydroLAKES database (http://wp.geog.mcgill.ca/hydrolab/, accessed on 1 May 2021). MOD and MYD datasets showed similar trends for the total, lake and land regions. Generally, we found warming of 0.  The annual averaged temperature datasets from MOD11A2 and MYD11A2 for the period of 2001-2018 were then used to compute the annual spatial patterns of temperature trend via linear regression. The spatial temperature trends of regression slopes over the study area, with the data passing the 95% confidence threshold, are illustrated in Figure  5. The annual temperature trends from the nighttime MOD11A2 and MYD11A2 are shown for both periods, 2001/2003-2018 (Figures 5a,b) and 2001/2003-2012 (Figures 5A,B). The boxplots for both datasets are also shown for 2001/2003-2018 in Figures 5(a1,b1). Additionally, the annual temperature trends of the six ground sites, corresponding to the spatial distribution, are also listed in Figures 5(a2,b2).
The values that failed the 95% confidence threshold are shown in white in Figure 5. The MYD nighttime datasets for both periods show much more apparent warming trends over the study area than that the MOD nighttime ones. As in the site analysis in Section 4.2, the warming trend for MYD nighttime datasets is much more apparent over the whole study area, including Lake Ontario and its surroundings during 2003-2012. Some areas warmed by 0.25 K year −1 according to MYD nighttime datasets ( Figure 5B).
For the boxplot analysis, data within the 95% confidence threshold are shown in Figures 5(a1,b1). Apparently, no trend exists for the lake's median value in either dataset over 2001/2003-2018. Meanwhile, a slight warming trend of around 0.1 K year −1 is shown for the land's median value.
The temperatures from the six ground sites' annual trends show peaks in 2012 for both MOD (Figure 5(a2)) and MYD ( Figure 5(b2)) datasets. As in the former results, during 2003-2012, the warming trends over the six ground sites were 0.18 K and 0.12 K year −1 for MYD and MOD datasets, respectively. During the whole study period of 2001-2018, a warming trend of 0.06 K year −1 was found by both datasets using the six ground sites.  The annual averaged temperature datasets from MOD11A2 and MYD11A2 for the period of 2001-2018 were then used to compute the annual spatial patterns of temperature trend via linear regression. The spatial temperature trends of regression slopes over the study area, with the data passing the 95% confidence threshold, are illustrated in Figure 5.  Figure 5(a1,b1). Additionally, the annual temperature trends of the six ground sites, corresponding to the spatial distribution, are also listed in Figure 5(a2,b2).
The values that failed the 95% confidence threshold are shown in white in Figure 5. The MYD nighttime datasets for both periods show much more apparent warming trends over the study area than that the MOD nighttime ones. As in the site analysis in Section 4.2, the warming trend for MYD nighttime datasets is much more apparent over the whole study area, including Lake Ontario and its surroundings during 2003-2012. Some areas warmed by 0.25 K year −1 according to MYD nighttime datasets ( Figure 5B).
For the boxplot analysis, data within the 95% confidence threshold are shown in Figure 5(a1,b1). Apparently, no trend exists for the lake's median value in either dataset over 2001/2003-2018. Meanwhile, a slight warming trend of around 0.1 K year −1 is shown for the land's median value.
The temperatures from the six ground sites' annual trends show peaks in 2012 for both MOD (Figure 5(a2)) and MYD ( Figure 5(b2)) datasets. As in the former results, during 2003-2012, the warming trends over the six ground sites were 0.18 K and 0.12 K year −1 for MYD and MOD datasets, respectively. During the whole study period of 2001-2018, a warming trend of 0.06 K year −1 was found by both datasets using the six ground sites.

Spatial Patterns of Monthly Temperature Trend via Linear Regression from MODIS Datasets
For the monthly linear regression trends of the data within the 95% confidence interval, apparent warming trends over MYD datasets in May, August, September and October were selected and shown in Figure 6. The corresponding ground site trends are also shown. May and August, no values were within the 95% confidence threshold for Lake Ontario itself, though a warming trend for the surroundings was found. On the contrary, the month of October showed an apparent warming trend over Lake Ontario, but no significant trend for its surroundings. For the ground sites, an apparent warming trend of 0.56 K year −1 was found for September during 2012−2018, and a declining trend of −0.08 K year −1 was found for September during 2003−2012. The declining trends here are contradictory to the former results of annual trend analysis.  The month of September (Figure 6c) had an apparent warming trend (mostly around 0.25 K year −1 ) over the whole of Lake Ontario and its surroundings. Furthermore, the center of the lake showed the most obvious warming trend in September. For the months of May and August, no values were within the 95% confidence threshold for Lake Ontario itself, though a warming trend for the surroundings was found. On the contrary, the month of October showed an apparent warming trend over Lake Ontario, but no significant trend for its surroundings.
For the ground sites, an apparent warming trend of 0.56 K year −1 was found for September during 2012−2018, and a declining trend of −0.08 K year −1 was found for September during 2003−2012. The declining trends here are contradictory to the former results of annual trend analysis.

The Relationship between LSWT and Air Temperature
The relationship between LST and air temperature have been widely discussed. The relationships between LSWT and air temperature in this study are shown in Figure 3 and Table 2. The air temperature measurements of EL and RL are of a lake, while the other four sites are situated on land. Therefore, both the relationships between LST/LSWT and air temperature are shown. However, there is not apparent difference between land and water surface temperatures in regard to air temperature. The RMSDs and correlations of EL and RL are not apparently better than those of the other four sites. Similarly to other studies of the relationships between LST and air temperature, the correlations of LSWT and air temperature are higher for nighttime temperatures than daytime temperatures.

The Similar Results from Other Lakes
Many similar analyses of LSWT, air temperature and their trends have been performed for both large and small lakes [40][41][42][43][44]. Some researchers focus on the Tibetan Plateau LSWT variations. Studies found that a global set of inland water bodies showed significant warming in seasonal nighttime LSWT between 1985 and 2009 [43].
Recently, many researchers conclude that the lakes in Tibetan Plateau partly show a warming trend [40,41], while others conclude a cooling trend [42]. The cooling trend in some way confirmed the global warming hiatus occurred during 1998 and 2012 [11,12].

The Possibilities for and Challenges of the MODIS Products in Lake Temperature Trend Analysis
In this part, pros and cons are summarized to illustrate the possibilities for and challenges of the MODIS-derived LSWT products in large lake temperature trend analysis: Pro 1: Satellite observed thermal imagery has the advantage of providing a timesynchronized dense grid of temperature data at local, regional or even global scales [20,45]. MODIS-derived LSWT products could give the continuous spatial and temporal temperature record over the large lake, which the traditional field observations could not achieve, especially for remote and large lake areas. Moreover, MODIS and some other satellite products could be obtained from online ways relatively easier and cheaper than the acquisition of in situ datasets.
Pro 2: MODIS LSWT is relatively easily obtained, and the derived LSWT product has high accuracy for the lake surface according to the validation results in Section 3.1.
Pro 3: MODIS can not only provide LSWT products, but also some surface parameters which are of great importance for large lake analysis, such as ice phenology.
Con 1: MODIS has a short time series-less than 20 years. For climate change analysis, the time period is often over 30 years [9]. However, some researchers also claimed that the length of 17 years' interval was required for identifying global warming trends [12]. This conclusion might give much more confidence in the satellite-derived temperature trend analysis.
Con 2: MODIS can only obtain records four times in a day. The clouds hinder the temporal acquisition greatly too. Even more, for small lakes, the LSWT is not highly correlated with the air temperature. However, some researchers have been making attempts at reconstructing daytime LST under cloud-covered conditions using MODIS land products and MSG geostationary satellite data [46]. Con 3: MODIS is unable to obtain a temperature depth profile of lake water. The above are the major pros and cons for the MODIS datasets. Obviously, MODIS could give accurate spatial and temporal LSWT datasets over the large lake's area. Although the duration may affect the trend analysis, the large number of spatial and temporal MODIS datasets could help the researchers to make up for the shortcomings when digging through the large amounts of information in satellite datasets. The combination of satellite observations, numerical simulations and ground measurements could allow much more accurate and consistent LSWT trend analysis of large lakes [16,22].

Conclusions
LSWT is of great importance in climate change studies, especially for large lakes. To evaluate the feasibility of MODIS LSWT products in warming trend analysis from satellite observations, this research first used in situ measurements collected by QUBS to compare them with satellite-derived daily LSWT. The results show that the satellite-derived LSWT correlated well with in situ air temperatures, especially at nighttime. We can clearly see that the air temperature is much closer to satellite LSWTs at nighttime due to the stable nighttime air conditions. The MODIS LSWTs still match well with in situ air temperatures of the daytime for the large lake's surface and dense vegetation surface in this research. More detailed analysis, including comparing freezing times using satellite data and air temperatures, should be carried out before the warming trend analysis in the future.
The 8-day MODIS datasets were used in the temperature trend analysis over the duration of nearly twenty years of 2001-2018. For the six ground sites, detailed temperature analysis was accomplished using the MK test method. A slight warming trend with 95% confidence was concluded over the ground sites from 2003-2012 for MYD11A1 night datasets. The warming trend for the whole region, including both lake and land, was about 0. Overall, our analysis indicates that the MODIS LSWT datasets have a powerful potential for large lake spatial-temporal temperature trend analysis. In this research, we only used the MODIS temperature products for the trend analysis. The combination and integration of various temperature products from past archives and upcoming satellites could overcome the shortcomings of single sensors in climate change analysis and should be the future trend. The satellite datasets could give relative high accuracy of temperature change over the observed time period. However, the systematic and other inner differences between different sensors have uncertainties. MODIS-based LSWT products cannot totally replace in situ observations; however, they could overcome some limitations of field measurements, such as human error and environmental contaminations.