Evaluation of Methods for Estimating Lake Surface Water Temperature Using Landsat 8

: Changes in lake water temperature, observed with the greatest intensity during the last two decades, may signiﬁcantly affect the functioning of these unique ecosystems. Currently, in situ studies in Poland are conducted only for 38 lakes using the single-point method. The aim of this study was to develop a method for remote sensing monitoring of lake water temperature in a spatio-temporal context based on Landsat 8 imagery. For this purpose, using data obtained for 28 lakes from the period 2013–2020, linear regression (LM) and random forest (RF) models were developed to estimate surface water temperature. In addition, analysis of Landsat Level-2 Surface Temperature Science Product (LST-L2) data provided by United States Geological Survey (USGS) and the National Aeronautics and Space Administration (NASA) was performed. The remaining 10 lakes not previously used in the model development stage were used to validate model performance. The results showed that the most accurate estimation is possible using the RF method for which RMSE = 1.83 ◦ C and R 2 = 0.89, while RMSE = 3.68 ◦ C and R 2 = 0.8 for the LST-L2 method. We found that LST-L2 contains a systematic error in the coastal zone, which can be corrected and eventually improve the quality of estimation. The satellite-based method makes it possible to determine water temperature for all lakes in Poland at different times and to understand the inﬂuence of climatic factors affecting temperature at the regional scale. On the other hand, spatial presentation of thermics within individual lakes enables understanding the inﬂuence of local factors and morphometric conditions.


Introduction
Water ecosystems are characterized by considerable dynamics of energy and matter exchange, and therefore reflect modern transformations of the natural environment. Among the key properties of rivers and lakes is water temperature. Its course and distribution strongly determine a number of other processes [1]. It is therefore not surprising that the first water temperature measurements in lakes date back to the 18th century, and the frequency and accuracy of such measurements have been successively improving since. Traditional methods of temperature records based on in situ measurements are of point character, and are often not representative for complex and dynamic ecosystems such as lakes. The methods also have logistic limitations, and are costly and time consuming [2]. Moreover, the measurements are performed at various time intervals and at different depths, leading to high heterogeneity of the available data [3]. A complex understanding of the functioning of lake ecosystems requires collecting credible data on the state of water at the global scale, over long time periods, with high spatial and temporal resolution [2,4]. Data regarding the surface temperature of lakes have been collected using remote sensing for many years [5]. Satellite image data have been recorded at the global scale since 1972, when the first civil Earth observation mission commenced under the name Landsat (formerly known as the Earth Resources Technological Satellite). For the first time in the Landsat program, a sensor with a thermal band was installed on the Landsat 3 satellite (launched in 1978), but it quickly failed. Each of the following satellites of this program was equipped with sensors that recorded thermal radiation emitted from the Earth's surface. An increase in water temperature in lakes shows a rate approximate to that of the growth of regional air temperature [6,7], although in the region of the Great Lakes and in northern Europe, a faster increase in lake water temperature was recorded than that of the surrounding air [7,8]. Satellite sensors used for the observation of the temperature of the surface water layer record emitted radiation at a range of 7-13 µm. The temperature of a very thin water layer is measured that can be referred to the temperature of the epilimnion in the conditions of lack of wind, although the external water layer is usually cooler than the 50 cm layer of sub-surface water [9]. Thermal scanners mounted on planes measure the temperature of the water layer on the interface of air and water with a thickness of 10-20 µm, but it can be several tenths of a degree warmer than temperature measured by a conventional sensor [10]. The difference between the temperature of the surface water layer and temperature in water depths is determined by factors such as the time of day, cloudiness, and wind speed [11,12]. To estimate the surface water temperature, data from one ("mono-window method") or two ("split-window method") thermal bands of the Landsat satellite can be used. The latter method is based on the difference in energy intensity in two thermal bands, which has the potential to improve temperature estimation [13]. In research determining the correlation between water temperature measured in situ and estimated from satellite data, thermal bands with low spatial resolution (i.e., Advanced Very High Resolution Radiometer) obtained very high values of determination coefficients (R 2 > 0.9) [14]. Similar results were obtained using data from Landsat satellites with higher spatial resolution, although with lower temporal resolution. Tavares et al. [15] recommends the use of thermal data from the Landsat 7 ETM+ sensor for the observation of small lakes because at a relatively high spatial resolution the temperature estimation error is very close to the error obtained from the MODIS sensor (RMSE 1.07 and 1.05 • C, respectively). The usefulness of other medium-resolution thermal satellite sensors for water surface temperature assessment, such as Landsat TM, ETM+, TIRS [16][17][18][19][20] and Terra ASTER [21][22][23], was also investigated. Lake temperature is primarily determined by meteorological factors (insolation, cloudiness, air temperature, and wind speed), and to a lower degree by geomorphometric factors (surface area and depth) [24]. The temperature of the surface layer of the lake increases with its insolation, although an increase in temperature also causes cloudiness, because the amount of long-wave radiation reaching the surface increases [25]. A negative correlation occurs between water temperature and mean lake depth, i.e., shallow lakes warm up faster and to higher temperatures [26,27].
In the case of the territory of Poland, a lot of research has been recently conducted on lake water temperature, in reference to both measurement techniques [28,29] and improvement of methods of its modelling [30]. Polish lakes (approximately 7000 natural lakes in total) with multiannual in situ datasets are suitable for searching solutions in both of the aforementioned aspects of research. Water temperature measurements in Polish lakes have a history of several decades, and systematic observations commenced at the end of the 1950s. Over the years, the population of lakes covered by field measurements has changed, and today it concerns several tens of objects. No comprehensive research with the application of satellite images in reference to lake water temperature in Poland has been undertaken to date, although studies in the scope dynamically develop in different regions of the globe [31][32][33][34]. Polish research so far covers three lakes, and the obtained results revealed high cohesion with in situ measurements [35], encouraging a broader Remote Sens. 2022, 14, 3839 3 of 21 approach to the issue. Polish lakes already currently show a considerable increase in water temperature [29,36] that will continue to progress in the future [37,38].
The objective of this article is the comparison of four methods of estimation of water surface temperature in lakes in Poland based on thermal images from Landsat 8. Other detailed objectives include identification of limitations of the applied methods, and determination of the spatial variability of water temperature within the selected lakes and in a broader regional approach.

