Seasonal Crop Water Balance Using Harmonized Landsat-8 and Sentinel-2 Time Series Data

Abstract

Efficient water management in agriculture requires a precise estimate of evapotranspiration ($ET$). Although local measurements can be used to estimate surface energy balance components, these values cannot be extrapolated to large areas due to the heterogeneity and complexity of agriculture environment. This extrapolation can be done using satellite images that provide information in visible and thermal infrared region of the electromagnetic spectrum; however, most current satellite sensors do not provide this end, but they do include a set of spectral bands that allow the radiometric behavior of vegetation that is highly correlated with the $ET$. In this context, our working hypothesis states that it is possible to generate a strategy of integration and harmonization of the Normalized Difference Vegetation Index ($NDVI$) obtained from Landsat-8 ($L8$) and Sentinel-2 ($S2$) sensors in order to obtain an $NDVI$ time series used to estimate $ET$ through fit equations specific to each crop type during an agricultural season (December 2017–March 2018). Based on the obtained results it was concluded that it is possible to estimate $ET$ using an $NDVI$ time series by integrating data from both sensors $L8$ and $S2$, which allowed to carry out an updated seasonal water balance over study site, improving the irrigation water management both at plot and water distribution system scale.

1. Introduction

Under the constant pressure of population growth, food consumption is increasing in nearly all regions of the world. The world population is expected to reach 9 billion by 2026 [1], which implies an increase in the worldwide area under irrigation of 300,000 km${}^2$, along with a 40% increase in water and energy demand over the next 20 years [2]. The sector with the greatest water use in Chile is agriculture, accounting for 78% of the total water availability [3], which supplies an irrigated area of 11,000 km${}^2$ [4]. However, in various regions of the country, water use rights exceed the actual availability of water resources, which has led to numerous regions being declared depleted in terms of both surface and groundwater [5]. In this scenario, it is clear that one of the main challenges of the 21st century is to increase agricultural water productivity [6]. Thus, precise information on agricultural water demands is crucial for efficient management of water and crop productivity [7,8]. This high-water-demand scenario has triggered a search for solutions to narrow the gap between irrigation water demand and availability in terms of quantity and quality through the use of new technologies [9,10,11].

Water resource management strategies must rely on the estimation of crop water demand. One of the most widely used methods to determine water demand is the estimation of evapotranspiration ($ET$). Several studies have been done to evaluate the effect of water applied in crop yield [12,13,14,15,16] to optimize the water management in agriculture. Nonetheless, this process is difficult to correctly quantify when dealing with large areas, as there is great spatiotemporal variability due to the complex interactions between the soil, vegetation, and climate [17]. Currently, $ET$ estimates are mainly based on observations from terrestrial weather stations. Several approaches have been proposed in the literature to estimate ET using weather station observations: (i) A two-step approach by multiplying the weather-based reference evapotranspiration $E{T}_r$ by crop coefficients ($K_c$) [18,19,20]. Crop coefficients are determined according to the in situ type of crop and the crop growth stage [19]. (ii) On the basis of the Penman–Monteith (P–M) equation [21], with crop to crop differences represented by the use of specific values of surface and aerodynamic resistances [22,23,24,25,26].

Presently it is possible to estimate $ET$ for different crops, providing spatial and temporarily distributed information over a wide area, using information gathered from aircraft or satellite platforms. Two main strategies for $ET$ estimation from remote sensing can be distinguished as follows: (i) Methods that use visible and near infrared sensors to extract a vegetation index ($VI$) and the radiative surface temperature to estimate the corresponding skin temperature [27,28]; (ii) Residual methods using the surface energy balance (SEB). These methods calculate $ET$ by subtracting sensible heat and soil heat fluxes from net radiation [29]. Among the best-known methods that belong to this category are the surface energy balance index ($SEBI$) [30], the two-source model ($TSM$) [31], the evapotranspiration mapping algorithm ($ETMA$) [32], and the Mapping Evapotranspiration at high Resolution using Internalized Calibration ($METRIC$) model [33], which is based almost entirely on the Surface Energy Balance Algorithm for Land ($SEBAL$) model developed by Bastiaanssen et al. [34]. $METRIC$ and $SEBAL$ estimate crop $ET$ by calculating the surface energy balance using spectral information from multispectral satellite images in the optical, near infrared, and thermal ranges. Only remote sensing imagery that provides spectral information in the thermal band may be used as input for these models. Unfortunately, most current satellite sensors do not provide this information, but they do include a set of spectral bands that allows the radiometric behavior of vegetation to be determined by focusing on the spectral contrast presented by plant cover in the red and near infrared bands [12]. Most $VI$ are based on this principle; one which allows multitemporal data series to be constructed, thus providing essential information on water consumption patterns in various crop types and helping keep information on different agricultural covers up to date, as well as for monitoring of the biophysical properties of plants, such as plant cover, vigor, and growth dynamics [35]. In particular, a strong relationship between $ET$ and the $VI$, known as the Normalized Vegetation Difference Index ($NDVI$), has been demonstrated [36]. In [6,37], the authors demonstrated that it is possible to estimate $ET$ using $NDVI$ using models trained previously with $ET$ maps obtained by the $METRIC$ model (considering visible, infrared, and thermal bands of Landsat 7 ETM+) and $NDVI$ maps.

