Exploring the Use of Vegetation Indices for Validating Crop Transpiration Fluxes Computed with the MOHID-Land Model. Application to Vineyard

Abstract

The success of an irrigation decision support system (DSS) much depends on the reliability of the information provided to farmers. Remote sensing data can expectably help validate that information at the field scale. In this study, the MOHID-Land model, the core engine of the IrrigaSys DSS, was used to simulate the soil water balance in an irrigated vineyard located in southern Portugal during three growing seasons. Modeled actual basal crop coefficients and transpiration rates were then compared with the corresponding estimates derived from the normalized difference vegetation index (NDVI) computed from Sentinel-2 imagery. On one hand, the hydrological model was able to successfully estimate the soil water balance during the monitored seasons, exposing the need for improved irrigation schedules to minimize percolation losses. On the other hand, remote sensing products found correspondence with model outputs despite the conceptual differences between both approaches. With the necessary precautions, those products can be used to complement the information provided to farmers for irrigation of vine crop, further contributing to the regular validation of model estimates in the absence of field datasets.

1. Introduction

Irrigation is fundamental to fulfill crop water requirements in many regions of the world. Yet, inefficient practices often lead to the degradation of soil and water resources by promoting nutrient leaching, surface runoff and soil erosion, salt accumulation in the rootzone, and the eutrophication of water bodies with associated biodiversity loss. There is thus the need for minimizing environmental risks through the accurate estimate of crop water requirements and the definition of irrigation schedules (irrigation timing, duration, and quantity) that maximize agricultural water productivity and farmers’ income [1].

Several decision support systems (DSS) were developed over the last few decades to aid farmers in the decision-making of irrigation. One example is the IrrigaSys DSS developed by academics and stakeholders to support farmers in the Sorraia valley irrigation district, in southern Portugal [2]. When running in operational mode this DSS computes the weekly soil water balance based on hindcast and forecast weather data as well as inputs from farmers. It then suggests an optimized irrigation schedule for the week that follows. Assuring the reliability of predictions has been fundamental for the success of the DSS. It has also been its most challenging issue considering that the core engine of the system is the MOHID-Land model [3], which is highly complex in terms of parametrization of soil and crop state variables. The current support of the DSS to 103 agricultural fields distributed throughout the region only further presses the need for correctly accounting for a wide variety of crop management and soil conditions to guarantee the quality of the service.

Remote sensing, with the ability to cover large remote areas at different spatial and temporal resolutions, compatible with the characterization of the state and dynamics of several meteorological, vegetation growth, and hydrological variables, has been a preferential approach for evaluating the performance of the DSS [4,5]. Ramos et al. [4] analyzed the impact of assimilation of leaf area index (LAI) data derived from Landsat 8 imagery on model simulations. The relevancy of this study was because MOHID-Land computes the soil water balance using mass and momentum conservation equations, with LAI being critical to the correct partition of crop evapotranspiration (ET_c) rates into crop transpiration and soil evaporation. Their main conclusion was that modeling vegetation growth and the partition of the ET_c components at the plot scale could not depend solely on inputs from LAI data assimilation because estimates of the soil water balance could diverge substantially from the reality, thus confirming the need to use a proper crop dataset for model calibration. On the other hand, Simionesei et al. [5] used LAI data derived from the same satellite sensor for simply adjusting crop growth parameters to calibrate/validate model results, providing a more feasible solution for rapidly developing a large database of crop variability at the regional scale. This would, however, depend on the availability of a relationship to be established between field LAI data and that estimated from the satellite sensor.

This paper now focuses on an alternative solution for evaluating the performance of the DSS, in this case by comparing MOHID-Land estimates of the transpiration component in ET_c with those derived from remote sensing-based vegetation indices (VI). Several procedures already exist that make use of spectral VI for estimating crop evapotranspiration fluxes and irrigation needs in agricultural fields [6]. Although in most of these procedures the soil water balance is computed according to the FAO56 approach [7], they can eventually be adapted for obtaining calibrated relations between the VI and fluxes from a complex mechanistic model such as MOHID-Land, providing a scalable solution for assessing model behavior in fields with similar management.

Hence, the objectives of this study are (i) to simulate soil water contents in an irrigated vineyard (Vitis vinifera L.) using the MOHID-Land model during the 2018–2020 growing seasons; (ii) to compute the soil water balance for the study period; (iii) to establish relationships between actual transpiration rates computed by the MOHID-Land model and those derived from the normalized difference vegetation index (NDVI) and Sentinel-2 satellite (European Space Agency, European Union) imagery. Results of this study can thus help define better management practices to be implemented in IrrigaSys and improve its performance in the Sorraia Valley region.

