A Comparison of Three Trapezoid Models Using Optical and Thermal Satellite Imagery for Water Table Depth Monitoring in Estonian Bogs

This study explored the potential of optical and thermal satellite imagery to monitor temporal and spatial changes in the position of the water table depth (WTD) in the peat layer of northern bogs. We evaluated three different trapezoid models that are proposed in the literature for soil moisture monitoring in regions with mineral soils. Due to the tight capillary connection between water table and surface soil moisture, we hypothesized that the soil moisture indices retrieved from these models would be correlated with WTD measured in situ. Two trapezoid models were based on optical and thermal imagery, also known as Thermal-Optical TRApezoid Models (TOTRAM), and one was based on optical imagery alone, also known as the OPtical TRApezoid Model (OPTRAM). The models were applied to Landsat imagery from 2008 to 2019 and the derived soil moisture indices were compared with in-situ WTD from eight locations in two Estonian bogs. Our results show that only the OPTRAM index was significantly (p-value < 0.05) correlated in time with WTD (average Pearson correlation coefficient of 0.41 and 0.37, for original and anomaly time series, respectively), while the two tested TOTRAM indices were not. The highest temporal correlation coefficients (up to 0.8) were observed for OPTRAM over treeless parts of the bogs. An assessment of the spatial correlation between soil moisture indices and WTD indicated that all three models did not capture the spatial variation in water table depth. Instead, the spatial patterns of the indices were primarily attributable to vegetation patterns.


