The Impact of Land Cover Change on Surface Water Temperature of Small Lakes in Eastern Ontario from 1985 to 2020

: Land Cover Change (LCC) has been shown to signiﬁcantly impact the magnitude and trend of Land Surface Temperature (LST). However, the inﬂuence of LCC near waterbodies outside of an urban environment remain less understood. Waterbodies serve as local climate moderators where nearby LCC has the potential to decrease their cooling ability. Altered water surface temperatures can lead to altered species migration and distribution in aquatic species depending on a given species thermal boundary. In this study, using remotely sensed land cover and surface temperature data, we investigate the role that LCC around small lakes (500 m) plays on the surface water temperature change of nine small lakes in the Cataraqui Region Conservation Authority’s watershed, located in Eastern Ontario, from 1985 to 2020. The Continuous Change Detection Classiﬁcation (CCDC) algorithm was used alongside the Statistical Mono-Window (SMW) algorithm to calculate LCC and LST, respectively. Results indicated a strong positive relationship (R2 = 0.81) between overall LCC and lake surface water temperature (LSWT) trends, where LSWT trends in all inland small lakes investigated were found to be positive. The land cover class sparse vegetation had a strong positive correlation with water temperature, whereas dense vegetation displayed a strong negative correlation. This 35-year study contributes to the broader understanding of the impact that LCC has on the surface water temperature trends of inland lakes.


Introduction
The impact of land cover change (LCC) on Land Surface Temperature (LST) has been well studied [1,2]. Land cover transitions are known to result in warmer temperatures when transitioning from vegetated, evaporative, pervious surfaces to less vegetated, impervious surfaces, and surfaces with less heat storage [3]. This is understood to occur on various scales and timespans across the Earth [4]. However, proximity to water bodies [3], parks [5], and wetlands [6] is known to cool nearby surface temperatures. Transitions from one land cover type to another can alter how water moves [7], sensible and latent heat, and thermal infrared radiation [8]. Water bodies have been noted to have a cooling effect by mitigating surface temperatures that are increasing nearby [9]. This phenomenon is generally known as 'cooling islands', with variations on the name depending on the LC type being investigated (urban cooling islands, water cooling islands) [10]. Waterbodies tend to have a high thermal inertia, high heat capacity, low thermal conductivity [11], and to provide evaporative cooling [3]. Consequently, the high heat capacity and evaporation of water allow it to create an 'oasis effect' by decreasing nearby surface and air temperatures [12]. Waterbodies can be seen as local climate regulators, and nearby LCC has the ability, through advection, to modify the rate at which water moderates the temperature of its surrounding neighborhood. However, the impact of LCC on water bodies is much less explored, resulting in a lack The CRCA watershed is characterized by a cool, humid climate, moderated by Lake Ontario. Nine individual small lakes within the CRCA were selected for detailed analysis of LCC and LST using the following criteria: (a) having a distance of a minimum of 500 m from another large body of water, and (b) not interconnecting or being a portion of another waterbody. Since all lakes investigated lie within the same study region, they experience the same climatic effects. To limit other factors that could contribute to surface water temperature trends, none of the selected lakes have incoming water flow from other lakes. Bodies of water exhibit a cooling effect on the land surrounding them, which, in turn, would affect the analysis of water temperature trends when nearby lakes are included in the analysis [3]. To determine the impact of LCC on the surface water temperature, the land cover change of the buffered areas within a distance of 500 m around the lakes was extracted. The buffer distance of 500 m was chosen based on the study by Cao et al. (2010) to investigate cooling islands [10]. Although a larger distance was found by Du et al. [24], their study included much larger waterbodies (1-171 km 2 ) than the ones used in this study (0.37-2.38 km 2 ). The lakes selected are displayed in Figure 1.

