Spatio-Temporal Assessment of Surface Moisture and Evapotranspiration Variability Using Remote Sensing Techniques

: The relationship between the physic features of the Earth’s surface and its temperature has been signiﬁcantly investigated for further soil moisture assessment. In this study, the spatiotemporal impacts of surface properties on land surface temperature ( LST ) were examined by using Landsat-8 Operational Land Imager (OLI) and Thermal Infrared Sensor (TIRS) and meteorological data. The signiﬁcant distinctions were observed during a crop growing season through the contrast in the correlation between different multi-spectral satellite indices and LST , in which the highest correlation of − 0.65 was found when the Normalized Difference Latent heat Index ( NDLI ) was used. A new index, named as Temperature-soil Moisture Dryness Index ( TMDI ), is accordingly proposed to assess surface moisture and evapotranspiration (ET) variability. It is based on a triangle space where NDLI is set as a reference basis for examining surface water availability and the variation of LST is an indicator as a consequence of the cooling effect by ET. TMDI was evaluated against ET derived from the commonly-used model, namely Surface Energy Balance Algorithm for Land (SEBAL), as well as compared to the performance of Temperature Vegetation Dryness Index ( TVDI ). This study was conducted over ﬁve-time points for the 2014 winter crop growing season in southern Taiwan. Results indicated that TMDI exhibits signiﬁcant sensitivity to surface moisture ﬂuctuation by showing a strong correlation with SEBAL-derived ET with the highest correlation of − 0.89 was found on 19 October. Moreover, TMDI revealed its superiority over TVDI in the response to a rapidly changing surface moisture due to water supply before the investigated time. It is suggested that TMDI is a proper and sensitive indicator to characterize the surface moisture and ET rate. Further exploitation of the usefulness of the TMDI in a variety of applications would be interesting.


Introduction
Surface moisture plays a significant role in land-surface processes by modulating the exchanges of water and energy between land and atmosphere [1][2][3][4]. Understanding surface moisture and evapotranspiration (ET) can provide valuable information for a variety of subsequent applications, such as hydro-meteorological predictions, drought monitoring, and agriculture management [5,6]. Various techniques have been proposed to assess surface moisture and quantify ET. Ground instruments have long been used. However, the limitation in spatial coverage of point-wise measurements, ground instruments cannot sufficiently provide complete coverage of surface signature over a large region or on a global scale [7]. ET estimation based on the surface energy balance theory using ground observations and remotely sensed data has been broadly demonstrated [8,9]. However, they still cannot be conveniently implemented globally due to the sparsely distributed ground observations and the complexity of integration models. In the recent decades, the approaches based on empirical parameterization of the LST and surface properties have been proposed for soil moisture monitoring. Temperature vegetation dryness index (TVDI) has been widely used for soil moisture assessment and drought monitoring [10,11]. This method is based on an interpretation of the simplified NDVI-LST space. That is, NDVI is sensitive to the appearance of vegetation and its influence on LST through transpiration and evaporation from plants and soil, respectively [12]. Many studies indicated that TVDI can reveal the soil moisture status and closely associated with ET [13,14]. However, TVDI is severely affected by the rapid water supply due to an immediate change in surface water availability and a delay in the response of vegetation. Moreover, no correlation between LST and NDVI was revealed when NDVI is less than 0.1 as presented by Du et al. [12], with open water and crop in the inundated regions excluded from the analysis.
Surface water has a great potential to moderate LST through absorbing solar radiation and converting it into latent heat in the form of ET [15]. The reduction of water coverage and decrease in surface water availability lead to increased sensible heat transfer as latent heat transfer is constrained, resulting in rising LST. Thus, the surface water availability indicators and their correlation with LST variability have the potential to reveal more complete information on surface moisture and ET rate from bare soil to highly vegetated cover. Thermal infrared satellites provide a straightforward and consistent method to derive the LST at regional scales with the different spatial and temporal resolutions [16][17][18]. Many studies have been proposed to estimate LST by using remotely sensed thermal infrared data with different methods to tackle the atmospheric effects, such as radiative equation-based method [19], split window algorithm [20], and single channel method [21]. On the other hand, multi-spectral remote sensing has been widely utilized to extract and analyze surface properties of top layer by taking the advantage in optimizing the spectral sensitivity on certain characteristic of surface materials [22]. Accordingly, various indices have been widely used for the purpose of drought and water status assessment in certain situations. Since the appropriate index selection is of particular importance, five frequently used drought indices are examined in this study as addressed hereafter. Normalized Difference Vegetation Index (NDVI) is widely used because of its sensitivity to the interaction between solar radiation and vegetation cover [23]; Normalized Difference Moisture Index (NDMI) is considered as a complementary index for remote sensing of vegetation water content because of its sensitivity to the parameter of interest with less influence by the atmospheric effect [23,24]; Normalize Difference Water Index (NDWI) is used to delineate open water features because it can enhance water signature due to nearly complete absorption of infrared light by water [24]; Bare Soil Index (BI) is used to collect the information of bare soil areas, fallow lands, and vegetation by using four different bands of satellite images [25]; and Normalized Difference Latent heat Index (NDLI) is proven to be an effective indicator in extracting surface latent heat information by examining water availability of the Earth's surface. The advantage of the NDLI is its uniqueness in optimizing the spectral sensitivity on biophysical land surface parameters by incorporating three commonly-used channels chosen for the satellite missions, including green, red, and shortwave-infrared [26].
This study aims to: (1) provide a reference for selecting the appropriate indicators for characterizing surface properties and their correlation with LST variability; and (2) propose a new and advantageous index for rapid surface moisture status and ET assessment by using remote sensing data, in which a potential multi-spectral index was set as a reference basis for examining surface water availability and the variation of LST was set as the parameter of ET cooling effect. It is expected to provide a potential indicator for further studies related to the large-scale assessment of Earth's surface properties, especially plant water stress and physical surface water scarcity, by using different spatial-temporal resolution images from the great number of available Earth observation satellites.

