Comparison of Methods to Estimate Lake-Surface-Water Temperature Using Landsat 7 ETM+ and MODIS Imagery: Case Study of a Large Shallow Subtropical Lake in Southern Brazil

Abstract

Water temperature regulates many processes in lakes; therefore, evaluating it is essential to understand its ecological status and functioning, and to comprehend the impact of climate change. Although few studies assessed the accuracy of individual sensors in estimating lake-surface-water temperature (LSWT), comparative analysis considering different sensors is still needed. This study evaluated the performance of two thermal sensors, MODIS and Landsat 7 ETM+, and used Landsat methods to estimate the SWT of a large subtropical lake. MODIS products MOD11 LST and MOD28 SST were used for comparison. For the Landsat images, the radiative transfer equation (RTE), using NASA’s Atmospheric Correction Parameter Calculator (AtmCorr) parameters, was compared with the single-channel algorithm in different approaches. Our results showed that MOD11 obtained the highest accuracy (RMSE of 1.05 ${}^{°}$C), and is the recommended product for LSWT studies. For Landsat-derived SWT, AtmCorr obtained the highest accuracy (RMSE of 1.07 ${}^{°}$C) and is the recommended method for small lakes. Sensitivity analysis showed that Landsat-derived LSWT using the RTE is very sensitive to atmospheric parameters and emissivity. A discussion of the main error sources was conducted. We recommend that similar tests be applied for Landsat imagery on different lakes, further studies on algorithms to correct the cool-skin effect in inland waters, and tests of different emissivity values to verify if it can compensate for this effect, in an effort to improve the accuracy of these estimates.

1. Introduction

Water temperature regulates many physical, chemical, and biological processes in lakes and reservoirs, such as reaction and metabolic rates, species distribution, nutrient cycling, and the concentration of dissolved gases [1]. Therefore, evaluating water temperature in lakes is essential to understanding their function and ecological status. In addition, lakes are considered sentinels of climate change due to their sensitivity to the environment [2], and lake-surface-water temperature (LSWT) is a proxy to analyse water-quality conditions and the impact of climate change on these systems [3,4]. A conventional way to assess LSWT is by using in situ sensors to measure temperatures at specific monitoring stations within a system. However, this approach requires the appropriate infrastructure and field work to install and maintain the sensors, which is expensive and time consuming. Moreover, for large aquatic ecosystems, this approach demands an extensive monitoring network to adequately measure the LSWT due to the spatiotemporal heterogeneity of these systems [5,6]. Unfortunately, technical and financial resources are not always available for such monitoring.

Data from satellite sensors may provide better information on LSWT variability than conventional field monitoring, since most modern sensors have improved capabilities with respect to spectral (e.g., spectral range of the thermal infrared), radiometric, temporal, and spatial resolutions [5,7]. Although satellite-derived water temperature is a description of the LSWT only in the top layer (i.e., approximately the upper 100 $\mathsf{\mu }$m, called “skin temperature”), it may provide important information about the patterns of water-temperature variations in lakes, and has been used in many studies, such as analysis of temperature patterns and heat budget [8,9], spatial distribution of water-quality variables [4], evaporation estimation [10], LSWT spatial gradients [6,11,12,13], LSWT temporal-variation patterns [14,15], and the assessment of the impact of El Niño/La Niña events [16] and climate change [3,17,18] on lakes.

Thermal infrared (TIR) bands are usually employed to retrieve accurate estimations of LSWT from satellite sensors, which are available only from a restricted number of these sensors. Two widely used thermal sensors for retrieving surface temperature (ST) are the Moderate Resolution Imaging Spectroradiometer (MODIS), aboard the Terra and Aqua satellites, and Thematic Mapper (TM), Enhanced Thematic Mapper Plus (ETM+) and Thermal Infrared Sensor (TIRS), aboard Landsat missions 4 to 8 (e.g., References [7,19,20,21]). Landsat sensors have a temporal resolution of 16 days, and a spatial resolution of 30 m, with the exception of their thermal bands, which have a spatial resolution of 60–120 m. MODIS thermal bands (31 and 32) have a spatial resolution of 1000 m and a daily overpass; however, its products have a resolution of 1–9 km [22]. The Terra and Landsat satellites have an identical orbit, with about 15 min between overpass times, which makes their sensor measurements comparable. For lakes and reservoirs, these sensors can be used for monitoring SWT, depending on the size of the water body as well as the purpose of the study (e.g., References [23,24]). Besides spatial and temporal resolution, radiometric precision differs among these sensors: MODIS TIR bands have a noise-equivalent change in temperature (NE$\Delta$T) of 0.05 ${}^{°}$C at 27 ${}^{°}$C [25], while Landsat 7 ETM+ thermal band has a NE$\Delta$T of 0.22 ${}^{°}$C at 7 ${}^{°}$C [26].

Most LSWT studies have evaluated the accuracy of individual sensors [27], such as MODIS (e.g., References [6,10,16,23]) and Landsat (e.g., References [12,24,28]), indicating the limitations of each sensor, namely, the spatial/temporal resolution trade-off, different Landsat methodologies, and sensor accuracy. Therefore, comparative analysis considering different sensors to estimate LSWT is still needed and may reveal important differences and limitations in the application of each sensor and direct studies [29]. Furthermore, while MODIS has SWT products that can be directly downloaded and used in these studies, for Landsat imagery, the SWT needs to be derived by the user, and it can be done using a few available methods [20]. Moreover, as noted by Allan et al. [12], these methods still have to be further validated using ground measurements, especially in the Southern Hemisphere. Additionally, few studies discuss the sensitivity of these methods to their input variables, such as water emissivity and atmospheric water vapour content, which can significantly impact the result [30,31,32,33]. Since few studies have assessed these methods [12,24,34,35] from different perspectives, it is also important to test and compare these multiple methods of temperature estimates.