Datasets
Google Earth Engine offers Landsat 4-8 TOA and surface reflectance (SR) (Tier 1, Collection 1) without the need to download images to an individual computer [42]. Landsat images used in this study were collected from Tier 1 images taken on their respective sensors. Images were sorted and analyzed from the Landsat satellite series archive, beginning in 1982 [43]. All images were then masked from clouds (< 20%), cloud shadows, and aerosols based on the quality assessment and metadata provided from the images [44]. Individual pixels were also filtered for clouds and cloud shadows based using the Band Quality Assessment (BQA) and pixel quality attributes (pixel_qa) bands for the TOA and SR images, respectively. A total of 695 images were collected and processed from 1982 to 2020. Figure 2 illustrates the number of images from 1982 to 2020 that were used in each year and the time of year they were captured. In total, there were 27 images from Landsat 4 TM, 281 images from Landsat 5 TM, 170 images from Landsat 7 ETM+, and 169 images from Landsat 8. The CRCA watershed is characterized by a cool, humid climate, moderated by Lake Ontario. Nine individual small lakes within the CRCA were selected for detailed analysis of LCC and LST using the following criteria: (a) having a distance of a minimum of 500 m from another large body of water, and (b) not interconnecting or being a portion of another waterbody. Since all lakes investigated lie within the same study region, they experience the same climatic effects. To limit other factors that could contribute to surface water temperature trends, none of the selected lakes have incoming water flow from other lakes. Bodies of water exhibit a cooling effect on the land surrounding them, which, in turn, would affect the analysis of water temperature trends when nearby lakes are included in the analysis [3]. To determine the impact of LCC on the surface water temperature, the land cover change of the buffered areas within a distance of 500 m around the lakes was extracted. The buffer distance of 500 m was chosen based on the study by Cao et al. (2010) to investigate cooling islands [10]. Although a larger distance was found by Du et al. [24], their study included much larger waterbodies (1-171 km 2 ) than the ones used in this study (0.37-2.38 km 2 ). The lakes selected are displayed in Figure 1.

Datasets
Google Earth Engine offers Landsat 4-8 TOA and surface reflectance (SR) (Tier 1, Collection 1) without the need to download images to an individual computer [42]. Landsat images used in this study were collected from Tier 1 images taken on their respective sensors. Images were sorted and analyzed from the Landsat satellite series archive, beginning in 1982 [43]. All images were then masked from clouds (< 20%), cloud shadows, and aerosols based on the quality assessment and metadata provided from the images [44]. Individual pixels were also filtered for clouds and cloud shadows based using the Band Quality Assessment (BQA) and pixel quality attributes (pixel_qa) bands for the TOA and SR images, respectively. A total of 695 images were collected and processed from 1982 to 2020. Figure 2 illustrates the number of images from 1982 to 2020 that were used in each year and the time of year they were captured. In total, there were 27 images from Landsat 4 TM, 281 images from Landsat 5 TM, 170 images from Landsat 7 ETM+, and 169 images from Landsat 8.
LC reference data was used from the Canadian Centre for Remote Sensing (CCRS) which has generated 30-m LC maps of Canada for the years 2010 and 2015 using Landsat images from TM, Enhanced TM (ETM), and Operational Land Imager (OLI) sensors [45]. The 2010 and 2015 CCRS land cover classifications were classified using different nomenclature, which required the classifications to be regrouped into suitable general LC types  (Table 1). Only consistent land cover data between 2010 and 2015 were used. The LC types were broadly grouped, as similar LCs would have a similar environmental effect. The LC types used were 'impervious surface', which represented all developed land; 'dense vegetation', which represented all forested areas; 'sparse vegetation', which represented all grasslands, shrublands and croplands, 'wetlands', and 'water'. LC reference data was used from the Canadian Centre for Remote Sensing (CCRS) which has generated 30-m LC maps of Canada for the years 2010 and 2015 using Landsat images from TM, Enhanced TM (ETM), and Operational Land Imager (OLI) sensors [45]. The 2010 and 2015 CCRS land cover classifications were classified using different nomenclature, which required the classifications to be regrouped into suitable general LC types for this study (Table 1). Only consistent land cover data between 2010 and 2015 were used. The LC types were broadly grouped, as similar LCs would have a similar environmental effect. The LC types used were 'impervious surface', which represented all developed land; 'dense vegetation', which represented all forested areas; 'sparse vegetation', which represented all grasslands, shrublands and croplands, 'wetlands', and 'water'.   Surface emissivity data are an important component of the LST calculation. Surface emissivity data were obtained from the Advanced Spaceborne Thermal Emission and Reflection Radiometer Global Emissivity Database (ASTER GEDv3) dataset. The ASTER GEDv3 dataset includes the emissivity for five emissivity bands in the TIR region from all clear-sky ASTER images from 2000-2008 at a spatial resolution of 100 m [46]. Previous research has found that this dataset has an error of approximately 1 K and an overall RMSE of 1.52 • C when calculating LST [46,47]. Since we are only interested in the emissivity of water, which remains stable throughout the study period, it is not necessary to account for emissivity values outside of the ASTER GEDv3 time span.
Atmospheric water vapor content data are necessary to more accurately account for atmospheric effects in the TIR region. Water vapor content data are available from the National Center for Environmental Predication (NCEP) and National Center for Atmospheric Research (NCAR), which generates environmental reanalyses of data [48]. The Total Column Water Vapor (TCWV) dataset from NCEP/NCAR reanalysis is available on GEE from 1948 to present, with a six-hour temporal resolution and a 2.5 arc degrees spatial resolution [48]. Due to the dataset's temporal resolution, the TCWV data were linearly interpolated to each Landsat image's acquisition time between noon and 6 pm local time [49,50].