Images from the Landsat Program have generally been used for water consumption analysis and crop yield estimation. This is because of their high-spatial resolution in the visible spectrum (30 m) and the inclusion of the thermal bands (60 m), which allow for the quantification of the evapotranspiration in large areas at a temporal resolution of approximately 16 days in a precise and consistent way [38]. This acquisition frequency, however, can become a problem, especially in areas where climate conditions do not always allow good quality data to be available. As one might expect, a higher temporal frequency would be better. Sentinel-2 from European Space Agency (ESA) provides good spatial (20–60 m) and temporal (five days) resolutions for the visible spectrum, but lacks the Thermal Infrared (TIR) band [38]. Recently, methods have been developed with the aim of taking advantage of the characteristics of both products (Landsat-8 and Sentinel-2), with efforts to generate harmonized time series of surface reflectances for land monitoring applications being especially relevant [39,40,41]. The combination of the products of both programs allows for an effective increase in spatial and temporal coverage, providing a greater availability of data to users [42].

In this context, our working hypothesis states that it is possible to prepare a seasonal water balance in crops using serial data from Landsat-8 and Sentinel-2, obtained through the integration and harmonization between the $NDVI$ of both sensors, known as $NDV{I}^{\prime }$. The expected result of this work is a $NDV{I}^{\prime }$ time series used to estimate $ET$ through specific adjustment equations for each type of crop, which allows one to continuously characterize the demand for water during an irrigation season.

2. Material and Methods

2.1. Study Site

The study was carried out in one of the characteristic zones of the Central Valley of Chile. It has an area of 70 km${}^2$ irrigated by the Convento Viejo reservoir, which has a water storage capacity of approximately 237 million m${}^3$. The area is located in the O’Higgins Region (252,626.83 E and 6,153,791.95 S, Zone 19, Datum $WGS84$) along the South Lolol Canal (Figure 1). It is characterized by a temperate Mediterranean climate (with winter rain), a condition that favors the development of various crops such as forestry, orchards, cereal, and table grapes and vine, with 85.8% of the region’s agroforestry area concentrated there [43]. The main agricultural land-use types identified in the study site were plums, olives, almonds, blueberries, table grapes and vines, industrial tomatoes, maize, wheat, cereals, and alfalfa. In this sense, for this work, the area was divided according to two types of cover: orchards (2.7 km${}^2$) and annual crops (1.1 km${}^2$) (Figure 2).

2.2. Satellite Imagery

For this study, satellite images obtained from the Landsat-8 ($L8$) (L1T processing level, radiometrically calibrated and orthorectified) and Sentinel-2 ($S2$) satellites (1C processing level, with upper atmosphere reflectance values and orthorectified) were used. In

Table 1. Technical characteristics of the Landsat-8 ($L8$) and Sentinel-2 ($S2$) satellite images.

Landsat-8 OLI			Sentinel-2 MSI
Temporal Resolution: 16 Days			Temporal Resolution: 5 Days
Bands ($\mathsf{\mu }$m)		Spatial Resol. (m)	Bands ($\mathsf{\mu }$m)		Spatial Resol. (m)
B1: Coastal	0.435–0.451	30	B1: Coastal	0.433–0.453	60
B2: Blue	0.452–0.511	30	B2: Blue	0.458–0.523	10
B3: Green	0.533–0.590	30	B3: Green	0.543–0.578	10
B4: Red	0.636–0.673	30	B4: Red	0.650–0.680	10
			B5: Vegetation Red Edge	0.698–0.713	20
			B6: Vegetation Red Edge	0.733–0.748	20
			B7: Vegetation Red Edge	0.773–0.793	20
			B8: NIR	0.767–0.908	10
