Improvement of the “Triangle Method” for Soil Moisture Retrieval Using ECOSTRESS and Sentinel-2: Results over a Heterogeneous Agricultural Field in Northern India

Abstract

For the purpose of deriving spatiotemporal estimates of soil moisture, the triangle method is one of the most widely used approaches today utilizing remote sensing data. Generally, those techniques are based on the physical relationships that exist when a satellite-derived land surface temperature (Ts) is plotted against a spectral vegetation index (VI). The present study proposes an improvement in the triangle method in retrieving soil moisture over heterogeneous areas. In particular, it proposes a new approach in robustly identifying the extreme points required for the technique’s implementation. Those extreme points are then used in calculating fractional vegetation cover (Fr) and scaled Ts. Furthermore, the study proposes a new approach for calculating the coefficients required to develop the relationships between surface soil moisture (SSM) and Fr/Ts, which is implemented using a model and field data. As a case study, an agricultural field in the Varanasi district in India has been used, on which the triangle method is implemented using ECOSTRESS and Sentinel-2 data. The much-improved spatial resolution satellite data of ~70 m from ECOSTRESS allowed deriving more vivid results of SSM spatial variability for the study area. Comparisons between field soil moisture calculated using the proposed method returned an RMSE of 0.03 and R² value of 0.84, which are considered very satisfactory. The methodology proposed herein and the results obtained are of significant value with regards to the triangle method, contributing to ongoing efforts at present examining its use for operational product development at a global scale.

1. Introduction

Water that is held between the pores of the soil particles is known as soil moisture. Several biological and physical phenomena that occur at the surface are affected by the level of water stored within the first layers of the soil, called Surface Soil Moisture (SSM) [1]. Soil moisture (SSM) has a significant role in the Earth’s water cycle, contributing to the change of weather as well as to the climate from regional to global scales [2,3]. SSM also exerts a major control over the exchange of heat and water between the interface of land surface and the atmosphere [4]. It is a variable playing a substantial role in seasonal climate and in simulating future projections of the global climate variables [5], in drought monitoring [6,7], weather forecasting [8], vegetation health index [9], risk assessment of forest fire [10], and irrigation management [11,12,13]. Soil moisture regulates runoff potential and the availability of water for plants and humans. As such, information on its temporal and spatial distribution at a variety of geographical scales is of prime importance [14].

Field-based methods consist of the traditional approach for retrieving information on soil moisture and there is a wide range of instruments available to facilitate this purpose [15]. However, soil moisture retrievals using those methods are generally point-based and cannot be used to interpret soil moisture for large areas [16]. The use of remote sensing data presents the most promising avenue to address this challenge. The ability of remote sensing techniques to provide repetitive, synoptic views in a spatially contiguous fashion without disturbing influence on the area to be surveyed and without site accessibility are the key advantages of this approach [17,18]. A plethora of methods have been proposed for deriving spatially distributed estimates of surface soil moisture utilizing spectral information acquired in all regions of the electromagnetic radiation spectrum [12,19,20]. Those methods, which vary from empirical to purely physically based, often combine remote sensing data with ancillary surface and atmospheric observations and are characterized by diverse strengths and limitations with regard to their practical implementation [21,22,23]. A recent discussion of the available methods and relevant satellite products can be found in [24,25].

A method that stands out from others in this group is the so-called “triangle method”. This approach uses a combination of optical and thermal infrared data to generate the variability of soil moisture over a given area [6]. Variants of this group of methods have been developed by different investigators, which include residual methods based on the energy balance equation and methods based on the information derived from the ‘scatterplot’ relationships between satellite-derived vegetation index (VI) and surface radiant temperature (Ts) measures [26]. Ts/VI methods have certain advantages compared to other EO-based modeling approaches employed in the retrieval of soil surface moisture. Some of those include their improved ability to deal with surface heterogeneity and to avoid dependence on external surface and meteorological parameters [14]. Additionally, key input data are relatively easy to obtain from space or airborne data over large areas, which make them ideal for operational implementation scenarios. However, most of these approaches have been applied with satellite data with a spatial resolution of ~100 m [27] or higher i.e., ~250 m [28], ~500 m, or ~1000 m [29]. Additionally, most approaches are characterized by a certain degree of subjectivity in regard to the estimation of the Ts/VI extreme points required in their implementation [30,31]. Thus, there is a clear requirement to investigate the use of those techniques using satellite data having a high spatial resolution (i.e., below 100 m) and also to automate the process of anchor/extreme points selection. The automation of the process including the coefficients retrieval will allow the user to easily calculate soil moisture of any region for any month of the year. The automation of coefficients retrieval will allow users to take a wide range of values between 0 and 1, i.e., to the exact decimal value the user intends. This will also enable the user to create time series of any large area without much effort.