Methods
The workflow used in this research is illustrated in Figure 3. This research has been completed via Google Earth Engine (GEE), a cloud-based platform for geospatial analysis that makes use of Google's computational capabilities [42]. Each lake investigated was first set to a buffer of 500 m in GEE using the buffer function, which uses the edge of each lake polygon to draw a distance of 500 m around the lake polygon. For LCC, the continuous change detection classification (CCDC) method was applied for the time series Landsat images. The LCC is calculated for the land between the lake itself and the 500 m lake boundary. The LCC percentage is calculated by dividing the LCC in the buffered region by the total area around the lake. The LCC percentage is the total amount of land within the 500 m boundary that has experienced a change, regardless of which type, within the study period. For a land surface temperature (LST) change, LST is calculated for every image using the Statistical Mono-Window (SMW) algorithm. Then a trend analysis is conducted based on the seasonal harmonic model to capture the gradual change of the waterbody temperature through the study period. The components of LC classification and LST are brought together to analyze the effects of LCC surrounding lakes in the study area.
RMSE of 1.52 °C when calculating LST [46,47]. Since we are only interested in the emissivity of water, which remains stable throughout the study period, it is not necessary to account for emissivity values outside of the ASTER GEDv3 time span.
Atmospheric water vapor content data are necessary to more accurately account for atmospheric effects in the TIR region. Water vapor content data are available from the National Center for Environmental Predication (NCEP) and National Center for Atmospheric Research (NCAR), which generates environmental reanalyses of data [48]. The Total Column Water Vapor (TCWV) dataset from NCEP/NCAR reanalysis is available on GEE from 1948 to present, with a six-hour temporal resolution and a 2.5 arc degrees spatial resolution [48]. Due to the dataset's temporal resolution, the TCWV data were linearly interpolated to each Landsat image's acquisition time between noon and 6 pm local time [49,50].

Methods
The workflow used in this research is illustrated in Figure 3. This research has been completed via Google Earth Engine (GEE), a cloud-based platform for geospatial analysis that makes use of Google's computational capabilities [42]. Each lake investigated was first set to a buffer of 500 m in GEE using the buffer function, which uses the edge of each lake polygon to draw a distance of 500 m around the lake polygon. For LCC, the continuous change detection classification (CCDC) method was applied for the time series Landsat images. The LCC is calculated for the land between the lake itself and the 500 m lake boundary. The LCC percentage is calculated by dividing the LCC in the buffered region by the total area around the lake. The LCC percentage is the total amount of land within the 500 m boundary that has experienced a change, regardless of which type, within the study period. For a land surface temperature (LST) change, LST is calculated for every image using the Statistical Mono-Window (SMW) algorithm. Then a trend analysis is conducted based on the seasonal harmonic model to capture the gradual change of the waterbody temperature through the study period. The components of LC classification and LST are brought together to analyze the effects of LCC surrounding lakes in the study area.

Continuous Change Detection Classification
LCC can be grouped into three distinct categories: intra-annual change, gradual inter-annual change, and abrupt change. Intra-annual change is caused by vegetation