B5: NIR	0.85–0.88	30	B8a: NIR	0.848–0.881	20
			B9: WV	0.931–0.958	60
B9: Cirrus	1.363–1.384	30	B10: Cirrus	1.338–1.414	60
B6: SWIR	1.567–1.651	30	B11: SWIR	1.539–1.681	20
B7: SWIR	2.107–2.294	30	B12: SWIR	2.072–2.312	20
B8: PAN	0.503–0.676	15
B10: TIRS	10.60–11.19	100
B11: TIRS	11.50–12.51	100

, the characteristics of the spectral bands of each satellite are presented. For the image selection, the absence of clouds over the study area during the agricultural season (December 2017 to March 2018) was considered. The images selected are shown in

Table 2. $L8$ and $S2$ satellite image dates (DOY: Day Of Year).

Date	Satellite	Date	Satellite	Date	Satellite
10 December 2017 (DOY 344)	$S2B$	19 January 2018 (DOY 19)	$S2B$	Mar/10/2018 (DOY 69)	$S2B$
20 December 2017 (DOY 354)	$S2B$	08 February 2018 (DOY 39)	$S2B$	Mar/15/2018 (DOY 74)	$S2A$
25 December 2017 (DOY 359)	$S2A$	13 February 2018 (DOY 44)	$S2A$	Mar/20/2018 (DOY 79)	$S2B$
30 December 2017 (DOY 364)	$S2B$	18 February 2018 (DOY 49)	$S2B$	Mar/25/2018 (DOY 84)	$L8$-$S2A$
04 January 2018 (DOY 4)	$L8$	21 February 2018 (DOY 52)	$L8$	Mar/30/2018 (DOY 89)	$S2B$
09 January 2018 (DOY 9)	$S2A$	05 March 2018 (DOY 64)	$S2A$

.

2.3. Methods

The workflow used to achieve the research objective is shown in Figure 3. At the start, the spectral bands from the $L8$ and $S2$ satellites were atmospherically corrected using the Semi-Automatic Classification plugin in the open-source software QGIS 2.8.9. Then, the co-registration of the $L8$ and $S2$ images was carried out with the AutoSync tool in Erdas IMAGINE v.2011. This tool algebraically recognizes the coordinates of common points between the two satellite images. To carry out this process, the $L8$ satellite image was chosen as a reference, such that the other images were geometrically fitted to it. Subsequently, the subsampling of $S2$ from 10 m to 30 m pixel size was carried out using the linear interpolation method.

To obtain the harmonization between the $NDVI$ maps of the $L8$ and $S2$ sensors, the first step is to generate the $NDVI$ maps during the season for both sensors, applying Equations (1) and (2), respectively; where $NDV{I}_{L8}$ and $NDV{I}_{S2}$ are the normalized difference vegetation indices for the $L8$ and $S2$ satellites, respectively.

$$NDV{I}_{L8}=\frac{B5-B4}{B5+B4}$$

(1)

$$NDV{I}_{S2}=\frac{B8-B4}{B8+B4}$$

(2)

For $NDV{I}_{L8}$ and $NDV{I}_{S2}$ harmonization, images from coinciding dates or with minimal temporal separation during the analyzed season were selected. For each pair of images, 618 sampling points in the annual crops and 638 in the orchard cover were selected. These points were chosen at random, and those out of the confidence interval were eliminated. The confidence interval is defined as absolute value difference between $NDV{I}_{L8}$ and $NDV{I}_{S2}$ for each pair of identified points. This value must be less or equal to the average value of difference. With the sampling points selected, a scatter plot for both agricultural covers was plotted, with the x-axis representing the $NDV{I}_{S2}$ variable and y-axis the $NDV{I}_{L8}$ variable. The fitting equation for the two variables was obtained from linear least squares regression, to which the 1:1 line represents perfect agreement between the $NDVIs$. In addition, for each graph, the Pearson correlation coefficient (r) was determined. The fitting equation that relates the harmonized vegetation index (${\widehat{NDVI}}_{L8}$) is presented in Equation (3).

$${\widehat{NDVI}}_{{L8}_i}=a_i^{NDV{I}_{L8}}·NDV{I}_{S{2}_i}+b_i^{NDV{I}_{L8}}$$