2. Materials and Methods

2.1. Description of the Study Area

This study was carried out at Companhia das Lezírias, Samora Correia, Portugal (38.808° N, 8.900° W, 45 m a.s.l.) from January 2018 to October 2020. The climate in the region is dry subhumid, with mild winters and hot, dry summers. The mean annual precipitation is 669 mm, mainly concentrated between October and May, while the mean annual temperature is 16.8 °C. The weather data (Figure 1) for the study area was obtained from the local weather station (Figueirinha) and included daily precipitation (P; mm); maximum (T_max; °C) and minimum (T_min; °C) surface air temperatures; maximum (RH_max; %), mean (RH_avg; %), and minimum (RH_min; %) relative humidity; solar radiation (Rs; MJ m⁻² d⁻¹); wind speed measured at 2 m height (u₂; m s⁻¹). This information was then used to compute the reference evapotranspiration (ET_o, mm) following the FAO Penman–Monteith method [7]. The soil was classified as a Haplic Fluvisol [8], with loamy-sand texture in the top 60 cm layer and sandy texture in the bottom 60–100 cm layer.

The selected field, planted in 2008, was relatively flat (slope < 2‰), and part of a larger (130 ha) vineyard area. The field was a drip-irrigated plot, 5 ha in size, planted with different varieties of wine grapes, with Touriga Nacional being dominant and Castelão, Moreto, and Alicante Bouschet in smaller proportions. The plants were grown on vertical shoot positioned trellis, with wood pruning during the dormant period. Plants presented a row distance of 1.0 m and a row spacing of 2.8 m, thus a plant density of approximately 3571 plants ha⁻¹, with an orientation in the east–west direction. Irrigation was delivered through a drip system, with management practices performed according to the standard practices in the region and decided by the farmer. Drippers were spaced 1 m apart, and the drip line was placed on the trellis 0.5 m above the soil surface. The total water applied through irrigation summed 470, 625, and 465 mm in 2018, 2019, and 2020 growing seasons, respectively. The application depth during irrigation events varied from 1 to 12 mm. Soil water contents were continuously monitored in two locations at depths of 10, 20, 30, 40, 50, 60, 70, and 80 cm using EnviroPro MT (MAIT Industries, Bayswater North, Australia) capacitance probes.

2.2. The MOHID-Land Model

2.2.1. Model Description

The MOHID-Land model [3] computes the variable-saturated one-dimensional water flow in porous medium using the Richards equation and following a finite volume approach. The unsaturated soil hydraulic properties are described according to the van Genuchten–Mualem functional relationships [9,10]:

S_e(h) = (θ(h) − θ_r)/(θ_s − θ_r) = (1 + |α h|^η) ^−m

(1)

K(h) = K_s S_e^ℓ [1 − (1 − S_e^1/m) ^m] ²

(2)

where S_e is the effective saturation (L³ L⁻³), θ_r and θ_s denote the residual and saturated water contents (L³ L⁻³), respectively, K_s is the saturated hydraulic conductivity (L T⁻¹), α (L⁻¹) and η (-) are empirical shape parameters, m = 1 − 1/η, ℓ is a pore connectivity/tortuosity parameter (-), and h is the soil pressure head (L).

Crop evapotranspiration (ET_c, L T⁻¹) is computed from the product of the reference evapotranspiration (ET_o, L T⁻¹) following the FAO Penman–Monteith method and a crop stage-dependent coefficient (K_c) [7], and then partitioned into potential soil evaporation (E_{s model}; L T⁻¹) and potential crop transpiration (T_{c model}; L T⁻¹) as a function of the simulated LAI [11]:

T_{c model} = ET_c(1 − e ^{(−λ LAI)})

(3)

E_{s model} = ET_c − T_{c model}

(4)

where λ is the extinction coefficient of radiation attenuation within the canopy (-). Potential root water uptake values, given by the T_{c model} rates, are then linearly distributed over the root zone (z), creating the function T_{c model}(z), which may be diminished due to water stress [12]. Actual transpiration (T_{c act model}; L T⁻¹) rates are obtained by limiting the potential values using the piecewise linear model proposed by Feddes et al. [12]. This approach considers that the water uptake is at the potential rate when the pressure head is between h₂ and h₃, drops off linearly when h > h₂ or h < h₃, and becomes zero when h < h₄ or h > h₁ (subscripts 1–4 denote different threshold pressure heads). E_{s model} is limited by a pressure head threshold value to obtain the actual soil evaporation rate (E_{s act model}; L T⁻¹) [13]. LAI evolution is simulated as a function of crop stage, potential heat units (PHU) for plants to reach maturity, and plant stress [14]. In this study, the subscript model identifies the variables computed using the MOHID-Land model.