Materials and Methods
To organize the presentation of the results obtained in the subsequent stages, the following subsections were distinguished in this paper. Section 2.1 describes the study lakes, while Section 2.2 presents methods for conducting daily water temperature measurements. Sections 2.3 and 2.4 describe the acquisition of Landsat 8 Collection 2 Level 1 (L8L1) data and Landsat Level-2 Surface Temperature Science Product data (LST-L2). Sections 2.5-2.7 describe the method of calculating lake water temperature using linear regression model (LM), random forest regression (RF), and Ermida et al. (2020) [39] method (LST), respectively. Section 2.8 describes the validation process. Finally, Section 2.9 shows the possibility of using results to present the spatial variability of water temperature at the scale of a single lake and at the regional scale for multiple lakes. The flowchart of this study is shown in Figure 1. obtained results revealed high cohesion with in situ measurements [35], encouraging a broader approach to the issue. Polish lakes already currently show a considerable increase in water temperature [29,36] that will continue to progress in the future [37,38]. The objective of this article is the comparison of four methods of estimation of water surface temperature in lakes in Poland based on thermal images from Landsat 8. Other detailed objectives include identification of limitations of the applied methods, and determination of the spatial variability of water temperature within the selected lakes and in a broader regional approach.

Materials and Methods
To organize the presentation of the results obtained in the subsequent stages, the following subsections were distinguished in this paper. Section 2.1 describes the study lakes, while Section 2.2 presents methods for conducting daily water temperature measurements. Sections 2.3 and 2.4 describe the acquisition of Landsat 8 Collection 2 Level 1 (L8L1) data and Landsat Level-2 Surface Temperature Science Product data (LST-L2). Sections 2.5-2.7 describe the method of calculating lake water temperature using linear regression model (LM), random forest regression (RF), and Ermida et al. (2020) [39] method (LST), respectively. Section 2.8 describes the validation process. Finally, Section 2.9 shows the possibility of using results to present the spatial variability of water temperature at the scale of a single lake and at the regional scale for multiple lakes. The flowchart of this study is shown in Figure 1.

Study Sites Description
The analysis of the quality of water temperature determination based on satellite images from Landsat 8 covered 38 lakes in north Poland ( Figure 2). All lakes selected for the analysis are natural water bodies. Their genesis is primarily related to the course of geomorphological processes during the Weichselian glaciation. The maximum extent of the ice during the Weichselian in north Poland is marked with a dotted line in Figure 2. It also constitutes the boundary of three main lakelands. The northern area of Poland is the richest in lakes, covering more than 500 lakes with a surface area larger than 1 km 2 . Lakes selected for the study have a variable surface area from 1.6 to 111.9 km 2 , and mean depths from 1.3 to 38.7 m. In genetic terms, postglacial lakes are largely dominant (channel, moraine, and kettle) related both to the erosional and accumulation activity of the Scandinavian ice sheet. Lakes developed in the Holocene also occur, namely Jamno (No. 28), Gardno (No. 29) and Łebsko (No. 30) which are also the northernmost study objects. Lake Sławskie (No. 1) is the southernmost one. Lake Morzycko (No. 2) is the westernmost lake, and Lake Studzieniczne (No. 27) has the easternmost location. Details of the morphometric parameters of lakes are presented in Supplementary Table S1. Mean annual temperatures of the analyzed lakes are at a level from 8.6 to 11.4 • C. In the regional approach, evident differences are observed in the period of winter-spring, when ice cover is still present on the lakes in the eastern part of Poland, but already absent in the lakes of west Poland. It is due to the features of transitional climate of the analyzed region, when a longer effect of cold continental air is recorded in the east.

Study Sites Description
The analysis of the quality of water temperature determination based on satellite images from Landsat 8 covered 38 lakes in north Poland ( Figure 2). All lakes selected for the analysis are natural water bodies. Their genesis is primarily related to the course of geomorphological processes during the Weichselian glaciation. The maximum extent of the ice during the Weichselian in north Poland is marked with a dotted line in Figure 2. It also constitutes the boundary of three main lakelands. The northern area of Poland is the richest in lakes, covering more than 500 lakes with a surface area larger than 1 km 2 . Lakes selected for the study have a variable surface area from 1.6 to 111.9 km 2 , and mean depths from 1.3 to 38.7 m. In genetic terms, postglacial lakes are largely dominant (channel, moraine, and kettle) related both to the erosional and accumulation activity of the Scandinavian ice sheet. Lakes developed in the Holocene also occur, namely Jamno (No. 28), Gardno (No. 29) and Łebsko (No. 30) which are also the northernmost study objects. Lake Sławskie (No. 1) is the southernmost one. Lake Morzycko (No. 2) is the westernmost lake, and Lake Studzieniczne (No. 27) has the easternmost location. Details of the morphometric parameters of lakes are presented in Supplementary Table S1. Mean annual temperatures of the analyzed lakes are at a level from 8.6 to 11.4 °C. In the regional approach, evident differences are observed in the period of winter-spring, when ice cover is still present on the lakes in the eastern part of Poland, but already absent in the lakes of west Poland. It is due to the features of transitional climate of the analyzed region, when a longer effect of cold continental air is recorded in the east.