Study Area
The study area is an alluvial plain located along the south-west coast of Taiwan, including the Yunlin and Chiayi counties, and Chiayi and Tainan cities as shown in Figure 1: Geographic locations of the study area. This area faces the Taiwan Strait on the west and lies to the west of foothills that extend from the Alishan Mountain. It is known as a largely flat landform created by alluvial deposits derived from the eastern foothill area through fluvial transportation [27]. It is of monsoon-influenced humid subtropical climate region with high humidity all year round, even during winter. It is once a place plenty of salt production, but most of it has been abandoned in the recent years. Additionally, many fish farms are located along the seacoast. Paddy rice is a major crop of the study region. The other primary crops include sugarcane, peanuts, corn, sweet potato, and some floricultural plants and vegetables. This study was conducted in the second growing season of 2014 since a significantly low amount of precipitation was recorded, especially in the southwest Taiwan. According to Taiwan's Central Weather Bureau (CWB), a lack of rainfall caused severe droughts and water scarcity for the needs of agricultural practice in the second half of 2014.

Study Area
The study area is an alluvial plain located along the south-west coast of Taiwan, including the Yunlin and Chiayi counties, and Chiayi and Tainan cities as shown in Figure  1: Geographic locations of the study area. This area faces the Taiwan Strait on the west and lies to the west of foothills that extend from the Alishan Mountain. It is known as a largely flat landform created by alluvial deposits derived from the eastern foothill area through fluvial transportation [27]. It is of monsoon-influenced humid subtropical climate region with high humidity all year round, even during winter. It is once a place plenty of salt production, but most of it has been abandoned in the recent years. Additionally, many fish farms are located along the seacoast. Paddy rice is a major crop of the study region. The other primary crops include sugarcane, peanuts, corn, sweet potato, and some floricultural plants and vegetables. This study was conducted in the second growing season of 2014 since a significantly low amount of precipitation was recorded, especially in the southwest Taiwan. According to Taiwan's Central Weather Bureau (CWB), a lack of rainfall caused severe droughts and water scarcity for the needs of agricultural practice in the second half of 2014.

Data and Image Pre-Processing
Remote sensing imageries data used in this study were acquired by the operational land imager (OLI) and thermal infrared sensor (TIRS) onboard the Landsat 8 satellite (Table 1). Satellite images were acquired from the United States Geographical Survey website