This study evaluated the performance of two thermal sensors, MODIS and Landsat 7 ETM+, and Landsat methods to estimate the SWT of a large subtropical lake in Southern Brazil.

2. Materials and Methods

2.1. Study Area and Field Measurements

The study subject is Lake Mangueira, a large shallow subtropical lake located in the coast of Rio Grande do Sul, close to the Brazil–Uruguay border. It has a surface area of 820 km${}^2$ and is elongated, with a length of 90 km and width of 3–10 km, mean depth of 2.6 m, and maximum depth of 6.9 m. Its trophic state ranges from oligotrophic to mesotrophic [36,37]. The local climate is subtropical (Cfa type in Köppen’s climate classification) [38], with a mean annual temperature of 16 ${}^{°}$C and annual precipitation of 1800–2200 mm. The water-temperature time series was measured over 18 years, between 2001 and 2018, at the water surface at three stations across Lake Mangueira, namely, TAMAN (Northern), TAMAC (Central), and TAMAS (Southern), as shown in Figure 1. All measurements were made during the morning, and temperature was considered steady during this time of day (i.e., there was no variation in time between temperature measurements and satellite overpass at around 10:30 a.m. local time), as there is little diel variation in temperature (around 0.5 ${}^{°}$C between morning and night) in large lakes [39].

To compare the water-temperature estimates, one Landsat sensor and two MODIS products were used. Landsat 7 ETM+ was chosen over other Landsat sensors because it has a time series that covers the entire period of the measured data, providing more matching dates of observed temperatures and satellite images, whereas the Landsat 5 mission ended in 2011, and Landsat 8 OLI was launched in 2013, resulting in fewer suitable images from these two sensors. In the case of MODIS, both MOD11 and MOD28, which have been used in previous studies of LSWT (e.g., References [8,10,16,23,41,42]), are available for Lake Mangueira, and are described in the following sections.

Landsat 7 ETM+ and MODIS images were selected in order to compare the estimated LSWT with the measured data. Cloud cover is of major concern in SWT estimation using the TIR bands, since water vapour both absorbs and emits radiance in this spectrum [19], limiting the SWT estimation to cloudfree pixels. Since the Landsat images coincided with water-temperature measurements on only a few days (cloud cover substantially reduced the number of available dates), in order to increase the available data, in this study we used a window of two days of difference between the measurement and the satellite image, as in large lakes the water temperature varies only slightly over this period of time [24,43]. When Landsat images were from a different date than the measurement, MODIS images from the same Landsat date were also used when possible, to better compare them. The dates used in this study, with the availability of images for each sensor and station, are shown in Table S1.

In total, 26 dates were used, spread over all seasons, with 68 in situ temperature measurements (considering that, in some field campaigns, water temperature was not measured at all three monitoring stations). Altogether, considering all three stations, there were 58 Landsat 7 ETM+ estimations, 46 MOD28 estimations, and 64 MOD11 estimations. In order to assess the performance of each method and sensor, we used bias, absolute error, root mean squared error (RMSE), and the correlation coefficient ($r^2$). A statistical procedure was done in order to verify the assumption of normality (Shapiro–Wilk test) and equal variance (Bartlett’s test) of the residuals (difference between the calculated and measured temperatures), as well as to assess whether the performance of the methods was statistically different (ANOVA followed by Tukey’s honestly significant difference test).

The procedure used in this work followed the flowchart in Figure 2. The following sections detail the methods used for each sensor.

2.2. MODIS Products

MODIS has two products of surface temperature available for download: MOD28, which is a product of sea-surface temperature (SST), and MOD11, which is a product of land-surface temperature (LST). In the case of the coastal Lake Mangueira, both products are available. MOD28 is available with two different spatial resolutions, 4 and 9 km, and based on two different wavelengths, in the TIR (∼11–12 $\mathsf{\mu }$m) and in the medium wavelength infrared (MWIR, ∼3.5–4.1 $\mathsf{\mu }$m) spectra [22]. For comparison with Landsat images, for which temperature estimates were made based on the TIR band, and because SST derived from the MWIR are only for night time, only MODIS TIR products were used.

MODIS has two thermal bands, band 31, with wavelengths in the range of 10.78–11.28 $\mathsf{\mu }$m, and band 32, in the range of 11.77–12.27 $\mathsf{\mu }$m. For the estimation of the ST, both products use a split-window algorithm, which bases the correction of the atmospheric effects on the variation of the radiance measured by each TIR band, although they have differences, which are briefly described next.

2.2.1. MOD11 LST

The MOD11 LST product uses other bands (e.g., red, near infrared) to determine land cover, and the ST is derived using generalised split-window algorithm of Wan and Dozier [44] for pixels with known emissivity, which is the case with water. It considers the difference of emissivity between bands 31 and 32, along with constants calibrated for the sensor, to estimate the ST. The main issue with its LSWT estimation is the assumption that emissivity is constant for water, whereas it varies with many components of water, such as salinity, for example, which can result in errors of 0.2–0.5 ${}^{°}$C [32,45].