In purview of the above, this study’s aim is twofold: (1) to evaluate the ability of the triangle method in deriving soil surface moisture by using ~70 m spatial resolution ECOSTRESS (ECOsystem Spaceborne Thermal Radiometer Experiment on Space Station) and 10m spatial resolution Sentinel-2 satellite data, and (2) to develop an automated approach for retrieving the extreme points required in the triangle method implementation and for calculating the coefficients of fractional vegetation (Fr) cover and scaled surface temperature (Ts), assisting in eliminating any subjectivity from this method in this regard. The soil moisture values obtained from the algorithm for soil moisture are then validated using Steven’s hydra probe’s soil moisture values.

2. Triangle Method

The triangle method uses optical and thermal data to calculate soil moisture [30] (Figure 1). This method does not require any initial surface condition or atmospheric correction to work with. Instead, it is based on the NDVI product from an optical sensor and LST from a thermal sensor to generate spatiotemporal estimates of soil moisture [32]. Carlson [33] demonstrated that the sensitivity of Ts for soil moisture is more for bare soil than for the vegetated soil. This also varies with the type of leaf and soil surface around the leaf. This is the reason for obtaining the triangle shape for Ts/Fr plot. Ridd [34] and Carlson [35] further illustrated different regions in the plot corresponding to the biophysical properties. Moran [36] and Gillies [37] defined the boundaries of the triangular shape and their physical constraints using spatial datasets and their use to determine surface energy flux and soil moisture. Initially, generation of certain soil moisture retrieval parameters are needed through ground measurements, but later we can solely rely on the satellite products, which are frequently available for a large region. This has been clearly demonstrated in a few studies, e.g., by Owen [38] and Chauhan [39].

The triangle shape is formed by plotting the scaled surface temperature and fractional vegetation cover in the form of a scatterplot. This plot resembles the shape of a triangle or trapezoid when there is a full variability in fractional vegetation cover over the area. Different regions of the plot depict different soil moisture conditions. (1) represents well-watered vegetation having high F_r and low T_s. (2) represents water-stressed vegetation having high F_r and relatively moderate T_s. (3) represents saturated bare soil; the region having low Fr and low Tr (the lower most portion on the left side of the plot). (4) represents dry bare soil; the area having low Fr but high Tr (the lower most portion on the right side of the plot). (5) represents wet edge; the straight line at a low Tr value goes all along the Fr. This line depicts the presence of moisture. The soil moisture values along this line are relatively high. (6) represents dry edge: The slope of graph towards the left (lower temperatures) with increasing vegetation fraction is a consequence of the fact that sunlit vegetation is generally cooler than sunlit bare soil. The vegetation in this region has less water content and the increase in temperature is due to a decrease in density of vegetation as we go down the dry edge. (7) represents full cover: The area having maximum values of F_r (the upper most portion of the plot) and values of Tr in the mid region show the region with high vegetation density. Crops in this region are considered to be in good health. (8) represents partial cover. This area has intermediate Fr values (the central portion of the plot) and moderate values of Ts. The vegetation lying in this region has moderate density. (9) represents bare soil: This area in the plot has less Fr (lower portion of the plot) and moderate Ts.

3. Materials and Methods

3.1. Satellite Data

To satisfy the study objectives, remote sensing data from Sentinel-2 and ECOSTRESS satellites were used herein (

Table 1. Properties of data sets.

Properties	Sentinel-2	Ecostress
Altitude	786 km	400 km
Swath	290 km (FOV:20.6°)	384 km
LTDN	10:30	05:00
Revisit	5 days	3 days
No. of Bands	10	6
Spatial Resolution	10, 20, 30 m	~60, ~70 m
Data extracted	NDVI	LST

). Sentinel-2 has an average height of 785 km and the presence of two satellites allow repeated surveys every 2 to 3 days. The satellite launched on 23 June 2015, has a lifespan of 7.25 years, and has been placed in sun-synchronous orbit. The data obtained are 100 × 100 km² ortho images in UTM/WGS 84 projection. The coverage limits are between 56° and 84° north. The orbit is 10 days at the equator with one satellite, and 5 days with 2 satellites under cloud-free conditions, which results in 2–3 days at mid-latitude. Sentinel-2 provides products with a high spatial resolution at every 5-day interval, which makes it most suitable for this purpose. The spatial resolution of Sentinel-2 is ~10 m. The reason for using the ECOSTRESS having a swath width of 402 km for the LST product is its 4-day revisit time along with the high spatial resolution of ~70 m. The selected image must have some bare soil pixels and some fully vegetated pixels for application of the triangle method [16]. These bare and vegetated pixels are used in the selection of extreme value pixels. The values of NDVI and LST for these extreme value pixels are used for the calculation of Fractional Vegetation Cover (Fr) and Scaled Surface Temperature (Ts). To make bare soil and vegetated land coexist, we use a large area into account.