The MOHID-Land model further includes an irrigation scheduling tool for automatizing irrigation using a system-dependent boundary condition that triggers the application of water when a certain threshold pressure head (h_t) is reached in different grid cells of the rootzone domain. Irrigation then ceases after a second target pressure head (h₀) is obtained in the same grid cells. Since the root zone domain is typically defined by a large and variable number of grid cells, MOHID-Land further includes a series of constraints that prevent the application of meaningless irrigation amounts and countless irrigation events, namely a maximum irrigation pulse (I_max) and a minimum irrigation interval (I_int). The model is thus automatized for triggering irrigation whenever h drops below h_t in different cells of the root zone domain, supplying them enough water to reach h₀ in those same cells based on a predefined irrigation strategy. Further details on the MOHID-Land model can be found in Ramos et al. [3].

2.2.2. Model Setup, Calibration, Validation

The simulation period covered the 2018–2020 growing seasons. The soil profile was represented by one vertical column discretized into 20 grid cells, 1 m wide, 1 m long, and with variable thickness (0.025 m on the top to 0.625 m at the bottom) considering the depth of the simulation domain, the root zone, and the measured soil moisture data. ET_c rates were computed by multiplying daily ET_o values with the respective crop coefficients (K_c) for the initial, mid-season, and late-season stages. The K_c values for wine grapes listed in Allen et al. [7] were used as default. The upper boundary condition was then determined by E_{s model} and T_{c model} rates, and net irrigation and precipitation fluxes. The soil hydraulic parameters were initially set according to the texture classes of the soil horizons/layers following the class pedotransfer functions in Ramos et al. [15,16]. The crop growth parameters in Neitsch’s et al. [14] crop database were also used as default settings. Root water uptake reductions were computed by considering the following parameters: h₁ = −10, h₂ = −25, h₃ = −1000, h₄ = −18,000 cm [17]. Free drainage was used as the bottom boundary condition.

Model calibration was performed during the 2018 growing season following Ramos et al. [3]. The trial-and-error procedure was first used to calibrate the crop growth parameters by making them vary within reasonable ranges until deviations between modeled and field LAI data were minimized. The Copernicus Global Land Service LAI dataset derived from the Sentinel-3 sensor at 300 m resolution was assumed here to represent field conditions [18]. Only pixels covering the vineyard area and not affected by adjacent vegetation were used to extract LAI data. LAI values were also compared to the existing literature to evaluate their adequacy [19,20]. LAI data corresponding to the non-growing season were simply ignored since it represented grass growth covering the vine’s interrow during the rainfall period, but which dynamics could not be considered in this application. The interrow crop had also little influence on soil moisture data as the capacitance probes were installed in the vine rows. Then, the same trial-and-error procedure was adopted for calibrating the K_c values for the different crop stages as well as the soil hydraulic parameters for the different soil layers. These parameters were also made to vary within reasonable ranges until deviations between observed and simulated soil water contents were minimized. The connectivity/tortuosity parameter ℓ was not adjusted, being set to 0.5 following Mualem [9]. The parameters θ_s, α, η, K_s, and the maximum value of LAI (LAI_max) are identified as the most sensitive parameters during model calibration [3]. The calibrated crop growth parameters, K_c values, and soil hydraulic parameters were then used to validate model simulations of crop growth and soil water contents during the 2019 and 2020 growing seasons. Seasons were run separately, with end results from one season being updated at the beginning of the following according to measured data of soil water contents.

The goodness-of-fit indicators adopted for comparing field and simulated LAI values and soil water contents were the coefficient of determination (R²), the root mean square error (RMSE), the normalized RMSE (NRMSE), the percentage bias (PBIAS), and the model efficiency (NSE). R² values close to 1 indicate that the model explains well the variance of observations. RMSE and NRMSE values close to zero indicate small estimation errors and good model predictions [21]. PBIAS values close to zero indicate that model simulations are accurate, while positive or negative values indicate under- or over-estimation bias, respectively. NSE values close to 1 indicate that the residuals’ variance is much smaller than the observed data variance, hence the model predictions are good. On the contrary, if NSE is less than zero the model-predicted values are worse than simply using the observed mean [22].

2.3. Data Processing of Sentinel-2 Imagery