Continuous Change Detection Classification
LCC can be grouped into three distinct categories: intra-annual change, gradual inter-annual change, and abrupt change. Intra-annual change is caused by vegetation phenology driven by seasonal patterns of environmental factors, such as temperature and precipitation, which occur within the same 12-month period or less. Gradual inter-annual change is caused by climate variability, vegetation growth, or gradual change in land management, and occurs over the period of a year or longer. Abrupt change is caused by deforestation, floods, fire, and urbanization, which are all high-impact and relatively quick types of changes. Due to these three types of LCC, it is important to use a time-series model that will account for and detect all types of changes. The Continuous Change Detection and Classification (CCDC) algorithm is a time-series model that has components of seasonality, trend, and segment breaks [51]. The model coefficients are estimated from the Ordinary Least Squares (OLS) method. The CCDC algorithm was applied to surface reflectance, brightness temperature, and NDVI of images that met the selection criteria for LC classification. The algorithm is able to identify LC changes through a time-series model consisting of the seasonal and trend components (Equation (1)).
where x is the Julian date, i is the Landsat Band being investigated, T is the number of days per year (set to 365), a 0,i is the coefficient for the overall value for a Landsat Band, a 1,i and b 1,i are the coefficients for intra-annual change for a Landsat band, c 1,i is the coefficient for inter-annual change for a Landsat band, T k is the number of break points, and p(i,x) OLS is the predicted value for a Landsat Band at Julian date x.
The basis of the CCDC algorithm is to compare model predictions with cloud-free satellite observations and to normalize their differences by three times the root mean square error (RMSE) (Equation (2)). LCC was identified if the normalized difference exceeded the predefined thresholds determined by Equation (2) below for three consecutive observations.
K is the dimension of the datasets used in the algorithm. The difference of three times the RMSE to detect change was decided upon, because LCC typically occurs when the spectral signature deviates by that amount from the modelled predicted value [31]. The model normalization period (initialization phase) was dependent on the number of clear-sky observations that were used in the model-fitting procedure. Due to the fact that LCC may occur during the first and last observations, or during the process of the model initialization for the first 12 observations, three conditions were set. The implementation of these conditions were used to identify irregular observations for the first and last observations of the model initialization and the observations between them (Equations (3)- (5)). This is an important step because when the model initialization is completed, it is used as the basis for the CCDC algorithm.
Condition for first observation: Condition for last observation: Condition for observation between the first and last observation: where x 1 is the Julian date for the first observation during model initialization, x n is the Julian date for the last observation during model initialization, and t model is the total time used for model initialization.
From here, the model is capable of classifying the given study area at any point in time after the initialization period. A Random Forest (RF) classifier is used with inputs of the coefficients derived from the time-series model. The main idea of the classification is that Land 2023, 12, 547 7 of 18 different LC will demonstrate different modeling characteristics for reflectance, NDVI, and brightness temperatures. The RF classifier [52] was used to train the reference data, from which the trained classifier was applied to the entire time series on a pixel-by-pixel basis to achieve LC maps at the beginning and end of the study period for each lake's surrounding land cover.
The classification accuracy for all lakes investigated was assessed within the time period of the reference data. Since stratified sampling can achieve better precision for small classes and each lake, a stratified random sample was performed with 50 samples per class to ensure each class had enough samples [53,54]. In total, 250 references points were selected around each lake, where possible.

Lake Surface Water Temperature (LSWT) Extraction
Lake Surface Water Temperature (LSWT) is calculated using the Statistical Mono-Window (SMW) algorithm [55]. The Statistical Mono-Window (SMW) algorithm is based on the relationship between the TOA brightness temperature of a TIR channel and LST and makes use of linear regression [9,27,56]. The SMW method corrects the TOA brightness temperature in a single TIR channel for atmospheric effects and surface emissivity [14]. The SMW algorithm is commonly used to calculate LST for sensors with a single TIR channel and is used for consistency across all sensors [55,57]. The SMW equation is given by Duguay-Tetzlaff et al. (2015) [55], as shown in Equation (6): where T b is the TOA brightness temperature in the TIR channel for a given pixel, ε is the surface emissivity in the same channel for a given pixel, and A, B, and C are algorithm coefficients. Following the procedure of Ermida et al. (2020) [57], the algorithm coefficients are determined from linear regressions of radiative transfer simulations performed for 10 classes of Total Column Water Vapor (TCWV), ranging from 0 cm to 6 cm of water vapor in increments of 0.6 cm. If values of TCWV are greater than 6 cm, they are assigned to the last class, as these values are rarely seen [58]. The look-up table for the TCWV values used was created by Ermida et al. (2020) [57], who generated the algorithm coefficients for the SMW algorithm. The look-up table was generated from a calibration database that was derived using a dataset of air temperature, ozone profiles, and water vapor compiled by Borbas et al. (2005) [59]. Further details on the process are provided in Ermida et al. (2020) [57].
This LWST retrieval method was validated with an overall accuracy of 0.5 K, 0.1 K, and 0.2 K for Landsat 5, Landsat 7, and Landsat 8, respectively [57]. The original code provided by Ermida et al. (2020) [57] in JavaScript was reorganized and reconfigured in GEE, as it was originally created for point LST analysis, whereas entire portions of the image were analyzed in this study.