3.2. Experimental Site and Ground Measurements

The study site consisted of an agricultural area, which lies under the Indo-Gangetic region in the south-eastern part of Uttar Pradesh, India. The extent of the region is from 25°10′3” and 82°40′50” N to 82°12′1” E. The river Ganga, which flows along the eastern border of the area from north to south, is the major source of irrigation to the agriculture prevailing in area. Soil moisture reference measurements were collected during a field campaign at the experimental site concurrently to the satellite overpass in the different patches of the agricultural area at the campaign site. Due to its rich importance in agriculture, the soil moisture products of the region are necessary to be calculated. This method uses reflectance from the land to calculate soil moisture, which comes as a downside when used on terrain with steep slopes. The steep-sloped terrain has some regions of shadow when capturing the image. This reduces the efficiency of the results. Thus, the area selected for applying the triangle method should have mild to moderate slope. These samples were taken from the field site using the Steven’s hydra probe, which is based on the principle of “dielectric impedance”. Its detailed mathematical and signal characterization of the dielectric spectrum helps factor out errors in the soil moisture measurement such as errors due to salinity, temperature effects, and soil type. Thus, it gives out results better than Babaeian [40] or Mananze [41]. The value of soil moisture varied between 0.11 (m³/m³) to 0.415 (m³/m³).

3.3. Data Pre-Processing

The process of SSM retrievals using the modified scheme proposed in this study is summarized in Figure 2. NDVI and LST were calculated as proposed by Chauhan [39]. NDVI estimation was based on band 4 and band 8 of Sentinel-2. For calculation of LST, a scale factor of 0.02 to the ECOSTRESS image was applied and temperature was also converted to Celsius units. NDVI of Sentinel was then resampled to match the spatial resolution LST i.e., ~70. As the minimum and maximum pixel values change with every image, it was deemed necessary to automate the range calculation for every image. This would have been easy taking the extreme pixel values, but it may not cover a considerable area, which can induce errors in calculations. Hence, the approach was to take into consideration the first 10% of pixels for minimum range and the last 10% of pixels as maximum range (after arranging pixel values in ascending order).

It is a known fact that barren soil will have high LST and low NDVI values which will give us LSTmax and NDVImin, while the most vegetated soil will have low LST and high NDVI values, which will give us LST_min and NDVI_max. For identification of these sites, a large enough area is needed which has homogeneous data. So, we chose a grid of 25 pixels which equates to 60 m × 60 m (which is 3600 m²). This grid of 25 pixels, which hovers across both NDVI and LST values, searches for values that are within the range for that image. If all 25 values come under the specified range, then the value is stored in the list for the barren or vegetated pixels list. This is done to ensure that barren or vegetated land values were chosen from a large enough area for reducing errors. After finding the list of values for barren and vegetated pixels, we needed to find the one with maximum difference between maximum and minimum values. Then, the percentage change of both NDVI and LST in barren as well as vegetated pixels was taken. The values with maximum difference were chosen as NDVImin and LSTmax from barren pixels and NDVImax and LSTmin from vegetated pixels. This step was implemented to ensure that maximum pixels were covered under fractional vegetation cover and scaled surface temperature.

3.4. Retrieval of Coefficients for Fractional Vegetation Cover, Scaled Surface Temperature, and SSM

F_r for every pixel is calculated using the NDVI_min and NDVI_max values and T_r for every pixel is calculated using the LST_min and LSTmax values. These values range between 0 and 1. For calculation of coefficients for A_i and A_j, 100 random values are generated between 0 and 1. The no. of random values can be changed depending upon the precision we need for the coefficient. These values are used along with field data to generate values for the right-hand side of the equation to match the left-hand side of the equation for every field point. This process eliminates the non-equal values of A_i and A_j, hence, giving out the most appropriate coefficient.

SM_field = 1 − ((A_i × T_s)/(1 − (A_j × F_r)))

(1)

Here, SM_field is the field soil moisture. A_i and A_j represent the coefficient for scaled surface temperature and fractional vegetation cover, respectively. F_r and T_s represent the fractional vegetation cover and scaled surface temperature, respectively.