T_{c act model} rates, and respective actual basal crop coefficients (K_{cb act model} = T_{c act model}/ET_o ratio), computed by the MOHID-Land model were compared with those derived using a second approach, where those variables were estimated from a vegetation index (VI). Satellite sensors have been extensively used for estimating crop evapotranspiration with real-time single (K_c) or basal crop coefficients (K_cb) estimated from VI data. Pôças et al. [6] provided a comprehensive listing of the many types of K_c-VI and K_cb-VI relationships developed for annual and perennial crops. Examples of applications in vineyards can be found in Campos et al. [23,24] and Er-Raki et al. [25]. Landsat imagery (National Aeronautics and Space Administration Agency, Washington, DC, USA) has been the most frequently used satellite data to generate VIs to estimate those crop coefficients. However, its revisiting period of 16 days, which is often extended due to high levels of cloud cover, and the relatively coarse resolution (30 m) of its multispectral imagery has limited the use of that sensor for irrigation water management at the field scale. In this sense, the launch of the Sentinel-2 mission (European Space Agency, European Union) in 2017 represented a major step forward on the use of satellite data for remotely monitoring crop irrigation needs by offering a revisiting time of 5 days under the same viewing angle, and multispectral imagery at 10 m (visible and broad near-infrared spectrum), 20 m (red edge and narrow and short-wave infrared), and 60 m (atmospheric bands) resolution.

In this study, the red (Red; band 4) and near infra-red (NIR, band 8) bands from Sentinel-2 image tiles having less than 10% cloud cover were downloaded from the Copernicus Open Access Hub [26] for the study period. The images were subjected to atmospheric correction of the downloaded scenes using the Sen2cor software [27], which is a processor for Sentinel-2 Level 2A product generation and formatting, performing the atmospheric, terrain and cirrus correction of top-of-atmosphere Level 1C input data, and creating bottom-of-atmosphere corrected reflectance images. A total of 65 images were available during the 2018–2020 growing seasons. The NDVI was then computed for all available images as follows [28]:

NDVI = (NIR − Red)/(NIR + Red)

(5)

This VI was chosen as it is the most used for establishing K_c-VI and K_cb-VI relationships for annual and perennial crops in the literature [6]. Following Campos et al. [23], the NDVI was calculated on a pixel-by-pixel basis and averaged for the area surrounding the location of soil moisture probes (100 m long × 60 m width), avoiding field edge pixels (Figure 2).

The NDVI data were used to derive values of the basal crop coefficients and actual transpiration rates for each available image date as follows [6,29]:

K_{cb act NDVI} = min [K_{cb max}, K_{cb max}/f_{c max} (NDVI − NDVI_min)/(NDVI_max − NDVI_min)]

(6)

T_{c act NDVI} = K_{cb act NDVI} ET_o

(7)

where K_{cb act NDVI} is the actual basal crop coefficient computed from the NDVI (-), T_{c act NDVI} is the actual transpiration rate computed from the NDVI (L T⁻¹), K_{cb max} is the maximum value of basal crop coefficient (-), f_{c max} is the fraction of ground cover corresponding to the maximum K_cb (-), and NDVI_min and NDVI_max represent the maximum and minimum values of NDVI corresponding to bare soil and effective full cover, respectively (-). Thus, in this study, the subscript NDVI identifies the variables computed from satellite imagery. The K_{cb act NDVI} represents primarily plant transpiration as well as a residual diffusive evaporation component supplied by soil water below the dry surface [7,30]. As the VI reflects the actual vegetation cover conditions, the estimated values represent actual rather than standard conditions [31]. For each growing season, the K_{cb max} (0.48) was defined according to Rallo et al. [32]. The f_{c max} was obtained from the maximum LAI value (LAI_max) and crop height following Allen and Pereira [33]. The NDVI_min and NDVI_max values were specified based on the evolution of the NDVI over the growing seasons, being set to 0.10 (assumed to represent bare soil conditions) and 0.58 (corresponding to a 10% increase of the maximum observed NDVI value), respectively, in each season.

The K_{cb act model} and T_{c act model} values computed using the MOHID-Land model were thus compared with the corresponding values derived from satellite imagery (K_{cb act NDVI} and T_{c act NDVI}) by linear regression analysis. The estimates from the MOHID-Land model were considered the dependent variables since they were the objective of validation. The resulting linear regression models were analyzed using a cross-validation technique, in which data from two years were used as the training subset and data from the remaining year was used as the validation subset. As data from only three growing seasons was used, this procedure was repeated three times to include all possible combinations as calibration and validation subsets. For each validation test, the performance of regression models was assessed using the same goodness-of-fit tests referred earlier except for the NSE.

3. Results and Discussion

3.1. Model Parametrization

Table 1. Calibrated crop growth parameters.