(3)

where i is the type of agricultural land-cover considered (orchards or annual crops), and $a_i^{NDV{I}_{L8}}$ and $b_i^{NDV{I}_{L8}}$ are the coefficients of linear fit between $NDV{I}_{L{8}_i}$ and $NDV{I}_{S{2}_i}$ for the analyzed season. Thus, the multimodal $NDVI$ time series is composed of a total of three $NDV{I}_{L8}$ and 14 ${\widehat{NDVI}}_{L8}$ maps, and is called $NDV{I}^{\prime }$.

Following the workflow outlined in Figure 3, linear fitting between $E{T}_{L8}$ and $NDV{I}_{L8}$ for each $L8$ image (Table 2) was carried out. The $E{T}_{L8}$ maps at $30m$ spatial resolution were generated from $L8$ images throughout the agricultural season using the $SEB$ methodology proposed in Allen et al. [33] and implemented and validated by Fonseca-Luengo et al. [44]. The fit equation for five defined type of annual crops and six types of orchard covers, considering all $L8$ images is as follows:

$$E{T}_{L{8}_{k,j}}=a_{k,j}^{E{T}_{L8}}·NDV{I}_{L{8}_{i,j}}+b_{k,j}^{E{T}_{L8}}$$

(4)

where k represents the different types of agricultural land use considered in Figure 2 (plums, olives, almonds, blueberries, table grapes and vines, industrial tomatoes, maize, wheat, cereals, and alfalfa), and j is the number L8 satellite images (DOY 4, DOY 52, and DOY 84, Table 2). Meanwhile, $a_{k,j}^{E{T}_{L8}}$ and $b_{k,j}^{E{T}_{L8}}$ are the coefficients of linear fit. Thus, for each map of $E{T}_{L8}$ and $NDV{I}_{L8}$ during the analyzed season, different fit equations for each orchard and annual crop type in the study area were obtained.

Finally, with Equation (5), it is possible to estimate evapotranspiration ($\widehat{ET}$) for different agricultural land uses considered, by means of the $NDV{I}^{\prime }$ time series.

$${\widehat{ET}}_{k,m}=a_{k,j}^{E{T}_{L8}}·NDV{I}_{i,m}^{\prime }+b_{k,j}^{E{T}_{L8}}$$

(5)

where m represents the total amount of $S2$ images that were captured 30 days before and 30 day after the day j when the $L8$ image were captured.

To determine daily evapotranspiration during the season, it was necessary to obtain the actual crop factor ($Fc$) from the ratio between ${\widehat{ET}}_{k,m}$ and reference evapotranspiration ($E{T}_r$), estimated using the Penman–Monteith model [19] at the time the images were captured. In Equation (6), the $Fc$ for each type of cover (k) and the whole series (m) is described:

$$F{c}_{k,m}=\frac{{\widehat{ET}}_{k,m}}{E{T}_{r,m}}$$

(6)

Thus, daily estimated evapotranspiration $\widehat{E{T}_c}$ was the product of $E{T}_r$ and the $Fc$ obtained in each of the satellite images and extrapolated to the analyzed agricultural season.

Whereas, potential evapotranspiration ($E{T}_p$) was determined from the potential crop factor ($Fcp$), obtained from the literature [8,45], associated with the percentage of coverage during plant development. However, in our proposed method, this factor was adjusted to the $NDVI$ values obtained from the satellite images during the season. This factor was called $NDVI$ potential crop factor ($Fc{p}^{NDVI}$) and the way to obtain it was through the linear relationship between the $Fcp$ and the $NDV{I}^{\prime }$ of each crop. For this, the extreme values of $Fcp$ were related to the $NDV{I}^{\prime }$ values for each crop, while for the orchards, the difference was made between young trees and adult trees. The $Fc{p}^{NDVI}$ is described in Equation (7).

$$Fc{p}_{k,m}^{NDVI}=a_{k,m}^{Fcp}·NDV{I}_{i,m}^{\prime }+b_{k,m}^{Fcp}$$

(7)

Thus, daily $E{T}_p$ was the product of $E{T}_r$ and adjusted $Fc$ during the agricultural season for each crop in the study area.

To determine the water balance during the agricultural season, the estimated water demand and the potential demand of each of the annual crop and orchard types in the study area were determined. The two water demands are described in Equations (8) and (9), respectively.