Data and Image Pre-Processing
Remote sensing imageries data used in this study were acquired by the operational land imager (OLI) and thermal infrared sensor (TIRS) onboard the Landsat 8 satellite ( were first captured with 100-m resolution, but are registered and delivered as the 30-m spatial resolution data product. Hourly average weather data were collected from the meteorological station over the entire Yun-Chia-Nan area. The station was set up by CWB with measurements taken at 1.25-2.00 m above the ground level with sensors sheltered from direct solar radiation (e-service.cwb.gov.tw, accessed on 19 October 2014, 4 November 2014, 20 November 2014, 6 December 2014, and 22 December 2014).

Land Use/Land Cover Mapping
As an agricultural region, the phenological stages of cultivation are associated with the changes of land cover characteristics. Several types of land covers, such as built-up, open water, and long-term vegetation, do not experience significant changes during a crop growing season. Croplands corresponding to high dynamic agriculture activities are related to crop rotation, including annual crops and semi perennial crops. A time series of spectral signatures from five Landsat images were analyzed to overcome single data weakness for land use/land cover (LULC) mapping in order to better capture the different phenology of land cover types over Yun-Chia-Nan Plain. The fieldworks, highresolution images from Google Earth, and expert knowledge were integrated to collect training samples by digitization in Landsat 8 OLI images. ENVI 5.3 software was used for image classification. Theory of supervised classification maximum likelihood method was applied to classify study area into five general classes, including water, built-up, forest, bare soil, and cropland, in which, the cropland class includes all main crops that were planted in one crop season.

Calculation of Spectral Indices
Multispectral remote sensing indices, including NDVI, NDWI, NDMI, BI, and NDLI were used to characterize the biophysical land surface parameters and quantify the relationships between land cover and thermal environment. These indices are expressed in the following equations: where ρ 2 , ρ 3 , ρ 4 , ρ 5 , and ρ 6

Land Surface Temperature Retrieval
The radiative transfer model was used to estimated LST from Landsat 8 TIRS [28,29]. Digital numbers of thermal band (10.62-12.51 µm) were firstly converted to Top of Atmosphere (TOA) radiance using radiometric rescaling coefficients provided in the product metadata file [30]. TOA radiance of thermal infrared band was then converted to LST by integrated with atmospheric simulation using radiative transfer equation. In this study, the atmospheric profile parameters, such as atmospheric transmittance, upwelling, and downwelling radiances, were required from the Moderate Resolution Atmospheric Transmission model for TIRS using the online calculator is available at https://atmcorr.gsfc.nasa.gov, accessed on 19 October 2014, 4 November 2014, 20 November 2014, 6 December 2014, and 22 December 2014 [31,32]. The equation can be expressed as: where C 1 is 14,387.7 µm·K; C 2 is 1.19104 × 10 8 W·µm 4 ·m −2 ·sr −1 ; L λ is top of atmosphere (TOA) radiance W/(m 2 ster µm); λ i is effective band wavelength for band i; I ↑ i is upwelling path radiance; I ↓ i is downwelling path radiance; τ i is atmospheric transmittance for channel i when view zenith angle is θ; and ε i is land surface emissivity for channel i.

Calculations of the Temperature-Surface Moisture Dryness Index
The spatial patterns of NDLI and LST and their interrelationships were used to characterize the surface properties and understand the response of surface water availability to the thermal environment. The scatter plot is established in the LST-NDLI triangle space to propose an advantageous method for assessing the surface moisture and ET through the deviation in LST with its minimum and maximum at the given NDLI intervals. The dry and wet edges were defined as the minimum and maximum LSTs, respectively, in the small NDLI interval of 0.02. TMDI is calculated using the approach of VIs-LST triangle as expressed below: where LST min is the minimum LST given NDLI along the wet edge; LST max is the maximum LST given NDLI along the dry edge; and LST i is the LST in any given pixel (K). The triangle space formed because LST decreases as ET increases at the given amount of water availability due to surface cooling by ET. This approach would enable us to understand the linkages between the surface water status and LST as the consequence of the cooling effect by ET at the earth's surface, and, thus, helping in the identification of the application potentials of TMDI in estimating ET rate and further monitoring the regional drought. The climate information and irrigation schedule were checked to make sure the fairness analysis throughout the entire process. In this study, ET derived from a simple surface energy balance model was used as reference data for comparison in order to demonstrate the ability of TMDI in terms of investigating surface moisture and ET [33,34]. The Surface Energy Balance Algorithm for Land (SEBAL) has been widely used and validated to estimate ET accurately using different satellite data. The advantages of SEBAL in terms of deriving ET was confirmed under the various conditions of different countries by ground-based measurements [35][36][37]. Following the SEBAL model, ET can be theoretically computed in terms of latent heat flux through a residual of land surface energy balance equation at the time of the satellite overpass as following equation: where R n is the net radiation (W m −2 ), G is the soil heat flux (W m −2 ), H is the sensible heat flux (W m −2 ), ET is the evapotranspiration (mm h −1 ), and λ is the latent heat of vaporization (MJ kg −1 ).

Land Cover Map and Validation
LULC of Yun-Chia-Nan Plain was aggregated into five major types, including built-up, vegetation, bare soil, water, and agriculture. Characteristics of LULC have been clearly shown in Figure 2. The visual interpretation shows that classification result is basically consistent with land surface features extracted from high spatial resolution image interpretation and information from field observations. Time-series data showed an additional source of information to identify land use types. This advantage is used to recognize the spectral signatures in accordance with the change in the surface canopy, which is not detectable in a single-date data. For instance, by using single-date data, the cultivated areas could be misidentified as built-up/dryland or water body in the land preparation stage and planting or after harvesting stage, respectively. Results indicated that the cropland dominates land cover, covering 59.86% of total area. Cropland consists of paddy field, dry land farming, and other agriculture purposes. Water bodies (lake, river, and aquaculture systems), mostly located along coast, are the second dominant land cover of this region, comprising approximately 16.62% of the land. Built-up, bare soil, and forest occupy 14.7%, 5.89%, and 2.92% of the land area, respectively. In order to make a quantitative analysis of land cover classification results, Google Earth Pro and expert knowledge of authors were used to understand the land surface characteristics as reference data for validation based on an independent set of training points randomly chosen from different sub-areas of the entire study area. Results show that time-series of Landsat 8 OLI can capture surface characteristics over different periods for LULC classification with overall accuracy of 88.6% and Kappa statistics of 0.83. where Rn is the net radiation (W m −2 ), G is the soil heat flux (W m −2 ), H is the sensible heat flux (W m −2 ), ET is the evapotranspiration (mm h −1 ), and λ is the latent heat of vaporization (MJ kg −1 ).

Land Cover Map and Validation
LULC of Yun-Chia-Nan Plain was aggregated into five major types, including builtup, vegetation, bare soil, water, and agriculture. Characteristics of LULC have been clearly shown in Figure 2. The visual interpretation shows that classification result is basically consistent with land surface features extracted from high spatial resolution image interpretation and information from field observations. Time-series data showed an additional source of information to identify land use types. This advantage is used to recognize the spectral signatures in accordance with the change in the surface canopy, which is not detectable in a single-date data. For instance, by using single-date data, the cultivated areas could be misidentified as built-up/dryland or water body in the land preparation stage and planting or after harvesting stage, respectively. Results indicated that the cropland dominates land cover, covering 59.86% of total area. Cropland consists of paddy field, dry land farming, and other agriculture purposes. Water bodies (lake, river, and aquaculture systems), mostly located along coast, are the second dominant land cover of this region, comprising approximately 16.62% of the land. Built-up, bare soil, and forest occupy 14.7%, 5.89%, and 2.92% of the land area, respectively. In order to make a quantitative analysis of land cover classification results, Google Earth Pro and expert knowledge of authors were used to understand the land surface characteristics as reference data for validation based on an independent set of training points randomly chosen from different sub-areas of the entire study area. Results show that time-series of Landsat 8 OLI can capture surface characteristics over different periods for LULC classification with overall accuracy of 88.6% and Kappa statistics of 0.83.

Spatio-Temporal Variation of LST
Time series variation of LST over Yun-Chia-Nan Plain was retrieved from Landsat 8 OLI/TIRS data. Figure 3 shows the spatial distribution of LST over five time-points, which

Spatio-Temporal Variation of LST
Time series variation of LST over Yun-Chia-Nan Plain was retrieved from Landsat 8 OLI/TIRS data. Figure 3 shows the spatial distribution of LST over five time-points, which were conducted during 2014 winter growing season. Results indicated that average LSTs  19 October, 4 November, 20 November, 6 December, and 22 December, respectively. LST values were then evaluated by the near-surface air temperature measured by CWB meteorological stations at 2-m height (T air ). Although there exist different physical meanings, the significant relationship between LST and T air has been reported in the literature [38,39]. Results indicated that the LST values are sufficiently accurate for further applications. LST and T air presented the strong correlation r 2 of 0.94 with RMSD of 5.13 Kelvins in the difference between LST and T air . The previous studies indicated that the difference between LST and T air can reach 20 (Kelvins) at noon under an extremely hot day [40].

Comparison of LST Retrievals with Spectral Indices
Five different multi-spectral indices, including NDVI, BI, NDMI, NDWI, and NDLI were derived from Landsat 8 OLI images. These indices were then used to establish and quantify their correlations with LST in order to understand the relationship between sur face properties and thermal environment variability. These correlations were assessed against an independent set of reference pixels randomly chosen from the entire study area To characterize the changes of cropland status and its driver on LST variability, we firstly explore the relationship between LST and multi-spectral indices within stable areas where land covers do not suffer significant change during a growing season, including built-up water, bare soil, and long-term vegetation. Afterwards, cropland was included in the re Results indicated that even the intensity of solar radiation varied with time, both day and season, the highest LSTs in all five time-points were mostly located in the northern region that covers some parts of Zhuoshui River basin with appearance of rock and sand. The relatively higher LST signatures mostly appeared in the built-up areas, such as the city centers of the two counties (Yunlin and Chiayi) and two cities (Chiayi and Tainan). Lower LST portions concentrated in the high vegetation coverage areas in the eastern study area. Water features, such as the lakes, rivers, and aquaculture systems in the western part of the study area, always correspond to the lowest LST in all five-time points. The significant Remote Sens. 2021, 13, 1667 8 of 15 variation of LST is recorded in cropland. Three areas of homogeneous agricultural zones are chosen and labeled as A, B, and C (as shown in Figure 2) for further examination. These three areas largely correspond to, respectively, maize, paddy, and sugarcane fields, which experience significant changes in the phonological stage during a growing season. The winter rice-growing season typically starts in early August and the harvest period is from late November to early December. While the major corn season is from September to February of the next year and the season of sugarcane starts in February and lasts for a year. Results indicated that (1) paddy rice field can be clearly seen as it has lower LST compared to the other croplands, such as corn and sugarcane, on 19 October and 4 November. This is in good accordance with the fact that, at these two time points, rice plants were at its vegetative growth stage, while corn maize field in the early stage. The paddy rice field often has a lower surface temperature than that of dry-soil planting areas due to evaporative cooling. (2) The similarity in the value of LST over three selected areas was observed on 20 November. That is, the growth and development of the corn plant enhance latent heat flux and decrease LST because vegetation often has a higher ET rate than dry land and absorbs more energy. While the paddy rice field is converted to the reproductive stage, which is associated with the decrease of surface ET rate and increase of LST. (3) On 6 and 22 December, LST values of paddy rice field were relatively higher than both corn and sugarcane fields. On these days, the paddy rice field is usually drained out of water before harvest, the surface water availability was significantly decreased, resulting in lower LST over paddy rice field compared to corn and sugarcane areas at their vegetative stage.

Comparison of LST Retrievals with Spectral Indices
Five different multi-spectral indices, including NDVI, BI, NDMI, NDWI, and NDLI, were derived from Landsat 8 OLI images. These indices were then used to establish and quantify their correlations with LST in order to understand the relationship between surface properties and thermal environment variability. These correlations were assessed against an independent set of reference pixels randomly chosen from the entire study area. To characterize the changes of cropland status and its driver on LST variability, we firstly explore the relationship between LST and multi-spectral indices within stable areas where land covers do not suffer significant change during a growing season, including built-up, water, bare soil, and long-term vegetation. Afterwards, cropland was included in the relationship analysis over a crop growing season. We found that the correlations between indices and LSTs over stable land covers are relatively higher than that over the entire study area. LST is influenced by soil moisture, crop types, and crop plant condition. The changes in cropland canopies lead to variation in solar energy absorption and significantly influence LST through ET, which can be suddenly changed between different periods of crop growth and agricultural transformation.
The different correlations between multi-spectral indices and LST are presented in Table 2. Results indicated that the properties of vegetation cover, open water, and bare soil represented, respectively, by NDVI, NDWI, and BI, exhibited insignificant impacts on near-surface thermal environment variability, as evidenced by their low correlations with LST. While the characteristics of vegetation water content and surface water availability represented, respectively, by NDMI and NDLI exhibited a significant impact on the temperature at the land-air interface, resulting in higher correlations with LST. This is in good accordance with the fact that the cooling effect as the consequence of ET is significantly controlled by the amount of water content in the surface canopies. The highest correlations between indices and LST at all five-time points are represented by NDLI with reliable negative correlation coefficients of −0.654, −0.608, −0.529, and −0.393 corresponding to 19 October, 4 November, 20 November, 6 December, and 22 December, respectively. It can be concluded that the characteristic of surface property related to both soil moisture and vegetation water content can be noticeably presented by NDLI as a sensitive and reliable indicator to determine the variation of near-surface temperature.

Integrating NDLI and LST for ET Examination
Integration between NDLI and LST provided an advantageous method for monitoring surface moisture and ET rate, in which NDLI was set as a reference basis for examining surface water availability, while the variation of LST is an indicator of the cooling effect by ET at the earth's surface. Figure 4 presents the scatter plot of the dry and wet edges in LST/NDLI triangle in order to determine the surface moisture and ET. The dry and wet borders represent the minimum and maximum of LSTs at the given NDLI intervals, respectively. The edges were obtained by linear functions of NDLI. We found that, although the LST and slopes of the edges at all 5-time points are different, they all tend to be closed as NDLI increases and can be overlapped as NDLI may saturate at later stages. Results show that, in the dry edge, LST max has a significant negative correlation with NDLI with correlation coefficients ranging from −0.88 to −0.96. All dry edges that have negative slopes with the driest points mostly appear close to the straight lines. That is, the linear model provided a decent fit to the data. In the wet edge, relatively small variation of minimum LST was exhibited at different NDLI values. The fluctuations are based on the differences in ET rate from land cover types. There exist some outlier points under special conditions of saturated air.
Based on the dry/wet edge fitting equations and the feature space of LST/NDLI, the TMDI was simulated for all of five days (as presented in Figure 5). Results indicated that temporal variation of TMDI do not change much in the areas of built-up and water surface due to the fact that water availability of these areas does not suffer a significant change during a growing season. At all five time points, relatively high TMDI appeared in the dry surface, such as built-up, bare land areas, and, especially, rock/sand of Zhuoshui River basin as dry season, while low and very low TMDIs were highly concentrated in vegetation and water areas, respectively. It fluctuated as the consequence of the cooling effect by ET due to the fact that ET of water and vegetative areas is considerably higher than the dry surface. Spatial and temporal variations of the TMDI in the cropland exhibited a distinctive change in the apparent property of crop during the growing process. The transformation of land-cover in the farmland during crop season leads to variation of surface water condition, resulting in difference of surface ET. Very low TMDIs were observed in the paddy rice fields, especially during the periods of growth, corresponding to very high value of ET compared to the other land cover types. Relative low TMDIs appear in the farmland during the vegetative period, such as maize, sugarcane, and pineapple, when canopy reaches its full coverage. The lower TMDI in these areas is a consequence of transpiration from crop plants. Meanwhile, high and very high TMDIs of croplands are mainly corresponding to dry farming or the areas at land preparation phase. We found that, over large open water surfaces, dry and wet edges tend to overlap since LST may be equal or lower than the upper-air temperature. Consequently, in some water areas of aquaculture systems, TMDI appears lower than the theoretical values. corresponding to dry farming or the areas at land preparation phase. We found that, over large open water surfaces, dry and wet edges tend to overlap since LST may be equal or lower than the upper-air temperature. Consequently, in some water areas of aquaculture systems, TMDI appears lower than the theoretical values   To assess the sensitiveness of TMDI in terms of determining the surface moisture and rate of surface ET, SEBAL-derived ET ( Figure 6) was used as reference data for comparison to evaluate the performance of TMDI as compared to TVDI. The comparisons were performed against an independent set of reference points chosen randomly. Results showed that both TVDI and TMDI have strong negative correlations with the SEBAL-derived ET (Table 3). Thus, both of them correspond well with surface moisture. Surface water scarcity leads to lower ET and reduces the cooling effect of surface, resulting in increase of TVDI and TMDI. In all five days, higher correlations with SEBAL-derived ET were presented by TMDI compared to TVDI. That is, TMDI more effectively reflects the surface water status and appropriately presents the rate of surface ET than TVDI in the southern parts of Taiwan. In addition, the rapid water supply, as the irrigation or rain occurred before the investigated time-points, can lead to sudden changes in surface water availability, surface ET rate, and LST. The time lag of the surface canopy's response to water supply leads to difficulties in surface property investigation by remote sensing. This study exhibits the superiority of TMDI in reflecting the rapid variation of surface water status. As compared to TVDI, the dramatic improvement in the correlation with SEBALderived ET was presented by TMDI in the time-points with the rain recorded in the few To assess the sensitiveness of TMDI in terms of determining the surface moisture and rate of surface ET, SEBAL-derived ET ( Figure 6) was used as reference data for comparison to evaluate the performance of TMDI as compared to TVDI. The comparisons were performed against an independent set of reference points chosen randomly. Results showed that both TVDI and TMDI have strong negative correlations with the SEBAL-derived ET (Table 3). Thus, both of them correspond well with surface moisture. Surface water scarcity leads to lower ET and reduces the cooling effect of surface, resulting in increase of TVDI and TMDI. In all five days, higher correlations with SEBAL-derived ET were presented by TMDI compared to TVDI. That is, TMDI more effectively reflects the surface water status and appropriately presents the rate of surface ET than TVDI in the southern parts of Taiwan. In addition, the rapid water supply, as the irrigation or rain occurred before the investigated time-points, can lead to sudden changes in surface water availability, surface ET rate, and LST. The time lag of the surface canopy's response to water supply leads to difficulties in surface property investigation by remote sensing. This study exhibits the superiority of TMDI in reflecting the rapid variation of surface water status. As compared to TVDI, the dramatic improvement in the correlation with SEBAL-derived ET was presented by TMDI in the time-points with the rain recorded in the few days before the days of concern.
Resulting in R-value's increase from −0.593 to −0.647, −0.734 to−0.88, and −0.8 to −0.852 is observed on the rain-affected days of 20 November, 6 December, and 22 December, respectively. Meanwhile, slighter increases in correlations with SEBAL-derived ET were presented by TMDI in the time-points with no water supply during the previous days. Resulting in R-values slightly improved from −0.856 to −0.89 and −0.8 to −0.823 is seen on 19 October and 4 November, respectively.

Conclusions
This study aims to characterize surface properties and its contribution to the LST variability in order to propose an advantageous indicator for assessing surface moisture and ET rate. The spatio-temporal LST maps of Yun-Chia-Nan area presented the significant variability in the agriculture areas during the second crop growing season of 2014. It corresponds to different characteristics of the canopies and surface water availability at different crop growing stages. The correlation of thermal variability and surface biophysical properties was subsequently determined via quantitative correlations between LST and multi-spectral satellite indices under investigation, i.e., NDVI, BI, NDMI, NDWI, and NDLI. The highest correlation coefficient was observed between LST and NDLI with the negative correlations of −0.617, −0.654, −0.608, −0.529, and −0.393 corresponding to 19 October, 4 November, 20 November, 6 December, and 22 December, respectively. This shows the superiority of NDLI over the other four indices in representing the variation of near-surface thermal signatures. The LST-NDLI triangle space-based, Temperature-surface Moisture Dryness Index (TMDI), was accordingly proposed for assessing the surface moisture and ET. It is based on the deviation of LST with its min/max LST at the given NDLI intervals. The TMDI exhibits good performance in representing surface moisture and ET by showing strong correlations with SEBAL-derived ET. TMDI also exhibited its improvement compared to TVDI in revealing the signatures of immediate surface water variability by presenting the slight and significant improvements in correlations with SEBAL-derived ET on the days no-water/water supply nearby, respectively. Correlations were found to be significantly increased in the investigated times with rains recorded in the previous days (R-values are increased from −0.593 to −0.647, −0.734 to −0.88, and −0.8 to −0.852 on 20 November, 6 December, and 22 December, respectively). The superiority of TMDI has been proven by its effective responses to the surface water status over the region of diverse land cover types, extreme conditions, and different time points. The use of TMDI gives a significant advantage in quickly investigating the thermal environment variability, taking into account the surface moisture and ET, and providing the information for fast decision-making, especially farming management decisions, such as precision irrigation. With diverse remotely-sensed data available, the TMDI is potentially helpful for a variety of practical applications related to the Earth's surface properties over a large scale since its formula is simple to be easily implemented and does not require any ground-based data. Data Availability Statement: In this study, satellite data were downloaded from United States Geological Survey (USGS) and weather data were collected from Central Weather Bureau (CWB), Taiwan.