Crop Parameter	Value
Optimal temperature for plant growth, Topt (°C)	20.0
Minimum temperature for plant growth, Tbase (°C)	8.0
Plant radiation-use efficiency, RUE [(kg ha−1) (MJ m−2)−1]	30
Total heat units required for plant maturity, PHU (°C)	3500
Fraction of PHU to reach the end of stage 1 (initial crop stage), frPHU,init (-)	0.05
Fraction of PHU to reach the end of stage 2 (canopy development stage), frPHU,dev (-)	0.15
Fraction of PHU after which LAI starts to decline, frPHU,sen (-)	0.55
Maximum leaf area index, LAImax (m2 m−2)	1.4
Fraction of LAImax at the end of stage 1 (initial crop stage), frLAImax,ini (-)	0.35
Fraction of LAImax at the end of stage 2 (canopy development stage), frLAImax,dev (-)	0.85
Maximum canopy height, hc,max (m)	1.5
Maximum root depth, Zroot,max (m)	0.8
Net radiation coefficient (-)	0.463
Photosynthetically active radiation coefficient (-)	0.650

presents the calibrated crop growth parameters for the vine. Like in previous applications of the MOHID-Land model [3,5], most default settings used for simulating crop growth needed to be modified to accurately describe field data. In this case study, LAI_max was set according to the maximum value extracted from the Copernicus Global Land Service LAI dataset [18] during the three growing seasons. The remaining parameters of the LAI curve were then modified using the data available in the calibration period. The maximum canopy height (h_c,max) was defined according to field observations while the maximum root depth (Z_root,max) was adjusted based on measured soil moisture profiles using the capacitance probes installed at the field plot. Lastly, the base temperature for vine growth (T_base), i.e., the minimum temperature required for crop development, was calibrated to 8 °C, which is slightly lower than the minimum threshold of 10 °C generally considered for vineyard in the literature [34]. The optimal temperature was set to 20 °C, in accordance with the existing literature.

Table 2. Calibrated soil hydraulic parameters.

Depth (m)	0–0.15	0.15–0.25	0.25–0.35	0.35–0.45	0.45–0.55	0.55–0.65	0.65–0.75	0.75–2.0
θr (m3 m−3)	0.057	0.057	0.057	0.057	0.057	0.057	0.037	0.037
θs (m3 m−3)	0.410	0.410	0.410	0.410	0.410	0.410	0.410	0.410
α (cm−1)	0.184	0.184	0.194	0.194	0.184	0.184	0.184	0.164
η (-)	2.0	1.9	1.8	1.8	1.8	1.8	1.8	1.8
ℓ (-)	0.5	0.5	0.5	0.5	0.5	0.5	0.5	0.5
Ks (cm d−1)	90.7	90.7	1339.2	907.2	565.9	565.9	907.2	1339.2

θ_r, residual water content; θ_s, saturated water content; α and η, empirical shape parameters; ℓ, pore connectivity/tortuosity parameter; K_s, saturated hydraulic conductivity.

presents the calibrated van Genuchten-Mualem parameters for different soil layers. Most parameters showed little variation with depth, with α and η reflecting as expected the characteristically relatively high values of coarse-textured soils [15,16]. The exception was the K_s, which values increased with depth on several orders of magnitude. On the other hand, the calibrated K_c values for the initial (K_{c ini} = 0.30), mid-season (K_{c mid} = 0.70), and late-season (K_{c end} = 0.45) crop stages agreed well with Allen et al. [7]. The K_c value for the initial stage was then adjusted for the frequency of the rainfall events and average infiltration depths, varying between 0.3 and 0.5. The K_c values for mid-season and late-season were adjusted for local climate conditions taking into consideration plant height, mean u₂, and mean RH_min for the period under consideration, varying from 0.69 to 0.71 and 0.44 to 0.46, respectively, during the different growing seasons.

3.2. Model Performance

Figure 3 shows the daily averages of the soil water contents measured at depths of 10, 40, and 80 cm during the 2018–2020 growing seasons and compares these values with the MOHID-Land simulations. Although measured and simulated data were compared for all monitored depths, results are presented graphically only for the depths mentioned above to limit the number of figures. Measured soil water contents increased sharply with precipitation to values close to saturation to then decrease also rapidly to lower levels due to the dominance of the gravitational gradient near saturation which promoted percolation, but also due to crop evapotranspiration. Irrigation was usually applied at small depths (1–12 mm) to maintain soil water contents relatively controlled during the growing seasons. Yet, large variations of soil moisture levels were still noticed during these periods, particularly at shallower depths.

Table 3. Statistical parameters for the agreement between model simulations and observed data.

Statistical Indicator	Soil Water Contents			LAI
	Validation 2018	Validation 2019	Validation 2020	Validation 2018	Validation 2019	Validation 2020
