Utilizing Optical Satellite Imagery to Monitor Temporal and Spatial Changes of Crop Water Stress: A Case Study in Alfalfa

: Since the 1980s, thermal imagery has been used to assess crop water stress. The increase in the temporal resolution of optical satellite sensors (in the range of 400–2500 nm) and the better spatial resolution compared to the thermal imagery call for the deﬁnition of a new way for crop water stress monitoring. Hence, we are suggesting a new method utilizing spectral indices from three subsequent images to address this challenge. This method predicts the current water stress with the two past images and compares it to the current stress: if the existing conditions are better than the predicted stress, the crop is not under stress and has sufﬁcient water for development. To evaluate the suggested method, we downloaded Sentinel-2 images and compared the stress found with that method to the leaf area index, leaf water potential, and yield from seven alfalfa growth cycles. The results outline the ability of the new optical stress index to depict spatial and temporal changes in the alfalfa water stress and especially illustrated the changes in the crop water stress over the growth cycle and after each irrigation. This new method needs to be validated with different crops and satellite sensors to verify its success.


Introduction
Soil water stress (SWS) is an important factor in irrigation strategy. It represents the availability, or limitation, of the water in the soil for the crop. The Food and Agriculture Organization paper No. 56 [1] refers to the SWS as Ks (measured by soil sensors), which ranges from 1.0 for a well-watered crop to 0.0 when the lack of available soil water limits transpiration. Furthermore, in the irrigation strategy, Ks is used for irrigation scheduling and thus to increase the water use efficiency (WUE). In addition, the water-stressed strategy aims to achieve better yield quality [2][3][4], decrease plant diseases [5], or even reduce soil moisture to facilitate the harvest with heavy machinery [6].
Alfalfa (Medicago sativa L.) is a perennial crop with a non-stress strategy. Several studies have suggested that the biomass and yield increase with irrigation [6,7] or decrease when water stress is increasing [8]. Further, the timing of the irrigation not only affects the yield quantity and quality, but it is also important for weed control [9]. Yet, the need to decrease soil moisture before the harvest may develop stress at the end of the growth cycle or the beginning of the next one. This may reduce the following growth cycle's biomass and yield.
In parallel to the SWS, the crop water stress (CWS) is highly affected by the local climate. For example, the vapor pressure deficit (VPD) above 3.0 kPa increases the crop water stress dramatically [10,11]. Unfortunately, the spatial variability of climate parameters is greater than the field scale and, hence, does not allow revealing differences within the field.
For in-field mapping of the CWS, remote sensing data are used via a combination of optical (400-2500 nm) and thermal (10,000-12,000 nm) imagery [12,13]. However, thermal application for the CWS has two main problems. (1) The thermal pixel size is usually courser for agriculture plot applications. For example, a Landsat-8 thermal pixel is 100 × 100 m, requiring the plot width to be greater than 200 m to obtain one pure pixel. It contradicts the precision agriculture practice of reducing the plot size for better management. (2) The availability of thermal data is mainly from specific satellite sensors such as Landsat, Moderate Resolution Imaging Spectroradiometer (MODIS), and Sentinel-3 while optical sensors are more common.
To address the need for in-field CWS mapping without thermal data, several optical indices were developed [14][15][16], which we refer to as the remotely sensed water stress index (RWSI). Yet, their values are specific for each image and thus do not permit depicting the water stress and its temporal changes, especially after irrigation. This effect is illustrated by the integration of thermal and optical imagery as presented by Moran [17]. She found that immediately after an irrigation event, there is a decrease in the CWS, and only after a day or two there is an increase in growth. Hence, we hypothesize that between irrigation events (meaning 2 days after irrigation and until the next event), the crop growth and the CWS linearly correlate and allow us to predict the CWS solely via the growth value. Further, after the irrigation, the real value of the CWS is better than this prediction, and the difference between the prediction and the current values is positive for non-stress situations and vice versa.
To test this hypothesis, we mapped the CWS via regression of the remotely sensed growth index (RGI, as the normalized difference vegetation index, NDVI) to the RSWI (as the normalized difference water index, NDWI). Repeated measurements of these indices for the same location generate the new optical stress index (OSI) by subtracting the current RSWI value from the predicted RSWI (which was calculated by the regression of the previous RSWI and RGI and the current RGI). This index ranges from 1 for a well-watered and non-stressed crop to −1 for a high-stress crop. This work aims to evaluate the performance of the new suggested index, OSI, as an indicator for the CWS. To do so, temporal and spatial changes in spectral data, in situ physiological measurements, meteorological information, and yield were collected, analyzed, and compared.