In Situ Data
In this study, the results of daily water temperature measurements conducted between 2013 and 2020 by the Institute of Meteorology and Water Management-National Research Institute (IMWM) were used. The period from April to October, when there is generally no ice cover on the analyzed lakes, was selected for the analysis. Water temperature measurements by the IMWM are performed with a frequency of once a day at 6 UTC. The measurement points are located in the coastal zone and the sensor is placed approximately 0.4 m under the water level. Temperature measurements are performed both automatically and manually (especially for previous years), and their accuracy is 0.1 • C. Based on earlier hourly measurements of lake water temperature in Poland [40], it was evidenced that the difference between water temperature measurements at 6 UTC (in accordance with the IMWM standard) and~9 UTC (overpass time of Landsat) in the summer half-year averages 0.1 • C, reflecting high thermal stability of water. Moreover, it should be emphasized that water temperature variability in the near-surface water layer (from the surface to a depth of 1 m) is inconsiderable, and as shown by earlier studies in the case of lakes in Poland [41], it averages 0.2 • C.

Landsat 8 TOA Data
In this paper, Landsat 8 data provided by the United States Geological Survey (USGS) was used. The data from April 2013 to October 2020 over north Poland were analyzed. The Landsat 8 Collection 2 Level 1 Tier 1 (L8L1) data were obtained using the Google Earth Engine (GEE; accessed on 24 February 2022) [42]. The L8L1 dataset contains calibrated top-of-atmosphere (TOA) reflectance derived from the data produced by the Operational Land Imager (OLI) and Thermal Infrared Sensor (TIRS). These L8L1 products contain 5 visible and near-infrared (VNIR) bands, 2 short-wave infrared (SWIR) bands with spatial resolution 30 m and 2 thermal infrared (TIR) bands resampled from 100 m to 30 m. Only scenes covering positions in the reference in situ data were considered ( Figure 3). Due to the fact that in situ lake water temperature measurement sites are located within the shoreline range, it is not possible to obtain L8L1 data for these locations specifically. That results in the pixels (resolution 30 m) that extend over land or shoreline vegetation as well. Therefore, L8L1 data were collected from a point located 60 m inside the lake in relation to the in situ lake water temperature measurement sites. The water difference in this neighborhood is negligible, and we have the confidence to exclude pixels that are land. For the analysis all scenes with cloud cover less than 60% were taken. Moreover, we excluded data marked as clouds and cloud shadows based on the CFMask algorithm [43]. CFMask derives from the Function of Mask (FMask) developed by Zhu and Woodcock [44] and Zhu et al. [45]. Additionally, we removed data with reflectance values below 0, temperature values below 0 • C in thermal band B10, and outliers based on the ultra-blue band (B1 values ≥ 0.14). Additionally, some observations were made on two scenes on the same day-we treated them as duplicates and removed them.
Due to cloud cover patterns in Poland, for each of the analyzed lakes a slightly different imagery dataset was obtained. The least, 38 L8L1 images were acquired for Białe Augustowskie Lake (No. 26), while 83 images were acquired for Niesłysz Lake (No. 3). The dates for Niesłysz Lake (No. 3) are shown in Figure 4. The red dots indicate imageries that were rejected from the analysis and did not meet previously defined criteria, while the green dots indicate imageries that were used for estimating water temperature in the lake. Only 43% of all available scenes were suitable for analysis. Supplementary Figure S1 shows the exact availability of satellite scenes by lake, month, year and water temperature. Remote Sens. 2022, 14, 3839 6 of 22 Due to cloud cover patterns in Poland, for each of the analyzed lakes a slightly different imagery dataset was obtained. The least, 38 L8L1 images were acquired for Białe Augustowskie Lake (No. 26), while 83 images were acquired for Niesłysz Lake (No. 3). The dates for Niesłysz Lake (No. 3) are shown in Figure 4. The red dots indicate imageries that were rejected from the analysis and did not meet previously defined criteria, while the green dots indicate imageries that were used for estimating water temperature in the lake. Only 43% of all available scenes were suitable for analysis. Supplementary Figure S1 shows the exact availability of satellite scenes by lake, month, year and water temperature.   Due to cloud cover patterns in Poland, for each of the analyzed lakes a slightly different imagery dataset was obtained. The least, 38 L8L1 images were acquired for Białe Augustowskie Lake (No. 26), while 83 images were acquired for Niesłysz Lake (No. 3). The dates for Niesłysz Lake (No. 3) are shown in Figure 4. The red dots indicate imageries that were rejected from the analysis and did not meet previously defined criteria, while the green dots indicate imageries that were used for estimating water temperature in the lake. Only 43% of all available scenes were suitable for analysis. Supplementary Figure S1 shows the exact availability of satellite scenes by lake, month, year and water temperature.