Soil moisture is calculated using the formula given:

SM_satellite = 1 − ((A_i × T_s)/(1 − (A_j × F_r)))

(2)

Here, here SM_satellite is the soil moisture obtained using the satellite images. A_i and A_j are the coefficient of scaled surface temperature (T_s) and coefficient of fractional vegetation cover (F_r), respectively. The F_r and T_s represent fractional vegetation cover and scaled surface temperature, respectively.

3.5. Validation Approach

Accuracy assessment of the triangle-derived SSM maps was based on point-by-point comparisons, which were performed between the predicted and measured parameters. For this purpose, a series of statistical scores were computed to quantify the agreement between predictions and observations, summarized in

Table 2. Quantitative measures used to assess the agreement between the “triangle”-derived estimates and in situ observations. Subscripts i = 1 … N denote the individual observations, P and O denote the predicted and observed values obtained from field measurements. The horizontal bar in Scatter/MSD ratio equation denotes the mean value.

Name	Description	Mathematical Definition
Bias/MBE	Bias (accuracy) or Mean Bias Error	$bias=MBE=\frac{1}{N}{\displaystyle \sum _{i=1}^N(P_i-O_i)}$
Scatter/SD	Scatter (precision) or Standard Deviation	$scatter=\frac{1}{(N-1)}{\displaystyle \sum _{i=1}^N\sqrt{{(P_i-O_i-\overline{(P_i-O_i)})}^2}}$
RMSE	Root Mean Square Error	$RMSD=\sqrt{bia{s}^2+scatte{r}^2}$

. The selected statistical metrics have already been used in analogous verification exercises.

4. Results

4.1. Evaluation of Variables for Triangle Method Implementation

The NDVI of the region varied from 0 to 0.70 (Figure 3a,c). The presence of blue color in the majority of the map shows presence of vegetation in the region. This is due to the prevalence of agriculture over the entire region. The LST values are lower in the regions with high NDVI values. The value of NDVI in the south-eastern corner is shows much less red color. This is due to the absence of vegetation in this region. The LST of the region varied from 20 °C to 35 °C in 2020 (Figure 3b), while in 2021 it was between 17 °C to 40 °C (Figure 3d). The value of NDVI in the south-eastern corner is very high, showing a blue color. This is due to the absence of vegetation in this region.

4.2. Extreme Points

The range for NDVI_min was between 0 and 0.27708, while NDVI_max was between 0.55416 and 0.6927. The range for LST_min was between 20.6865 °C and 23.49402 °C, while for LST_max it was 27.7053 °C and 34.7241 °C. Using these range values, extreme points were calculated. The extreme points (EP) model gives out 15536 combinations of barren pixels when NDVI was taken between 0 and 0.27708 and LST was taken between 27.7053 °C and 34.7241 °C. In these combinations, the value with maximum difference was taken as barren pixels having an NDVI_min value of 0.277031004 and LST_max value of 33.64904404. The number of combinations for vegetated pixels is 11241. These values were generated using NDVI values between 0.55416 and 0.6927 and LST values between 20.6865 °C and 26.30154 °C. In these combinations, the value with minimum difference was taken as vegetated pixels having an NDVI_max value of 0.554162562 and LST_min value of 21.3511467.

4.3. Optimized Parameters and Coefficients

The model randomly generates 100 normally distributed values between 0 and 1. The combination of these values is made for assigning these 100 values to each of the coefficients. This makes the total set of values equal to 100 × 100, i.e., 10,000. These values are then swapped on the soil moisture results to find the best suited value. The coefficient for fractional vegetation cover (A_i) is 0.99 and the coefficient for scaled surface temperature is 0.74. Using these coefficients, satellite soil moisture was calculated for all field data points. The scatter plot between field soil moisture and satellite soil is given in Figure 4. The root mean square error for the soil moisture results was 0.03 and the R² value was 0.835.

4.4. Triangle Construction and Soil Moisture Retrieval

The scatter plot between the pixels of Tr and Fr resembles a triangle shape. The horizontal and vertical axis depicts the variation of temperature and vegetation, respectively. The shape of the plot in this study is deformed at the top due to the presence of water-stressed vegetation (point 2). This scatter plot shows the density of pixels from the study area under different conditions.