Introduction
Knowledge of soil moisture is important for a wide range of applications in hydrology, meteorology, climatology, and biogeochemistry [1][2][3]. Soil moisture is pivotal to understand dynamics of net ecosystem exchange and ecosystem feedbacks to global climate change [4]. This is especially the case for peatlands that are defined by a surface soil layer of peat that is at least 30 cm thick and that act as large terrestrial carbon pools [5][6][7]. The soil moisture in surface peat layers and the position of the water table in the peat layer relative to the ground surface, i.e., the water table depth (WTD), are important features of peatlands, as they control the interface between aerobic and anaerobic conditions. This, in turn, causes a shift in fluxes of CO 2 , N 2 O, and CH 4 [8][9][10].
Mineral and peat soils differ from each other primarily in the proportion of solid inorganic and organic particles, leading to differences in their hydraulic properties [11,12]. Within a region with mineral soils, the spatial differences in water retention curves can be attributed to a mineral texture (from sand to clay) [13,14], while, in peat soils, differences can be attributed to the vegetation and the related plant debris, and its degree of decomposition over time [15]. The temporary unsaturated surface layer of the peat profile (acrotelm) is characterized by a high organic matter content that is poorly decomposed. This structure of the acrotelm leads to a high water storage capacity and high hydraulic conductivity [16][17][18].
Shallow WTD conditions and the typical hydraulic properties of the acrotelm cause a tight connection between WTD and soil moisture content in peatlands [7,[19][20][21][22][23][24]. In several studies, the tight connection is illustrated by scatter plots of temporal data of soil moisture and WTD at individual locations. In these plots, points typically align in narrow bands, which indicates that a new soil water equilibrium is quickly established after a precipitation event. Figure 1 illustrates how this tight connection can be used to derive information about WTD dynamics in time from the surface characteristics sensed with optical and thermal satellite imagery. Note that the soil moisture in peatlands is highly variable at the scale of a few meters due to the presence of small-scale elevation differences also known as hummock and hollow microtopography. The typical optical and thermal satellite imagery has a spatial resolution of tens of meters and integrates over this small-scale variability.
Remote Sens. 2020, 12, x FOR PEER REVIEW 2 of 24 Mineral and peat soils differ from each other primarily in the proportion of solid inorganic and organic particles, leading to differences in their hydraulic properties [11,12]. Within a region with mineral soils, the spatial differences in water retention curves can be attributed to a mineral texture (from sand to clay) [13,14], while, in peat soils, differences can be attributed to the vegetation and the related plant debris, and its degree of decomposition over time [15]. The temporary unsaturated surface layer of the peat profile (acrotelm) is characterized by a high organic matter content that is poorly decomposed. This structure of the acrotelm leads to a high water storage capacity and high hydraulic conductivity [16][17][18].
Shallow WTD conditions and the typical hydraulic properties of the acrotelm cause a tight connection between WTD and soil moisture content in peatlands [7,[19][20][21][22][23][24]. In several studies, the tight connection is illustrated by scatter plots of temporal data of soil moisture and WTD at individual locations. In these plots, points typically align in narrow bands, which indicates that a new soil water equilibrium is quickly established after a precipitation event. Figure 1 illustrates how this tight connection can be used to derive information about WTD dynamics in time from the surface characteristics sensed with optical and thermal satellite imagery. Note that the soil moisture in peatlands is highly variable at the scale of a few meters due to the presence of small-scale elevation differences also known as hummock and hollow microtopography. The typical optical and thermal satellite imagery has a spatial resolution of tens of meters and integrates over this small-scale variability. Different types of satellite observations have been used to estimate soil moisture. The large difference between the dielectric constants of water (~80) and soil particles (~4) permits the monitoring of soil moisture with passive and active microwave remote sensing [25][26][27][28]. Microwave observations have the advantage that they are directly sensitive to soil moisture, although the sensing depth is typically limited to the top few centimeters of soils. The longer wavelength (L band) of the passive microwave satellite missions that are currently dedicated to soil moisture monitoring are less sensitive to vegetation than the operational active microwave missions. However, the passive observations have a coarse resolution (~40 km) [25,26], whereas some active missions allow high spatial resolution mapping (<100 m) [29,30]. A few studies have analyzed and used passive [31][32][33] and active [34,35] microwave measurements for the purpose of soil moisture monitoring in peatlands. One of the major constraints of applicability of active microwave observations over peatlands is a dependency of radar signal on land surface properties [36]. The microwave backscatter signal over natural peatlands showed a low sensitivity to WTD dynamics, which might be related to strong confounding effects that occur due to complex scattering in vegetation and the fibric surface peat layer [37][38][39][40].
Soil moisture can also be estimated using optical and thermal imagery. One of the most wellknown approaches for soil moisture estimation is the "trapezoid", or "triangle", model concept [41,42]. The concept takes advantage of the existence of a relationship between optical and thermal Figure 1. Sketch illustrating the concept of the link between water table depth (WTD) and remotely sensed parameters via the capillary connection between WTD, moisture in the soil, and vegetation. The remotely sensed parameters used in this study: normalized difference vegetation index (NDVI), fractional vegetation cover (FVC), land surface temperature (LST), and shortwave infrared transformed reflectance (STR). Different types of satellite observations have been used to estimate soil moisture. The large difference between the dielectric constants of water (~80) and soil particles (~4) permits the monitoring of soil moisture with passive and active microwave remote sensing [25][26][27][28]. Microwave observations have the advantage that they are directly sensitive to soil moisture, although the sensing depth is typically limited to the top few centimeters of soils. The longer wavelength (L band) of the passive microwave satellite missions that are currently dedicated to soil moisture monitoring are less sensitive to vegetation than the operational active microwave missions. However, the passive observations have a coarse resolution (~40 km) [25,26], whereas some active missions allow high spatial resolution mapping (<100 m) [29,30]. A few studies have analyzed and used passive [31][32][33] and active [34,35] microwave measurements for the purpose of soil moisture monitoring in peatlands. One of the major constraints of applicability of active microwave observations over peatlands is a dependency of radar signal on land surface properties [36]. The microwave backscatter signal over natural peatlands showed a low sensitivity to WTD dynamics, which might be related to strong confounding effects that occur due to complex scattering in vegetation and the fibric surface peat layer [37][38][39][40].
Soil moisture can also be estimated using optical and thermal imagery. One of the most well-known approaches for soil moisture estimation is the "trapezoid", or "triangle", model concept [41,42]. The concept takes advantage of the existence of a relationship between optical and thermal land surface properties and soil moisture [43,44]. Due to the high sensitivity to vegetation properties, such as the moisture status in the plants, it is proposed that derived moisture estimates also reflect the root-zone moisture [44].
Two different types of the trapezoid concept exist and they mainly differ in the used moisture-sensitive signal. One of the types uses the land surface temperature (LST) and the physical principle of evaporative cooling (wetter = cooler); the other one uses the strong absorption of shortwave infrared (SWIR) radiation of water and calculates the shortwave infrared transformed reflectance (STR) from SWIR [44]. In both types, the moisture-sensitive signal is combined with a vegetation index (VI), e.g., the normalized difference vegetation index (NDVI) or the fractional vegetation cover (FVC), to construct a trapezoid space [44][45][46][47][48][49][50]. The lowest LST or highest STR along the vegetation cover gradient defines the wettest soil moisture, whereas the highest LST or lowest STR defines the driest soil moisture in the landscape. The trapezoid space is then used to linearly scale soil moisture. The type of trapezoid model that utilizes LST is further called the Thermal-Optical TRApezoid Model (TOTRAM) [42,[51][52][53][54], and the type of trapezoid model that utilizes STR is further called the OPtical TRApezoid Model (OPTRAM) [44].
There exist several different algorithms of the TOTRAM concept. One type of algorithm derives the dry edge directly from the LST over the driest areas in the scene. In another type of algorithms, additional weather data are used to theoretically estimate the temperature of the driest land fraction, which then defines the dry edge of the trapezoid space [55,56].
Previous studies evaluated the applicability of TOTRAM and OPTRAM mainly in mineral soils [44,48,[57][58][59]. Very little attention has been paid to the models' applicability in organic soils [60,61]. In this study, we explored the feasibility of applying the TOTRAM and OPTRAM soil moisture indices for the WTD monitoring in two ombrotrophic peatlands, i.e., bogs, in Estonia. The specific objectives of this study were: (i) to estimate the soil moisture indices based on two different TOTRAM and one OPTRAM concepts proposed in the literature, using Landsat 5, 7, and 8 observations for the growing seasons of 2008-2019; and (ii) to evaluate the temporal and spatial correlation of the indices with in-situ measured WTD of the same period.

TOTRAM: Thermal-Optical Trapezoid Model
TOTRAM is one of the most widely used approaches to estimate surface moisture from remote sensing data over mineral soils [42,53,54,62,63]. Initially, this method was developed to interpret the trapezoidal-shaped distribution of pixel values in the space defined by LST and a vegetation index (VI), such as FVC, NDVI, or leaf area index (LAI) [47,48,64]. The highest and lowest values of LST along the vegetation cover gradient represent the so-called dry and wet edges ( Figure 2). The wet and dry edges are isopleths and indicate the uniform soil moisture availability for two extreme surfaces: wet and dry [42,46]. Evapotranspiration on the wet edge is maximum and not limited by water availability, whereas evapotranspiration at the dry edge is minimum because of extreme water limitation [42]. In addition to the wet and dry edges, more isopleths could be drawn indicating transitional moisture condition [65]. Consequently, the TOTRAM soil moisture index, W TOTRAM , of each pixel i is calculated on the basis of its position relative to the wet and dry edges within LST-VI space at a single time step or scene, The accuracy of the TOTRAM depends on several criteria: • Negligible spatial variability in weather conditions. The study area should be limited and have a minimal topographic variation to guarantee that the location of the isopleths within LST-VI space is determined by the water availability and not by the difference in atmospheric conditions. For this reason, TOTRAM also demands an individual parameterization of the trapezoid space for each observation scene.

•
Negative temporal correlation between LST and VI at a given location. The temporal applicability of the TOTRAM concept requires that LST decreases not only in space but also in time with increasing VI [66,67]. We applied two different TOTRAM scenarios from the literature that differ in the treatment of these application criteria in the context of the dry edge determination.

Scenario 1: Observed dry edge
In Scenario 1, the trapezoid was solely based on observed pixel values. This scenario requires a sufficient range in vegetation and moisture conditions. The study area should be large enough to include sufficient land surface variability for a precise determination of the dry and wet edges [46]. The area size needs be balanced against the TOTRAM requirement of homogenous weather conditions. The dry and wet edges were determined from the observed highest and lowest LST values along the FVC gradient following the algorithm in [41]. This algorithm splits the LST pixels based on FVC values in LST-FVC space into 20 intervals and, after that, each interval into 5 subintervals. For each interval, the minimum and maximum LST values are calculated. Thereafter, the wet and dry edges are modeled as the linear fit to the mentioned minimum and maximum LST values. Scenario 2: Modeled dry edge During wet periods following widespread precipitation, the dry edge is difficult to define based on observed pixel values alone. To cope with this limitation, it has been proposed to derive the dry edge theoretically, i.e. based on models [55,68,69]. In Scenario 2, the modeled dry edge was LST is the land surface temperature, FVC is fractional vegetation cover, NDVI is the normalized difference vegetation index, and STR is the shortwave infrared transformed reflectance. For TOTRAM, both the observed (Scenario 1) and the modeled (Scenario 2) dry edges are presented. The dry edges are indicated by points LSTs max and LSTc max for TOTRAM Scenario 1, LSTs max " and LSTc max " for TOTRAM Scenario 2, and STRs min and STRc min for OPTRAM. The wet edges are indicated by points LSTs min and LSTc min for both TOTRAM scenarios, and STRs max and STRc max for OPTRAM. The color gradient shows the soil moisture availability from blue (wet edge) to red (dry edge). Point i is a surface with LST i , FVC i , STR i , and NDVI i . For i within LST-FVC space, the temperature of wet edge is LST min,i , observed dry edge is LST max,i , and modeled dry edge is LST max,i ". For i within STR-NDVI space, the STR value for the wet and dry edge are STR max,i and STR min,i , respectively.
The accuracy of the TOTRAM depends on several criteria:

•
Negligible spatial variability in weather conditions. The study area should be limited and have a minimal topographic variation to guarantee that the location of the isopleths within LST-VI space is determined by the water availability and not by the difference in atmospheric conditions. For this reason, TOTRAM also demands an individual parameterization of the trapezoid space for each observation scene.

•
Negative temporal correlation between LST and VI at a given location. The temporal applicability of the TOTRAM concept requires that LST decreases not only in space but also in time with increasing VI [66,67].
We applied two different TOTRAM scenarios from the literature that differ in the treatment of these application criteria in the context of the dry edge determination.

Scenario 1: Observed Dry Edge
In Scenario 1, the trapezoid was solely based on observed pixel values. This scenario requires a sufficient range in vegetation and moisture conditions. The study area should be large enough to include sufficient land surface variability for a precise determination of the dry and wet edges [46]. The area size needs be balanced against the TOTRAM requirement of homogenous weather conditions. The dry and wet edges were determined from the observed highest and lowest LST values along the FVC gradient following the algorithm in [41]. This algorithm splits the LST pixels based on FVC values in LST-FVC space into 20 intervals and, after that, each interval into 5 subintervals. For each interval, the minimum and maximum LST values are calculated. Thereafter, the wet and dry edges are modeled as the linear fit to the mentioned minimum and maximum LST values. During wet periods following widespread precipitation, the dry edge is difficult to define based on observed pixel values alone. To cope with this limitation, it has been proposed to derive the dry edge theoretically, i.e., based on models [55,68,69]. In Scenario 2, the modeled dry edge was determined by two modeled temperatures: the temperature of the driest bare surface (LSTs max ") and the driest fully vegetated surface (LSTc max ") ( Figure 2a). Because the in-situ meteorological measurements needed for modeling the dry edge were only available starting from June 2009, Scenario 2 was applied for a shorter period (2009-2019). We followed the algorithm for retrieving LSTs max " and LSTc max " developed by Long et al. [55], which is based on solving for the system of physically-based non-linear equations in an iterative manner (Appendices A and B). Knowing these two extreme temperatures, it becomes possible to model the theoretical dry edge with linear regression approach, where (LSTc max "-LSTs max ") is a slope and LSTs max " is an intercept.

OPTRAM: Optical Trapezoid Model
Sadeghi et al. [44] introduced the concept of OPTRAM which uses STR as soil moisture sensitive parameter. Figure 2b illustrates the concept of OPTRAM. In OPTRAM, the trapezoid is formed by NDVI as a measure of vegetation cover and STR as a measure of moisture content [70][71][72]. In contrast to LST, STR is not assumed to depend on the specific weather conditions, and the trapezoid can therefore be determined based on observations across all observed scenes. To determine the wet and dry edges, firstly, we randomly sampled 10% (approximately four million in space and time) of the pixels from the total dataset. Secondly, the full range of sampled NDVI values was divided into 100 intervals. Within each interval, the median, standard deviation, maximum, and minimum STR values were obtained. The dry edge was distinct and clearly observable, and it was modeled as a linear fit to the minimum STR values. Due to several obvious outliers caused by open water surfaces and remaining shadow effects, the wet edge was less explicitly defined than the dry one and it was modeled as a linear fit to the "median + standard deviation" values of the 100 intervals.
As noticed by Sadeghi et al. [44], OPTRAM approach is sensitive to oversaturated pixels. Therefore, we excluded pixels with negative NDVI values which correspond to water, as suggested by Babaeian et al. [58]. Additionally, all the pixels located above the wet edge were excluded from further statistical analysis.
The OPTRAM soil moisture index, W OPTRAM , is then derived by a linear scaling between the maximum and minimum STR for a specific value of NDVI, where STR i is the STR value of the pixel i, while the STR max,i and STR min,i are the STR values of the dry and wet edge at the NDVI of pixel i (Figure 2b).

Study Area and Field Measurements
Our study focused on the two bogs, Linnusaare and Männikjärve, because of their available long-term records of WTD, which were measured since 1950 [73][74][75]. These bogs belong to the Endla mire system and they are located in eastern central Estonia ( Figure 3). The climate of the study region is temperate transitional between maritime and continental. The annual precipitation is about 660 mm and the average annual temperature is 5.2 • C. The bogs are of limnogenic origin and located within the East-Baltic Bog Province [76]. The tree layer of both bogs consists mainly of sparse Scots pines (Pinus sylvestris). In the inner-bog area, dwarf pines grow on ridges and they mark the treed ridge-hollow and treed ridge-pool complexes. The pine forest grows only at the bog's edges. The grass and dwarf shrub layers consist of Calluna vulgaris, Eriophorum vaginatum, Chamaedaphne calyculata, Andromeda polifolia, Rhynchospora alba, Ledum palustre, Oxycoccus microcarpus, and Oxycoccus palustris [77]. The typical moss species are Sphagnum fuscum, Sphagnum balticum, Sphagnum magellanicum, Remote Sens. 2020, 12, 1980 6 of 24 and Sphagnum rubellum [77]. Nowadays, Linnusaare and Männikjärve bogs have a conservation status and are on the list of Ramsar Sites [78].
Remote Sens. 2020, 12, x FOR PEER REVIEW 6 of 24 Air temperature and air pressure were measured at hourly resolution at Tooma weather station beginning from June 2009 ( Figure 3) [80]. Eight WTD monitoring wells were installed in hollows along a transect ( Figure 3). All the WTD measurements were referenced to the average height of the hollow surface. WTD was measured in sampling wells of 150 cm depth and 20 cm diameter; wells were permeable for water over the full length. In this study, we used daily measured WTD data from the growing seasons (May to September) from 2008 to 2019. The temporal dynamics in WTD were similar in different parts of the bogs. We observed a strong temporal correlation between WTD data from eight wells (Pearson correlation coefficients vary from 0.85 to 0.96). In addition, over short distances of a few meters, previous studies have shown only very small hydraulic gradients between hummock and hollow locations that can be neglected for our type of study [81]. However, the amplitudes of WTD fluctuation did differ: fluctuations were lower for the central part and higher for the bog edges, which is a peculiarity of peatland hydrology [82].

ERA 5 Data
ERA 5 is the atmospheric reanalysis of the global climate produced by the European Centre for Medium-Range Weather Forecasts [83]. ERA5 provides hourly atmospheric parameters as regular latitude-longitude grids at 0.25° spatial resolution. Data of friction velocity and total column water vapor were used to perform the modeling of the dry edge required for TOTRAM Scenario 2 (Appendix A and B). We selected ERA 5 data which in terms of the time were the nearest to Landsat acquisition time.

Remote Sensing Sources and Ancillary Data
Remotely sensed Landsat and MODerate-resolution Imaging Spectroradiometer (MODIS) images together with National Centers for Environmental Prediction and the National Center for Air temperature and air pressure were measured at hourly resolution at Tooma weather station beginning from June 2009 ( Figure 3) [80]. Eight WTD monitoring wells were installed in hollows along a transect ( Figure 3). All the WTD measurements were referenced to the average height of the hollow surface. WTD was measured in sampling wells of 150 cm depth and 20 cm diameter; wells were permeable for water over the full length. In this study, we used daily measured WTD data from the growing seasons (May to September) from 2008 to 2019. The temporal dynamics in WTD were similar in different parts of the bogs. We observed a strong temporal correlation between WTD data from eight wells (Pearson correlation coefficients vary from 0.85 to 0.96). In addition, over short distances of a few meters, previous studies have shown only very small hydraulic gradients between hummock and hollow locations that can be neglected for our type of study [81]. However, the amplitudes of WTD fluctuation did differ: fluctuations were lower for the central part and higher for the bog edges, which is a peculiarity of peatland hydrology [82].

ERA 5 Data
ERA 5 is the atmospheric reanalysis of the global climate produced by the European Centre for Medium-Range Weather Forecasts [83]. ERA5 provides hourly atmospheric parameters as regular latitude-longitude grids at 0.25 • spatial resolution. Data of friction velocity and total column water vapor were used to perform the modeling of the dry edge required for TOTRAM Scenario 2 (Appendices A and B). We selected ERA 5 data which in terms of the time were the nearest to Landsat acquisition time.

Remote Sensing Sources and Ancillary Data
Remotely sensed Landsat and MODerate-resolution Imaging Spectroradiometer (MODIS) images together with National Centers for Environmental Prediction and the National Center for Atmospheric Research (NCEP/NCAR) reanalysis data were used for the growing seasons from 2008 to 2019. Data processing was performed on the Google Earth Engine (GEE) online platform.
We processed Landsat 5 data for the time period between 2008 and 2012 since Landsat 5 was decommissioned in 2013. The resolution of visible and near-infrared (VNIR) bands is 30 m. Additionally, Landsat 5 had one thermal infrared band originally sensed with a resolution of 120 m and later resampled to 30 m. Landsat 7 data, unlike Landsat 5 and 8 data, cover the entire period of our study. Its data are recorded in 8 spectral bands; VNIR data have a spatial resolution of 30 m, but thermal data are recorded with 60 m spatial resolution and later resampled to 30 m [84]. Due to the Landsat 7 Scan Line Corrector failure, the study area was partly covered with stripes of missing data. We did not apply any blending techniques to fill in stripes of missing data. Landsat 8 (launched in 2013) data cover only the last six years of our study period. It provides data with the same spatial resolution as Landsat 7, except for the two thermal bands, which are recorded with 100 m spatial resolution and resampled to 30 m [85].
Each Landsat satellite overflew the study area every 8 days since the study area overlaps with two Landsat paths: 186 and 187 ( Figure 3). For the next processing steps, we used only those Landsat images that included at least 400 of the 10,549 bog pixels after filtering the pixels affected by shadows, clouds, and Landsat 7 scan-line-off stripes. This minimum threshold was subjectively defined to include a sufficient amount of bog pixels during the definition of the trapezoid space. In total, 168 Landsat images were selected for further calculation and analysis. We used the following Landsat products: at-sensor radiance, brightness temperature, cloud mask, and metadata information about solar zenith angle and relative Earth-Sun distance.
MODIS aboard Terra provides MOD11A1 V6 product, which includes LST and emissivity data. Both have 1 km spatial resolution and were daily sensed over the study area. We used the MODIS emissivity product sensed in two bands, namely 31 and 32, which in terms of acquisition time was the nearest to the Landsat acquisition time.
The data of atmospheric water vapor content produced by NCEP/NCAR were used. This dataset is available in GEE at 2.5 • spatial and six hourly temporal resolutions (00:00, 06:00, 12:00, and 18:00 UTC).

Variable Derivation
The basic concept of data preparation to derive moisture indices using either of the two TOTRAM scenarios or OPTRAM is presented in Figure 4. The figure illustrates the sequence of data analysis steps with inputs, intermediate variables, and outputs. As mentioned above, in Scenario 1, the TOTRAM estimation was solely based on observed pixel values. The remotely sensed LST and FVC data were used to construct trapezoid and indicate the dry and wet edges. In Scenario 2, in addition to LST and FVC, auxiliary data were used to model LSTs max " and LSTc max ", including friction velocity and total column water vapor data from ERA5, land surface albedo, and in-situ meteorological data. The OPTRAM soil moisture index was calculated based on remotely sensed STR and NDVI data. Remote Sens. 2020, 12, x FOR PEER REVIEW 8 of 24 The single-channel algorithm was used to calculate Land Surface Temperature (LST) for Landsat 5, 7, and 8 data [86]. We implemented the single-channel methodology and calculated LST in GEE following the algorithm described in [87].
The single-channel algorithm requires only one thermal channel to retrieve LST (K), where ε is the land surface emissivity; ψ1, ψ2, and ψ3 are atmospheric functions which can be obtained as a function of the total atmospheric water vapor content w; Lsen is a thermal radiance at sensor level in W m -2 sr -1 μm -1 ; and γ and δ are two parameters which are dependent on the Planck's function and were calculated using: where Tsen is at-sensor brightness temperature (K) and bγ is 1277 K for Landsat 5 and 7 and 1324 K for the 10th band of Landsat 8. In Equation (3), the emissivity ε was considered as the average value of bands 31 and 32 provided by MODIS MOD11A1 V6 product. We chose MODIS derived emissivity since it has several advantages for our study. Firstly, the MODIS dataset fully covers the study period. To be precise, GEE provides MOD11A1 dataset available started from 2000 until the present time. Secondly, the algorithm for MODIS emissivity calculation considers the land cover types [88]. The sixth collection of MODIS emissivity product, which is used in this study, includes "Permanent wetlands" as a class of land cover, which likely increases the accuracy of the subsequent LST estimates. Thirdly, previous study has demonstrated the high accuracy of MODIS derived emissivity for LST calculation [87].
GEE gives access to the Landsat products of different processing levels. The dataset with raw scenes contain orthorectified, scaled, and calibrated at-sensor radiance images. The thermal infrared radiance-at-sensor band Lsen was taken from this collection. Another Landsat dataset is the surface reflectance, which is an atmospherically corrected and orthorectified collection of images. This product includes a brightness temperature band Tsen, which was used in LST calculation. In addition, The single-channel algorithm was used to calculate Land Surface Temperature (LST) for Landsat 5, 7, and 8 data [86]. We implemented the single-channel methodology and calculated LST in GEE following the algorithm described in [87].
The single-channel algorithm requires only one thermal channel to retrieve LST (K), where ε is the land surface emissivity; ψ 1 , ψ 2 , and ψ 3 are atmospheric functions which can be obtained as a function of the total atmospheric water vapor content w; L sen is a thermal radiance at sensor level in W m −2 sr −1 µm −1 ; and γ and δ are two parameters which are dependent on the Planck's function and were calculated using: where T sen is at-sensor brightness temperature (K) and b γ is 1277 K for Landsat 5 and 7 and 1324 K for the 10th band of Landsat 8. In Equation (3), the emissivity ε was considered as the average value of bands 31 and 32 provided by MODIS MOD11A1 V6 product. We chose MODIS derived emissivity since it has several advantages for our study. Firstly, the MODIS dataset fully covers the study period. To be precise, GEE provides MOD11A1 dataset available started from 2000 until the present time. Secondly, the algorithm for MODIS emissivity calculation considers the land cover types [88]. The sixth collection of MODIS emissivity product, which is used in this study, includes "Permanent wetlands" as a class of land cover, which likely increases the accuracy of the subsequent LST estimates. Thirdly, previous study has demonstrated the high accuracy of MODIS derived emissivity for LST calculation [87].
GEE gives access to the Landsat products of different processing levels. The dataset with raw scenes contain orthorectified, scaled, and calibrated at-sensor radiance images. The thermal infrared radiance-at-sensor band L sen was taken from this collection. Another Landsat dataset is the surface reflectance, which is an atmospherically corrected and orthorectified collection of images. This product includes a brightness temperature band T sen , which was used in LST calculation. In addition, Landsat surface reflectance dataset contains cloud, shadow, and water mask called CFMASK, which we applied for per-pixel quality estimation.
To eliminate the atmospheric effect on LST, Equation (3) includes atmospheric functions (ψ 1 (w), ψ 2 (w), and ψ 3 (w)) and, thus, indirectly accounts for different atmospheric water vapor content (w) across the study time. For atmospheric water vapor content, we used reanalysis data produced by NCEP/NCAR. To cope with time misalignment between Landsat and NCEP/NCAR datasets, we followed the approach presented in [87] and took the value of column water content as the averaged one between two water content values, which were the closest in time to Landsat acquisition time.
To calculate atmospheric function for Landsat 5 and 7 we used the following coefficients obtained in [86]: Another matrix of coefficients was used for the 10th band (10.9 µm) of Landsat 8 [89]: Figure 5 shows the time series of LST along with other variables and indices over the studied bog area. For the whole study period, LST values were higher than the air temperature, which is in agreement with the previous study of this bog [80]. LST and air temperature follow similar seasonal curves with maximum values in July.
Remote Sens. 2020, 12, x FOR PEER REVIEW 9 of 24 Landsat surface reflectance dataset contains cloud, shadow, and water mask called CFMASK, which we applied for per-pixel quality estimation.
To eliminate the atmospheric effect on LST, Eq. (3) includes atmospheric functions (ψ1(w), ψ2(w), and ψ3(w)) and, thus, indirectly accounts for different atmospheric water vapor content (w) across the study time. For atmospheric water vapor content, we used reanalysis data produced by NCEP/NCAR. To cope with time misalignment between Landsat and NCEP/NCAR datasets, we followed the approach presented in [87] and took the value of column water content as the averaged one between two water content values, which were the closest in time to Landsat acquisition time. To calculate atmospheric function for Landsat 5 and 7 we used the following coefficients obtained in [ Figure 5 shows the time series of LST along with other variables and indices over the studied bog area. For the whole study period, LST values were higher than the air temperature, which is in agreement with the previous study of this bog [80]. LST and air temperature follow similar seasonal curves with maximum values in July.

Normalized Difference Vegetation Index (NDVI)
NDVI was proposed in the 1970s and since that time it has been widely used as an indicator of vegetation photosynthetic activity [90]. NDVI is calculated using near-infrared (NIR) and red reflectivity, ρNIR and ρRed, since chlorophyll in canopy highly reflects in the NIR range and absorbs red wavelengths, Figure 5. Time series of daily air temperature measured at 9 am UTC (12 am local time) at the weather station, and LST, STR, and NDVI as median values over the Linnusaare and Männikjärve bogs. The gaps in the air temperature line mark missing data due to the technical problems.

Normalized Difference Vegetation Index (NDVI)
NDVI was proposed in the 1970s and since that time it has been widely used as an indicator of vegetation photosynthetic activity [90]. NDVI is calculated using near-infrared (NIR) and red reflectivity, ρ NIR and ρ Red , since chlorophyll in canopy highly reflects in the NIR range and absorbs red wavelengths, Figure 5 shows time series of NDVI derived from the Landsat surface reflectance dataset in GEE. CFMASK was used for cloud and shadow masking. The spatial median of the NDVI values varied in time from 0.46 to 0.71. The highest NDVI values were observed in July and August.

Fractional Vegetation Cover (FVC)
FVC was obtained from NDVI through the equation [68,91], where NDVI max and NDVI min are the maximum and minimum values of NDVI corresponding to the values of dense vegetation and bare soil.

Shortwave Infrared Transformed Reflectance (STR)
STR is a remotely sensed variable used to estimate soil water content [72]. The STR is calculated from SWIR reflectance ρ SWIR , STR was calculated using the SWIR channel (2210 µm) taken from Landsat surface reflectance dataset. Similar to NDVI, we used CFMASK to mask clouds and shadows. The spatial median of the STR values varied in time between 2.93 and 7.77 ( Figure 5).

Broadband Albedo of Vegetated and Bare Surfaces
The broadband albedo was calculated for Landsat data according to the technique proposed by Liang [92], where α is a broadband albedo and ρ is surface reflectance measured in blue, red, NIR, and two SWIR spectra. This conversion from narrow-to broadband albedo was originally developed for Landsat 5 and 7, but it is commonly applied also to Landsat 8 bands [93][94][95]. The calculated albedo for a pixel is assumed to be the weighted sum of two components: albedo of canopy α c and albedo of soil α s [55]. To estimate both albedo components, Long et al. [55] suggested fitting a linear model to the maximum α values which lie within the interval "mean plus standard deviation" along FVC axis. This approach is reasonable when a linear relationship is observed between α and FVC. However, many studies revealed a non-linear relationship between α and vegetation cover [96][97][98][99]. We observed in the data of our study area a non-linear relationship between α and FVC and, therefore, the method suggested in [55] was modified. We fitted two linear models: the first one for albedo data with FVC range between 0 and 0.5 and the second one for albedo data with FVC between 0.8 and 1. In such a way, the first model was fitted to the pixels with dominant bare soil, the second one with dense healthy vegetation. Knowing the slope and intercept of the linear model, it became possible to estimate α for dry bare and vegetated surfaces using LSTs max " and LSTc max " (dry edge values, see Figure 2a), considering that their FVC values are 1 and 0.

Correlation Analysis
We determined the temporal Pearson correlation coefficients (R) between soil moisture indices of the three different trapezoid models and in situ WTD measurements. Anomaly Pearson correlation coefficients (anomR) were calculated to evaluate the capability of the trapezoid models to monitor the short-term and interannual variability of WTD. Anomaly time series of in-situ and soil moisture indices were obtained by removing the multi-year one-month-smoothed average from the original time series values. A limited spatial analysis on the correspondence between in situ WTD patterns and retrieved soil moisture indices was also included. time series values. A limited spatial analysis on the correspondence between in situ WTD patterns and retrieved soil moisture indices was also included. Figure 6 shows box plots (over eight wells) of the temporal correlation coefficients between WTD and soil moisture index estimated with two TOTRAM scenarios and OPTRAM as average values across four pixels closest to the wells. Both TOTRAM scenarios resulted in negative average R and anomR values of -0.19 and -0.23, respectively, for Scenario 1, and -0.16 and -0.08, respectively, for Scenario 2. In contrast, OPTRAM resulted in a weak positive R and anomR with average values of 0.41 and 0.37, respectively. The variability of R and anomR between the wells was also the highest for OPTRAM, with both low and high correlation values for the wells in the treed parts of the bog.   Figure 7b) as an average across four pixels closest to the wells. The WTD data highlights a fundamental relationship between WTD and vegetation type in bogs. Tree cover is an indicator of the deeper long-term WTD and higher temporal WTD variability [82].

Temporal Correlation of Soil Moisture Indices with WTD
As already indicated by the correlation coefficients, Figure 7 illustrates that the temporal variation of soil moisture indices estimated from TOTRAM Scenarios 1 and 2 does not agree with temporal changes in WTD. In contrast, the OPTRAM soil moisture index follows reasonably well the WTD dynamics, in particular in the treeless bog. However, some obvious outliers are still present.   Figure 7b) as an average across four pixels closest to the wells. The WTD data highlights a fundamental relationship between WTD and vegetation type in bogs. Tree cover is an indicator of the deeper long-term WTD and higher temporal WTD variability [82].
As already indicated by the correlation coefficients, Figure 7 illustrates that the temporal variation of soil moisture indices estimated from TOTRAM Scenarios 1 and 2 does not agree with temporal changes in WTD. In contrast, the OPTRAM soil moisture index follows reasonably well the WTD dynamics, in particular in the treeless bog. However, some obvious outliers are still present.
As mentioned above, we observed high spatial correlation among the water levels measured in different monitoring wells. For the next analysis, we therefore assumed coherent fluctuations of WTD over the whole bog by calculating the mean WTD of the eight monitoring wells. We then analyzed the temporal correlation coefficient R between the mean WTD and the three different soil moisture indices for each individual Landsat pixel. As mentioned above, we observed high spatial correlation among the water levels measured in different monitoring wells. For the next analysis, we therefore assumed coherent fluctuations of WTD over the whole bog by calculating the mean WTD of the eight monitoring wells. We then analyzed the temporal correlation coefficient R between the mean WTD and the three different soil moisture indices for each individual Landsat pixel.
The resulting spatial differences in R were assumed to be dominated by the local specifics of the temporal relationship between the WTD and soil moisture index. Figure 8 shows that, for the two TOTRAM scenarios, the R values are close to zero. In contrast, R values of OPTRAM are positive throughout the whole bog area. The highest R values (up to 0.8) can be observed for treeless areas, where no pools and dwarf pines are present.  Figure 9 shows scatterplots between the moisture indices and WTD across all times for two exemplary monitoring wells (same as in Figure 7): well 225 located in the treed part and well 323 The resulting spatial differences in R were assumed to be dominated by the local specifics of the temporal relationship between the WTD and soil moisture index. Figure 8 shows that, for the two TOTRAM scenarios, the R values are close to zero. In contrast, R values of OPTRAM are positive throughout the whole bog area. The highest R values (up to 0.8) can be observed for treeless areas, where no pools and dwarf pines are present. As mentioned above, we observed high spatial correlation among the water levels measured in different monitoring wells. For the next analysis, we therefore assumed coherent fluctuations of WTD over the whole bog by calculating the mean WTD of the eight monitoring wells. We then analyzed the temporal correlation coefficient R between the mean WTD and the three different soil moisture indices for each individual Landsat pixel.

Spatial Variability of Soil Moisture Indices and WTD
The resulting spatial differences in R were assumed to be dominated by the local specifics of the temporal relationship between the WTD and soil moisture index. Figure 8 shows that, for the two TOTRAM scenarios, the R values are close to zero. In contrast, R values of OPTRAM are positive throughout the whole bog area. The highest R values (up to 0.8) can be observed for treeless areas, where no pools and dwarf pines are present.  Figure 9 shows scatterplots between the moisture indices and WTD across all times for two exemplary monitoring wells (same as in Figure 7): well 225 located in the treed part and well 323  Figure 9 shows scatterplots between the moisture indices and WTD across all times for two exemplary monitoring wells (same as in Figure 7): well 225 located in the treed part and well 323 located in the treeless part of the bogs. Figure 9a,b indicates that the soil moisture indices from both TOTRAM scenarios do not positively correlate with WTD. Only the OPTRAM index shows a positive correlation over the whole WTD range (Figure 9c). Furthermore, Figure 9c shows that OPTRAM values are systematically higher over treed parts of the bog and lower over treeless parts of the bog. located in the treeless part of the bogs. Figure 9a,b indicates that the soil moisture indices from both TOTRAM scenarios do not positively correlate with WTD. Only the OPTRAM index shows a positive correlation over the whole WTD range (Figure 9c). Furthermore, Figure 9c shows that OPTRAM values are systematically higher over treed parts of the bog and lower over treeless parts of the bog.  Figure 10 presents spatial patterns of soil moisture index generated by the two TOTRAM scenarios and OPTRAM for an exemplary wet (5 May 2016) and dry (23 August 2018) conditions. The TOTRAM Scenario 1 and OPTRAM indices indicate a basic agreement with in-situ WTD in terms of spatial averages of the two contrasting conditions. This is not the case for the TOTRAM Scenario 2 in which a wetter condition was indicated for the date with the dry condition. For this date, this was caused by a modeled dry edge that was too warm and therefore led to an underestimation of the TOTRAM Scenario 2 index.

Spatial Variability of Soil Moisture Indices and WTD
None of the three methods was able to indicate plausible soil moisture variability. Under the dry conditions of the second date, the bogs typically show drier soil moisture and deeper WTD towards the margins. The opposite can be observed for OPTRAM, while TOTRAM does show very little variability. Furthermore, in-situ data of the second date indicated deepest WTD at wells 217 and 222, which are installed in the treed parts of the bogs. The OPTRAM index, however, did show relatively high values around those two wells.  Figure 10 presents spatial patterns of soil moisture index generated by the two TOTRAM scenarios and OPTRAM for an exemplary wet (5 May 2016) and dry (23 August 2018) conditions. The TOTRAM Scenario 1 and OPTRAM indices indicate a basic agreement with in-situ WTD in terms of spatial averages of the two contrasting conditions. This is not the case for the TOTRAM Scenario 2 in which a wetter condition was indicated for the date with the dry condition. For this date, this was caused by a modeled dry edge that was too warm and therefore led to an underestimation of the TOTRAM Scenario 2 index.  Figure 10 presents spatial patterns of soil moisture index generated by the two TOTRAM scenarios and OPTRAM for an exemplary wet (5 May 2016) and dry (23 August 2018) conditions. The TOTRAM Scenario 1 and OPTRAM indices indicate a basic agreement with in-situ WTD in terms of spatial averages of the two contrasting conditions. This is not the case for the TOTRAM Scenario 2 in which a wetter condition was indicated for the date with the dry condition. For this date, this was caused by a modeled dry edge that was too warm and therefore led to an underestimation of the TOTRAM Scenario 2 index.
None of the three methods was able to indicate plausible soil moisture variability. Under the dry conditions of the second date, the bogs typically show drier soil moisture and deeper WTD towards the margins. The opposite can be observed for OPTRAM, while TOTRAM does show very little variability. Furthermore, in-situ data of the second date indicated deepest WTD at wells 217 and 222, which are installed in the treed parts of the bogs. The OPTRAM index, however, did show relatively high values around those two wells. None of the three methods was able to indicate plausible soil moisture variability. Under the dry conditions of the second date, the bogs typically show drier soil moisture and deeper WTD towards the margins. The opposite can be observed for OPTRAM, while TOTRAM does show very little variability. Furthermore, in-situ data of the second date indicated deepest WTD at wells 217 and 222, which are installed in the treed parts of the bogs. The OPTRAM index, however, did show relatively high values around those two wells. Figure 11 shows maps of the anomaly time series of the soil moisture indices for four exemplary dates and the corresponding soil moisture anomaly averaged over all monitoring wells. In terms of spatial averages, there is again a good agreement between anomalies in WTD and OPTRAM index. In contrast, the TOTRAM scenarios did not yield changes that could be related to WTD.
Remote Sens. 2020, 12, x FOR PEER REVIEW 14 of 24 Figure 10. Maps of soil moisture index generated with TOTRAM Scenario 1, TOTRAM Scenario 2, and OPTRAM together with indicated water table depth (WTD) for each monitoring well. The white areas within the bogs represent missing data resulting from Landsat 7 Scan Line Corrector failure (striped pattern on 05/05/2016) or methodological constraints (filtering of oversaturated pixels in OPTRAM). The black hatched pattern indicates the treed bog areas. Figure 11 shows maps of the anomaly time series of the soil moisture indices for four exemplary dates and the corresponding soil moisture anomaly averaged over all monitoring wells. In terms of spatial averages, there is again a good agreement between anomalies in WTD and OPTRAM index. In contrast, the TOTRAM scenarios did not yield changes that could be related to WTD.

Lack of Correlation between TOTRAM Index and WTD
Given the strong linkage between soil moisture and WTD present in peatlands, a positive correlation between the TOTRAM soil moisture index and WTD was expected. However, such a positive correlation could not be observed for TOTRAM Scenario 1 or TOTRAM Scenario 2. Table 1 summarizes the limitations of TOTRAM approaches. Here, we emphasize two aspects of TOTRAM that hampered its applicability to our study area; both aspects are related to the generally wetter conditions in northern latitude ecosystems.
The most important limitation is that the TOTRAM indices are computed per scene. This limits the ability to sufficiently constrain the trapezoid space. The generally wetter conditions in northern

Lack of Correlation between TOTRAM Index and WTD
Given the strong linkage between soil moisture and WTD present in peatlands, a positive correlation between the TOTRAM soil moisture index and WTD was expected. However, such a positive correlation could not be observed for TOTRAM Scenario 1 or TOTRAM Scenario 2. Table 1 summarizes the limitations of TOTRAM approaches. Here, we emphasize two aspects of TOTRAM that hampered its applicability to our study area; both aspects are related to the generally wetter conditions in northern latitude ecosystems. The most important limitation is that the TOTRAM indices are computed per scene. This limits the ability to sufficiently constrain the trapezoid space. The generally wetter conditions in northern latitudes may not sufficiently cover the range of moisture variability to properly constrain the trapezoid space. Our attempt to overcome this limitation by modeling a dry edge based on weather conditions (Scenario 2) did not improve the performance of TOTRAM.
A further possible reason for the lack (or even negative) correlation between WTD and TOTRAM indices might be that evapotranspiration and vegetation growth in northern latitudes are much more limited by energy than by moisture [100]. This likely leads to the violation of the assumption of TOTRAM that lower soil moistures lead to an increase in LST. Karnieli et al. [101] showed that LST and VI can have a positive association in the area where vegetation growth is limited by energy. For example, in high latitude regions, warming is accompanied by vegetation growth and, thus, an LST increase cannot be a reliable indicator of vegetation stress. Garcia et al. [102] demonstrated that this is a limitation of TOTRAM and suggested to apply it only to the regions where water and not energy is limiting. Further research should investigate the applicability of TOTRAM for the WTD monitoring in peatlands of climatic zones for which energy is not a limiting factor of vegetation growth (Table 1).

Potential and Challenges of Using OPTRAM for WTD Monitoring
Our results demonstrate a high potential of the OPTRAM index to monitor temporal changes of WTD dynamics in northern bogs. The average correlation coefficients were as high as the ones obtained for radar backscatter coefficient over natural peatlands in an earlier study [38]. Especially over treeless areas, correlation coefficients achieved with OPTRAM index were very promising and reached values up to 0.8. The good temporal correlation statistics suggests that OPTRAM index might be a new type of observation suitable for assimilation into peatland-specific hydrological models [33].
However, a closer investigation of the local time series showed some obvious outliers. The identification of the reasons for these outliers and methods to filter such data will be critical to increase the applicability of the OPTRAM index for WTD monitoring in peatlands. Most importantly, it is possible that the OPTRAM index represents more or less the vegetation status or the soil moisture at various depths at different times. As OPTRAM relies on the dependency of the SWIR reflectance on plant moisture content, the soil depth dependency needs to be developed for different vegetation types and their typical rooting depths.
Similar to the microwave backscatter signal and coherence, the spatial variation of the OPTRAM index is highly influenced by the vegetation patterns [38,103]. For our study area, the simple differentiation into treed (covered with dwarf pines) and treeless (covered mainly with Sphagnum mosses) bog areas showed that vegetation variability was the dominating factor controlling spatial variability of the OPTRAM index. The importance of vegetation types and leaf internal structure of different species likely goes beyond such a coarse differentiation of vegetation types [44,70] as indicated by the large variance of the OPTRAM index within the treeless and treed areas (Figure 9). This variance is likely attributed to the differences in SWIR reflectance of species within one genus. For example, differently structured Sphagnum species have uneven ability to retain water under the same environmental conditions; thus, their behavior in the SWIR channel varies [104][105][106][107][108][109].
The high OPTRAM dependency on vegetation cover warrants a note of caution for its use in long-term monitoring studies (Table 1). Vegetation types and specific species communities can vary a lot in time in reaction to long-term changes of WTD [110]. This limitation is an important issue for future research, because the long Landsat record would be of particular interest for such a long-term monitoring of moisture changes in peatlands due to climate change or human disturbance. By now, it has been suggested combining SWIR channel data with NIR data to minimize the uncertainty in the estimation of soil moisture for different vegetation cover [44]. In further studies, it should be tried to integrate the NIR channel into the OPTRAM approach.
Finally, it should be reiterated that the WTD estimation with OPTRAM is successful here, mainly because of the strong vertical coupling of surface (and root-zone) soil moisture with the typical shallow WTD over bogs. It is unlikely (but not tested) that the method would perform as well in areas with deep WTD or mineral soil properties.

Conclusions
We evaluated the temporal and spatial relationships between in-situ measured WTD and soil moisture indices from three different trapezoid models over two northern bogs in Estonia. We compared the performance of two trapezoid models based on optical and thermal imagery (TOTRAM Scenarios 1 and 2) and one solely based on optical imagery (OPTRAM). Our results demonstrate: 1.
a general inapplicability of the TOTRAM index for the spatial and temporal WTD monitoring in our study area; 2.
a high potential of OPTRAM index for monitoring temporal changes in WTD with average temporal Pearson correlation coefficients of 0.41 for original and 0.37 for anomaly time series; 3.
the highest temporal correlation coefficients (0.8) for the OPTRAM index over treeless bog areas with little or no surface water (no bog pools); and 4. a high sensitivity of OPTRAM index to the vegetation type. Together with unknown spatial variability of the soil moisture-WTD relationship, this strongly limits the spatial interpretability and probably also the long-term temporal interpretability of the OPTRAM index for WTD monitoring under progressive changes of vegetation and peat properties.
The present study provides the first comprehensive assessment of the potential of TOTRAM and OPTRAM to monitor WTD dynamics over two selected northern bogs. Further studies are necessary to assess the performance over other peatlands in northern as well as tropical regions. Further research should be undertaken to improve the accuracy of OPTRAM and minimize the uncertainties caused by the differences in vegetation cover. Funding: This research was funded by the European Social Fund's Dora Plus Programme, grant number 36.9-6.1/1222, the Ministry of Education and Science of Estonia, grants numbers IUT2-16 and PRG352, and the EU through European Regional Development Fund for the Centre of Excellence "Ecology of Global Change: Natural and Managed Ecosystems" (EcolChange), grant number is not applicable. M. Bechtold thanks the Alexander von Humboldt Foundation for a Feodor Lynen Fellowship.

Acknowledgments:
We are grateful to the staff of the Tooma mire research station, who provided hydrometeorological data for this study.

Conflicts of Interest:
The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

Appendix A
Estimation of the temperature of the driest bare surface (LSTs max ) Following the two-source model described in [55], LSTs max could be implicit in radiation budget and energy balance equations for the driest bare soil surface which does not evaporate: where R s is the net radiation of soil surface (W m −2 ); α s is the albedo for the soil surface (unitless); S d is the downwelling shortwave radiation (W m −2 ); L d and L u are the downwelling and upwelling longwave radiation (W m −2 ); ε s is the broadband emissivity of the soil surface (unitless); ε a is the emissivity of the atmosphere (unitless); ρ is the air density (kg m −3 ); c p is the air specific heat at the constant pressure taken as 1.013 × 10 3 (J kg −1 K −1 ) [111]; r a,s is the aerodynamic resistance of the bare soil surface (s m −1 ); σ is the Stefan-Boltzmann constant 5.67 × 10 −8 (kg s −3 K −4 ); c is the calibrated proportionality coefficient equal to 0.35 (unitless); T a is the air temperature (K); and LE s is the evaporation from soil surface (W m −2 ). The aerodynamic resistance for the bare soil surface could be calculated in the following way: where u 1m is the wind velocity at the height of 1 m above the bare soil surface (m s −1 ), which can be estimated from the logarithmic wind profile function: where u * is the friction velocity (m s −1 ); k is the von Karman's constant equal to 0.41 (unitless); z m is the reference height of the wind velocity measurements, which is taken as 1 m; d is the zero plane displacement taken as 0 m; z om is the roughness length for momentum transfer assumed to be 0.005 m [112]; and ϕ m(z m ) and ϕ m(z om ) are the factors used for the stability correction at 1 m and 0.005 m height (unitless), which depend on the Monin-Obukhov length L (m) and can be expressed for the stable conditions (L > 0) as follows: where L is calculated as follows: where g is gravitational acceleration taken as 9.8 (m s −2 ). Finally, LSTs max can be obtained iteratively using Equations (A1)-(A4) and (A7).

Appendix B
Estimation of the temperature of the driest fully vegetated surface (LSTc max ) In a similar way, Long et al. [55] suggested calculating LSTc max from the radiative budget equation: where R c is the net radiation of fully vegetated surface (W m −2 ); α c is the albedo for vegetation (unitless); ε c is the broadband emissivity of vegetation (unitless); LE c is the evaporation from vegetation (W m −2 ); and r a,c is the aerodynamic resistance above the canopy surface (s m −1 ) expressed as follows: where u * is calculated following Equation (A5), but ϕ m(z m ) and ϕ m(z om ) for the stable condition (L > 0) are calculated in the following way: For the unstable condition (L < 0), ϕ m(z m ) and ϕ m(z om ) are: ϕ m(z m ) = 2 ln 1 + x (z m ) ϕ m(z om ) = 2 ln 1 + x (z om ) where z T is the reference height for temperature measurements assumed to be 2 m and the vegetation height is h c = 1 m; and d, z om , and z oh are taken as: LSTc max can be derived by solving for the system of Equations (A13)-(A15) and (A20).