Landsat Level-2 Surface Temperature Science Product
In this paper, the results of Landsat Level-2 Surface Temperature Science Product provided by the USGS were used [46]. The LST-L2 product contains the temperature of the Earth's surface in Kelvin (K). This product is generated mainly from the top-of-atmosphere thermal infrared data using a single-channel algorithm. Nevertheless, the algorithm is very extensive and uses auxiliary data, including data from the ASTER satellite and atmospheric data from reanalysis. The product is available in a resolution of 30 m as georeferenced rasters (.tif format) in the Universal Transverse Mercator coordinate system. Information on the applied techniques related to atmospheric compensation, calibration methodology, and validation are widely presented and discussed by Cook [47], Cook et al. [48], Malakar et al. [18], Schaeffer et al. [49]. In the case of LST-L2, the data availability differed for all analyzed lakes due to cloud cover pattern in Poland. From among all analyzed lakes, the least number of LST-L2 images, i.e., 38, was acquired for Biale Augustowskie Lake (No. 26), while the most number of 88 was acquired for Jamno Lake (No. 28).

Simple Linear Model
We used multiple linear regression as the first choice method because of its simplicity. The purpose of this model is to determine the relationship between multiple explanatory variables and the modelled (response) variable. In other words, we want to predict the unknown values of one variable based on the known values of other variables using the following formula: where f (x)-is the response variable; b 1 , b 2 , b n -are regression coefficients calculated for individual explanatory variables; x 1 , x 2 , x n -are the values of the explanatory variables; a-is the intercept.
Therefore, our modelled variable is water temperature and our explanatory variables are the two thermal bands B10 and B11. The main assumption of linear regression is the linear relationship between the modelled variable and the predictors, which is obviously not clear in the case of complex natural processes. However, this is not the only assumptionfurther limitations are related to the lack of multicollinearity, independence and normal distribution of errors, and a limited number of explanatory variables. For this reason, a more robust machine learning model may be a better solution.

Random Forest Model
Random Forest (RF) is a method that enables the estimation of searched values based on a large set of uncorrelated and random trees [50]. In our work, we used it for regression analysis. RF is based on bagging and random subspace methods. Bagging is related to the creation of regression trees, which are then combined into ensembles to obtain an overall prediction. To design regression trees, a number of independent bootstrap samples are created from the original training dataset. Each bootstrap sample (D b ) is created by randomly sampling n subsamples from the original training data D, containing N subsamples. The bootstrap sample (D b ) usually consists of 2/3 of D and contains no duplicate subsamples. Then, K independent regression trees are created for bootstrap samples with input vector x. In the regression analyses, the average prediction of the K regression trees, h k (x), is calculated to obtain an overall prediction [51,52] as follow, Since the regression trees have high variance, the bagging is designed to reduce the variance and to prevent overfitting the complex RF model. Therefore, the learning trees cannot be correlated. Samples derived from the dataset that were not selected to train the k-th regression tree during the bagging process are compiled into an out-of-bag (OOB) subdataset. The OOB subdataset contains the remainder of D dataset ( 1 / 3 samples). Based on OOB subdataset, the performance of the k-th regression tree is calculated by the mean squared error (MSE OOB ) as, where y i and y i are the prediction and the mean of the i-th prediction from all regression trees. The determination coefficient R 2 OOB of the OOB subdataset, can be calculated based on MSE OOB and the total variance Var y , of the output parameter, using the following formula, The input data during the RF analyses were TOA reflectance values (B1, B2, B3, B4, B5, B6, B7, B10, and B11). Moreover, the spectral indices the Normalized Difference Vegetation Index (NDVI) [53] and the Normalized Difference Water Index (NDWI) [54], the month index as auxiliary variables to maximize prediction performance were used. All calculations related to the RF were accomplished using the ranger [55] package within the statistical software R 4.1.2 [56].

Land Surface Temperature Model
Ermida et al. (2020) [39] provides a code on GEE to derive LST from Landsat Collection 1 Level-1 thermal infrared bands. The method, referred to hereafter as LST, is based on the Statistical Mono-Window (SMW) algorithm developed by the Climate Monitoring Satellite Application Facility (CM-SAF) for deriving LST climate data records from Meteosat First and Second Generation [57]. The algorithm uses a single thermal infrared band (band 10 in the case of the Landsat 8) for consistency across all Landsat series. This method also uses the ASTER Global Emissivity Dataset (GED) database and applies vegetation cover correction based on Landsat NDVI. Since not all lakes considered in this study were correctly identified by the ASTER GED dataset, we have updated the code to use the Landsat water mask (available through the quality flag). For water bodies, the emissivity was set to 0.99.

Model Validation Procedure
To assess the accuracy of the LM, RF, LST and LST-L2 model predictions, we randomly selected 10 of the 38 lakes studied. In total, we used 538 in situ lake water temperature measurements to validate the temperature prediction results, which is approximately 23% of the input dataset. The root mean square error (RMSE), mean bias error (MBE), coefficient of Pearson correlation (R), between measured and predicted values were used for validation. Additionally, we calculated the standard deviation (SD) of the predicted values.

Model Application
Results obtained from the best model were used for the presentation of the spatial distribution of water temperature in lakes in the studied area at the local and regional scale. The objective of the presentation of the spatial variability of water temperature in single lakes was to show the role of local factors, i.e., the morphometry of the lake basin, range of impact of river waters, groundwater supply, and effect of wind on the distribution of water temperatures in the lake. The presentation of the thermal regime of lake waters at the regional scale aimed at illustrating the effect of climatic factors on the spatial distribution of water temperature in the lakes of north Poland in different time periods from April to October.
We tested the applicability of the developed RF model for the entire example scene of 28 April 2021 (ID: LC08_191023_20210428) and made the temperature prediction for 1568 lakes located in the north-western part of Poland with an area larger than 3 ha (approximately 30 pixels). On a local scale, we chose 4 lakes on different dates-Gopło, Łebsko, Drawsko, Ełckie (IDs: LC08_190023_20210608, LC08_191022_20150818, LC08_192023_20180529, and LC08_187022_20211025). We used the terra [58] package to process spatial data (both raster and vector) and perform spatial predictions.