MOD11A1 LST images are daily and have a spatial resolution of 1 km, and are available directly via NASA’s LAADS DAAC page or on demand (available on https://ladsweb.modaps.eosdis.nasa.gov/search/), where users can select the products, and the images can be reprojected to the desired projection/UTM zone. MODIS Collection 6 LST products for the dates shown in Table S1 were downloaded, and the ST was extracted for each date and for each station using ArcGIS 10.3.

2.2.2. MOD28 SST

MOD28 SST is derived using a split-window algorithm, which is a modified version of the nonlinear split-window algorithm of Walton et al. [46]. It depends on brightness temperatures calculated by each thermal band and an a priori (first-approximation) estimate of the ST, and uses empirical coefficients for each month of the year, derived from in situ buoy and satellite observations and different atmospheric profiles based on latitude. The main issues with to the SST estimate and its application to lakes, besides its coarse spatial resolution, are the effect of wind on the relationship between skin and subsurface temperatures, and the fact that it was designed for sea water, so emissivity can be slightly overestimated. For more information on MOD28 characteristics, please refer to Kilpatrick et al. [22].

MOD28 SST is available directly via NASA’s PODAAC/JPL FTP (available on ftp://podaac-ftp.jpl.nasa.gov/OceanTemperature/modis/L3/terra/11um/v2014.0/4km/daily/). As was done with MOD11, the SST products for the dates shown in Table S1 were downloaded, and the ST was extracted for each date and for each station using ArcGIS 10.3.

2.3. Landsat 7 ETM+ Images

2.3.1. Methods to Estimate LSWT

While MODIS has specific products of ST, with standard algorithms designed and applied by the MODIS team (e.g., Reference [44]), the user needs to estimate the ST from Landsat images, using the radiative transfer equation (RTE) or algorithms developed for this purpose [20].

The radiance measured by sensors aboard satellites is different from that leaving the surface, since the atmosphere, mainly water vapour, both absorbs and emits radiance in this spectrum [31]. For a precise estimate of the ST, it is necessary to remove these atmospheric effects considering three atmospheric parameters: atmospheric transmittance ($\tau$), and the emitted ($L^{\mathrm{up}}$) and absorbed ($L^{\mathrm{down}}$) radiance. Their relationship to the radiance measured by the sensor is described by the radiative transfer equation, given by Reference [47]:

$$L_{\mathrm{sensor}}=[\epsilon B+(1-\epsilon )L^{\mathrm{down}}]\tau +L^{\mathrm{up}}$$

(1)

where $L_{\mathrm{sensor}}$ is the top-of-atmosphere (TOA) radiance measured by the sensor, $\epsilon$ is the surface emissivity, and B is the spectral radiance emitted by a black body with temperature T.

While MODIS has two thermal bands and its own algorithm to correct for these effects, Landsat TM and ETM+ sensors have only one thermal band, which makes this correction using its own band measurements unfeasible [20,26]. These correction parameters must be derived from additional data provided from other sources, and are most accurately calculated using radiative transfer codes, with in situ radiosonde data and weather models to provide the input data for the atmospheric profiles [12,31]. When no local atmospheric data are available to apply this equation, as is the case of most lake applications, algorithms can be used. For sensors with a single thermal band, Jiménez-Muñoz and Sobrino [19] developed a single-channel algorithm to estimate surface temperatures, focusing on Landsat 5 TM, replacing the atmospheric parameters by functions of the total columnar atmospheric water vapour content and then adjusting second-degree polynomials to it, considering a number of atmospheric profiles and using the MODTRAN radiative transfer code. Jiménez-Muñoz et al. [48] updated the algorithm to include Landsat 4 TM and Landsat 7 ETM+. In their algorithm, surface temperature is given by:

$$\mathrm{ST}=\gamma [{\epsilon }^{-1}({\psi }_1{L}_{\mathrm{sensor}}+{\psi }_2)+{\psi }_3-L_{\mathrm{sensor}}]+T_{\mathrm{sensor}}$$

(2)

where $T_{\mathrm{sensor}}$ is the brightness temperature, ${\psi }_1$, ${\psi }_2$, and ${\psi }_3$ are atmospheric functions, and $\gamma$ is a function calculated by:

$$\gamma ={\left[\frac{c_2{L}_{\mathrm{sensor}}}{T_{\mathrm{sensor}}^2}\left(\frac{{\lambda }_{\mathrm{eff}}^4}{c_1}{L}_{\mathrm{sensor}}+{\lambda }_{\mathrm{eff}}^{-1}\right)\right]}^{-1}$$

(3)

where $c_1=1.191\times {10}^8$ W $\mathsf{\mu }$m${}^4/$m${}^2$ sr and $c_2=1.439\times {10}^4$ K$\mathsf{\mu }$m are the Planck’s constants, and ${\lambda }_{\mathrm{eff}}$ is the effective wavelength of the thermal band, which is 11.269 $\mathsf{\mu }$m for Landsat 7 ETM+ [49]. The atmospheric functions are related to atmospheric parameters as:

$$\begin{array}{c}{\psi }_1=\frac{1}{\tau }\hfill \\ {\psi }_2=-L^{\mathrm{down}}-\frac{L^{\mathrm{up}}}{\tau }\hfill \\ {\psi }_3=L^{\mathrm{down}}\hfill \end{array}$$

(4)

and calculated by:

$${\psi }_n=c_n{w}^2+b_{n}w+a_n$$

(5)

where $n=1,2,3$; w is the atmospheric water vapour content, and a, b, and c are constants given by Jiménez-Muñoz et al. [48], according to the chosen atmospheric profile database. In this study, we used two databases: SAFREE, which is a compilation of maritime profiles, and thus should be useful for coastal zones; and TIGR3, which has a wide range of atmospheric profiles, and showed better results in our preliminary tests.

For the atmospheric water vapour content, we used the MODIS product MOD07L2, which is daily and has a spatial resolution of 5 km, considering that this parameter does not have a large spatial variability in comparison to its spatial resolution; thus, it could be used here [31,50], and because it showed good results in the study of atmospheric parameter sources by Jiménez-Muñoz et al. [31].

However, as these constants are generalized for a number of atmospheric profiles, it is expected that estimates are not quite precise, especially when atmospheric water vapour content is high [19]. Because of this, in this study we also used NASA’s Atmospheric Correction Parameter Calculator (AtmCorr) [51,52], a freely available online platform (available on https://atmcorr.gsfc.nasa.gov/) that provides these parameters by applying MODTRAN 4.0 to atmospheric profiles modeled by the National Centers for Environmental Prediction (NCEP/NOAA), generated every 6 h for each node of meridians and parallels (1${}^{°}$ × 1${}^{°}$ grid), interpolating the terms of atmospheric correction for Landsat thermal band for the time, date and location desired. With these data, the radiance measured by the thermal band can be corrected from the atmospheric effects ($L_{\mathrm{corr}}$), and the estimated ST is given by brightness temperature. At-sensor radiance was corrected using the following equation [47]:

$$L_{\mathrm{corr}}=\frac{L_{\mathrm{sensor}}-L^{\mathrm{up}}-\tau L^{\mathrm{down}}(1-\epsilon )}{\epsilon \tau }$$

(6)

Emissivity strongly influences the estimated temperatures but, as it is difficult to estimate [53], a fixed value is commonly employed (e.g., References [15,24,28]). Here, water emissivity was fixed at 0.99, as was done in Díaz-Delgado et al. [50], Okwen et al. [54], and Deng and Wu [55], for example, and found in MODIS emissivity product MOD11 (0.992 for band 31 and 0.988 for band 32) for Lake Mangueira. Emissivity is discussed further in Section 2.4.

Altogether, three different methods were used to estimate Lake Mangueira’s SWT using Landsat 7 ETM+ images. The nomenclature for each method is described in

Table 1. Nomenclature of methods used for the estimation of LSWT using Landsat 7 ETM+ imagery. Note: AtmCorr, NASA’s Atmospheric Correction Parameter Calculator.

Method	Description
AtmCorr	Radiative transfer equation applied with atmospheric correction parameters from AtmCorr
SC-SAFREE	Single-channel algorithm using SAFREE database and MOD07L2’s water vapour
SC-TIGR	Single-channel algorithm using TIGR3 database and MOD07L2’s water vapour

.

2.3.2. Processing Landsat Images

Only Landsat 7 ETM+ T1 images (Tier 1 images, the highest-quality data available, with geometric and radiometric correction) were used. We used the radiance measured by Landsat 7 $B{6}_2$ (high gain) thermal band, which had higher radiometric precision than $B{6}_1$ (low gain), although it worked in a narrower temperature range [56]. First, images were cropped to a polygon of 1 × 1 km around each station, in order to provide more pixels for the estimation, given the limiting factor of cloud cover in the region and the known gaps in the images due to a malfunction of the satellite’s scan-line corrector. Low-quality pixels (due to cloud cover and shading, for example) were removed using the pixel quality band (BQA) with a 30 m buffer around them. Then, the DN (Digital Number) values in the image were transformed into spectral radiance (in W/m${}^2$ sr $\mathsf{\mu }$m), using values in the metadata file, and then to brightness temperature, in K, by inverting Planck’s law, according to Chander et al. [56]:

$$T=\frac{K2}{\ln \left(\frac{K1}{L_{\mathrm{sensor}}}+1\right)}$$

(7)

where K1 and K2 are calibration constants, and for Landsat 7 ETM+ they are 666.09 W/(m${}^2$ sr $\mathsf{\mu }$m) and 1282.71 K, respectively.

The LSWT was then calculated for each Landsat image and for each polygon; then, the calculated temperature for each station was assumed to be the mean temperature of the pixels in each polygon. The algorithm for all these processes was written in R [57], with the help of packages maptools [58], sp [59], raster [60], rgdal [61], and rgeos [62].

2.4. Sensitivity Analysis

Besides the spatial, spectral, and radiometric resolution of sensors, the accuracy of the temperature estimate depends mainly on two variables: the emissivity and the atmospheric correction parameters [28,30]. Therefore, an assessment of the sensitivity of the temperature estimates to these variables can improve the understanding of the errors resulting from remote sensing-derived temperatures [33,63]. Because in MODIS products they are calculated using internal algorithms developed by the MODIS team, we only analyzed the sensitivity of the estimates for Landsat images.

2.4.1. Emissivity

Emissivity is difficult to estimate from remote-sensing data. A few methods to estimate surface emissivity have been developed [20,64,65], most of them fixing values for different land covers. MOD11, for example, provides a product of surface emissivity for each of its thermal bands, and uses fixed values of emissivity for water, which depend on the characteristics of the water-surface images in the red and near infrared bands. Water emissivity varies with temperature and dissolved materials, for example, but in the TIR region, it approaches a blackbody [66]. A range of emissivity values has been used in the literature (e.g., 0.9925 in Reference [28] and 0.98 in Reference [67]), but the most common values are 0.9885 [15,24,35], which is the mean of the mean water-emissivity values for MODIS bands 31 and 32 [68], and 0.99 [50,54,55], which is an approximate value derived by Masuda et al. [69] for distilled water for a range of wavelength and sensor viewing angles.

For sensitivity analysis, we estimated the LSWT using the RTE with AtmCorr parameters, varying emissivity from 0.97 to 1.0 by increments of 0.005.

2.4.2. Atmospheric Correction Parameters

There are a few sources of atmospheric parameters that can be used when no local data are available. Jiménez-Muñoz et al. [31] compared three different sources of atmospheric profiles, AtmCorr, MOD07L2, and MAPRI (developed by the authors from reanalysis products from the NCEP), with local radiosonde data for ST estimation. In their study, however, the authors used only one Landsat 5 TM image, limiting its conclusions. The authors found errors of up to 3.5% for $\tau$ and up to 20% for $L^{\mathrm{up}}$ and $L^{\mathrm{down}}$ when using MOD07L2 for ASTER thermal bands, and 8% for $\tau$ and 30% for $L^{\mathrm{up}}$ and $L^{\mathrm{down}}$ when using AtmCorr, and 10% for $\tau$ and 60% for $L^{\mathrm{up}}$ and $L^{\mathrm{down}}$ when using MAPRI for the Landsat 5 TM thermal band. Coll et al. [70] ran a similar analysis, finding RMSE of 0.9, 1.0, and 1.1 K when using MOD07, AtmCorr, and radiosonde data, respectivelly.

For sensitivity analysis, we only used the RTE with the atmospheric parameters given by AtmCorr. The Reanalysis product was not used due to its coarse spatial resolution (2.5${}^{°}$ × 2.5${}^{°}$). Since the results calculated by Jiménez-Muñoz et al. [31] were only for a single date, and are therefore difficult to extrapolate, here we considered a variation of −5% to +5% with increments of 2.5% for $\tau$ and −10% to +10% with increments of 5% for $L^{\mathrm{up}}$ and $L^{\mathrm{down}}$, which are more realistic in terms of the variations of these parameters. Although the variation was consistent for all parameters, that is, all parameters varied positively or negatively, we note that, as atmospheric water vapour increases, $\tau$ tends to decrease (the atmosphere interacts more with the radiance emitted and reflected by the surface) and $L^{\mathrm{up}}$ and $L^{\mathrm{down}}$ tend to increase (the atmosphere both emits and absorbs more radiance), which we would expect to compensate each other.

We also assessed the sensitivity of the single-channel algorithm to atmospheric water vapour content. Since AtmCorr also provides an estimation of atmospheric water vapour content, it was compared with MOD07L2. As it has been reported, the algorithm is not very sensitive to the water vapour content (according to Prats et al. [15], a 50% decrease in the value of w produced a decrease in the estimated temperature of 0.2–0.3 ${}^{°}$C), we varied the water vapour from −50% to +50% to assess its sensitivity.

3. Results

3.1. Temperature Pattern in the Lake

The largest temperature difference between the three stations (Figure 3) was observed on 15 March 2013 (day 74), during late summer, when the temperature at the Central station was 21.79 ${}^{°}$C, and at the Southern station was 18.83 ${}^{°}$C. However, most of the time, the difference in temperature between the three stations was ≤1.5 ${}^{°}$C.

The lake shows a clear seasonal pattern, with temperatures in the range of 20–25 ${}^{°}$C during summer, and reaching less than 8 ${}^{°}$C during winter. The Central and Southern stations showed the lowest temperatures in almost the same number of measurements each, while the Northern station showed the highest temperatures in most measurements, as previously reported [36,37].

3.2. Estimated Temperatures

All three methods using Landsat images showed good performance (Figure 4), with most estimates within a window of +2 ${}^{°}$C/−3 ${}^{°}$C of the measured temperatures. In most cases, LSWT was underestimated, which is supported by the negative bias found for all methods (

Table 2. Metrics calculated for each method of estimating Lake Mangueira’s SWT using Landsat 7 ETM+ imagery and MODIS products.

Method	Absolute Error (${}^{°}$C)	Bias (${}^{°}$C)	RMSE (${}^{°}$C)	${\mathit{r}}^2$
AtmCorr	0.880	−0.465	1.07	0.964
SC-SAFREE	1.19	−0.990	1.50	0.936
SC-TIGR	1.00	−0.560	1.32	0.938
MOD11	0.825	−0.523	1.05	0.962
MOD28	0.989	0.329	1.25	0.906

). As ANOVA showed significant statistical difference between all methods, we also ran Tukey’s post hoc test to pairwise analyze the differences between the methods (

Table 3. p-value of the Tukey’s test between the residuals of each method.

Method	AtmCorr	SC-TIGR	SC-SAFREE	MOD11	MOD28
AtmCorr	NA	0.967	<0.05	0.990	<0.05
SC-TIGR	0.967	NA	0.198	0.999	<0.05
SC-SAFREE	<0.05	0.198	NA	0.113	<0.05
MOD11	0.990	0.999	0.113	NA	<0.05
MOD28	<0.05	<0.05	<0.05	<0.05	NA

). AtmCorr estimated LSWT with the highest accuracy, showing the best values for all calculated metrics, with values of RMSE and $r^2$ of 1.07 ${}^{°}$C and 0.964, respectively. Among the estimates using the single-channel algorithm, the TIGR3 database functioned better than SAFREE for this region, reducing absolute error and RMSE by around 0.2 ${}^{°}$C, and bias by almost 0.5 ${}^{°}$C, although no significant statistical difference between the two methods was found (p-value = 0.198, Table 3). Tukey’s post hoc test also showed that AtmCorr and TIGR3 residuals (Figure 5) had no significant statistical difference (p-value = 0.967).

Both MODIS products also showed good performance (Figure 4), but MOD11 performed better than MOD28, with values of RMSE and $r^2$ of 1.05 ${}^{°}$C and 0.962, respectively (Table 2), and statistically similar to the residual results found for AtmCorr (Table 3). Despite showing the smallest bias (in terms of magnitude) among all estimators, MOD28 showed higher inconsistencies, in some cases overestimating the temperatures by more than 3 ${}^{°}$C (Figure 5), with a significant statistical difference of residuals from all methods (Table 3). It also consistently overestimated temperatures, with a positive bias of 0.329 ${}^{°}$C, which was not found for the other estimators. This is also reflected in the low values of the correlation coefficients between the MOD28 and Landsat methods (

Table 4. Correlation coefficient between the residuals of each method.

Method	AtmCorr	SC-TIGR	SC-SAFREE	MOD11	MOD28
AtmCorr	NA	0.768	0.730	0.606	0.355
SC-TIGR	0.768	NA	0.962	0.670	0.324
SC-SAFREE	0.730	0.962	NA	0.704	0.425
MOD11	0.670	0.606	0.496	NA	0.644
MOD28	0.324	0.355	0.425	0.644	NA

).

The accuracy of the estimates using Landsat images is directly affected by atmospheric water vapour content (Figure 6), as seen in the regression lines (

Table 5. p-value of the linear regression between the residuals and w for each method.

Method	Slope	Intercept	p-Value
AtmCorr	0.587	−1.41	<0.01
SC-TIGR	0.620	−1.53	<0.01
MOD11	−0.041	−0.413	0.813
MOD28	−0.522	1.30	0.0513

), while MODIS estimations do not seem to follow this pattern (p-value > 0.05).

3.3. Sensitivity Analysis

Sensitivity analysis (Figure 7) showed that, for emissivity, there is a considerable variation in error (

Table 6. Metrics calculated for sensitivity analysis of the temperature estimates to emissivity.

Emissivity	Bias (${}^{°}$C)	RMSE (${}^{°}$C)	Maximum Error (${}^{°}$C)	${\mathit{r}}^2$
1.00	−0.916	1.34	−3.34	0.965
0.995	−0.691	1.19	−3.08	0.965
0.990 *	−0.465	1.07	−2.82	0.965
0.985	−0.237	0.981	−2.55	0.965
0.980	−0.0068	0.943	−2.28	0.965
0.975	0.225	0.962	2.17	0.965
0.970	0.459	1.04	2.40	0.965

* value of emissivity used in this study.

), with RMSE values varying by almost 12% for emissivities in the 0.97–0.99 range. Bias and the absolute maximum error showed an almost linear increase, with minimum values found for emissivities of 0.975 and 0.980, respectively. In the case of the atmospheric correction parameters, the temperature was very sensitive to changes (when applying the radiative transfer equation): even a small change of 5% in these parameters (2.5% for atmospheric transmissivity) increased the RMSE and the maximum error by 2 ${}^{°}$C (

Table 7. Metrics calculated for sensitivity analysis of the temperature estimates to atmospheric parameters.

Atmospheric Parameters	Bias (${}^{°}$C)	RMSE (${}^{°}$C)	Maximum Error (${}^{°}$C)	${\mathit{r}}^2$
AtmCorr	−0.465	1.07	−2.82	0.965
−10%	−5.39	5.51	−8.28	0.960
−5%	−3.27	3.40	−4.85	0.945
+5%	2.44	2.75	4.64	0.963
+10%	4.74	5.03	8.12	0.957

). The increases in errors were sharper when the values of the atmospheric parameters are decreased.

On the other hand, the single-channel algorithm is only slightly sensitive to the atmospheric parameters (simplified as atmospheric water vapour content). The algorithm showed a very similar performance when using MOD07L2 or AtmCorr as the source of water vapour (

Table 8. Metrics calculated for sensitivity analysis of the temperature estimates to atmospheric water vapour content.

Water Vapour Source	Bias (${}^{°}$C)	RMSE (${}^{°}$C)	Maximum Error (${}^{°}$C)	${\mathit{r}}^2$
AtmCorr	−0.483	1.26	−2.82	0.944
MOD07L2 *	−0.594	1.31	−2.99	0.938
MOD07L2 −50%	−0.937	1.43	−3.12	0.944
MOD07L2 +50%	−0.195	1.80	5.46	0.912

* source of water vapour content used in this study.

), despite the significant changes in the values given by each source, with values of RMSE of 0.39 g/cm${}^2$ (23.5%), maximum difference of −1.26 g/cm${}^2$ (64%), and $r^2$ of 0.675. An even more pronounced difference in water vapour content, increasing or decreasing its value by 50%, showed a variation in the results considerably smaller than in the case of the radiative transfer equation, considering its order of magnitude. The same pattern was observed when varying AtmCorr’s water vapour content by −50% and +50%, with quite similar values of bias, RMSE, and maximum absolute error.

4. Discussion

4.1. Comparison of Landsat 7 ETM+ and MODIS Performance

When using TIGR3 database constants, AtmCorr and the single-channel algorithm both showed good performance (with no statistical difference between their residual data), corroborating with similar studies found in the literature [15,24,33]. TIGR3 outperforming SAFREE (Table 3) for coastal-lake SWT estimation was previously reported [35], suggesting that the SAFREE database should not be used in such studies.

The main limiting factor for LSWT estimates using Landsat imagery is the correction of the atmospheric effects. Using commercial radiative transfer codes to estimate these parameters and applying the RTE, Schneider and Mauser [28] found an RMSE of 0.73 ${}^{°}$C when applying Lowtran with Landsat 5 images, and Allan et al. [12] found an RMSE between 0.48 ${}^{°}$C and 0.94 ${}^{°}$C when applying MODTRAN with Landsat 7 ETM+. These studies showed how specific, more accurate atmospheric correction parameters can improve LSWT estimation. While AtmCorr is the best alternative (to the authors’ knowledge) for freely obtaining these parameters, it provides only a point estimate, which is interpolated from the closest nodes of parallels and meridians (1${}^{°}$ × 1${}^{°}$ grid). For large lakes, this can be a problem when estimating the temperature for its whole surface and, considering the operational limitation of AtmCorr, MOD07L2 may be the best option when using Landsat imagery in these cases, since it provides water vapour content with spatial variability. The main issue with the estimates using MOD07L2 with the single-channel algorithm in our study is the lower accuracy in some SWT estimates, reaching out to a $\Delta$T of ±3 ${}^{°}$C in some cases. Many of these low-quality estimates are consistent with the results found by Jiménez-Muñoz et al. [48], in which the error of the ST estimates increased considerably for atmospheric water vapour content higher than 3.0 g/cm${}^2$ (Figure 6).

The performance found for MODIS also compared well with the literature, with both products showing performance similar to what was found in most studies [6,10,11,14,41]. Reinart and Reinhold [23] used MOD28 to estimate the LSWT of two large lakes in Sweden, and found a mean absolute difference of 0.41 ${}^{°}$C and $r^2$ of 0.993 for both lakes. In spite of the good performance found by the authors, LSWT estimates were compared with measured bulk temperatures, and so this may not be the case for the ST, since many of the data were from the summer months, when both lakes may have been thermally stratified.

Our findings showed higher similarity between LSWT estimations based on MOD11, AtmCorr, and SC-TIGR, different from LSWT estimations based on MOD28. Despite the similar performance shown by MOD28 and SC-TIGR (in terms of absolute error and RMSE), MOD28 showed more erratic behavior, with a smaller correlation coefficient and different linear regression (Figure 4). Besides that, MOD28 showed a positive bias, while all other estimators showed a negative bias. Additionally, good agreement between the estimates using Landsat and MOD11 (Table 3 and Table 4) indicates that they are consistent estimators over this lake, while MOD28 is not, as there is a high correlation between the residuals for the Landsat methods and MOD11, but not between the former and MOD28. Even though high correlation was found between MOD11 and MOD28, it may be due to the fact that both MODIS products similarly overestimated small-to-medium temperatures and underestimated higher temperatures (despite MOD28 showing a positive bias), which was not observed in the Landsat estimates.

The lower performance obtained for MOD28 may be explained by its low spatial resolution of 4 km, which might be too coarse for Lake Mangueira due to its narrowness (5 km wide over the Northern station, and less than 4 km wide over the Southern station). In fact, errors in the estimates were considerably higher for the Southern station than for the Central and Northern stations, pointing to this. Considering that temperature varies more widely near the shores [6,8,24], it can be inferred that this product is not recommended for applications to lakes with narrow widths (in comparison to its spatial resolution), as is the case for Lake Mangueira.

4.2. Sensitivity to Parameters and Source of Errors

As seen in the literature, error in lake-water temperatures derived from remote-sensing imagery can occur mainly due to poor correction of atmospheric effects [30,71], water-surface effects [72], and due to the misrepresentation of water emissivity [30,45].

There is difficulty associated with estimating and removing atmospheric effects, which can promote substantial errors [31,70]. In fact, estimates tend to have lower accuracy when the water vapour content in the atmosphere is higher, since this increases the uncertainty regarding the effect of the atmosphere on the radiance measured by the sensor and, thus, in the correction parameters. Jiménez-Muñoz et al. [48] showed that the error in ST estimates considerably increased for atmospheric water vapour content higher than 3.0 g/cm${}^2$. Here, Landsat-derived estimates using AtmCorr showed a similar (although less pronounced) pattern. As atmospheric water vapour content (from the MOD07L2 water vapour product) increases, so does the error in the estimates of both SC-TIGR and AtmCorr, although this is not the case for MOD11 and MOD28 (Figure 6, Table 5). Although in other studies the accuracy of MODIS products decreased with higher temperatures (e.g., Reference [11]), which is related to higher water vapour content, this was not evident here (p-value > 0.05 for both products, Table 5). This may be due to the refined quality of the MODIS algorithm to remove these effects using its own band measurements [31], and to directly apply it in its ST products.

Temperature estimates using the RTE are very sensitive to atmospheric parameters (Figure 7), particularly to atmospheric transmissivity. For instance, a small change of 2.5% in atmospheric-transmissivity value implied an increase of about 2 ${}^{°}$C in the RMSE. The increase in errors is more pronounced when decreasing the values of the atmospheric parameters, showing that $\tau$ (which decreases as atmospheric water vapour increases) has a more important effect on $L_{\mathrm{corr}}$ than $L^{\mathrm{up}}$ and $L^{\mathrm{down}}$ (both increase as atmospheric water vapour increases), which is expected considering Equation (6). In this case, negative bias also contributes to this sharp increase in errors, as it increased in magnitude when the atmospheric parameters negatively varied. This result highlights the importance of an accurate, precise source of estimation of the influence of the atmosphere on the radiance measured by the sensor, as has been shown in the literature by good results in estimation of LSWT when using in situ radiosonde data and radiative transfer models (e.g., Reference [12]).

The single-channel algorithm is much less sensitive to changes in atmospheric parameters, i.e., atmospheric water vapour content. Prats et al. [15] reported that a 50% decrease in water vapour content produced a decrease of 0.2–0.3 ${}^{°}$C in the estimated temperatures; in our sensitivity analysis, this decrease produced a (negative) increase of 0.3 ${}^{°}$C in the bias, and only marginally increased the RMSE (0.07 ${}^{°}$C). On the other hand, for increased w values, the RMSE and maximum error increased by 0.5 ${}^{°}$C and 1 ${}^{°}$C, respectively, which is consistent with the report of Jiménez-Muñoz et al. [48]. Considering the small size (in comparison to emissivity and the radiance measured by the sensor) of the contribution of the water vapour content to the temperature estimates in this algorithm, MOD07L2 can be used for retrievals of LSWT using Landsat images, since data are measured daily and are spatially distributed, and this product is ready to use.

In our study, only a few estimates were made when water vapour content was high (>3 g/cm${}^2$, Figure 6), which limits further understanding of the performance of these methods under such conditions. However, as seen in these few estimates and in the literature (e.g., References [11,48]), studies estimating LSWT in areas with high atmospheric water vapour, such as tropical areas, should keep in mind that estimation is bound to produce higher errors, especially when using Landsat imagery (Table 8). For the Landsat methods, other sources of atmospheric parameters, such as reanalysis products [15,31], should also be further studied in order to be used with the radiative transfer equation in an effort to increase the accuracy of the estimates.

The cool bias found in our study was also found in most LSWT studies (e.g., References [6,10,11,15,73]). This can be a result of the cool-skin effect, when the few millimetres at the top of the water column are colder than the water underneath. Usually, this is not measured by in situ sensors but affects estimates made using remote sensing [74]. The cool-skin effect reduces water temperature over a range of 0.1–0.6 ${}^{°}$C, occurring as a result of heat loss from water to air due to evaporation, and is affected mainly by wind speed [72,75]. A warm-layer effect can also be observed during calm sunny days in deeper lakes, where thermal stratification can occur, and can increase the ST up to several degrees [76], but as observed here and in other studies (e.g., Reference [36]), this is not the case for the shallow Lake Mangueira.

While this effect is different for each lake and depends on in situ lake-temperature measurements to be understood, techniques are available to address the cool-skin effect in sea water, which assume this to be a function of wind speed [75]. Although we did not attempt to assess and remove the cool-skin effect in the estimations here, our study shows that it is important in LSWT studies using remote sensing. Prats et al. [15] estimated the LSWT for French water bodies for a long time series using Landsat 5 TM and 7 ETM+ images, applying an adjusted sea-water algorithm to remove the effects of the cool-skin and warm-layer effects, and testing these approaches on a reservoir with in situ data to validate it. The estimated cool-skin effect reduced the LSWT by about 0.3 to 0.6 ${}^{°}$C for most of the data, and even though the authors concluded that the algorithm required further adaptations for use on continental waters, it demonstrated how techniques to remove these effects can be useful to improve the accuracy of remote sensing-derived SWT in lake studies.

In the case of emissivity, considering that MODIS products use a fixed value of emissivity for water, as was also used in our study for Landsat 7 ETM+ images (as well as in most studies of water-temperature estimations using remote sensing, e.g., References [15,24]), it can also be a source of considerable error. Emissivity varies with several parameters, such as salinity and turbidity, which can result in errors of estimation ranging from 0.2 to 0.5 ${}^{°}$C [45]. Hulley et al. [32] used an algorithm to improve the estimation of water emissivity and the coefficients for the split-window algorithm for MOD11, in order to improve its estimation of inland SWT. They found differences of emissivity in up to 0.01 for MODIS, which resulted in errors of 0.3–0.7 ${}^{°}$C.

As shown in the sensitivity analysis, emissivity has a significant impact on the estimated temperatures, and showed an optimal value of 0.98, reducing the RMSE and maximum absolute error by 10% and 20%, respectively, compared to our standard value of 0.99. This, however, can be a result of the balancing of the cool-skin effect: as the latter reduces water temperature at the surface, a lower value of emissivity (i.e., further from its physical, more-likely value) makes the surface appear warmer than it is, partially compensating for the cool-skin effect. Although this 0.98 emissivity value produced the best results, we assumed (taking into account that it was not measured) that the closest value, from 0.97 to 1, to its real value is 0.99, which was calculated by Masuda et al. [69] and also found in MODIS product MOD11 (0.992 for band 31 and 0.988 for band 32) for Lake Mangueira. Due to the limitation of water emissivity estimation [66], most approaches still use a fixed, tabulated value of emissivity, found in reflectivity and emissivity data available in spectral libraries such as MODIS’ [77]. As emissivity and the removal of the cool-skin effect are limiting factors to be improved in LSWT estimations, the use of lower emissivity values can be tested in future studies to reduce errors resulting from the latter in an effort to improve the accuracy of LSWT estimates.

5. Conclusions

In this study, we compared different methods for estimating lake-surface-water temperatures using Landsat 7 ETM+ and MODIS products. Considering the results from our study and the literature, MOD11 is the recommended product to be used in LSWT studies due to its accuracy (values of bias, RMSE, and $r^2$ of −0.523 ${}^{°}$C, 1.05 ${}^{°}$C and 0.962, respectively), temporal resolution, and ease of use, and should be preferred over MOD28, which gave less accurate results. Landsat 7 ETM+ should be used in small lakes, due to its higher spatial resolution, and we recommend the application of the radiative transfer equation with atmospheric correction parameters when possible. Despite its high sensitivity to these parameters, as shown in our work, which can result in gross errors in temperature estimates, this equation was successfully applied in Lake Mangueira with AtmCorr, and showed quite similar performance to MOD11 (values of bias, RMSE, and $r^2$ of −0.465 ${}^{°}$C, 1.07 ${}^{°}$C, and 0.964, respectively). The limitations of this application are difficulties in operating this platform, and that it only provides correction parameters for a point (i.e., it is not spatially distributed). Although not as accurate, the single-channel algorithm can also be used with Landsat imagery using the TIGR3 database, and most sources of atmospheric water-vapour content, since it is not very sensitive to this parameter. We recommend MOD07L2 to be used with this algorithm, because the data are measured daily and spatially distributed, and it is ready to use. However, users should take into account that the accuracy of this algorithm is considerably reduced in cases when atmospheric water vapour content is high.

We recommend that similar tests be done for Landsat imagery (especially on RTE sensitivity to different sources of atmospheric parameters) on other lakes with different geographical and morphological characteristics to support our findings. We also recommend the further development and adaptation of techniques to assess and remove the cool-skin effect in SWT estimations in inland waters to reduce the cool bias seen in most studies on this subject, and tests of different emissivity values, to verify if they can compensate, to some degree, for the cool-skin effect, in an effort to improve the accuracy of these estimates.