R2	0.671	0.683	0.676	0.680	0.445	0.842
RMSE (L L−1)	0.014	0.012	0.015	0.155	0.462	0.191
NRMSE	0.102	0.097	0.115	0.186	0.794	0.246
PBIAS (%)	0.508	−1.550	−0.614	−6.035	−67.403	14.667
NSE	0.653	0.615	0.658	0.639	−9.311	0.513

R², coefficient of determination; RMSE, root mean square error; NRMSE, normalized RMSE; PBIAS, percentage bias; NSE, modeling efficiency.

presents the statistical indicators used for evaluating the level of agreement between measured and simulated values. The MOHID-Land model performed reasonably well when simulating soil water contents during the 2018 calibration period. The value of R² was relatively high (0.671), showing that the model could explain most of the variability of the observed data. The errors of the estimates were quite small, resulting in a RMSE value of 0.014 m³ m⁻³ and a NRMSE value of 0.102. The PBIAS value was 0.51%, indicating no under or overestimation trend when simulating the measured data. The NSE value was also high (0.653), indicating that the residual variance was much smaller than the measured data variance. The parameters calibrated in 2018 were then validated during the 2019 and 2020 growing seasons, producing similar goodness-of-fit indicators. These were also within the range of values reported in the literature for soil water content simulations using the MOHID-Land model [3,5]. As such, the model was considered adequate to simulate soil water dynamics during the three growing seasons.

Figure 4 presents the simulated LAI using the MOHID-Land model as well as the Copernicus Global Land Service dataset used for representing field LAI data during the growing seasons [18]. As referred earlier, only satellite data corresponding to the vine growing seasons is shown, with data from the non-growing period being ignored since the development of the interrow plants could not be considered in this application. Table 3 describes the goodness-of-fit indicators obtained when comparing the simulated and field datasets during the three growing seasons. The correspondence between those datasets was quite satisfactory during the 2018 calibration period, resulting in relatively high R² (0.680) and NSE (0.639) values, and relatively low RMSE (0.155 m² m⁻²) and NRMSE (0.186) values. However, model simulations failed to reproduce crop growth during the 2019 validation season, with the main stressor affecting crop development during that year not being identified. This may have been the fire occurrence that affected the area west of the study vineyard at the end of the 2018 growing season, with effects on the 300 m resolution LAI product being then particularly noticed during the following year. Nevertheless, as satellite data in 2019 revealed an uncharacteristic trend of the LAI curve the model was also considered to be calibrated for simulating vine growth in the study area since the goodness-of-fit indicators were again satisfactory in 2020.

3.3. Soil Water Balance

Table 4. Soil water balance in the studied vineyard.

	Farmer’s Schedule			Optimized Schedule
	2018	2019	2020	2018	2019	2020
Inputs (mm):
P	512	241	399	512	241	399
I	470	625	465	55	98	105
CRmodel	0	0	0	0	0	0
ΔSSmodel	−1	9	1	−2	9	1
Outputs (mm):
Tc act model	82	105	117	115	139	146
Tc model	130	164	177	130	164	177
1-Tc act model/Tc model	0.37	0.36	0.34	0.11	0.15	0.17
Es act model	210	210	203	68	44	49
Es model	315	333	378	315	333	379
DPmodel	688	560	544	382	166	309
Runoffmodel	0	0	0	0	0	0
Error (%)	0	0	0	0	0	0

P, precipitation; I, irrigation; CR_model, capillary rise; ΔSS_model, soil water storage variation; T_{a act model}, actual transpiration; T_{c model}, potential transpiration; E_{a act model}, actual soil evaporation; E_{s model}, potential soil evaporation; DP_model, deep percolation. The subscript model corresponds to results computed with the MOHID-Land model. Model error = 100 (∑inputs − ∑outputs)/∑inputs.

shows the components of the soil water balance for each growing season. Figure 5 also presents the daily fluxes of ET_c, T_{c model}, T_{c act model}, E_{s act model} computed with the MOHID-Land as well as the daily irrigation and precipitation depths during the 2018–2020 growing seasons. The seasonal ET_c values ranged from 445 mm in 2018 to 555 mm in 2020. The seasonal ET_{c act model} values were smaller due to the impact of water stress on T_{c model}, varying from 293 to 320 mm in the same years. Although these values depended on the seasonal atmospheric demand as well as soil moisture conditions, they were found to be comparable to estimates of 321 mm in Campos et al. [23], 274–354 mm in Cancela et al. [35], and even 395–567 mm in Fandiño et al. [36] for different regions in Spain, 239–382 mm in Phogat et al. [37] for South Australia, and 320–480 mm in Wilson et al. [38] for California, USA.