Comparison of Methods
At the first stage of work, based on a training sample covering measurement data from 28 lakes (1828 measurements) and the corresponding data from Landsat 8 Collection 2 Level 1 (L8L1) meeting the criteria described in the methodology, a multiple linear regression model (LM) was developed, as well as a random forest model (RF). The following equation parameters were obtained for the linear model: Water temperature = 2.9 × B10 − 2.07 × B11 + 48.48 (5) where bands B10 and B11 represent Top of Atmosphere Reflectance in Kelvins. Considerably more credible from the point of view of further application of the model are results of the validation of statistical models on independent data not used at the stage of model development. The results are obviously worse than those obtained at the stage of model development. In this study, for LM, the obtained values R and RMSE were lower than in the case of RF model (Table 1). In the case of the LST model, RMSE values were almost twice higher than for results obtained from the RF model. The worst results of the estimation of lake water temperature were obtained for Landsat Level-2 Surface Temperature Science Product (LST-L2). In this case, MBE and RMSE values were 2.55 and 3.68 • C, respectively. Considering the fact that next to values of global statistics obtained for the entire dataset, the selection of the appropriate model for the estimation of water temperature in lakes requires broader analysis, at least in the temporal approach-of particular months. In the period from April to October, a substantial change in water temperature occurs in particular lakes (from 0.3 • C to 28.3 • C), and a thermal gradient occurs from the east to the west of Poland modified between lakes by the impact of the Baltic Sea. The analysis of results obtained from the validation of LM and RF models evidently shows that RMSE results are lower than those obtained from finished products LST-L2 and LST (Table 2). An advantage of LST-L2 and LST models is their global cover and no requirement for measurement data for their calibration. The comparison of RMSE values between LM and RF models using the Wilcoxon signed-rank test showed a significant difference between them (p-value < 0.05). It suggested that RMSE values are lower for RF model than for LM. Differences between the estimation of water temperature are particularly evident in April (more than 1 • C), and the smallest differences occur in August (0.13 • C). Considering all models, the lowest RMSE values occur in October, and the highest in May. The analysis of the values of correlation coefficients for all methods showed that they are at a similar level, and the occurring differences are not statistically significant. Results of the validation of the LM, RF, LST, and LST-L2 models are presented in Figure 5. The greatest dispersion of points presenting the comparison of measured and estimated values evidently occurs for April and May. According to Skowron [59], in that period water temperature shows high dynamics due to water mixing related to the development of stratification. Moreover, water temperatures estimated by means of the LST and LST-L2 models are evidently higher than the measured ones, suggesting the occurrence of a systematic error. For models RF and LM, the dispersion of points focuses within a straight line of y = x, suggesting the occurrence of only incidental/random errors. Differences between measured water temperatures and those estimated by means of particular models are presented in Figure 6. Distributions of deviations obtained for the LST and LST-L2 models are evidently shifted right from the value of 0, pointing to the occurrence of systematic errors of 2.0 and 2.6 °C, respectively (for a testing sample of 10 lakes, 538 measurements). Due to this, improvement of the estimation of water temperatures can be relatively easily obtained by means of LST and LST-L2 by introducing a correction factor. Results of such analyses are presented in the next sub-chapter. In the Differences between measured water temperatures and those estimated by means of particular models are presented in Figure 6. Distributions of deviations obtained for the LST and LST-L2 models are evidently shifted right from the value of 0, pointing to the occurrence of systematic errors of 2.0 and 2.6 • C, respectively (for a testing sample of 10 lakes, 538 measurements). Due to this, improvement of the estimation of water temperatures can be relatively easily obtained by means of LST and LST-L2 by introducing a correction factor. Results of such analyses are presented in the next sub-chapter. In the case of the remaining LM and RF models, the average deviation between the measured and forecasted value is 0.0 and 0.1 • C, respectively. This points to the lack of occurrence of systematic errors. Histograms between actual and predicted values by all models indicate that they are normally distributed, although this is not apparent from the Shapiro-Wilk tests (p-value < 0.05). This is mainly due to single outliers. Figure 5. Scatterplots for the models tested. The black dashed line indicates the ideal prediction. A systematic error can be noticed for the LST and LST-L2 models. Differences between measured water temperatures and those estimated by means of particular models are presented in Figure 6. Distributions of deviations obtained for the LST and LST-L2 models are evidently shifted right from the value of 0, pointing to the occurrence of systematic errors of 2.0 and 2.6 °C, respectively (for a testing sample of 10 lakes, 538 measurements). Due to this, improvement of the estimation of water temperatures can be relatively easily obtained by means of LST and LST-L2 by introducing a correction factor. Results of such analyses are presented in the next sub-chapter. In the case of the remaining LM and RF models, the average deviation between the measured and forecasted value is 0.0 and 0.1 °C, respectively. This points to the lack of occurrence of systematic errors. Histograms between actual and predicted values by all models indicate that they are normally distributed, although this is not apparent from the Shapiro-Wilk tests (p-value < 0.05). This is mainly due to single outliers. The results obtained both globally and for individual months suggest that the best estimate of lake water temperature is possible using the RF model. This probably results, among others, from the fact that the collection of explanatory (independent) variables covers as many as 12 parameters (B1, B2, B3, B4, B5, B6, B7, B10, B11, the NDVI, the NDWI and the month index), and in LM only 2 parameters are used (B10 and B11). Therefore, an analysis of the suitability of the explanatory variables for estimating water temperature using the impurity algorithm from the random forest model was conducted (Figure 7). We found that B10 was more significant than B11, which is consistent with other publications. This is also confirmed by the result obtained for the linear model-a slightly higher Pearson correlation coefficient was noted for band B10 (R = 0.88) than for band B11 (R = 0.84). The results of the study by Barsi et al. [60] showed that due to the occurrence of out-of-field stray light, the calibration system of the TIR sensor is unstable, and the recorded thermal data may carry a substantial error. The effect of dispersed light on data recorded in band B10 is relatively small and therefore should be preferred. Considerably greater errors can occur in the case of data in band B11, they should be subject to additional external calibration. Another important variable is the month index-this is related to the seasonality of water temperature. Specifying the month significantly reduces errors resulting from incorrect satellite measurements. For example, in the warm months, we expect statistically higher temperatures than in the cold months. Thus, excluding weather anomalies, it can be expected that the water temperature range in June will be approximately from 18 • C (lower quartile) to 21.6 • C (upper quartile). The other variables are of little importance. Of the remaining spectral bands, the green band (B3) is the most significant. Green radiation that carries a lot of energy easily penetrates water. The depth of such penetration determines the magnitude of changes in water temperature, and depends on, among others, water transparency. Therefore, the amount of radiation reflected from water recorded in the green band can be related to water temperature [61][62][63]. months, we expect statistically higher temperatures than in the cold months. Thus, excluding weather anomalies, it can be expected that the water temperature range in June will be approximately from 18 °C (lower quartile) to 21.6 °C (upper quartile). The other variables are of little importance. Of the remaining spectral bands, the green band (B3) is the most significant. Green radiation that carries a lot of energy easily penetrates water. The depth of such penetration determines the magnitude of changes in water temperature, and depends on, among others, water transparency. Therefore, the amount of radiation reflected from water recorded in the green band can be related to water temperature [61][62][63].
We also evaluated the usefulness of spectral indices, the NDVI (used to calculate emissivity) and the NDWI (used to distinguish water from land), but they did not significantly affect the results. We also evaluated the usefulness of spectral indices, the NDVI (used to calculate emissivity) and the NDWI (used to distinguish water from land), but they did not significantly affect the results.
Although the LM provides somewhat worse prediction results than the RF model, it has its evident advantages. Predictions of RF model employ 12 independent variables, which in the case of the analysis of the entire Landsat scene requires a total of approximately 5.4 GB RAM. In the case of LM, it will be only approximately 900 MB (2 variables). Next to the number of variables, the complexity of the model is also important-the linear model is much simpler, allowing for considerably faster prediction in terms of calculations. The last difference is related to the application of the month index-in the case of RF, it substantially improves prediction results, although it results in limiting the applicability of the model only to the geographic area of Poland (approximate latitudes and longitudes related to climate). In LM, no such assumptions were adopted; therefore, the potential of its applications is greater. Nonetheless, part of predictors can also be excluded post hoc in RF with no loss of effectiveness.
The spatial analysis of values of RMSE characteristics for the best RF model shows better adjustment of the model to the training sample of lakes (orange color) than for the test sample (points marked with red color) (Figure 8). The spatial analysis of RMSE values within the training sample showed no occurrence of spatial autocorrelation, i.e., grouping of lakes with evidently lower or higher RMSE values. Among the lakes used for testing the model, the lowest RMSE values were obtained for Lake Ełckie (No. ues within the training sample showed no occurrence of spatial autocorrelation, i.e., grouping of lakes with evidently lower or higher RMSE values. Among the lakes used for testing the model, the lowest RMSE values were obtained for Lake Ełckie (No. 23): 1.40 °C. The highest RMSE values were obtained for Lake Hańcza (No. 37): 2.83 °C. Both Lakes Hańcza and Ełckie are located in east Poland. The correlation analysis between RMSE values and lakes' areas, depths and volumes showed no significant relationship for both the training set and the test set.

LST-L2 Calibration
Because LST-L2 data are a finished product that can be practically applied by a large group of scientists, including hydrologists, limnologists, and water ecologists with no knowledge of specialized tools for obtaining and processing satellite data or development of statistical models (such as LM or RF), a decision was made to develop a calibration system for Poland for the estimation of water temperatures in lakes based on LST-L2 data provided by USGS. The analysis employed an identical method of determination of calibration parameters (based on 1828 measurements) and validation (538 measure-

LST-L2 Calibration
Because LST-L2 data are a finished product that can be practically applied by a large group of scientists, including hydrologists, limnologists, and water ecologists with no knowledge of specialized tools for obtaining and processing satellite data or development of statistical models (such as LM or RF), a decision was made to develop a calibration system for Poland for the estimation of water temperatures in lakes based on LST-L2 data provided by USGS. The analysis employed an identical method of determination of calibration parameters (based on 1828 measurements) and validation (538 measurements) as described in Section 2. The comparison of LST-L2 values with values of temperature measured in situ showed the occurrence of a systematic error of approximately 2.42 • C (Figure 9).  (Figure 9). We tested three ways of correcting it using a linear equation (with and without an intercept) and adding the mean value of this error. We obtained successively RMSE values on the test set: 2.59, 2.86, and 2.88 °C (the initial RMSE value of LST-L2 was 3.68 °C). We strongly recommend calibrating the LST-L2 data for estimating the water surface temperature by using the formula below: LST-L2 corrected = 0.806 × LST-L2 + 54.37 (6) Figure 9. The difference between in situ temperature and LST-L2 in the scatter plot (first column) and histogram (second column). The first row represents uncalibrated training data, while the second row represents calibrated testing data. The blue dashed line indicates the mean systematic error of approximately 2.42 K. After the calibration, the systematic error was corrected.

Discussion
Satellite images are an important tool in the case of research on water ecosystems, and find application in reference to many issues concerning both biotic [64][65][66] and abiotic processes [67][68][69]. A number of papers refer to the elementary feature, namely water temperature, and the applied methods analyze the issue at various temporal and spatial scales [15,[70][71][72][73][74]. As mentioned in the introduction to this paper, research on lake water Figure 9. The difference between in situ temperature and LST-L2 in the scatter plot (first column) and histogram (second column). The first row represents uncalibrated training data, while the second row represents calibrated testing data. The blue dashed line indicates the mean systematic error of approximately 2.42 K. After the calibration, the systematic error was corrected.
We tested three ways of correcting it using a linear equation (with and without an intercept) and adding the mean value of this error. We obtained successively RMSE values on the test set: 2.59, 2.86, and 2.88 • C (the initial RMSE value of LST-L2 was 3.68 • C). We strongly recommend calibrating the LST-L2 data for estimating the water surface temperature by using the formula below: LST-L2 corrected = 0.806 × LST-L2 + 54.37 (6)

Discussion
Satellite images are an important tool in the case of research on water ecosystems, and find application in reference to many issues concerning both biotic [64][65][66] and abiotic processes [67][68][69]. A number of papers refer to the elementary feature, namely water temperature, and the applied methods analyze the issue at various temporal and spatial scales [15,[70][71][72][73][74]. As mentioned in the introduction to this paper, research on lake water temperature in Poland with the application of satellite images has found no broader application to date, as confirmed by scarce papers in the scope. Ptak et al. [35] compared in situ water measurements with Landsat images (4, 5 and 7) for three coastal lakes, obtaining high coefficients of determination R 2 (0.87-0.95). Based on in situ measurements and Landsat 8 data regarding surface water temperature in Lake Raduńskie Górne, Szkwarek and Wochna [75] obtained credible results allowing for detailed presentation of its spatial distribution. Monitoring programs are often based on a limited spatial and temporal range. Remote sensing offers promising tools for large-scale observation, improving our possibilities of comprehensive research of the indices of lake properties [32]. Therefore, in the context of dynamically progressing remote sensing technology, its limited scope of applicability in the analysis of lake water temperature in Poland constitutes a research gap that is attempted to be filled based on Landsat 8.
The application of Landsat images in this paper corresponds with the intensively addressed global research trend regarding lake water temperature. The applied methodical approach showed an RMSE error at a level of 1.8 • C, a result approximate to that in other studies using Landsat images. Simon et al. (2014), analyzing two freshwater lakes in France, determined the relationship between Landsat thermal data and water temperature by measuring the temperature in situ at various depths, from 0 to 0.5 m [20]. The authors concluded that the obtained results (R 2 ranging from 0.90 to 0.94 and RMSE ranging from 1.753 to 2.397 • C) were satisfactory and coherent with literature data, where the values of mean square error are within a range from 1 to 2 • C. Giardino et al. [17] presented the Landsat TM-derived surface temperature of the sub-alpine lake with the RMSE of 0.328 • C; however, the procedure they used was not based on empirical relationships between satellite data and in situ measurements. For several hundred lakes analyzed based on data from the period 1999-2016 (Landsat 5 and 7), the total RMSE of temperature measurements obtained from the satellite was approximately 1.2 • C [19].
Mean absolute error for the analyzed lakes averaging from 1.38 to 2.39 • C (depending on the applied method) was approximate to the results of research conducted by Schaeffer et al. (2018) for lake waters in the United States, and reached 1.34 and 4.89 • C [49]. Results obtained in the study were also very approximate to those obtained by Tavares et al. (2019), who compared different methods of temperature estimation based on data from Landsat 7 and MODIS satellite [15]. In the estimation of water temperature (LSWT), the authors applied the single-channel algorithm method and two databases, SAFREE and TIGR3, obtaining similar R 2 values (0.936 and 0.938, respectively) as in this study. A weaker correlation with in situ water temperature was also shown by MOD28 data from the MODIS satellite (R 2 = 0.906), and a stronger one by MOD11 data (R 2 = 0.962). The most useful for the estimation of temperature (R 2 = 0.964) proved to be data from Landsat 7 satellite with the application of radiative transfer equation applied with atmospheric correction parameters from AtmCorr.
Considering the lack of possibility of conducting regular field observations, satellite remote sensing was recognized as a cost-effective way of monitoring water surface temperature at large spatial and temporal scales [76]. This statement is confirmed by results obtained in this paper, as exemplified by Figure 10, presenting a set of thermal data regarding three hundred lakes-including only several monitored through field observations. It should be emphasized that it is the first elaboration of the type for a fragment of the Polish territory, providing a promising basis for conducting monitoring of water temperature in lakes for a population of objects not achievable so far. Because water temperature is an elementary property of lake ecosystems [77,78] directly determining a number of processes and phenomena (biotic and abiotic), detailed information in the scope is necessary in reference to each lake. The obtained spatial image presented in Figure 10 points to the thermal variability of the analyzed lakes. The distribution in the orientation from the north (middle fragment) to the south is generally relatively clear, with occurrence of, respectively, cooler and warmer lakes. The situation illustrates the effect of a combination of environmental factors responsible for the thermal conditions of lakes. They include the surface area of the lake, land use structure around the lake, length, location (coordinates), and elevation above sea level [36]. According to the conducted research, multiple regression results showed that the recorded spatial distribution was significantly determined by, among others, the latitude, surface area, distance from the sea, presence of forests, and elevation above sea level. Similar research conducted for lakes on the Arctic coastal plain showed that, among others, parameters such as geographic location and morphometry of the lakes control their temperature at the regional scale [79]. The analyzed lakes with cooler water temperature occurred with greater density in the catchments of the Gwda and Drawa Rivers, on postglacial outwash plains built of highly permeable formations, i.e., gravels and sands. The situation is evident in the case The analyzed lakes with cooler water temperature occurred with greater density in the catchments of the Gwda and Drawa Rivers, on postglacial outwash plains built of highly permeable formations, i.e., gravels and sands. The situation is evident in the case of runoff, where the groundwater component is dominant in both of the aforementioned catchments [80], causing greater supply of cooler waters to the lakes.
Detailed distribution of water temperature is presented for two lakes with variable morphometric parameters (Figure 11). In both cases, certain patterns were recorded also observed in other studies of the type, i.e., distribution of water temperature referring to the depth of particular sectors of the lake and changes in water temperature in estuarial zones of rivers. According to Cao et al. [81], temperature in the center of the East Lake (Wuhan, China) is higher than in the surrounding area. Maps of surface temperature obtained by means of Landsat 8 satellite for three lakes (Maggiore, Lugano, and Como) showed variability reaching two or three degrees between different areas, i.e., between centers of lakes and coastal zones where inlets of rivers are located [82]. Based on data from Landsat 8 satellite, Sener and Sener [83] determined surface water temperature for Lake Beysehir (Turkey), recording a range of water temperature from 20.1 to 26.8 • C. They also determined that it was the highest in shallower sectors of the analyzed lake. Aïtelghazi et al. [84] emphasise that important parameters to be considered in the context of thermal changes include depth and suspended particles transported by the river. In the case of dimictic Lake Drawsko (Figure 11a), the distribution of water temperature in reference to depth is evident-the warmest in the coastal zone, and coldest in the deepest areas. A different situation is recorded in the case of polymictic Lake Łebsko (Figure 11b), where water temperature is uniform throughout the water column, and its variability on the surface is determined by inflows of surface waters-the coldest zone occurs in the southeastern part (mouth of the Łeba River with average water flow), and in the north-eastern part (inflow of marine water), as well as in the north-western part, where a cooler water mass was pushed to an isolated bay with limited circulation. north-eastern part (inflow of marine water), as well as in the north-western part, where a cooler water mass was pushed to an isolated bay with limited circulation. The potential of Landsat images for monitoring small water bodies can be used for forecasting their evolution, and will be useful in the future [85]. It is important in the context of the progressing climate change and the related threat in reference to, among The potential of Landsat images for monitoring small water bodies can be used for forecasting their evolution, and will be useful in the future [85]. It is important in the context of the progressing climate change and the related threat in reference to, among others, water quality or species composition of the fauna and flora of lake ecosystems.

Conclusions
This study concerned the comparison of the effectiveness of four methods of estimation of lake temperature based on satellite data. Two empirical models were developed (linear regression and random forest), and compared with two state of the art approaches employing atmospheric calibration. The comparative analysis shows the best performance of the random forest method, providing the highest correlation coefficient and the lowest RMSE.
Many papers recommend the application of atmospheric correction before using satellite data for the estimation of lake water temperature. It should be remembered, however, that the accuracy of these correction methods can be considerably reduced under certain conditions, e.g., with high content of water vapor in the atmosphere, and atmospheric models themselves show substantial limitations in spatial and temporal resolution. Therefore, continuous work on the improvement of correction algorithms should be accompanied by research on the possibility of application of new methods of data analysis (e.g., machine learning) and consideration of empirical data in the improvement of prediction results. This paper proposes an appropriate calibration correction for data from Landsat Level-2 Surface Temperature Science Product (LST-L2) that considerably improves results (root mean square error reduced by 30%) of estimation of lake temperature in the territory of Poland. Nevertheless, despite the correction, the random forest model provided the best results.
Future works should focus on determining the pixel based temperature uncertainty, analyzing the potential of using low-resolution data from the MODIS sensor for small surface lakes in Poland and determining the effectiveness of other water and vegetation indices. As of 31 October 2021, imagery is available from the next Landsat series satellite with the improved TIRS-2 sensor. Along with the increase in data availability, it should also be used for analysis to maintain continuous monitoring.

Supplementary Materials:
The following supporting information can be downloaded at: https:// www.mdpi.com/article/10.3390/rs14153839/s1, Table S1: Morphometric parameters of the analyzed lakes. Table S2: Linear and random forest model performance per lake on train and test sets. Test lakes are marked with "x". Figure S1: Availability of satellite scenes depending on (a) lake; (b) month; (c) year; (d) water temperature (as intervals). Reference [86] is cited in the Supplementary Table S1.
Author Contributions: K.D.: methodology, software, validation, formal analysis, data curation, writing-original draft preparation, review, editing, and visualization. S.E.: software and writingreview, editing. M.P.: writing-original draft preparation, review, and editing. J.P.: writing-original draft preparation, review, and editing. M.S.: formal analysis, writing-original draft preparation, review, editing, and visualization. All authors have read and agreed to the published version of the manuscript.
Funding: This research received no external funding.

Data Availability Statement:
In this study, we used publicly available datasets. Hydrologic data, satellite measurements, random forest regression model, and statistical analysis results are available in this repository: https://github.com/kadyb/lakes_temp (accessed on 1 May 2022). This repository also includes fully functional and reproducible codes under the MIT license.