$$EstimatedDeman{d}_{k,m}=\sum _{d=1}^n{\left({\widehat{ET}}_c\right)}_{k,m}^n$$

(8)

$$PotentialDeman{d}_{k,m}=\sum _{d=1}^n{\left(E{T}_p\right)}_{k,m}^n$$

(9)

where d is the start day and n the final day of the analyzed agricultural season. Thus, the seasonal water balance is composed of a comparison between estimated water demand and potential demand of each of the crops in the studied area.

The methodology evaluation was carried out by comparing the estimated water demand with the volume of water applied during the analyzed period. To this end, annual crops and orchards representative of the study site were selected. Meanwhile, information on the volumes applied was compiled from the records of the volume meters installed on the different farms. To compare the two variables, statistical indicators used root mean square error ($RMSE$), which measures the variation of the estimated values with respect to the observed values, and the bias indicator ($BIAS$), which provides information on the tendency of the model to over- or underestimate a variable.

3. Results

In

Table 3. Image dates used to fit the harmonization between $L8$ and $S2$.

Landsat 8		Sentinel 2
04 January 2018 (DOY 4)	⟶	09 January 2018 (DOY 9)
21 February 2018 (DOY 52)	⟶	18 February 2018 (DOY 49)
25 March 2018 (DOY 84)	⟶	25 March 2018 (DOY 84)

, the image dates used to fit the harmonization between $L8$ and $S2$ are shown.

The scatter graph between $NDV{I}_{S2}$ and $NDV{I}_{L8}$ for the three pairs of images for the orchards and annual crops are shown in Figure 4.

Using Equation (3), the harmonized relation for $NDV{I}_{L8}$ and $NDV{I}_{S2}$ to orchards is

$${\widehat{NDVI}}_{L8}=0.8288·NDV{I}_{S2}+0.1172.$$

(10)

In addition, the harmonized relation for the annual crops is

$${\widehat{NDVI}}_{L8}=0.8588·NDV{I}_{S2}+0.0997.$$

(11)

In Figure 4, it is possible to observe that the relationship between $NDV{I}_{L8}$ and $NDV{I}_{S2}$ is close to 1:1 for both cover types. In addition, the Pearson coefficient values (r) were 0.88 for the orchards and 0.94 for the annual crop cover. The difference observed between orchards and the annual crops could be due to the partial cover of the orchards, a topic that has been investigated by Souto et al. [46].

In addition, in Figure 4, it can be seen for both cover types that $L8$ tends to give large values of $NDVI$ compared to $S2$, with values almost up to 0.6. Although Tello et al. [47], found a correlation coefficient (r) of 0.99 for the $NDVI$ in agricultural crops, this relationship was determined without distinguishing between different cover types.

The available $NDV{I}^{\prime }$ maps during the growing season are shown in Figure 5. The time series was composed by $NDV{I}_{L8}$ and ${\widehat{NDVI}}_{L8}$, making a total of 17 images during the agricultural season.

Using this multimodal $NDV{I}^{\prime }$ time series, the $\widehat{ET}$ was estimated for each orchard and annual crop in the study site.

The results show that there is a significant correlation between $E{T}_{L8}$ and $NDV{I}_{L8}$ in each of the orchards and annual crops on the three analyzed days. This is reflected by the Pearson correlation coefficient (r) values, which overall are above 0.8. Upon comparing this correlation to what appears in the literature, it can be seen that the r values are very similar in [48,49,50]; this being dependent on the type of cover under study. Thus, a total of 18 and 15 fitting equations for the orchards and annual crop agricultural covers, respectively, were obtained. Using each of these fitting equations, it was possible to estimate $\widehat{ET}$ for the season which extends from December 2017 to March 2018, specifically for each type of orchard and annual crop type in the study area. Nonetheless, $\widehat{ET}$ estimation could have been improved if more images from the season had been obtained for the $E{T}_{L8}$–$NDV{I}_{L8}$ fitting, which in our study area was impeded by the presence of clouds.

By replacing the coefficients from