In terms of ET_c partitioning, T_{c model} and E_{s model} rates varied from 130 to 177 mm and 315–378 mm, respectively, throughout the different seasons, being mostly dependent of the quality of the adjustment of the LAI curve to the LAI dataset derived from Sentinel-3. On the other hand, T_{c act model} and E_{s act model} rates ranged from 82 to 117 mm and 203–210 mm, respectively. Thus, the T_{c act model} accounted for only 28–37% of the ET_{c act model}. The T_{c act model} values were also comparably smaller than the range of 137–278 mm estimated in Fandiño et al. [36] or to the 183–263 mm reported in Cancela et al. [35] for vines subjected to different irrigation treatments. To refer that besides the obvious differences in climate conditions during different growing seasons, irrigation management, cultivars, soil textures, and soil water storage capabilities, the soil water balance in these latter studies was computed using the FAO56 dual-K_c approach [7,30], which estimation of the individual components of the ET_c generally results more precise as commonly recognized in the literature [31,39]. Nevertheless, the small T_{c act model} values estimated by the MOHID-Land model find some explanation in the coarse textures (loamy-sand to sandy) and limited water holding capacity of the study soil as well as irrigation management. The reduction of T_{c model} values was relatively large due to water stress (34–37%) and was usually observed from mid-June onwards (Figure 5) as the farmer attempted to improve berry and wine quality while controlling shoot vigor and ripening through deficit irrigation.

While T_{c act model} values were relatively small and crop stress was significant, the amount of water applied through irrigation was relatively high, varying from 465 mm in 2020 to 625 mm in 2019 (Table 4). Irrigation was used to compensate for the deficit in precipitation during the dry period, with seasonal values increasing when annual precipitation was less. However, due to limited soil water storage capacity percolation losses were quite high, ranging from 544 mm in 2020 to 688 mm in 2018. From these, 39–61% occurred during the irrigation period, corresponding to 45–58% of the total irrigation water applied. The occurrence of high percolation losses with the simultaneous observation of crop water stress only evidenced the imperative need for better irrigation protocols that duly consider the physical characteristics of the study soil.

Table 4 further shows the soil water balance in the studied vineyard following an optimized irrigation schedule using the MOHID-Land model. The threshold pressure head (h_t) for triggering irrigation was set at −1200 cm, i.e., slightly below the h₃ value in the Feddes et al. [12] model to induce some water stress to the plant. The target pressure head (h₀) was set at −100 cm, here assumed to represent field capacity. The maximum irrigation pulse (I_max) was set to 5 mm, with a minimum irrigation interval (I_int) of 1 day. Following these settings, the model proposed the net application of only 55, 98, and 105 mm in 2018, 2019, and 2020, respectively, which compared to farmer’s inputs are considerably lower. This inevitably led to lower percolation, with values ranging now from 166 mm in 2019 to 382 mm in 2018, resulting mostly from precipitation events. Additionally, the E_{s act model} decreased considerably since the soil surface was less moistened with the absence of successive irrigation events. On the other hand, T_{c act model} values increased between 25 and 40%, corresponding to a less pronounced, more controlled water stress (11–17%). Hence, this exercise exposed the advantages of using a modeling tool for optimizing irrigation schedules, helping to save substantial amounts of water in the process. Those low irrigation depths were only possible because the model was not subjected to constraints that usually occur in the decision-making of irrigation, being hard to match in field conditions. Additionally, model estimates much depended on how well soil hydraulic properties were able to represent actual flow conditions in the vineyard soil, how well the three-dimensional flow from the drip irrigation system was represented by a simple one-dimensional modeling approach as the one used here, how representative were the Feddes et al. [12] pressure head threshold values for describing the response of this particular variety to water stress, and how reliable was the partitioning of the ET_c and the computation of the soil water balance based on LAI evolution. Most of these modeling approach errors were already discussed in Ramos et al. [3,40,41] and will not be further extended here.

3.4. MOHID-Land vs. Remote Sensing

Table 5. Regression models between the MOHID-Land model and satellite sensor estimates of crop coefficients and transpiration fluxes (values in brackets correspond to the standard deviation of errors) ^a.

Model	Equation	R2 (-)	RMSE (b)	NRMSE (-)	PBIAS (%)
1	Kcb act model = 1.006 NDVI − 0.258	0.534 (0.050)	0.050 (0.008)	0.421 (0.113)	−0.555 (39.96)
2	Kcb act model = 0.875 Kcb act NDVI − 0.168	0.550 (0.054)	0.049 (0.007)	0.414 (0.097)	0.041 (39.84)
3	Tc act model = 2.216 NDVI − 0.612	0.365 (0.066)	0.432 (0.102)	0.839 (0.182)	57.057 (10.55)
4	Tc act model = 0.718 Tc act NDVI − 0.443	0.782 (0.069)	0.174 (0.022)	0.339 (0.119)	0.065 (10.07)