Trend Analysis
From the SMW LST method, the LST was generated for a given pixel in an area of interest over all of the images that meet our image criteria. Subsequently, we are interested in what the LST overall trend from all investigated images can tell us. However, these images only provide us with the LST values at the time of image acquisition, whereas we are interested in approximating the overall trend throughout the study period outside of the image acquisition. Due to the Landsat satellite series' temporal resolution and image contamination, there are gaps in the data where it becomes necessary to interpolate LST values. This accounts for seasonality, as opposed to simply using the individual LST values of each image. In order to accomplish this, a harmonic model was adopted for its ability to manage the unevenly distributed remote sensing images and to characterize periodic patterns throughout the study period. A linear change of the harmonic model was assumed to capture the gradual change of the waterbody temperature throughout the study period [60]. For each individual pixel investigated, a seasonal harmonic model is applied, from which the linear trend is then taken. The seasonal harmonic model is generated from the equation by Shumway & Stoffer (2014) [61], which is shown below: where e t is random error, B 0 is the constant, B 1 t is the linear trend, A is amplitude, ω is the frequency, and ϕ is the phase. This process can be performed in GEE by first adding the harmonic components to each individual image, fitting the model with the linear trend, and then plugging the coefficients into the equation to arrive at a harmonic time series of fitted values. From the harmonic model, the next step is to create a general linear trend in order to capture the overall surface water temperature trend of each individual pixel. The equation for a linear trend is shown below: where β 0 is the intercept, β 1 is the slope, and e t is the random error. This step of determining the overall linear trend was performed in GEE with the use of the ee.Reducer.linearRegression function. The trending rate of LST change in the study period was generated for each pixel.

CCDC Classification Maps and LCC around Lakes
Land cover classified maps were generated from the CCDC method for each year. Figures 4 and 5 illustrate the classified LC maps from the CCDC classification for the nine lakes' buffered areas in 1985 and 2020 which are the beginning and end of the study period, respectively. The land cover percentages of each class in 1985 and 2020 are summarized in Table 2 and their land cover changes as percentages are listed in Table 3. Due to the lack of reference data in 1985 and 2020, we only conducted the accuracy assessment for the LC classified maps generated from CCDC in 2010 and 2015 using stratified sampling. The overall accuracies are reported in Table 3, with the most accurate buffer classification being Moulton Lake (96.61%), the least accurate being Cedar Lake (74.64%), and an average overall classification for all lakes of 84.30%.
From Figures 4 and 5, and Table 2, it is evident that most LC classes around nine lakes are sparse and dense vegetation, and impervious surfaces only cover a very small proportion of the buffered areas of nine lakes. Cedar Lake has the largest cover of impervious surface, at 11.96%, while Moulton Lake has almost no impervious surface within its 500 m surrounding buffer region, indicating no or little developed land around these lakes.
Among the nine lakes, Horseshoe Lake, Cedar Lake, and Moulton Lake have less than 10% of LCC in 35 years (1985-2020). The largest LCC percentage, at 18.02%, occurred in the buffer area of Elodia Lake, followed by Lambs Lake, with 17.26% LCC. The buffered region of both Elodia Lake and Lambs Lake contains a mix of all five land cover classes. Notably, Elodia Lake has wetlands on the eastern shore, a road on the western side, and sparse vegetation scattered throughout. The LCC around Elodia is highly concentrated on the northwestern side, with much less LCC taking place in the southeastern section. From 1985 to 2020, there was a notable decrease in the land cover of water and a similar increase in wetland cover around Elodia. Compared with Elodia Lake, the roads and the wetland to the northeast of Lambs Lake were stable throughout the study period. Contrarily, change can be observed from dense to sparse vegetation in the east of the buffer and an increase of sparse vegetation and wetland in the west portion of the buffer. Although LCCs occurred within the buffered region, they are distributed sparsely in general and LC classes remained relatively stable in most of the buffer regions during the study period.