Table 4. Coefficients of linear fit between evapotranspiration ($E{T}_{L8}$) and $NDV{I}_{L8}$ and Pearson correlation coefficient (r) for orchards and annual crops on the three analyzed days.

		04 January 2018			21 February 2018			25 March 2018
		DOY 4			DOY 52			DOY 84
		${\mathit{a}}^{\widehat{{\mathit{NDVI}}_{\mathit{L}\mathbf{8}}}}$	${\mathit{b}}^{\widehat{{\mathit{NDVI}}_{\mathit{L}\mathbf{8}}}}$	$\mathit{r}$	${\mathit{a}}^{\widehat{{\mathit{NDVI}}_{\mathit{L}\mathbf{8}}}}$	${\mathit{b}}^{\widehat{{\mathit{NDVI}}_{\mathit{L}\mathbf{8}}}}$	$\mathit{r}$	${\mathit{a}}^{\widehat{{\mathit{NDVI}}_{\mathit{L}\mathbf{8}}}}$	${\mathit{b}}^{\widehat{{\mathit{NDVI}}_{\mathit{L}\mathbf{8}}}}$	$\mathit{r}$
Orchards	Almonds	9.6331	−2.2964	0.95	9.5884	−2.6277	0.89	5.5133	−0.4960	0.79
	Blueberries	4.970	−0.0216	0.84	3.5310	0.2247	0.84	2.6192	0.7221	0.80
	Plums	7.6880	−1.3594	0.77	6.3071	−1.1207	0.84	5.0152	−0.5987	0.79
	Olives	8.0158	−1.8227	0.89	9.8622	−2.7053	0.77	7.5290	−1.3265	0.77
	Table grapes	9.5418	−3.0590	0.94	3.1905	0.9724	0.55	3.3327	0.3498	0.86
	Vine	9.2916	−2.5368	0.96	6.0001	−1.0156	0.93	3.1834	0.9047	0.84
Annual Crops	Corn	8.1881	−1.5872	0.97	5.0615	−0.4432	0.90	2.2771	1.3574	0.70
	Alfalfa	5.3963	−0.4837	0.80	4.4676	−0.6641	0.82	3.1753	0.0803	0.81
	Cereals	8.1596	−1.7149	0.98	6.6392	−1.3811	0.93	6.3676	−0.8058	0.82
	Wheat	7.8890	−1.5880	0.95	5.5146	−0.7694	0.93	3.5682	0.2759	0.63
	Tomato	8.9797	−2.1934	0.98	5.8506	−1.0944	0.94	4.2877	−0.1040	0.89

in Equation (3), the $\widehat{ET}$ maps were obtained for each crop of the study site (Figure 6). In Figure 6, it is possible to observe the temporal and spatial distribution of $\widehat{ET}$ represented by monthly images, which was estimated using the different $E{T}_{L8}$–$NDV{I}_{L8}$ fitting equations for each orchard and annual crop type in the study area.

The overall results show that $\widehat{ET}$ presents the highest values on 20 December and 19 January, in the range of 5 and 6 mm·day${}^{-1}$. Meanwhile, on 13 February and 10 March the values are between 4 and 5 mm·day${}^{-1}$ This lessening of $\widehat{ET}$ starting on 13 February can be attributed to the senescence period, particularly for the annual crops, which present different phenological periods than the orchards.

$\widehat{ET}$ of industrial tomatoes varies from 3 mm·day${}^{-1}$ in December to up to 6 mm·day${}^{-1}$ in January, a variation that coincides with the development of the crop as it reaches maturity in late February and early March. However, the difference in $\widehat{ET}$ of plums is lower than that of industrial tomatoes, with values between 3 and 5 mm·day${}^{-1}$.

In general, orchards presented a relatively homogeneous water demand over the irrigation season because the canopy cover does not change during the season. It is important to point out that $ET$ is closely related to the cover percentage of orchards, as well as annual crops [8,45].

The selected crops are identified in Figure 7 to show the behavior of the daily $\widehat{E{T}_c}$ during the analyzed period, which were obtained from the $\widehat{ET}$ maps that were generated during the agricultural season following the proposed methodology.

The results shown in Figure 8 show the variability in water consumption through $\widehat{E{T}_c}$ of the different orchards and annual crops, showing that in general the analyzed orchards present a similar $\widehat{E{T}_c}$ behavior during the season. However, the behavior of table grapes presents higher values than those of the other orchards, reaching 7 mm·day${}^{-1}$ in December. Additionally, it is observed that during the first days of December orchards obtained the highest $\widehat{E{T}_c}$ values recorded during the season; attributed to better irrigation management during this period or possibly the contribution of soil evaporation rather than the transpiration of the trees themselves.