In regard to the warm edge (pixel 6), this sharp edge is clearly visible in both the plots. The wet edge is visible in both the images at the left side of the triangle (point 5). The shape is sometimes deformed at the top in the form of a plateau region, as in this case (point 2). This shows the presence of water-stressed vegetation due to irrigation in this region. At the top-most tip of the triangle, there is well-watered soil (point 1). Saturated bare soil (point 3) is abundant in the region in contrast to dry bare soil (point 4). There are some pixels in the lower middle portion showing the presence of barren land (point 9). In the center of the triangle, partially covered soil is found (point 8). The soil moisture values of the agricultural land varied between dry soil, i.e., 0.0 (m³/m³), to saturated soil, i.e., 0.6 (m³/m³). The map for the same is given in Figure 5. The soil moisture of the region is on the higher side due to irrigation in the month of February. The presence of higher values can also be seen in the triangle in Figure 5 along with soil moisture maps. The presence of bare dry soil is much lower, which was seen in the triangle and is clearly visible at point 4.

Field measurements show that there was high variability in soil moisture during the month of February 2020. It is the period of irrigation in this region. The results can be seen with some places having high soil moisture. This can be verified by the scatter plot between F_r and T_s. The optimization for the study area was conducted using these points for the month of February. If we can obtain more field data points for different months of the year, we can retrieve coefficients for year-long data. These can be used for generating time series data for the region using this method. The soil moisture estimated through this study gives results of R² = 0.835 and RMSE = 0.032, which are better than some soil moisture results using microwave data [32,33].

5. Discussion

This study improves the potential of visible/infrared sensors in determination of soil moisture suggested by Fuzzo [31] and Carlson [32]. Many studies are carried out applying these methods but most of them were limited to determination of TVDI (temperature vegetation dryness index). Some studies explored the relation of evapotranspiration and soil moisture using vegetation and temperature indices [42,43,44,45]. The spatial resolution of these studies [27,28,29] resorting to similar models was ~100m (in best case), which has been improved in this study to ~70m. The automation of extreme values, which were previously done manually [16], will reduce time spent and produce unbiased results. Moreover, the study provides a profound dry edge and wet edge, which can be implemented for calculation of evaporative fraction. Results obtained in the study yielded satisfactory results. The drawback of the study is that it cannot be applied in the regions with steep slopes and we need cloud-free images every time it is required to estimate soil moisture. However, these can be overcome by applying this method on areas with flat terrain and using these satellites’ frequent visits as an advantage. The annual optimization of this model using field measurements to retrieve coefficients for different times of the year will enable us to create time series data for soil moisture estimation for that particular region.

The approach used in this study can produce non-linear results of soil moisture from large datasets. The process of generation of series data is quite easy as we have automated the whole process and we need to establish the results only once to produce results. This can be done in later work, where a database can be created for soil moisture data for the whole world by initializing the initial results and creating high-resolution global soil moisture data. The major advantage for this lies in its use of pixel distribution for soil moisture generation, as this process results in not needing any ambient atmospheric or surface parameters for modeling. The universal triangle and coefficient for every area are supposed to be calculated once and we can accommodate successive images, thus, creating time series data. Different regions in a triangle show different soil moisture values, these can be used as an alternative for land use classification, as the barren regions will remain dry throughout the year, while the forest area will have soil moisture values throughout the year and crop area will have dry and wet months. The major drawback of this method is that it cannot be used in regions with steep slopes. It will fail in regions with dense forest or in semi-arid regions, as there will not be enough values of barren or vegetated land to calculate extreme points. There are many applications of the soil moisture data, and each of them requires data with different levels of precision and spatial extent. To provide data with a large spatial extent with more accurate results, this method can be a viable option. With the use of this study, we can automate the process which is used to calculate the coefficient of fractional vegetation cover and scaled surface temperature. This can further be used to extract soil moisture of the region in any time frame without any hassle. Furthermore, we can use this improved method on any satellite data sets having thermal and multispectral bands.

6. Conclusions

A modification of the triangle approach has been provided by combining optical and thermal infrared data to generate the variability of soil moisture over a given area. The ability of the triangle method has been demonstrated in deriving the soil surface moisture by using the new ~70 m spatial resolution ECOSTRESS and 10m spatial resolution Sentinel-2 satellite data and an automated approach for retrieving the extreme points required has been developed for the triangle method’s implementation. The soil moisture values obtained from the algorithm for soil moisture are then validated using Steven’s hydra probe’s soil moisture values extensively collected in the region. The soil moisture estimated through this study, when compared with the ground observations, yielded a square of correlation (R² = 0.835) and low error (RMSE = 0.032), which are quite satisfactory for different applications. The overall analysis showed that the method is estimating quite well the soil moisture values and could be used for large-scale implementation. The method has a potential to be used operationally and efforts are ongoing towards this direction and outcomes remain to be seen.