Study Area
The alfalfa field is located in the northeastern part of Israel (coordinates: 35.558 E, 32.486 N). The crop was seeded in November 2019 and the plot was divided into the Northern (14.5 ha) and Southern (12.9 ha) management zones ( Figure 1). Center pivot irrigation at night started at the southern corner, three times in every growing cycle (25-28 days). Due to technical constraints, the irrigation amounts were similar for each cycle, regardless of reference evapotranspiration, an average of around 8 mm day −1 . At the beginning of September, the center pivot got stuck in the middle of the northern part and continued to work only five days later. Therefore, only two irrigations were applied in that cycle, and the yield was recorded only in the Southern management zone (Table 1). Weather data were collected from the Khavat-Eden station (https://www.meteo.co.il/home. accessed on 23 May 2022) 7 km south-west of the field. It included the daily values of the minimum (Tn) and maximum (Tx) air temperature ( • C), the minimum (VPDn) and maximum (VPDx) VPD (kPa), and the daily sum of the Penman-Monteith evapotranspiration (ETo, mm day −1 ).

Crop Measurements
During the 2021 growing season (21 March to 12 September), we measured the leaf area index (LAI), the leaf water potential (LWP), and the topsoil water content (SWC) at nine sites ( Figure 1). The specific location of each measurement site was determined to represent the inner (sites 1 and 5) and outer sides of the center pivot (sites 2-4 and 6-10). We measured the LAI, the LWP, and the SWC in a radius of 10 m around each site and recorded the time after irrigation (TAI) in hours for each site on each measurement day.
The LAI was averaged from 25 readings of a plant canopy analyzer (SunScan model SS1, Delta-T Devices, Cambridge, UK). The LWP was averaged from 6-8 young fully developed leaves with a pressure chamber (model 655, PMS Instrument Company, Albany, OR, USA). The SWC was averaged from 12-15 readings of topsoil (10 cm) with a handheld sensor (SWM5000, STEP Systems GmbH, Nuremberg, Germany). The instrument was calibrated at each site, and its values range from 100% for full-field capacity to 0% for fully dry soil.

Name Equation Source
Remotely sensed growth index Remotely sensed water stress index 4 Modified tasseled cap wetness To calculate the optical stress index (OSI), we calculated the above indices for each site, for each image ( Figure 2A). Next, for each image, we utilized the previous two images' RWSIs and RGIs for regression calculation. This regression was calculated between any of the RWSIs and any of the RGIs to predict the current RWSI value ( Figure 2B): Then, this value was subtracted from the actual RWSI ( Figure 2C): In order to keep the OSI values between −1 and 1, if the values were greater than 1.0 or below −1.0, they were changed to 1.0 or −1.0, respectively. The OSI serves as the stress indicator: positive or negative values of the OSI indicate that the crop is well-watered or is under stress, respectively. To define the best way to generate the OSI, we calculated it on the basis of 12 combinations (three RGIs multiplied by four RWSIs) and compared the results with the crop measurements.

Statistical Analyses
To identify the best OSI combination and report if it better maps the crop measurements (LWP, LAI, and yield) compared to the above spectral indices, we calculated the Pearson correlation [23].
where x and y are two arrays, n is the number of observations, Xi is the i th value of observation in the x array and Yi is the i th value of observation in the y array. We selected this statistical indicator to illustrate not only the magnitude of the relationship (weak or strong), but also its direction (negative or positive).
We also deployed a one-way analysis of variance (ANOVA) to investigate the effect of the time after irrigation (TAI) on the LAI and LWP. HSD Tukey's test (α = 0.05) was used to compare the TAI timeframe where plants experienced the highest and lowest water stress and the statistical assumptions for normal distribution and equal variances were tested. All the data analyses were performed using Microsoft Excel (Microsoft Corporation, Redmond, WA, USA), and the Python programming language.

Results
The results were divided into three sections: the first section defines the best OSI combinations by comparing the spatial variability of the OSI combinations and the spectral indices to the spatial variability of the crop measurements. The following Section 3.2 outlines the ability of the best OSI and two reference spectral indices (NDVI and NDWI) to track temporal changes in the alfalfa CWS at different growth cycles. Lastly, Section 3.3 focuses on the relationships between the yield, weather, spectral data, and crop measurements to explore how spectral indices and the OSI can predict the yield.

Spatial Variability of the Spectral Data: Alfalfa LWP and LAI Measurements
We examined the spatial variability at seven specific dates, where we had simultaneously crop and imagery data (only on 12 September 2021 the crop measurements included only the LWP, while the rest included the LAI and the LWP). A summary of the correlations between the LWP and the spectral data ( Figure 3) reveals that only the NDWI and the OSI combinations of ndvi2ndwi and ndvi2wetness had positive values. Although the wetness index had only one correlation value of −0.02, the rest of the spectral data had values greater than −0.15. Thus, only the maps that were based on the above three indices showed correctly the location of high stress and full crop moisture.
In contrast, the correlation values between the LAI and the spectral data were mostly positive and higher (Figure 4).
Not only the spectral indices were highly correlated with the LAI, but they were also highly correlated among themselves (Table 3). Further, the relationships between the selected OSI combinations ndvi2wet and ndvi2ndwi were also high, but the relationship with the other indices was not strong (r < 0.4). Thus, the temporal variability is demonstrated with the NDVI, the NDWI, and OSI ndvi2ndwi.

Temporal Variability of the OSI and the Crop Parameters
The ability of the OSI to track temporal changes is demonstrated with three growth cycles and four sites (two in each management zone). In these periods, we had a lot of images and crop measurements for presentation and comparison. The selected sites represent the inner sites in the center pivot (sites 1 and 5 in the Southern and Northern zones, respectively) and the outer locations (sites 2 and 7, in the Southern and Northern zones, respectively).
Sites 1 and 2, the second growth cycle (14 April to 9 May 2021; Figure 5). The general structure of the OSI values formed a semicircular curve, with stress conditions at the beginning and end of the cycle. It is associated with the shape of the LWP measurements. Yet, these measurements also showed the direct effect of an irrigation event, with the peak LWP one or two days after that event (the effect of irrigation on the crop water stress is present at the end of this section). Although the OSI, in general, does not show this irrigation influence, possibly because its dataset is not daily, it is not associated with the SWC ( Figure 5F). As such, the OSI relationships to the soil water stress are less pronounced and are stronger than the CWS.
The LAI measurements, as well as the NDVI and NDWI values, present a linear growth regardless of the irrigation timing. Both spectral indices, and especially the NDWI, mapped the LAI difference between sites 1 and 2. Yet, this continuance difference between these sites was not observed by the OSI or the LWP. Sites 5 and 7, the second growth cycle (17 April to 11 May, 2021; Figure 6). The temporal variability at these sites reveals the same trend, with one exception. As the third irrigation was only 3 days before harvest, the OSI and the LWP values after the 11 May harvest were higher than after the 16 April harvest.
Sites 1 and 2, the fourth growth cycle (4 to 29 June 2021; Figure 7). The main characteristic is the arc curve of the OSI and the LWP data, in contrast to the LAI, NDVI, and NDWI shapes. The increased numbers of soil moisture and LWP measurements outline the difference in the irrigation, where the peak crop LWP does not immediately follow the irrigation, but rather after a delay of more than a day, especially after the first irrigation.
Sites 5 and 7, the fourth growth cycle (6 June to 1 July 2021; Figure 8)-the OSI reveals site 5 had a lower stress than site 7 only at the beginning of the growth cycle, an observation that confirms the LWP measurements. After the first irrigation, there were not differences between these sites, as outlined by the spectral and the crop data.   Sites 1, 2, 5, and 7, the seventh growth cycle (16 August to 14 September 2021; Figure 9). The last cycle was mainly characterized by irrigation malfunctioning and the grower's decision to harvest without the third irrigation. As a result, the OSI and the LWP present a high CWS at sites 1 and 2 even before the harvest ( Figure 9A,B) while no stress was present at sites 5 and 7 ( Figure 9D,E). It is also reflected in the unusual NDVI curves, especially at sites 5 and 7 ( Figure 9F), with a decrease just after the second irrigation, and then it increases again.
To test if the OSI can detect water stress changes, we grouped all the LWP measurements and the spectral indices as presented above (Figures 5-9) and calculated the correlation coefficients. The results between the LWP and the OSI values are higher (r = 0.52) than the correlation of the LWP and the NDVI and NDWI (0.12 and 0.04, respectively). These results (see also Appendix C) highlight the greater efficiency of the OSI to detect changes in the water stress over time, while the NDVI and the NDWI could not be used for this proposal. Furthermore, maps of the three spectral indices from the 12 September 2021 image, just two days before the harvest, reveal how the OSI illustrates the CWS compared to the NDVI and the NDWI (Figure 10). In each map, the values for sites 5 and 7 were better than for sites 1 and 2, yet only the OSI map showed the high stress at site 2, which corresponds to an LWP of −28.4 kPa. The −28 value for alfalfa LWP was indicated by [24] as representing high stress in this crop. As such, these three maps illustrated the problem with the NDVI or the NDWI, as their values need interpretation to make a decision. For example, what is the meaning of the NDVI equal to 0.5? In contrast, an OSI value below 0.0 represents stress and makes it easier for understanding. The disadvantage of the OSI is presented in the western side of the map, where clouds in the previous image require removing these pixels from the analysis.  To examine the abovementioned irrigation effect on the peak LWP lag in time, we grouped these measurements per hour after irrigation: 0-12, 12-24, 24-36, 36-48, 48-72, and 72-96 ( Figure 11). It confirms, together with the ANOVA results (Appendix B), that the lowest LWP values (i.e., the highest level of crop moisture) does not occur instantly after the irrigation, but rather 12-36 h later.

Temporal Variability of the Alfalfa Yield
The growth cycle yield amounts (for each management zone, Table 1) were negatively correlated to the weather parameters (Table 4), highlighting the limitation of the irrigation amounts to overcome the extreme weather. The positive and highly correlated values between the yield and the LWP, especially the maximum LWP values representing the peak crop moisture, also emphasize the effect of the crop water stress on the yield. To represent each growth cycle with the spectral data, we could not use the average or median value because sometimes there were two images between the harvest and the first irrigation ( Figure 5) while other times there was only one ( Figure 6). Therefore, we utilized the minimum and maximum values, and the latter achieved high positive correlation values to the yield (Table 5). Five spectral indices, the NDVI, NIRv, EVI, wetness, and NDWI, have better and more linear relationships than the weather parameters as illustrated in Figure 12. Yet, the OSI values do not correspond well to the yield, in contrast to the LWP parameters, probably because each value of this method analyzes data collected over 10 or 15 days. With 24-27 days for each alfalfa growth cycle, it is not sufficient time for this kind of self-calibration method.

Discussion
Most of the spectral indices tested in this study were autocorrelated and highly correlated to the LAI and the yield regardless of their original target: growth or water stress. This issue was reported in previous studies. For example, the VARI, which was originally developed for mapping the vegetation cover [22], was found to have better relationships to the vegetation water content than the NDVI [25]. Another example is the NDWI, originally proposed to estimate the vegetation moisture, but found to have the best relationships with the LAI in corn, soybean, and almond [26,27]. However, none of these indices are sensitive to temporal changes in the crop water stress as illustrated in Figures 5-9. Hence, the suggested OSI (composed via ndvi2ndwi and ndvi2wetness) can be used more accurately. Utilizing it to depict the changes after an irrigation event is comparable to the water deficit index [12,17], yet it does not require thermal imagery. Its disadvantage, however, is the requirement of past data (see Figure 10), which suggests testing it via integration with other satellite sensors to achieve more frequent data. These tests should be in a controlled environment, to eliminate cloud cover problems or differences in sensor calibration.
The OSI values below −0.25 for the maximum alfalfa CWS correspond to the LWP values below −25, which are also reported as high water stress for this crop [24,28]. Several publications mentioned values of −10 to −15 LWP for alfalfa in non-stress conditions [7,11,29], which parallel the values we measured at the peak biomass ( Figures 5-8). Yet in the seventh cycle, the LWP values ranged from −16 to −28 (Figure 9), resulting in a very low yield compared to non-stress conditions ( Figure 11).
Another feature this method uncovers is the repeated alfalfa CWS cycles: high CWS immediately after the harvest; decreases with growth and irrigation; and increases again toward the next harvest. These CWS cycles, also seen by the LWP measurements ( Figures 5-8) and reported elsewhere [17], were broken when irrigation scheduling failed in the seventh growth cycle (Figure 9). The impact of irrigation on the alfalfa CWS has a lag in time, with the maximum LWP values 12-36 h after irrigation ( Figure 11). This difference in time between the irrigation and the peak crop moisture was also reported in cotton [30] and almond [31]. Furthermore, this lag in time could be used to explain the differences between the soil water stress (as measured by soil sensors and tensiometers), or Ks as defined by Allen et al. [1], and the crop water stress. As seen especially in Figures 7 and 8, the OSI is more associated with the CWS as measured by the LWP than the topsoil water stress (expressed by the SWC).
The reduced CWS 12-36 h after irrigation with no effect of the irrigation timing on growth (LAI) supports our motivation to develop the OSI. We hypothesized (after Moran (1994)) that irrigation causes a positive difference between the actual value of spectral CWS and the predicted value calculated from the spectral growth index. Our results for OSI, and especially the interpretation of the maps in Figure 10, highlight how the OSI better detects the CWS compared to the known spectral indices.
The OSI does not correlate well with the yield, although the yield is highly correlated with the LWP and weather parameters (Table 4). Nevertheless, the better interactions between the spectral indices and the yield and their linear relationships (Table 5 and Figure 12), reveal the importance of crop monitoring by remote sensing imagery [32]. As illustrated in Figure 9, the failure of the irrigation system in the last cycle empowers the information from the imagery. As so, any more knowledge we can obtain from the satellite images, such as the suggested OSI, will be helpful for growers on a global scale.
From a practical point of view, the OSI (in the form of ndvi2ndwi) should be used in two aspects: (1) spatially-to depict the spatial pattern of the crop water stress and (2) temporally-to outline deviation from standard practice, i.e., when time after irrigation is too long (Figures 8A and 9B) or when the soil was not dry out before harvest ( Figure 6A,B). The usefulness of the spectral indices should also consider the spatial and temporal aspects: except GB, all spectral indices outline the spatial and temporal changes in growth. Further, their maximum value for each cycle can be used to predict deviation from the standard yield. However, these indices' disadvantage is their insensitivity to the temporal changes of the water stress, which highlights the importance of the OSI.
Lastly, the crop measurements, and especially the LWP, are good indicators of the final yield. This practice, and the remote sensing data, are more effective than any weather parameter because of management practices, as shown in the system failure at the end of the 7th cycle. For implementation, the best time to define the maximum LWP, is not immediately after the irrigation, but rather about 24 h later.

Conclusions
This study suggests a new stress index termed optical stress index (OSI). The OSI is calculated using the NDVI and the NDWI and its main goal is to map the crop water stress on the pixel level. We collected alfalfa measurements such as LAI, LWP, and yield over 15 growth cycles to test it.
The main finding from this research is that the OSI has the potential to depict the temporal change in the alfalfa water stress from high stress after harvest, reduce by irrigation and growth and increase again toward the next harvest. Because the revisit time of Sentinel-2 is 5-days or more, the OSI failed to graphically illustrate the gap between the irrigation event and the peak of crop moisture, which occurs between 12 to 36 h after the irrigation. As so, the next studies should focus on increasing the temporal resolution and test water-stressed crops (such as cotton) to identify deviations from the irrigation strategy.

Acknowledgments:
We would like to thank the managers and workers of the Field Crops branch in Kibbutz Maoz-Haim.

Conflicts of Interest:
The authors declare no conflict of interest. Table A1. Sentinel-2 imagery (for 2021) utilized in this study.