Modeling Vegetation Water Stress over the Forest from Space: Temperature Vegetation Water Stress Index (TVWSI)

: The conventional Land Surface Temperature (LST)–Normalized Difference Vegetation Index (NDVI) trapezoid model has been widely used to retrieve vegetation water stress. However, it has two inherent limitations: (1) its complex and computationally intensive parameterization for multi-temporal observations and (2) deﬁciency in canopy water content information. We tested the hypothesis that an improved water stress index could be constructed by the representation of canopy water content information to the LST–NDVI trapezoid model. Therefore, this study proposes a new index that combines three indicators associated with vegetation water stress: canopy temperature through LST, canopy water content through Surface Water Content Index (SWCI), and canopy fractional cover through NDVI in one temporally transferrable index. Firstly, a new optical space of SWCI–NDVI was conceptualized based on the linear physical relationship between shortwave infrared (SWIR) and soil moisture. Secondly, the SWCI–NDVI feature space was parameterized, and an index d(SWCI, NDVI) was computed based on the distribution of the observations in the SWCI–NDVI spectral space. Finally, standardized LST (LST/long term mean of LST) was combined to d(SWCI, NDVI) to give a new water stress index, Temperature Vegetation Water Stress Index (TVWSI). The modeled soil moisture from the Australian Water Resource Assessment—Landscape (AWRA-L) and Soil Water Fraction (SWF) from four FLUXNET sites across Victoria and New South Wales were used to evaluate TVWSI. The index TVWSI exhibited a high correlation with AWRA-L soil moisture (R 2 of 0.71 with p < 0.001) and the ground-based SWF (R 2 of 0.25–0.51 with p < 0.001). TVWSI predicted soil moisture more accurately with RMSE of 21.82 mm (AWRA-L) and 0.02–0.04 (SWF) compared to the RMSE ranging 28.98–36.68 mm (AWRA-L) and 0.03–0.05 (SWF) were obtained for some widely used water stress indices. The TVWSI could also be a useful input parameter for other environmental models.