In the case of the annual crops, the behavior of the $\widehat{E{T}_c}$ is extremely variable among them, reflecting the natural differences in the water needs of each crop, with corn and the industrial tomato registering the highest values of $\widehat{E{T}_c}$, in a range varying between 3 to 6 mm·day${}^{-1}$.

The global results of water demand for the crops estimated with our proposed methodology (Figure 9) show that the applied volume of water to the crops was lower than the potential demand for all the season. This indicates that the water applied to different crops was insufficient or inadequate during the season.

In the particular cases of olives, table grapes, and maize (Figure 9, B, D and H), the estimated water demand was practically identical to the potential demand, showing that the amount of water applied during the season was the required one. This scenario is the least frequent in the studied area, even though the availability of the resource in the sector does not show restrictions. However, the rest of the crops showed an estimated water deficit, ranging from 6% to 48% for olive trees (B) and wheat (I), respectively.

The water deficit obtained by the different crops during the analyzed period can generate considerable damage to its yield. In the case of wheat, a water deficit causes a decrease in grain yield, especially when water scarcity occurs during the sensitive stages of growth (the flowering and grain filling stage) [51]. In blueberries, water deficit in any period of its phenological stage affects the vegetative growth, substantially reducing the potential yield, and in addition, it is a determining factor for the production in the following season [45,52]. Almond trees are considered drought tolerant, but irrigation scarcity is critical to the production of high yields with high quality nuts, mainly in the preharvest stage [53]. In plum trees, water stress in the final stages of fruit growth significantly decreases their size, but accelerates ripening, and increases the concentration of sugars [54]. For vines, water deficit in the stages after grapes veraison can improve the quality of the wine with almost no reduction in yield [55,56]. On the other hand, the imposition of severe drought stress in inappropriate phenological stages can cause a significant decrease in yield and even a decrease in quality in extreme cases [57]. As we have seen, the response of annual crops and orchards to a lack of water is different in each case and the decisions that must be taken in irrigation management must be appropriate for each one of them.

The results shown in Figure 10 show that the estimated volume of water applied for maize, plums, and olive trees (A, B, C, and D) satisfies the estimated water requirements through the potential demand. In the case of the adult plum tree (B), the water applied during the season is less than the amount of water required by the orchard. In addition, it is shown that the potential demand obtained for the adult plum tree was lower than that of the young plum, which can be attributed to the lack of water applied during the season or inadequate water management.

Olive trees in general are resistant to a lack of water and show a high recovery capacity after prolonged periods of drought, but adequate irrigation management is required to produce economical yields. However, a certain degree of stress improves oil quality [58,59].

For industrial tomatoes (F), the applied volume exceeded the estimated volume by 1300 m${}^3·$ha${}^{-1}$, but this is almost identical to the potential volume. This shows that the plants only evapotranspire the amount of water that is able to be retained in the soil at the root zone and not the excess applied water. Excessive irrigation in tomato plants can cause excessive leaf growth, plants with high vegetative vigor tend to produce low quality fruit, and even in some varieties, large variations in moisture levels in the soil during fruit ripening can cause fruit cracking, spots, rot, and variation in size and shape [60].

Statistically, the total estimated water demand through the proposed methodology based on the volume applied had an RMSE of 0.6 mm·day${}^{-1}$ and a BIAS of −0.4 mm·day${}^{-1}$.

Although the methodology implemented is able to effectively reproduce the amount of water applied, the considerable lack of water during the vital stages of plant development is demonstrated. Both water scarcity and inadequate applications can be combined with an uneven distribution of water in the field, often due to the low uniformity of irrigation systems. This especially occurs when irrigation systems are not properly designed, operated, and maintained. This can cause a substantial decrease in crop yield and economic income [61].

4. Conclusions

On the basis of the obtained results, we conclude that it is possible to construct an $NDVI$ time series by integrating $L8$ and $S2$ data, allowing the estimation of $\widehat{ET}$ during an agricultural season and suitably demonstrating the responses of crops to irrigation excess and shortage problems associated with water management. This is based on the Pearson correlation coefficient ($r>0.8$) values that were obtained for the orchards and annual crop cover. In addition, a comparison between estimated water demand and information on potential requirements meant that the water balance was kept up to date, according to the water needs of different agricultural cover types. This can encourage better management of applied water levels and production in different stages of the irrigation season. Further research on this topic, focusing on potential water demands for orchards and annual crops, will allow for better and more precise application of the proposed method.