LCC and Lake Surface Water Temperature (LSWT) Change
The LCC and rate of LSWT change were explored for all nine lakes investigated. For each individual pixel within a given lake, a time series model of the LSWT was generated of which a harmonic model was created. From the harmonic model for a given pixel, the line of best fit was taken which illustrated a given pixel's average trend for the study period. Figure 6 depicts the LCC and the rate of LSWT change for nine lakes. The minimum, maximum, and average rate of changes around each lake are listed in Table 3.
From Figure 6, the rate of LSWT change is not consistent among nine lakes. On average, Lambs Lake has the highest rate of change in WST and Horseshoe Lake has the lowest. The percentage of LCC is 17.26% for Lambs Lake and 5.2% for Horseshoe Lake. Elodia Lake, Lambs Lake, Wiltse Lake and Grippen Lake have higher average rates of change than other lakes, with higher percentages of LCC in their surrounding buffer regions. The two lakes with the highest maximum rates of change are Elodia and Lambs Lake, which also have the highest percentages of LCC, indicating a positive relationship between the LCC and water surface change.      The rate change of LSWT is also not evenly distributed within each lake. The distribution of rate change appears to be impacted by the LCC and its type. For example, the LCC around Elodia is highly concentrated on the north side and much less LCC is taking place in the south section. Notably, the south section of Elodia Lake contains a less positive temperature trend than the rest of the lake (Figure 6b). Similarly, the southeastern part of Grippen Lake with more LCC experiences a stronger positive rate of temperature change than the middle of the lake with the strongest maximum positive rates of temperature change with a rate of 0.139 × 10 −4 • Kelvin per year. The change rate may also be impacted by lake size, shape, lake depth, watershed flow direction and terrain in the buffer area. For example, the water surface near the shoreline experiences a higher change rate of LSWT than water near the lake center. This is more evident in Grippen Lake, Horseshoe Lake and Little Lake (Figure 6C,D,G). The largest lake, South Lake, contains an island in the middle of it (Figure 6h). The eastern portion of South Lake displays a relatively moderate positive trend while the western portion of the lake features diverging trends. The western edge of South Lake contains the lowest temperature trend increase, with small, fragmented areas of LCC nearby, while the shore north of the island contains the greatest positive temperature trends within the lake. There are larger extents of LCC in this north area which could help to explain this phenomenon. Figure 7 lists the scatter plots of the percentage of LCC and the minimum, maximum, and average rates of water temperature change of nine small lakes. The relationship between LCC and average rate of water surface change fits strongly positively with a R 2 value of 0.80 (Figure 7c). This relationship demonstrates that as the amount of LCC increases, as does the rate of water temperature warming for the lakes investigated. In Figure 6, there is a cluster of lakes around 14-18% LCC and a water temperature rate of 0.1 × 10 −4 -0.14 × 10 −4 • Kelvin per day. Although a water temperature rate of 0.1 × 10 −4 -0.14 × 10 −4 • Kelvin per day is relatively minimal, it translates to 0.0365 and 0.0511 • Kelvin per decade, which presents substantial long-term change. From Figure 6, the rate of LSWT change is not consistent among nine lakes. On average, Lambs Lake has the highest rate of change in WST and Horseshoe Lake has the lowest. The percentage of LCC is 17.26% for Lambs Lake and 5.2% for Horseshoe Lake. Elodia Lake, Lambs Lake, Wiltse Lake and Grippen Lake have higher average rates of change than other lakes, with higher percentages of LCC in their surrounding buffer regions. The two lakes with the highest maximum rates of change are Elodia and Lambs Lake, which also have the highest percentages of LCC, indicating a positive relationship between the LCC and water surface change.
The rate change of LSWT is also not evenly distributed within each lake. The distribution of rate change appears to be impacted by the LCC and its type. For example, the LCC around Elodia is highly concentrated on the north side and much less LCC is taking place in the south section. Notably, the south section of Elodia Lake contains a less positive temperature trend than the rest of the lake (Figure 6b). Similarly, the southeastern part of Grippen Lake with more LCC experiences a stronger positive rate of temperature change than the middle of the lake with the strongest maximum positive rates of temperature change with a rate of 0.139 × 10 −4 °Kelvin per year. The change rate may also be impacted by lake size, shape, lake depth, watershed flow direction and terrain in the buffer area. For example, the water surface near the shoreline experiences a higher change rate of LSWT than water near the lake center. This is more evident in Grippen Lake, Horseshoe Lake and Little Lake (Figure 6C,D,G). The largest lake, South Lake, contains an island in the middle of it (Figure 6h). The eastern portion of South Lake displays a relatively The western edge of South Lake contains the lowest temperature trend increase, with small, fragmented areas of LCC nearby, while the shore north of the island contains the greatest positive temperature trends within the lake. There are larger extents of LCC in this north area which could help to explain this phenomenon. Figure 7 lists the scatter plots of the percentage of LCC and the minimum, maximum, and average rates of water temperature change of nine small lakes. The relationship between LCC and average rate of water surface change fits strongly positively with a R 2 value of 0.80 (Figure 7c). This relationship demonstrates that as the amount of LCC increases, as does the rate of water temperature warming for the lakes investigated. In Figure 6, there is a cluster of lakes around 14-18% LCC and a water temperature rate of 0.1 × 10 −4 -0.14 × 10 −4 °Kelvin per day. Although a water temperature rate of 0.1 × 10 −4 -0.14 × 10 −4 °Kelvin per day is relatively minimal, it translates to 0.0365 and 0.0511 °Kelvin per decade, which presents substantial long-term change.  Figure 7a illustrates the relationship between land cover change percentage and the minimum rate of water temperature change, which has a R 2 value of 0.72. This indicates that as the amount of LCC increases in the buffered region, so does the minimum rate of water temperature change. Even the least positive water temperature trend within a given lake in the study region experiences a response to increasing the LCC. A stronger relationship was shown between the maximum water temperature trend, which had an R 2 value of 0.79 (Figure 7b). This demonstrates that the maximum water temperature rate for each  that as the amount of LCC increases in the buffered region, so does the minimum rate of water temperature change. Even the least positive water temperature trend within a given lake in the study region experiences a response to increasing the LCC. A stronger relationship was shown between the maximum water temperature trend, which had an R 2 value of 0.79 (Figure 7b). This demonstrates that the maximum water temperature rate for each individual lake increases as the amount of overall LCC increases. Both relationships between the minimum and maximum water temperature trend are moderately positive, demonstrating the clear relationship between the two variables. All trends demonstrated illustrate a strong positive water temperature trend and LCC relationships which can result in the modification of biochemical compositions of some algal species [62].
The relationship between LCC and surface temperature is well understood; however, the influence of LCC near waterbodies outside of an urban environment remains less understood [63]. Contrary to previous findings that impervious surfaces had played a warming role on land surface warming, impervious surface change does not show any positive correlation with the lake surface temperature change (Table 4). This is likely caused by the small proportion of impervious surface. Surrounding these inland lake regions investigated, there were no major impervious surfaces within the study area except for local roads and small settlements. If there were significant settlements, such as a suburb or a town near within the buffered region, a higher rate of water temperature change would be expected. Previously, Yang, Yu, & Luo (2020) investigated the water temperature trend of lakes but only had hypothesized that impervious surface had played a role [64]. In Table 4, the most notable positive correlation was between the LC class of sparse vegetation surface, which had a correlation of 0.57. A sparse vegetation surface increase would be related to residential development of clearing dense forests and increasing residential grassland areas. It should also be noted that parks which have similar LC properties to those of sparse vegetation have been known to exhibit similar cooling effects [10]. These findings consistently demonstrate that vegetation decrease has a strong positive impact on lake surface water temperature increase. The strongest negative correlation was associated with dense vegetation (Table 4), as dense vegetation would be expected to be the most dominant LC for cooling. An interesting finding is the positive correlation associated between the water temperature trend and wetlands, which had a positive correlation of 0.52. This finding could be due to wetlands being close to shores and small vegetation within them contributing to warmer water temperatures. This is consistent with the findings by Xue et al. (2019) [25], who found that wetlands exhibited a strong positive relationship between the cooling capability of urban wetlands and the area, shape, and hydrologic connectivity of wetlands.