Introduction
Large variations in rainfall and temperature may push forests into a non-recoverable stage of water stress [1][2][3]. In Australia, several studies documented the severe effects of long droughts or heat waves on the ecosystem in the form of canopy decline and mortality [4][5][6][7][8][9].
Eucalypts are the dominant forest type in Australia. Although many eucalypt species are adapted to seasonal water limitation, extreme droughts impose water stress in many areas, and hot summers can cause seasonal water deficits in non-drought years. Vegetation water stress is a significant factor when defining ecosystem-level health and has been used as a critical indicator of environmental quality [10,11]. This indicator can vary significantly Despite the wide acceptance of the LST-NDVI feature space in estimating soil moisture, it has some limitations. Improvements in parameterization were proposed to enhance LST-NDVI capability. For example, Rahimzadeh-Bajgiran et al. [17] considered a non-linear relationship between LST and soil moisture content rather than considering it linear, which was the original concept of LST-NDVI. Sun [42] proposed a two-stage trapezoid because the LST of a vegetative surface is less responsive to surface soil moisture than the bare surface, as trees can extract deep soil water, making the LST-soil moisture relationship non-linear. Tao et al. [43] identified a problem in the parameterization of the dry edge due to a mismatch of the spatial scale between visible and thermal bands. They proposed a new index by substituting LST by Modified Perpendicular Dryness Index (MPDI) in the LST-NDVI feature space. Ghulam et al. [44] substituted LST with surface broadband albedo in the near-infrared band to propose a vegetation condition albedo drought index (VCADI). However, one of the significant limitations of the LST-NDVI feature space is its complex and computationally intensive parameterization for temporal observations. It must be parameterized separately for each day due to LST's dependency on the region-specific properties such as canopy structure and the ambient atmospheric conditions such as near-surface air temperature, relative humidity, and wind speed, which are temporally variable at a location [32].
Unlike LST, the surface reflectance does not significantly vary with the ambient atmospheric condition. Therefore, the optical reflectance can be employed with temporal observations without parameterizing the model separately for each date. Studies have used optical spectral spaces to estimate soil moisture and drought conditions (eg. Ghulam et al. [45] and Ghulam et al. [46]. These authors developed the Perpendicular Dryness Index (PDI) and Modified PDI (MPDI), respectively. Amani et al. [14] proposed the Triangular Soil Moisture Index (TSMI) using Near Infrared (NIR)-Red triangular space. Alternatively, some authors have used SWIR instead of Red in NIR-Red spectral space, such as the Shortwave-infrared Perpendicular Water Stress Index (SPSI) proposed by Ghulam et al. [47].
Although LST's dependency on the ambient atmospheric condition makes it computationally expensive to calculate model parameters, it is a subtle measure of vegetation water stress regulated by soil moisture dynamics, which cannot be realized by optical reflection. One alternative is to construct a feature space using optical reflectance, retrieve model parameters out of this feature space, and then integrate standardized LST as was Unlike LST, the surface reflectance does not significantly vary with the ambient atmospheric condition. Therefore, the optical reflectance can be employed with temporal observations without parameterizing the model separately for each date. Studies have used optical spectral spaces to estimate soil moisture and drought conditions (e.g., Ghulam et al. [45] and Ghulam et al. [46]. These authors developed the Perpendicular Dryness Index (PDI) and Modified PDI (MPDI), respectively. Amani et al. [14] proposed the Triangular Soil Moisture Index (TSMI) using Near Infrared (NIR)-Red triangular space. Alternatively, some authors have used SWIR instead of Red in NIR-Red spectral space, such as the Shortwave-infrared Perpendicular Water Stress Index (SPSI) proposed by Ghulam et al. [47].
Although LST's dependency on the ambient atmospheric condition makes it computationally expensive to calculate model parameters, it is a subtle measure of vegetation water stress regulated by soil moisture dynamics, which cannot be realized by optical reflection. One alternative is to construct a feature space using optical reflectance, retrieve model parameters out of this feature space, and then integrate standardized LST as was performed by Wu and Lu [37] in modeling of Modified Vegetation Water Supply Index (MVWSI). The standardized LST is the ratio of the daily LST to its long-term mean for that time and region. The literature suggests that using standardized LST in the model of MVWSI has eliminated much of the differences caused by time and space; hence, we assume an index with mathematical similarity with MVWSI will also be independent of ambient atmospheric conditions. The Short-Wave Infrared (SWIR) region of the electromagnetic (EM) spectrum is sensitive to leaf water content and has been widely used to estimate canopy water content [48][49][50][51][52]. In this study, we aim to improve the soil moisture estimation over vegetative soil by combining all the three significant indications of water stress: vegetation transpiration status, vegetation water content, and fractional vegetation cover, into one temporally transferrable index. Previous models did not include these three significant indicators together, and they were computationally expensive to parameterize for temporal observations.
The primary objectives of this work are: (1) Test the hypothesis that an improved vegetation water stress index could be constructed by the representation of canopy water Remote Sens. 2021, 13, 4635 4 of 23 content information to the LST-NDVI trapezoid model, thus augmenting the conventional LST-NDVI feature space to include SWIR that accounts for leaf/canopy water status; (2) Derive a new temporally transferrable vegetation water stress index from the newly constructed feature space of SWIR-NDVI-LST; (3) Evaluate the performance of the proposed index against the widely used vegetation water stress indices such as PDI, MPDI, MVWSI, TVDI, and Crop Water Stress Index (CWSI).

Study Area
The state of Victoria lies in the Southeastern part of Australia and has a variable climate from extremely drier western regions to wet eastern regions. The annual rainfall and temperatures range from 400 mm, 20 • C in western parts, to >1500 mm, 5 • C in eastern parts of the state [53]. The Koppen-Geiger climate type in Victoria is mainly classified as temperate and arid. The eastern and coastal regions represent a temperate climate, while the North-Western parts represent an arid climate [54]. From November to March, the temperature remains high, while from May to September, temperatures are low, and the rainfall during this period accounts for 80% of the rainfall during the whole year.
The study area encompasses a mix of eucalypt forest types, from wet montane forests in the east to sparse dry woodlands in northwest Victoria. The Ecological Vegetation Classes (EVC) is widely used as a landscape-scale classification scheme for Victoria's vegetation types, classifying vegetation into 20 simplified native vegetation groups and 35 sub-groups [55]. These vegetation groups span a wide range of forest structural formations from low woodland to open forest and tall closed forest [56]. For the selection of the study forest patches, hereafter study sampling windows, the whole of the forest region of Victoria was sampled into spatial windows measuring (2 km × 2 km). The study sampling windows can accommodate 16 (4 × 4) pixels of MODIS-optical and 4 (2 × 2) pixels of MODIS-LST, as shown in Figure 2. The reason for selecting a study sampling window size of (2 km × 2 km) is three-fold: (1) To remove the artifacts caused by the spatial discrepancy between MODIS-optical (available at 500 m) and MODIS-LST (available at 1000 m); (2) To make spatial resolution comparable, as the modeled soil moisture data from AWRA-L, along with other meteorological data used in this study such as air temperature and vapor pressure deficit are available at 5 km spatial resolution; (3) The FLUXNET spatial location was mostly falling in the corners of the MODIS-LST/MODIS-optical pixels, so the averaging of the neighboring pixels were required.
A series of spatial analyses were executed over these study sampling windows to remove sampling windows that may have confounding data, based on the following criteria: (1) fire history (not burnt after 1990), (2) logging history (not logged after 1990), (3) continuity of forest patch (all pixels inside the sampling window should fully overlap the forest area), (4) distance from the coast (more than 30 km from the coast to avoid the continuous humid climate and probably unlimited moisture accessibility), and (5) proximity from rain gauge stations (ten or more Australian Bureau of Meteorology (BoM) rain gauge station in a radius of 30 km). As shown in Figure 2, finally, sixty study sampling windows were found suitable for historical time series analysis at the ecological level, covering large variations in aridity and vegetation densities from different species of Eucalypt.  [57]. This multispectral dataset is produced daily using a 16-day retrieval period available at a 500 m spatial resolution which combines data from Terra and Aqua. The LST data were extracted from MODIS dataset MYD11A1.006 [58], which provides daily Land Surface Temperature and Emissivity with 1000 m spatial resolution. The temperature values are derived from the MOD11_L2 swath product. It provides both daytime (passes around 1:30 pm local time) and nighttime (passes around 1:30 am in Australia) values accompanied by a quality layer for each. This study used daytime LST values.

Datasets
the study [57]. This multispectral dataset is produced daily using a 16-day retrieval period available at a 500 m spatial resolution which combines data from Terra and Aqua. The LST data were extracted from MODIS dataset MYD11A1.006 [58], which provides daily Land Surface Temperature and Emissivity with 1000 m spatial resolution. The temperature values are derived from the MOD11_L2 swath product. It provides both daytime (passes around 1:30 pm local time) and nighttime (passes around 1:30 am in Australia) values accompanied by a quality layer for each. This study used daytime LST values.
The datasets are publicly available and can be downloaded through NASA Land Processes Distributed Active Archive Center (LP DAAC, see https://lpdaac.usgs.gov, accessed on 29 August 2021). The data were downloaded in bulk from NASA's FTP server (LP DAAC FTP server, see https://e4ftl01.cr.usgs.gov/, accessed on 29 August 2021).  The datasets are publicly available and can be downloaded through NASA Land Processes Distributed Active Archive Center (LP DAAC, see https://lpdaac.usgs.gov, accessed on 29 August 2021). The data were downloaded in bulk from NASA's FTP server (LP DAAC FTP server, see https://e4ftl01.cr.usgs.gov/, accessed on 29 August 2021).

Meteorological and Model Soil Moisture Data
The modeled daily soil moisture and meteorological datasets (temperature and rainfall) and Vapour Pressure Deficit (VPD) grids at a spatial resolution of 0.05 • by 0.05 • were obtained from the Australian Bureau of Meteorology (BoM) for a period of 17 years (2002-2018). The daily modeled root zone soil moisture (top 1 m) from the AWRA-L model was used to develop an index and perform meta ground-truthing. AWRA-L (Australian Water Resource Assessment-Landscape) is a landscape model developed by BoM and Commonwealth Scientific and Industrial Research Organisation (CSIRO). The model runs on a daily timestep and simulates the landscape water balance for Australia from 1911 to yesterday [59]. The key outputs from the model include surface runoff, soil moisture, evapotranspiration, and deep drainage [60]. The model is validated and calibrated at distributed spatial locations using ground observed data from point-scale soil moisture probe data, flux tower estimates, and ground recharge estimates. The root zone soil moisture from BoM-AWRA-L is the sum of water in upper and lower soil layers. It represents the percentage of available water content in the top 1 m of the soil profile. Each data cell is independent of its neighbors and represents different water compartments [59,61]. The AWRA-L version 6 gives improved results, particularly for soil moisture and runoff compared to its previous version 5.

In Situ Soil Moisture Data
The study also used measured soil moisture data for the top 85 cm of soil from four OzFlux sites located at Wombat, Wallaby, and Whroo in Victoria, and Tumbarumba in New South Wales, Australia, to validate the proposed index [67,68]. The Level 3 datasets were downloaded from the OzFlux portal (OzFlux, see http://data.ozflux.org.au/portal/, accessed on 29 August 2021).
The University of Melbourne maintains the Wombat monitoring site, which has been operational since late January 2010. The site is in a cool temperate dry sclerophyllous eucalypt forest in Central Victoria. The annual rainfall is approximately 871 mm (142year long term average), and the temperature ranges (1-30 • C) with a mean annual air temperature (2001-2012) of 12.1 • C [69]. Wallaby flux station is located in Kinglake National Park, approximately 85 km north of Melbourne. The site's elevation is 720 m, and the site is dominated by Mountain ash (Eucalyptus regnans). In 2009, the site was damaged by a bushfire and was re-established again in 2010. The climate is cool and temperate with an annual rainfall ranging 1100-2000 mm, and the data are available from 2005 onwards [70]. The Whroo site is located approximately 160 km north of Melbourne at the lowest elevation among all the sites, with 165 m. The mean annual precipitation as recorded by the Bureau of Meteorology at the site is 558 mm. The site is classified as the box woodland, and the data are available from 2011 to the present [71]. The Tumbarumba flux station is situated in the Bago state forest site in southeastern New South Wales. CSIRO Land and water maintain this site, and the data are available from 2001 to the present. The forest is wet sclerophyll with a mean annual rainfall of 1159 mm [67,72].

Conceptual Model: The Connection between Canopy Temperature, Canopy Water Content, and Canopy Fraction during Soil Moisture Deficits
LST can be used to characterize the increase in canopy temperatures following stomatal closure in response to water supply limitations. The magnitude of change in the canopy temperature captured by satellite sensors is also dependent on the fractional ground cover (which recognizes the proportion of evaporation and transpiration fluxes, photosynthetic activities, and canopy resistance). A key assumption is that canopy water content and canopy temperature are inversely related.
Our proposed index is presumably tested for Eucalypts water stress reactions due to limited access to soil moisture in their root zones. We assume that a decline in rootzone soil moisture will alter the stomata functioning to save water loss, accompanied by a decrease in canopy water content; further, a persistent reduction in soil moisture can alter the canopy fraction. Thus, the proposed index is designed to measure these three physiological responses from the plant under soil moisture deficit conditions in the root zone. Please note that the primary source of stress is from the root-zone soil moisture.
The proposed index (TVWSI) is a function of canopy reflection from thermal (to infer transpiration status), Visible Near Infrared (to infer vegetation fraction/photosynthetic activity), and Short-Wave Infrared (to infer vegetation water content) regions of the EM spectrum. Figure 3 shows the conceptual diagram of the roles of VisNIR, thermal, and SWIR in identifying water stress in the vegetation.
matal closure in response to water supply limitations. The magnitude of change in the canopy temperature captured by satellite sensors is also dependent on the fractional ground cover (which recognizes the proportion of evaporation and transpiration fluxes, photosynthetic activities, and canopy resistance). A key assumption is that canopy water content and canopy temperature are inversely related.
Our proposed index is presumably tested for Eucalypts water stress reactions due to limited access to soil moisture in their root zones. We assume that a decline in root-zone soil moisture will alter the stomata functioning to save water loss, accompanied by a decrease in canopy water content; further, a persistent reduction in soil moisture can alter the canopy fraction. Thus, the proposed index is designed to measure these three physiological responses from the plant under soil moisture deficit conditions in the root zone. Please note that the primary source of stress is from the root-zone soil moisture.
The proposed index (TVWSI) is a function of canopy reflection from thermal (to infer transpiration status), Visible Near Infrared (to infer vegetation fraction/photosynthetic activity), and Short-Wave Infrared (to infer vegetation water content) regions of the EM spectrum. Figure 3 shows the conceptual diagram of the roles of VisNIR, thermal, and SWIR in identifying water stress in the vegetation.

Temperature Vegetation Water Stress Index (TVWSI)
The modeling of TVWSI is described under two sections: Section 3.2.1 describes the basis for selecting the SWIR band index Surface Water Content Index (SWCI) and conceptualizing a new optical spectral space constructed using SWCI-NDVI. In addition, SWCI-NDVI was parameterized, and a new intermediate index, d(SWCI, NDVI), was devised using NDVI and SWCI spectral space. Further to this, Section 3.2.2 describes the modeling

Temperature Vegetation Water Stress Index (TVWSI)
The modeling of TVWSI is described under two sections: Section 3.2.1 describes the basis for selecting the SWIR band index Surface Water Content Index (SWCI) and conceptualizing a new optical spectral space constructed using SWCI-NDVI. In addition, SWCI-NDVI was parameterized, and a new intermediate index, d(SWCI, NDVI), was devised using NDVI and SWCI spectral space. Further to this, Section 3.2.2 describes the modeling of the TVWSI, using the modeling concept of MVWSI, by replacing NDVI with d(SWCI, NDVI).

SWCI-NDVI Spectral Space and Derivation of d(SWCI, NDVI)
In developing our index, atmospherically stable multi-SWIR channel indices were given preference over a single SWIR channel to deal with errors from atmospheric water vapor [73,74]. SWCI was selected as the candidate SWIR index to augment NDVI. SWCI was developed to estimate vegetation canopy water content dynamics regulated by changes in soil moisture content, and it was used by numerous studies to estimate surface soil moisture [75,76]. The metric is the normalized ratio between SWIR (B 6 ) and SWIR (B 7 ) of MODIS, where B n is MODIS band n and R m represents the reflectance of waveband with a wavelength m in µm.
SWCI uses bands near the absorption valley and peak of the water spectra; therefore, the difference of these bands should increase with vegetation water content. As both bands have the same atmospheric scattering and radiation values, taking the normalized ratio reduces the atmosphere's effect on their individual bands [77,78].
The Normalized Difference Vegetation Index (NDVI) uses Red and NIR channels from the EM spectrum [79]. It is widely used to estimate vegetation growth parameters such as LAI, greenness, photosynthetic activeness, vigor, and phenology. MODIS Band (B 1 ) and Band (B 2 ) from the MODIS were used to calculate NDVI, R represents the radiance in reflectance units.
We assume a linear relationship between root zone soil moisture and vegetation canopy water content. This assumption conforms to previous studies that have used remotely sensed vegetation indices to quantify plant vigor and relate it to root zone soil moisture [80][81][82][83][84][85]. As for LST-NDVI, a feature space of SWCI-NDVI can be constructed, as shown in Figure 4. This figure is adopted from Sadeghi et al. [86], who used a similar feature space constructed using shortwave infrared transformed reflectance (STR) -NDVI to estimate root-zone soil moisture. It is to be noted that the reflection in the SWIR region is sensitive to canopy water content and the pooled water on the surface. The meaningful relationship between root-zone soil moisture and SWCI is valid for partially or fully saturated soil but not beyond this state (i.e., when water accumulates on the soil surface). Although SWCI will continue to increase with standing water, that does not mean that soil moisture is increasing, as it is bound to a certain maximum limit (saturation). This issue was pointed out in defining the wet edge in optical spectral space; hence, there must be some assumptions for the upper SWCI values corresponding to the maximum limit of soil moisture, shown by (?) in Figure 4 (oversaturated edge). Therefore, we only used the dry edge to remove this ambiguity while defining the new index.  Let us suppose the equation of the dry edge is given by where m is slope and c is the intercept of the line. Then a perpendicular line to dry edge, passing through A ( , ) can be given as the perpendicular distance AB from a point A ( , ) to dry edge can be derived as Therefore, the perpendicular distance corresponding to the ith observation of NDVI and SWCI, denoted as d(SWCI , NDVI ), was calculated as SWCI − NDVI − Figure 4. Trapezoidal spectral space, assumed to be between SWCI and NDVI. The figure is adopted from Sadeghi et al. [86].
Suppose a hypothetical dry line can be drawn in the spectral space of SWCI-NDVI, as shown in Figure 4. In that case, the perpendicular distance of observations to this line can represent soil moisture content, i.e., the distance of observation to this line will increase with soil water content. The dryline can be estimated by fitting a linear regression line on the lowest SWCI from NDVI bins across the x-axis. The number of bins in the NDVI range can be estimated using the method proposed by Sturges [87]. The method states that for a series of range R with N items, the optimal class interval (bins) may be estimated from the formula: class interval (bins) = R 1 + 3.322 log 10 N Let us suppose the equation of the dry edge is given by where m is slope and c is the intercept of the line. Then a perpendicular line to dry edge, passing through A (x i , y i ) can be given as the perpendicular distance AB from a point A (x i , y i ) to dry edge can be derived as Therefore, the perpendicular distance corresponding to the ith observation of NDVI and SWCI, denoted as d(SWCI i , NDVI i ), was calculated as This perpendicular distance may account for the effect of vegetation fraction (due to NDVI) and canopy water content (due to SWCI) concurrently and can be linked more efficiently than NDVI and SWCI to the ground soil moisture conditions.

Derivation of TVWSI
TVWSI is based on the theoretical principle of the modified version of the simple ratio index VWSI (NDVI/LST) [21,35,36]. It is the most straightforward derivation of the slope from the LST-NDVI feature space. Numerous studies have found that LST-NDVI slope dynamics are sensitive to plant-available soil moisture and plant water stress [17,19,21,25,88], while other studies have found a significant relation between VWSI and vegetation water supply, which can directly link to soil moisture conditions, transpiration regimes, and vegetation water stress [36]. A study by McVicar and Bierwirth [41] over Papua New Guinea suggested that the ratio of NDVI and LST is a direct measure of the drought.
The major limitation in using LST on spatio-temporal data is its dependency on the region-specific properties such as canopy structure and ambient atmospheric conditions such as air temperature and wind speed. However, a study by Wu and Lu [37] found that the relative LST, i.e., the real-time LST divided by long term mean of LST corresponding to that time and region, can solve this limitation. Consequently, we employed relative LST instead of LST to construct a spatially transferrable vegetation water stress index.
The standardized LST or relative LST denoted by RLST is the ratio of LST with its long-term mean (LST). This mathematical conversion of LST makes its application feasible in time and space.
TVWSI is devised by employing d(SWCI, NDVI) in place of NDVI in the model of MVWSI and is given by: Retrieving changes in soil moisture using TVWSI will depend on the dynamics of d(SWCI, NDVI) and LST during varying soil moisture conditions. Under normal soil moisture conditions, both d(SWCI, NDVI) and LST are expected to vary in some reasonable range; d(SWCI, NDVI) should have relatively higher values, and LST should have relatively smaller values. Overall, TVWSI will have a relatively high value for high soil moisture conditions. In contrast, under water-limited conditions, LST increases due to less evaporative cooling, and correspondingly, d(SWCI, NDVI) decreases due to a reduction in water content and vegetation growth, resulting in a decrease in the slope of SWCI-NDVI spectral space.

Evaluation of TVWSI
A comparative analysis was performed between TVWSI and existing vegetation water stress indices, based on their skills to predict AWRA-L soil moisture from sixty sampling study windows and observed Soil Water Fraction (SWF) from four FLUXNET sites. First-, second-, and third-order polynomial models were fitted between the indices and AWRA-L soil moisture/SWF. Then, Root Mean Square Error (RMSE) (predicted AWRA-L soil moisture/SWF vs. actual AWRA-L soil moisture/SWF), Mean Absolute Percentage Error (MAPE) (predicted AWRA-L soil moisture/SWF vs. actual AWRA-L soil moisture/SWF), and R 2 (predicted AWRA-L soil moisture/SWF vs. actual AWRA-L soil moisture/SWF) were used to compare the skills of the water stress metrics in estimating soil moisture. Table 1 shows the existing water stress indices against which TVWSI was evaluated.

Indices
References The dataflow diagram describing the construction and evaluation of TVWSI is shown in Figure 5. The dataflow diagram describing the construction and evaluation of TVWSI is shown in Figure 5.

Results
The proposed indices were evaluated using 17-year (2002-2018) weekly AWRA-L soil moisture from sixty study sampling windows and Soil Water Fraction (SWF) from the four FLUXNET stations. The daily data were converted into a weekly temporal resolution to deal with outliers or errors in the daily datasets. The weekly maximum NDVI, SWCI and mean LST, AWRA-L soil moisture, and SWF were considered for the analysis. The analysis was confined to the summer season (December to February) as these are the months are when vegetation may be exposed to soil moisture deficits and more likely to experience water stress. We assume that root zone soil moisture (moisture in the top 1 m of soil) is the primary driver of vegetation water stress. Separate regression analyses were performed between the proposed and existing water stress indices (that we consider are

Results
The proposed indices were evaluated using 17-year (2002-2018) weekly AWRA-L soil moisture from sixty study sampling windows and Soil Water Fraction (SWF) from the four FLUXNET stations. The daily data were converted into a weekly temporal resolution to deal with outliers or errors in the daily datasets. The weekly maximum NDVI, SWCI and mean LST, AWRA-L soil moisture, and SWF were considered for the analysis. The analysis was confined to the summer season (December to February) as these are the months are when vegetation may be exposed to soil moisture deficits and more likely to experience water stress. We assume that root zone soil moisture (moisture in the top 1 m of soil) is the primary driver of vegetation water stress. Separate regression analyses were performed between the proposed and existing water stress indices (that we consider are the remotely sensed indicator of vegetation water stress) and AWRA-L soil moisture. The comparison of proposed indices was made with the existing vegetation water stress indices, presented, and discussed in the following sections.

Estimation of d(SWCI, NDVI) from SWCI-NDVI Spectral Space and Analysis of Correlation with AWRA-L Soil Moisture
To calculate the dry edge in SWCI-NDVI spectral space, the NDVI were classified into different bins across the x-axis. The number of bins was calculated using the formula proposed by Sturges [87]. Then, in each bin of the NDVI, the lowest SWCI was identified, and a linear regression line fitted to these points ( Figure 6). Finally, the distance d(SWCI, NDVI) of each observation to this lowest slope line was calculated by the method described in Section 3.2. proposed by Sturges [87]. Then, in each bin of the NDVI, the lowest SWCI was identified, and a linear regression line fitted to these points ( Figure 6). Finally, the distance d(SWCI, NDVI) of each observation to this lowest slope line was calculated by the method described in Section 3.2.  Figure 7 shows the density scatter plots of NDVI, SWCI, and d(SWCI, NDVI) vs. AWRA-L soil moisture with their corresponding R 2 value. The shape of the plot and R 2 score suggest that NDVI, SWCI, and d(SWCI, NDVI) have a positive correlation with AWRA-L soil moisture, and d(SWCI, NDVI) has significantly improved the correlation and the linearity of the relationship compared with NDVI and SWCI. In addition, Figure  7c confirms our assumption that the distance of an observation relative to the lowest slope line in NDVI-SWCI spectral space increases with soil moisture.  Figure 7 shows the density scatter plots of NDVI, SWCI, and d(SWCI, NDVI) vs. AWRA-L soil moisture with their corresponding R 2 value. The shape of the plot and R 2 score suggest that NDVI, SWCI, and d(SWCI, NDVI) have a positive correlation with AWRA-L soil moisture, and d(SWCI, NDVI) has significantly improved the correlation and the linearity of the relationship compared with NDVI and SWCI. In addition, Figure 7c confirms our assumption that the distance of an observation relative to the lowest slope line in NDVI-SWCI spectral space increases with soil moisture. the color represents the number of observations per pixel, the red color points are the lowest observed SWCI in the bins of NDVI. The line (Dry edge) is the best-fitted regression line to these lowest SWCI points. The equation of the line is shown in the figure. Figure 7 shows the density scatter plots of NDVI, SWCI, and d(SWCI, NDVI) vs. AWRA-L soil moisture with their corresponding R 2 value. The shape of the plot and R 2 score suggest that NDVI, SWCI, and d(SWCI, NDVI) have a positive correlation with AWRA-L soil moisture, and d(SWCI, NDVI) has significantly improved the correlation and the linearity of the relationship compared with NDVI and SWCI. In addition, Figure  7c confirms our assumption that the distance of an observation relative to the lowest slope line in NDVI-SWCI spectral space increases with soil moisture.  In order to analyze the improvements d(SWCI, NDVI) bought over NDVI and SWCI, we classified the data from sixty study sampling windows into five classes representing 0-20, 20-40, 40-60, 60-80, and 80-100 soil moisture percentile. Then, the ranges of each of the indices NDVI, SWCI, and d(SWCI, NDVI) was estimated for each of these soil moisture percentile classes, shown by the boxplots in Figure 8. These boxplots suggest d(SWCI, NDVI) improved the skills of NDVI and SWCI in separating different ranges of soil moisture, especially at the end members of soil moisture, i.e., at lower (0-20) and upper (60-100) ranges (Figure 8c). These results suggest the distance metric may have higher skill for retrieving soil moisture dynamics in drier and wetter zones where other vegetation indices such as NDVI and SWCI generally saturates or is contaminated by background noise [90][91][92]. In order to analyze the improvements d(SWCI, NDVI) bought over NDVI and SWCI, we classified the data from sixty study sampling windows into five classes representing 0-20, 20-40, 40-60, 60-80, and 80-100 soil moisture percentile. Then, the ranges of each of the indices NDVI, SWCI, and d(SWCI, NDVI) was estimated for each of these soil moisture percentile classes, shown by the boxplots in Figure 8. These boxplots suggest d(SWCI, NDVI) improved the skills of NDVI and SWCI in separating different ranges of soil moisture, especially at the end members of soil moisture, i.e., at lower (0-20) and upper (60-100) ranges (Figure 8c). These results suggest the distance metric may have higher skill for retrieving soil moisture dynamics in drier and wetter zones where other vegetation indices such as NDVI and SWCI generally saturates or is contaminated by background noise [90][91][92].

Estimation of TVWSI and Analysis of Correlation with AWRA-L Soil Moisture
The construction of TVWSI is based on the LST-NDVI slope relationship with soil moisture. The simpler model of the LST-NDVI slope is their ratio (NDVI/LST). The principle behind the working of this ratio is the energy partitioning between vegetative and non-vegetative surfaces, which NDVI can identify for LST. In order to overcome the limitation of LST dependency on ambient atmospheric conditions, we employed the method proposed by Wu and Lu [37], who used standardized LST instead of LST. Further, we replaced NDVI with d (SWCI, NDVI) in the model of MVWSI to construct TVWSI.

Estimation of TVWSI and Analysis of Correlation with AWRA-L Soil Moisture
The construction of TVWSI is based on the LST-NDVI slope relationship with soil moisture. The simpler model of the LST-NDVI slope is their ratio (NDVI/LST). The principle behind the working of this ratio is the energy partitioning between vegetative and non-vegetative surfaces, which NDVI can identify for LST. In order to overcome the limitation of LST dependency on ambient atmospheric conditions, we employed the method proposed by Wu and Lu [37], who used standardized LST instead of LST. Further, we replaced NDVI with d (SWCI, NDVI) in the model of MVWSI to construct TVWSI. TVWSI showed the highest correlation with AWRA-L soil moisture than its input parameters: NDVI, SWCI, d(SWCI, NDVI), and LST. The coefficient of determination (R 2 ) of AWRA-L soil moisture vs. TVWSI, NDVI, SWCI, d(SWCI, NDVI), and LST are 0.71, 0.49, 0.42, 0.62, and 0.51, respectively. Figure 9 shows the density scatter plot of TVWSI vs. AWRA-L soil moisture. The locally estimated scatterplot smoothing (LOESS) shown in Figure 9 suggests that TVWSI has the sensitivity to estimate both low and high soil moisture events, with a slight tendency to underestimate soil moisture after around 75 mm. In comparison, the interannual variability in TVWSI also suggests its capability to identify these extremes, i.e., lowest and highest rainfall years. For example, the summer of 2006 was the driest year during the millennium drought (1997)(1998)(1999)(2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009) in SE Australia; TVWSI is relatively low in this time and has identified this driest year efficiently ( Figure S1). In contrast, the millennium drought ended with massive flood events in 2010; TVWSI has also identified this wet event efficiently, which is reflected by the highest value of TVWSI in 2010 ( Figure S1).

Evaluation of TVWSI in Predicting AWRA-L Soil Moisture, and Comparison with Existing Water Stress Metrics
This section presents a comparative analysis between TVWSI and other existing water stress indices in predicting AWRA-L soil moisture. The first-, second-, and thirdorder polynomial models were developed for TVWSI, PDI, MPDI, MVWSI, TVDI, and CWSI, to predict AWRA-L soil moisture. R 2 and RMSE between actual vs. predicted AWRA-L soil moisture were used to compare the performances of these indices. The seventeen-year weekly data from sixty study sampling windows were randomly divided into 70% and 30% sample sizes. The sample size comprising 70% of data was used to train the model, while the remaining 30% was used to predict AWRA-L soil moisture. This procedure was repeated 10,000 times on the different sets of randomly selected training and testing data samples. The R 2 and RMSE reported are the averages of the outcomes of these iterations. Boxplots show the range of R 2 and RMSE obtained by running 10,000 iterations in Figure S2. Figure 10 shows the density scatter plot (Actual AWRA-L soil moisture vs. Predicted AWRA-L soil moisture) with the model equation, R 2 and RMSE, corresponding to the first-order polynomial model for each index. In comparison, the model performance measures (R 2 and RMSE) from the second-and third-order polyno- The locally estimated scatterplot smoothing (LOESS) shown in Figure 9 suggests that TVWSI has the sensitivity to estimate both low and high soil moisture events, with a slight tendency to underestimate soil moisture after around 75 mm. In comparison, the interannual variability in TVWSI also suggests its capability to identify these extremes, i.e., lowest and highest rainfall years. For example, the summer of 2006 was the driest year during the millennium drought (1997-2009) in SE Australia; TVWSI is relatively low in this time and has identified this driest year efficiently ( Figure S1). In contrast, the millennium drought ended with massive flood events in 2010; TVWSI has also identified this wet event efficiently, which is reflected by the highest value of TVWSI in 2010 ( Figure S1).

Evaluation of TVWSI in Predicting AWRA-L Soil Moisture, and Comparison with Existing Water Stress Metrics
This section presents a comparative analysis between TVWSI and other existing water stress indices in predicting AWRA-L soil moisture. The first-, second-, and third-order polynomial models were developed for TVWSI, PDI, MPDI, MVWSI, TVDI, and CWSI, to predict AWRA-L soil moisture. R 2 and RMSE between actual vs. predicted AWRA-L soil moisture were used to compare the performances of these indices. The seventeen-year weekly data from sixty study sampling windows were randomly divided into 70% and 30% sample sizes. The sample size comprising 70% of data was used to train the model, while the remaining 30% was used to predict AWRA-L soil moisture. This procedure was repeated 10,000 times on the different sets of randomly selected training and testing data samples. The R 2 and RMSE reported are the averages of the outcomes of these iterations. Boxplots show the range of R 2 and RMSE obtained by running 10,000 iterations in Figure S2. Figure 10 shows the density scatter plot (Actual AWRA-L soil moisture vs. Predicted AWRA-L soil moisture) with the model equation, R 2 and RMSE, corresponding to the first-order polynomial model for each index. In comparison, the model performance measures (R 2 and RMSE) from the second-and third-order polynomial models are shown in Table 2.
Remote Sens. 2021, 13, x FOR PEER REVIEW 15 of 24 The soil line in the NIR-Red and the dry edge in the LST-NDVI feature space were estimated to calculate PDI, MPDI, and TVDI from NIR-Red and LST-NDVI feature spaces.First, the values of Red-reflectance were divided into the class intervals to generate a soil line in the NIR-Red spectral space. The number of intervals was defined by the method proposed in Sturges [87]. Similarly, NDVI values were divided into classes to generate the dry edge in the LST-NDVI feature space; then, PDI and MPDI were extracted from NIR-Red and TVDI were computed using the LST-NDVI feature space.  Table 2. R 2 (index vs. AWRA-L soil moisture) and RMSE (Actual vs. Predicted AWRA-L soil moisture) corresponding to the first-, second-, and third-order polynomial model fitted to AWRA-L soil moisture vs. TVWSI, MVWSI, CWSI, PDI, MPDI, and TVDI.   The soil line in the NIR-Red and the dry edge in the LST-NDVI feature space were estimated to calculate PDI, MPDI, and TVDI from NIR-Red and LST-NDVI feature spaces.First, the values of Red-reflectance were divided into the class intervals to generate a soil line in the NIR-Red spectral space. The number of intervals was defined by the method proposed in Sturges [87]. Similarly, NDVI values were divided into classes to generate the dry edge in the LST-NDVI feature space; then, PDI and MPDI were extracted from NIR-Red and TVDI were computed using the LST-NDVI feature space. Figure 10 and Table 2 clearly show that TVWSI (R 2 = 0.71, RMSE = 21.82) has the highest linear correlation with AWRA-L soil moisture. TVDI performed worst in predicting soil moisture, followed by MPDI and PDI. These three indices lost sensitivity at high soil moisture, which can be observed from the flattening of scatter plots at the higher soil moisture, as shown in Figure 10d-f. Apart from MVWSI, none of the indices has shown any noticeable enhancement in predicting AWRA-L soil moisture with higher-order polynomial models.

VI's Actual vs. Predicted AWRA-L Soil
The prediction abilities of indices to AWRA-L soil moisture were also estimated across each of the sixty individual study sampling windows. The R 2 and mean absolute percentage error (MAPE) ranges corresponding to each water stress index are shown by the boxplot in Figure 11. Notably, for PDI, MPDI, TVDI, MVWSI, and CWSI, many sampling windows showed weak or no correlation with R 2 as low as 0 for PDI, MPDI, TVDI, CWSI, and 0.1 for MVWSI. This weakening of relationship at individual study sampling window suggests that large observed data points are required to parameterize these indices. On the other hand, the MAPE is as high as greater than 300% for PDI, MPDI, and CWSI and greater than 250% for MVWSI and TVDI. The performance of TVWSI was relatively better across the sampling study windows, with the lowest R 2 of 0.31 and the highest MAPE of 145%, which is significantly less than other indices.
Remote Sens. 2021, 13, x FOR PEER REVIEW 16 of 24 Figure 10 and Table 2 clearly show that TVWSI (R 2 = 0.71, RMSE = 21.82) has the highest linear correlation with AWRA-L soil moisture. TVDI performed worst in predicting soil moisture, followed by MPDI and PDI. These three indices lost sensitivity at high soil moisture, which can be observed from the flattening of scatter plots at the higher soil moisture, as shown in Figure 10d-f. Apart from MVWSI, none of the indices has shown any noticeable enhancement in predicting AWRA-L soil moisture with higher-order polynomial models.
The prediction abilities of indices to AWRA-L soil moisture were also estimated across each of the sixty individual study sampling windows. The R 2 and mean absolute percentage error (MAPE) ranges corresponding to each water stress index are shown by the boxplot in Figure 11. Notably, for PDI, MPDI, TVDI, MVWSI, and CWSI, many sampling windows showed weak or no correlation with R 2 as low as 0 for PDI, MPDI, TVDI, CWSI, and 0.1 for MVWSI. This weakening of relationship at individual study sampling window suggests that large observed data points are required to parameterize these indices. On the other hand, the MAPE is as high as greater than 300% for PDI, MPDI, and CWSI and greater than 250% for MVWSI and TVDI. The performance of TVWSI was relatively better across the sampling study windows, with the lowest R 2 of 0.31 and the highest MAPE of 145%, which is significantly less than other indices.

Evaluation of TVWSI in Predicting Soil Water Fraction (SWF) from Flux Tower Data, and Comparison with Existing Water Stress Metrics
The proposed water stress index was evaluated against the existing indices in predicting soil water fraction from four contrasting FLUXNET sites across Victoria and New South Wales. The first-, second-, and third-order polynomial functions were fitted using soil water fraction as the dependent variable and the indices as an independent variable. The scatter plots in Figure 12 show the observed vs. predicted soil water fraction from four FLUXNET sites using the first-order polynomial model, along with R 2 , RMSE, and model equation. R 2 and RMSE are the averages of the outcomes obtained after 10,000 iterations on the randomly selected training and testing datasets. CWSI was not calculated over FLUXNET stations due to the unavailability of vapor pressure deficit.

Evaluation of TVWSI in Predicting Soil Water Fraction (SWF) from Flux Tower Data, and Comparison with Existing Water Stress Metrics
The proposed water stress index was evaluated against the existing indices in predicting soil water fraction from four contrasting FLUXNET sites across Victoria and New South Wales. The first-, second-, and third-order polynomial functions were fitted using soil water fraction as the dependent variable and the indices as an independent variable. The scatter plots in Figure 12 show the observed vs. predicted soil water fraction from four FLUXNET sites using the first-order polynomial model, along with R 2 , RMSE, and model equation. R 2 and RMSE are the averages of the outcomes obtained after 10,000 iterations on the randomly selected training and testing datasets. CWSI was not calculated over FLUXNET stations due to the unavailability of vapor pressure deficit. TVWSI showed improved prediction skills for SWF than other indices at Wombat, Wallaby, and Tumbarumba FLUXNET sites and slightly better at Whroo (Table 3). Tumbarumba and Wallaby sites are the wettest sites where the forest type is classified as a temperate evergreen wet sclerophyll forest [70,72]. These energy-limited systems with continuous access to soil moisture show minimal deviations in the greenness. This is evident from the NDVI range over these forests, which ranges between 0.65 and 0.8. The low performance of the indices over Tumbarumba might have been due to the saturation of red reflectance and thus NDVI at high biomass [93,94] or its considerable lag with soil TVWSI showed improved prediction skills for SWF than other indices at Wombat, Wallaby, and Tumbarumba FLUXNET sites and slightly better at Whroo (Table 3). Tumbarumba and Wallaby sites are the wettest sites where the forest type is classified as a temperate evergreen wet sclerophyll forest [70,72]. These energy-limited systems with continuous access to soil moisture show minimal deviations in the greenness. This is evident from the NDVI range over these forests, which ranges between 0.65 and 0.8. The low performance of the indices over Tumbarumba might have been due to the saturation of red reflectance and thus NDVI at high biomass [93,94] or its considerable lag with soil moisture [90,95]. The improved performance of TVWSI in predicting SWF could be due to the presence of SWCI in its model, which can perceive changes in canopy water content even from higher vigor canopy. The poor performance of PDI and MPDI at Whroo could be due to the absence of LST in its model. In contrast, the relatively better performance of the indices that include LST in their model may be related to the high correlation of LST and surface soil moisture content in the driest regions with open-canopy forests [96][97][98]. The relationship of LST and soil moisture content generally gets decoupled and becomes more complex at the wetter sites such as Tumbarumba and Wallaby because dense closed forest canopies usually characterize the wetter sites. A considerable amount of variations in the canopy temperature come from the canopy structure and background soil viewed to the sensor [35]. Therefore, to efficiently estimate soil moisture content, there should be a clear partitioning between temperature sensed from the canopy and the background soil, which is still a significant limitation of the satellite remote sensing [99].

Discussion
NDVI and SWCI are the indirect indicators of soil moisture over the vegetative surface. Their functioning depends on the biophysical links between soil moisture dynamics and its feedback with the vegetation canopy. NDVI correlates with soil moisture because of high vegetation fractions in wetter regions and vice versa. In comparison, the SWCI's association with soil moisture is based on a linear relationship between soil and vegetation water content. Situations could exist in a forest ecosystem when the vegetation fraction at two different locations is the same, but soil moisture is different. These scenarios can be well estimated by combining the indicator of vegetation fraction and vegetation water content. Replacing LST by SWCI in the LST-NDVI trapezoid has fulfilled this inability of NDVI.
The index d(SWCI, NDVI) is more sensitive than NDVI and SWCI at the lower and higher end of the soil moisture spectrum. The boxplots shown in Figure 8 demonstrates that d(SWCI, NDVI) has improved retrieval of soil moisture dynamics at the lower and higher ranges of soil moisture. This capability of retrieving soil moisture at extreme ends is a significant improvement over NDVI, which generally saturates and tends to represent a general nature of lagged canopy response to root zone soil moisture [90][91][92]95,100]. The improved correlation of soil moisture vs. d(SWCI, NDVI) at the higher and lower moisture regions could also be due to SWCI in its model, which can perceive water dynamics from the canopy and the soil [101].
TVWSI was constructed using d(SWCI, NDVI) and LST, and it is based on the theoretical principle of MVWSI, which is a simple derivation of slope from LST-NDVI and is sensitive to plant available water [21,35,37,101]. However, LST's dependence on ambient atmospheric conditions is the major limitation of parameterizing LST-NDVI spectral space on temporal observations. That is why we have extracted the parameters using SWCI-NDVI spectral space and then combined the extracted intermediate index with the standardized LST. TVWSI has shown significant improvement in R 2 and linearity with AWRA-L soil moisture over d(SWCI, NDVI), as shown in Figures 7c and 9. The reason could be that the recent rainfall abruptly increases soil moisture in the root zone and decreases the canopy temperature, which LST can perceive efficiently compared to the optical reflection. TVWSI has also shown sensitivity to the wettest and driest events during the millennium drought, as shown in Figure S1. This capability of TVWSI makes it a reliable indicator of root zone soil moisture, and we contend that TVWSI can be employed to estimate drought intensity.
The improved performance of TVWSI in predicting AWRA-L soil moisture and soil water fraction shows its superiority over other the existing water stress indices compared in this paper. The improvement over PDI and MPDI demonstrates that optical reflectance on its own is not a good predictor of soil moisture dynamics in the root zone. LST is a subtle measure of plant water stress, and LST's sensitivity to transpiration cannot be realized by optical reflectance. The improvement over TVDI can be explained by the fact that LST depends on the ambient atmospheric condition. That is why the spectral space of LST-NDVI cannot be globally parameterized for the temporal observations. In comparison, the improvement of TVWSI over CWSI and MVWSI could be due to the inclusion of canopy water content, which is missing in the models of MVWSI and CWSI.
TVWSI's performance decreased when applied to the prediction of point-scale SWF values from ground monitoring sites when compared to AWRA-L soil moisture. This may be due to the large discrepancy in their support scales. It had a superior skill in root-zone soil moisture prediction compared with the other indices. This study has considered an area of (2 km × 2 km) for the study sampling windows, which could be spatially comparable with AWRA-L (approximately 5 km × 5 km), but their spatial representation is far greater than the "point" scale ground soil moisture measurements. Furthermore, studies have shown that the soil moisture features very high spatial variability, changing even at the scale of meters for the surface soil moisture [102]. Therefore, it could be assumed that the point measurement is not representative of a 2 km × 2 km footprint. The other potential cause could be that SWF was measured at variable soil depth across the year, and also, the method to calculate SWF across different sites are different [67]. The data availability was also not continuous across the sites; for example, the availability of data was Wombat

Conclusions
This study investigated the hypothesis that an improved water stress index could be constructed by representing canopy water content information to the LST-NDVI trapezoid model. Moreover, the LST-NDVI feature space cannot be used to parameterize temporal observations due to LST's dependency on ambient atmospheric conditions. The following conclusions can be made from the study:

1.
The intermediate vegetation index d(SWCI, NDVI) is modeled using NDVI and SWCI to simultaneously consider the effects of vegetation growth and vegetation water content. d(SWCI, NDVI) showed an improved correlation with AWRA-L soil moisture than NDVI and SWCI. Moreover, d(SWCI, NDVI) improved retrieving soil moisture dynamics at the lower and higher ranges of soil moisture, which is a significant improvement over the saturation nature of NDVI over high-density vegetation and its inability on open canopy regions; 2.
Compared to PDI, MPDI, MVWSI, TVDI, and CWSI, the newly devised TVWSI improved prediction ability for the modeled AWRA-L root zone soil moisture from sixty study sampling windows and soil water fraction (SWF) from the four contrasting FLUXNET sites across Victoria and New South Wales; 3.
As the optical spectral space of SWCI-NDVI and standardized LST is used to model TVWSI, it can be employed to estimate soil moisture dynamics on the temporal data.