^a K_{cb act model} and K_{cb act NDVI}, actual basal crop coefficient computed, respectively, from the MOHID-Land and satellite imagery; NDVI, normalized difference vegetation index; T_{a act model} and T_{a act NDVI,} actual transpiration computed, respectively, from the MOHID-Land and satellite imagery; R², coefficient of determination; RMSE, root mean square error; NRMSE, normalized RMSE; PBIAS, percentage bias. ^b Units are the same as the variable units.

presents the regression models obtained after relating the MOHID-Land outputs with the corresponding ones derived from Sentinel-2 imagery. Figure 6 shows the scatterplots of those relations. Considering the mismatch between the simulated LAI curve and the field LAI dataset derived from the Sentinel-3 sensor during the 2019 growing season (Figure 4), the simulated LAI values were first correlated to the NDVI values computed from Sentinel-2 imagery to evaluate the quality of both datasets and whether the same disagreement was observed here. This was not verified as the comparison resulted in Pearson correlation coefficient (r) values of 0.80, 0.84, and 0.92 for 2018, 2019, and 2020 growing seasons (r = 0.65 for total data), with the correlation significant at the 0.01 level.

The direct relationships between the NDVI and the estimated K_{cb act model} or T_{c act model} values generally resulted in poor regression models (models 1 and 3). On the other hand, the scaling of the NDVI values between minimum and maximum values defined based on observed data improved the agreement between the MOHID-Land outputs and those derived from the satellite sensor. This improvement was only minor for K_{cb act model} (model 2), with the resulting R² value (0.55) remaining relatively low. This was attributed to the fact that the lowest NDVI value (0.232), computed at the end of the 2020 growing season (11/10/2020), corresponded to a K_{cb act NDVI} value of 0.16, which is close to the minimum K_cb value (0.15) to be expected for bare soils in the FAO56 procedure [7]. However, in mechanistic models such as the MOHID-Land model where LAI is used for the partition of ET_c rates, this does not occur. At the beginning or close to the end of the crop growing season, when LAI values are null or very small, the corresponding T_{c model}/ET_o ratio (or K_cb) is equally null or very small. For that, the K_{cb act model} values varied from 0 to 0.27 in MOHID-Land simulations while the range of variation for the K_{cb act NDVI} was from 0.16 to 0.48. This conceptual difference also affected the relationship between T_{c act model} and the T_{c act NDVI}, with the latter values being slightly higher than the former. Yet, the regression model obtained between these two parameters (model 4) was quite good, with the R² value (0.782) showing the ability of the model to explain most of the variability observed in the MOHID-Land dataset while the RMSE of 0.174 mm d⁻¹ showed the low error of the estimate. On the other hand, the large variation in the PBIAS revealed the conceptual differences in the two approaches that need to be considered in the IrrigaSys decision support system. Nevertheless, the relatively low RMSE suggested that the satellite approach if duly calibrated could be considered as a reliable approach for validating transpiration fluxes from the MOHID-Land model and assuring the reliability of the weekly recommendations issued by IrrigaSys to farmers. The assimilation of these independent predictions into IrrigaSys may, however, be dependent on the plot location as the number of Sentinel-2 images available in this case study was much lower compared to Ramos et al. [42], who had twice as many images for a study carried out 15 km north during the 2017–2019 growing seasons.

4. Conclusions

In this study, the MOHID-Land model was able to successfully estimate the soil water balance in an irrigated vineyard located in southern Portugal. Estimated actual transpiration rates were relatively low, revealing considerable water stress during the three monitored growing seasons. At the same time, percolation losses were high, showing the need for improved agricultural water use. Precise irrigation scheduling tools such as the one in MOHID-Land may thus contribute to ensuring optimum soil moisture levels for vine growth. Yet, the particularity of vine irrigation management is difficult to scale up in a decision support system such as IrrigaSys, much influenced by the characteristics of crop variety and fruit quality. Expert knowledge on these issues seems thus fundamental for setting up the model correctly and ensuring the accuracy of the irrigation schedules provided by the DSS. While model support can be validated by remote sensing data, the established relationships between model outputs and the NDVI were much dependent on experimental conditions, including irrigation management and climate conditions. This should be always considered when generalizing these products to other vineyards in the region. Nonetheless, taking the necessary precautions, Sentinel-2 data can well provide important information to validate model outputs during the decision-making of irrigation on a regular basis. In the absence of calibration datasets, the regression models developed in this study can be further helpful to rapidly identify fields covered by the DSS where the hydrological model may be acting poorly.