Conclusions
This research has illustrated the feasibility of the CCDC algorithm to accurately capture LCC and calculate LST trends with Google Earth Engine in rural small mid-latitude lakes. When compared with previous studies, this research investigated previously unexplored inland lake water temperature trends, LCC impacts on lake surface water temperatures, and made contributions to the analysis of LCC on water bodies by investigating long-term trends of lake water surface temperature. Individual lake buffers (500 m) were classified into five LC classes, while average LSWT trends from 1985 to 2020 were calculated for each waterbody. All LSWTs of the nine lakes investigated demonstrated positive trends with varying degrees of magnitude. Spatially, positive stronger water temperature trends were The selection of the studied waterbodies revolved around similar climatic effects over the study period, distance from another waterbody, and the inhibition of another waterbody in order to limit additional cooling effects. Hence, there was a limited selection of waterbodies within the same watershed. Although strong conclusions were made regarding the previous results, it is challenging to claim that they are all-encompassing, given the small sample size.
All dataset and workflow used in this research were directly conducted in GEE. GEE's datasets were not required to be downloaded, Google servers were used with GEE as opposed to an individual's computer, and all data was saved to the cloud. While GEE is a useful tool for geospatial analysis, it does restrict the user to functions within its interface. The CCDC algorithm was demonstrated to be useful at accomplishing accurate land classification, as was the SMW method for LST with GEE.
There were some uncertainties and limitations that should be clarified. First, the limitation of remote sensing studies was inevitable due to the number of clear images. The CCDC algorithm takes advantage of all clear available imagery throughout the study period to classify LC. It is important to note that water temperature trends analyzed in this study are a glimpse into the past based on image availability and their image acquisition time. Second, the land cover classification accuracy assessment was conducted only between 2010 and 2015 due to the availability of CCRS land cover data. The average 84.30% overall accuracy of CCDC LC classifications were achieved for all lakes for those two years based on a limited number of stratified 50 samples for each class. The accuracy of land cover classification maps of other years is assumed to be similar to those two years due to the consistency of the data and methods. It is possible that the CCRS land cover data contains errors and the accuracy assessment has sampling biases. Since the land cover changes are relatively small around the small inland lakes, the potential land cover classification errors might impact the assessment of the land cover change results.
Third, the 500 m buffers around the lakes were used as a region where land cover would have a direct impact on the lake water temperature. This selection was based on the conclusion from a previous study and the relatively flat terrain in the watershed study. This assumption may not hold true for all lakes due to the size and landscape variations. Future studies could explore how different buffer sizes could impact the results. Besides land cover, many other factors can contribute to the lake surface temperature changes. Due to global warming, the atmospheric temperature near the Earth's surface has a warming trend. This has impacted lake water and atmosphere interaction and lake water temperature. Similar to land cover classification errors, the remote sensing temperature products contain errors which could have impacted the LSWT trend values obtained through this research. The actual change rates of LSWT may be different from those values listed.
Effective accounting of the impact of LCC on water surface temperature changes can provide the knowledge for better management where needed. Mismanagement of LCC near waterbodies could lead towards warmer and unsuitable surface water temperatures for aquatic species and ecosystems [65][66][67]. Future research could expand the findings of this study by investigating additional lakes in different watersheds. Other factors to consider include the incorporation of bathymetric data and additional variables such as wind speed, surface pressure, and cloud cover. Investigating waterbodies in different watersheds could lead towards a better understanding of local inputs on water temperature trends. Exploring additional watersheds could provide the opportunity to determine the influence of regional variables such as rainfall, wind, and ice cover on water temperature trends.
The role that bathymetry plays on water temperature trends in small lakes remains unknown. It is understood that lake depth contributes to the relationship between lake depth and surface water temperature trends [68]. However, a challenge exists for obtaining lake depth for relatively small lakes, as only one lake (Grippen Lake) out of the nine lakes investigated in this study had existing bathymetric data. Additional variables such as wind speed, surface pressure, cloud cover, and Secchi depth can be included in future studies. These variables have previously been investigated, but only on much larger lakes (>25 km 2 ) than the ones in this study [34]. Without the completion of field work, data availability is a challenge when investigating smaller water bodies, especially readily available additional data apart from images on GEE.

Data Availability Statement:
The data used for this study are openly available through GEE. The code for processing the data is available on request from the corresponding author.