Analyzing the role of soil and vegetation spatial variability in modelling hydrological processes for irrigation optimization at large scale

Abstract

The main purpose of this paper was to study the effect of spatial variability of soil hydraulic properties (HP) and vegetation parameters (VP) (e.g., leaf-area index, LAI, and crop coefficient, K_c) on modelling agro-hydrological processes and optimising irrigation volumes at large scale. Based on this analysis, the effect of partly overlooking the spatial variability of soil HP and/or VP inputs was verified on a 140 ha irrigation sector in “Sinistra Ofanto” irrigation system in Apulia Region, Southern Italy. Five soil profiles were excavated and the HP were measured in all the soil horizons. Additionally, measurements of soil HP were taken in the surface soil layer in ninety sites distributed over the whole irrigation sector. All the HP measurements were carried out using tension infiltrometer. Remote sensing applications were used to obtain LAI and K_c using European Space Agency’s (ESA) Sentinel-2 images with 10 m resolutions. First, distributed (on ninety polygons with an average area of about 1.5 ha) optimal irrigation volumes and related deep percolation volumes at a depth of 80 cm, were computed using an agro-hydrological model and accounting for the actual observed variability of soil HP and VP inputs. The sector scale irrigation and deep percolation volumes were obtained by aggregating the distributed irrigation volumes. This was considered as the reference scenario (hereafter DVS—Detailed Variability Scenario). Then, reduced variability scenarios (hereafter RVS—Reduced Variability Scenario) were considered, where the information on the actual spatial variability of the soil HP and VP was gradually overlooked to find the minimum data set needed to still have sector scale irrigation volumes and related deep percolation volumes comparable to those obtained under the DVS. Results showed that overlooking VP (RVS-VP) variability did not significantly change the optimal irrigation volumes and the deep percolation fluxes. By contrast, neglecting the HP variability (RVS-HP) showed significant effects on both the irrigation and percolation volumes compared to the DVS. The main practical finding was that, at least for the area investigated in this study, hydraulic characterization of one soil profile in an area of approximately 30 ha provides sector scale irrigation volumes and percolation fluxes comparable to those obtained under the DVS, thus by accounting for all the observed local variability.

Introduction

Several well-known reasons related to issues such as water and environmental quality, energy costs, economic policies, impose a more parsimonious and effective use of water resources. As the agricultural sector is responsible for most of water consumptions (more than 70% in the Mediterranean area, for example), a more efficient irrigation water management is now indispensable (Coppola et al. 2019b). Actually, there is a general perception that agricultural water use is often wasteful, frequently with losses exceeding 50% (Hsiao et al. 2007), coming partly from dilapidated water distribution networks, but more frequently coming from inefficient irrigation techniques at farm scale and scarce management capability at the larger irrigation district scale. Moreover, water losses are generally accompanied by agro-chemicals (for example, nitrates) losses, with additional serious impacts on quality of surface and groundwater bodies. These losses also impact on the energy consumption when the irrigation network is served by a pumping system.

Optimising irrigation water use may be achieved at both farm and irrigation network scales. At farm scale, sensors to reliably measure soil water content and/or water potential, such as sensors based on Time Domain Reflectometry (TDR) and tensiometers, are now available allowing to establish irrigation timing and volumes, even if their use may be limited to the field scale. Other approaches are based on measuring the upward fluxes above the canopy to estimate evapotranspiration (Paciolla et al. 2021). Recently, proximal sensors based on Electromagnetic Induction (EMI) technique are increasingly used to estimate water storage in the root zone, along with the solute concentration. As these sensors may be moved at soil surface without inserting probes in the soil profile, they allow a rapid monitoring of relatively large areas (farm scale), so to also establish farm scale variability of water available in the root zone and its salinity (see, for example, Kachanoski et al. 1988; Hedley et al. 2013; Coppola et al. 2016; Dragonetti et al. 2018, 2022; Farzamian et al. 2021). Airborne sensors have also been used in an attempt to estimate water stored in the soil thickness explored by roots by measuring the surface soil thermal inertia, which is strictly related to the soil water content (see, for example, Coppola et al. 2007; Minacapilli et al. 2012).

At the irrigation network scale, the issue is to provide irrigation water resources managers with tools aiming to optimize the scheduling in terms of irrigation timing and discharges to be delivered to the different irrigation sectors in an irrigation district. This may be desirable for reasons related to the hydraulics and the pressure head distribution in the irrigation network (Zaccaria et al. 2006; Coppola et al. 2019b), but also to optimize water allocation to different irrigation sectors, characterized by different soils and crop distribution, for example, under conditions of water scarcity.

In this perspective, agro-hydrological models can be used as decision support tools to quantify and schedule optimal irrigation volumes at large scale, while minimizing deep percolation fluxes of water and nutrients, using soil, vegetation and meteorological parameters as inputs. In their more evolute versions, agro-hydrological models use mechanistic approaches to calculate water flow and solute transport in the soil–plant–atmosphere system (Jarvis et al. 1991; van Dam et al. 1997; Abrahamsen and Hansen 2000; Šimůnek et al. 2008; Coppola et al. 2011; Coppola et al. 2015; Kroes et al. 2017). Water flow is generally described by the Richards equation (RE). In distributed approaches, the single field is considered as the elementary simulation unit of an irrigation sector. Richards’ equation is thus applied to describe 1D vertical water flow processes in each of the fields in the sector, to estimate field scale irrigation volumes and related deep percolation fluxes. The latter are then aggregated to give the sector scale irrigation volumes and percolation fluxes (see Coppola et al. 2019b, as an example). This distributed approach allows to obtain sector scale irrigation volumes by accounting for the actual spatial variability of the soil HP and VP in the sector. However, it is generally considered to be impracticable as it requires measuring HP and VP in several sites and carrying out a high number of simulations. Nowadays, obtaining distributed VP is not an insurmountable problem, as they may be quite reliably estimated by remote sensing techniques (D’Urso 2010; D’Urso et al. 2010; Calera et al. 2017). To the contrary, soil HP are known to vary considerably in space (Comegna et al. 2010, 2013; Severino et al. 2010; Terribile et al. 2011; Severino et al., 2017) and frequently require cumbersome techniques. Consequently, finding distributed datasets for modelling hydrological processes at large scale remains one of the most frustrating issues and has frequently been eluded by replacing measured HP with pedotransfer functions (PTFs), with all the complications related to their physical significance and realistic modelling (Hassan et al. 2022).

An additional issue in using distributed models is to identify the simulations grid size (the size of the elementary simulation area) able to effectively catch the effects of variability in the simulations results. The grid size must be small enough to include most of the variability effects, large enough to minimize the number and time of simulations. The grid size and the number of distributed input data for simulations are obviously inversely related.

An alternative way to reduce the number of simulations and solve the issue of the discretization of an irrigation sector area could be to solve the Richards’ equation directly at the sector scale with effective HP coming from some averaging schemes (Zhu and Mohanty 2002; Zhu et al. 2007; Severino and Coppola 2011) or ensemble HP deduced from the statistical distributions of the hydraulic parameters (Coppola et al. 2009). However, even in these cases, the problem remains of establishing the number of soil measurements to obtain robust effective property estimations. Also, effective properties are known to change according to hydrological process to be described (Coppola et al. 2009; Zhu and Sun 2009).

To partly solve some of the issues discussed above about distributed approaches, this paper proposes a methodology aiming at identifying the optimal size of the elementary area for simulations within a 140-hectare irrigation sector in Apulia Region, Southern Italy, so to minimize the number of soil measurements, as well as of simulations, without significantly affecting the simulation results at large scale (the irrigation sector scale, in our case). This was done by firstly simulating the sector scale irrigation volumes and deep percolation fluxes by including all the available detailed information on the actual spatial variability of soil HP and VP. Then, this information was gradually reduced, to find the minimum data set needed to still have sector scale comparable to those obtained by accounting for the detailed spatial variability. An agro-hydrological model called FLOWS (Coppola et al. 2019a, b) was used to simulate soil–water processes using the registered irrigation fluxes of the year 2020 as an input.

Materials and Methods

FLOWS agro-hydrological model

FLOWS is a dynamic physically based model to simulate water flow and solute transport in the soil–plant-atmosphere system (Coppola et al. 2019b). The numerical code is written in MATLAB and is based on a standard finite difference scheme. One-dimensional vertical transient water flow is simulated by numerically solving the 1D form of the Richards equation (RE):

$$C(h)\frac{\partial h}{\partial t}=\frac{\partial }{\partial \mathrm{z}}\left(K\left(h\right)\frac{\partial h}{\partial z}-K\left(\mathrm{h}\right)\right)-{S}_{w}(h)$$

(1)

using an implicit, backward, finite differences scheme with explicit linearization. In Eq. 1, C(h) = dθ/dh (cm⁻¹) is the soil water capacity, θ is the soil water content (cm/cm), h (cm) is the soil water pressure head, t (d) is time, z (cm) is the vertical coordinate being positive upward, K(h) (cm d⁻¹) is the hydraulic conductivity and S_w(h) (d⁻¹) is a sink term describing water uptake by plant roots and drainage.

FLOWS model provides information along the soil profile about the daily values of water contents, pressure heads, root uptake and actual evapotranspiration, water and solute fluxes.

For the processes related to the issues dealt with in the paper, the model requires the following inputs:

• 1. Soil hydraulic parameters: soil water retention, θ(h), and hydraulic conductivity parameters, K(h);

• 2. Vegetation parameters: Leaf-area index, LAI, and crop coefficient, K_c;

• 3. Reference evapotranspiration (ET_r);

• 4. Irrigation management parameters: irrigation periods, irrigation management depth (the depth at which irrigation fluxes should keep the pressure head above the critical threshold) and the critical pressure head threshold allowed in the irrigation management depth; and

• 5. Registered irrigation volumes.

Richards equation is solved numerically node by node in several small discrete thicknesses (Dz). In this paper, the solution was calculated for 100 Dz up to a maximum depth of 150 cm. Thus, each Dz is 1.5 cm in thickness.

The sink term, S_w, in Eq. (1) consists of two parts: the root water uptake per unit length of the root depth, S_r, and the artificial drainage, S_d. In this paper, the sink term represents only the root water uptake, which is calculated using a macroscopic approach (Feddes et al. 1976). The potential transpiration, T_p, is obtained from potential evapotranspiration, ET_p (which equals K_c ET_r), by partitioning ET_p into potential evaporation, E_p, and potential transpiration, T_p on the basis of Beer’s law (Ritchie 1972):

$$\begin{gathered} E_{p} = {\text{E}}T_{{\text{p}}} \times {\text{exp}}\left( { - k_{{e{\text{x}}}} {\text{LAI}}} \right) \hfill \\ T_{p} = ET_{p} - E_{p} \hfill \\ \end{gathered}$$

(2)

The potential root water uptake, S_r,p, is calculated by distributing the potential transpiration, T_p, over the root depth in proportion to the root density distribution:

$${S}_{r,\mathrm{p}}=g\left(z\right){T}_{p}$$

(3)

$$\int_{0}^{{{\text{Dr}}}} {g\left( z \right){\text{d}}z} = 1$$

(4)

where g(z) is the root-density distribution function (cm/cm³).

In general, the normalized g(z) may be obtained by normalizing the root length density distribution R_ld (cm cm⁻³) by its integral across the rooting depth, D_r:

$$g(z)=\frac{{R}_{l\mathrm{d}}\left(z\right)}{{\int }_{0}^{{\mathrm{D}}_{\mathrm{r}}}{R}_{l\mathrm{d}}\left(z\right)dz}$$

(5)

R_ld may be calculated as the ratio of the total length, L_r(z), of roots in a sample to the sample volume. g(z)dz is the fraction of roots located between z and z + dz.

Several root density distributions, g(z), may be selected in FLOWS model for simulating the sink term in Eq. (1), assuming root distributions to be either homogeneous (Feddes et al. 1978) or variable with depth (Raats 1974; Prasad 1988; Vrugt et al. 2001), the latter accounting for the fact that in a moist soil the roots can mainly extract water from the upper root zone layers.

FLOWS consider several root density distribution functions. A uniform distribution can simply be:

$$g\left(z\right)=\frac{1}{{D}_{r}}$$

(6)

where D_r is the maximum root length (cm).

A Prasad-type triangular distribution (Prasad 1988), with root water uptake at the bottom of the root zone, D_r, equal to zero is:

$$g\left(z\right)=\frac{2\left({D}_{\mathrm{r}}-\mathrm{z}\right)}{{D}_{\mathrm{r}}^{2}}$$

(7)

A Vrugt-type distribution (Vrugt et al. 2001) allows for calculating the dimensionless spatial root distribution, ω(z):

$$\omega \left(z\right)=\left(1-\frac{z}{{D}_{\mathrm{r}}}\right)\mathrm{exp}\left(\frac{{\mathrm{p}}_{\mathrm{z}}}{{\mathrm{D}}_{\mathrm{r}}}\left|{z}^{*}-z\right|\right)$$

(8)

where p_z (dimensionless) and z^* are empirical parameters. The model allows for non-symmetric root water uptake. The non-symmetry with soil depth is determined by the ratio of the p_z value for z ≤ z^* and the p_z value for z > z^*. To reduce the number of parameters to be optimized, p_z may be set to unity for z > z^*, whereas it is fitted for z ≤ z^*.

For Vrugt’s model, the g(z) distribution is:

$$g\left(z\right)=\frac{\omega \left(z\right)}{{\int }_{0}^{{\mathrm{D}}_{\mathrm{r}}}\omega \left(z\right)\mathrm{dz}}$$

(9)

A Logistic distribution gives the cumulative root density distribution:

$$\int_{0}^{{{\text{z}} = {\text{Dr}}}} {g\left( z \right){\text{dz}}} = \frac{1}{{\left[ {1 + a\;exp\left[ { - b\left( {\frac{z}{{D{\text{r}}}} - c} \right)} \right]} \right]^{{\frac{1}{a}}} }}$$

(10)

so that

$$g\left(z\right)=\frac{b\mathrm{exp}\left[b\left(c-\frac{\mathrm{z}}{{D}_{\mathrm{r}}}\right)\right]}{{D}_{\mathrm{r}}{\left\{1+\mathrm{aexp}\left[b\left(c-\frac{\mathrm{z}}{{D}_{\mathrm{r}}}\right)\right]\right\}}^{\frac{1}{a}+1}}$$

(11)

where a, b and c are the coefficients of the logistic function.

In this paper, due to the lack of information on the root density distribution in the study area, g(z) was considered a uniform distribution following Eq. (6) (Feddes et al. 1978).

Low water contents and/or the presence of soluble salts in the soil lower the total hydraulic head and may reduce the water fluxes to the roots, thus reducing root activity and water uptake. Reduction coefficients to decrease the maximum water uptake according to the water and osmotic stresses may be calculated independently and multiplied to calculate the actual root uptake, S_r, as:

$${S}_{\mathrm{r}}={\mathrm{\alpha }}_{\mathrm{w}}\left(h\right){\mathrm{\alpha }}_{\mathrm{s}}\left({\mathrm{h}}_{\mathrm{os}}\right){S}_{\mathrm{p}}={\mathrm{\alpha }}_{\mathrm{w}}\left(\mathrm{h}\right){\mathrm{\alpha }}_{\mathrm{s}}\left({\mathrm{h}}_{\mathrm{os}}\right)g\left(z\right){T}_{\mathrm{p}}$$

(12)

with α_w and α_s being reduction factors depending on the local (at a given z) water pressure head, h, and osmotic head, h_os, respectively. Accordingly, the actual transpiration, T_a, can be calculated as:

$$T_{a} = \int_{0}^{{D{\text{r}}}} {S_{{\text{r}}} \left( z \right){\text{dz}}}$$

(13)

In this paper, only water reduction due to water pressure head (α_w) was accounted for.

Irrigation criterion used in FLOWS

The model uses two configurations to simulate irrigation: (1) Irrigation provided by the user (IU), allowing to introduce irrigation volumes as an input, or (2) Irrigation computed by the model (IM), where the model calculates optimal irrigation fluxes based on a pressure head based criterium and accounting for the farmers’ behaviour.

The criterion used by the model to calculate the time for irrigation and the irrigation volume is schematically summarized in Fig. 1a, b. First, the model calculates a pressure head, h_av, averaged over the root depth, D_r.

Fig. 1

Graphical view of the criterion used by the model to calculate the time for irrigation and the irrigation volume (Coppola et al. 2019b). a h_av higher than h_crit, no irrigation is required; b h_av lower than h_crit, irrigation is required to bring the pressure head at the field capacity, h_fc

If h_av remains higher than a critical pressure head, h_crit, inducing some undesired stress to the crop (in terms of yield, product quality, etc.), this means that on average the pressure head h lies above the stress condition and no irrigation is required (Fig. 1a). Irrigation starts any time h_av becomes lower than h_crit (Fig. 1b). In other words, it is assumed that stress starts when the average pressure head in the root zone becomes lower than the threshold for water stress.

Irrigation aims to bring the actual water content at each depth to the water content at field capacity, h_fc. At the beginning of each simulation input time, the model calculates the gross irrigation height to be supplied on that day as the difference of water storage between the water content at the field capacity and the actual water content in the irrigation depth. The irrigation depth, hereafter z_c, may be set at any value even if a reasonable value should be the maximum root depth D_r. It corresponds to the depth over which the model computes h_av, to be compared to h_crit. Finally, the net irrigation height, q_irr,net, is calculated as the difference between the gross irrigation and the rain height eventually falling in the same day. It is also possible to select an on-farm irrigation efficiency, IE, of the irrigation system used for the crop considered, so that the actual irrigation amount, q_irr,act, may be obtained as (q_irr,net/IE).

Soil hydraulic properties models

In this paper, the van Genuchten (1980) soil water retention model was adopted in the case of unimodal porous system:

$$\begin{array}{*{20}l} {S_{\text{e}} {\mkern 1mu} = {\mkern 1mu} \frac{{\theta {\mkern 1mu} - {\mkern 1mu} \theta_{\text{r}} }}{{\theta_0 {\mkern 1mu} - {\mkern 1mu} \theta_{\text{r}} }}{\mkern 1mu} = {\mkern 1mu} \left[ {1{\mkern 1mu} + {\mkern 1mu} \left| {\alpha_{{\text{VG}}} {\mkern 1mu} h} \right|^n } \right]^{ - m} {\mkern 1mu} } & {h{\mkern 1mu} < {\mkern 1mu} 0} \\ {\theta {\mkern 1mu} = \theta_0 {\mkern 1mu} } & {h{\mkern 1mu} \ge 0} \\ \end{array}$$

(14)

where h is the pressure head, S_e is effective saturation and θ is the water content (θ₀ and θ_r are the water contents at h = 0 and for h → ∞, respectively). α_VG (cm⁻¹), n and m = 1–1/n are shape parameters. The effective saturation, S_e, may be considered as a cumulative distribution function of pore size with a density function f(h) which may be expressed by:

$$f\left( h \right) = \frac{{{\text{d}}S_e }}{{{\text{d}}h}}$$

(15)

In the case of bimodal porous medium, Durner (1994) water retention model was used:

$${S}_{e}=\sum {\upbeta }_{i}{\left[\frac{1}{1+{\left({\mathrm{\alpha }}_{\mathrm{VG},i}\mathrm{h}\right)}^{{\mathrm{n}}_{i}}}\right]}^{{\mathrm{m}}_{i}}0<{\upbeta }_{i}<1\text{ and }\sum {\upbeta }_{i}=1 \, i=\text{1,} \, {2}$$

(16)

in which β₁ and β₂ are the weighting of the total pore space fraction to be attributed, respectively, to inter-aggregate pores (the macropores) and to the intra-aggregate pores (the micropores or matrix), and α_VG,i, n_i and m_i still represent the shape parameters for each of the two partial curves.

The relative unsaturated hydraulic conductivity, K_r, was described using Mualem model (Mualem 1976). In case of unimodal van Genuchten water retention, K_r becomes:

$${\text{K}}_{\text{r}} \left( {S_e } \right){ = }\,\frac{{{\text{K}}\left( {S_e } \right)}}{{{\text{K}}_{0} }}S_{\text{e}}^{\uptau } \left[ {{1} - \left( {{1} - S_e^{{\raise0.7ex\hbox{${1}$} \!\mathord{\left/ {\vphantom {{1} {\text{m}}}}\right.\kern-0pt}\!\lower0.7ex\hbox{${\text{m}}$}}} } \right)^{\text{m}} } \right]^{2}$$

(17)

where K₀ (cm/d) is the hydraulic conductivity at h = 0.

In case of bimodal Durner’s water retention, K_r becomes:

$${\mathrm{K}}_{\mathrm{r}}\left({S}_{e}\right)=\frac{\mathrm{K}\left({S}_{e}\right)}{{\mathrm{K}}_{0}}={\left({\sum }_{i=1}^{2}{\upbeta }_{i}{\mathrm{S}}_{\mathrm{e},i}\right)}^{\uptau }{\left\{\frac{{\sum }_{i=1}^{2}{\upbeta }_{i}{\mathrm{\alpha }}_{i}\left[1-{\left(1-{\mathrm{S}}_{e,i}^{1/{m}_{i}}\right)}^{{\mathrm{m}}_{i}}\right]}{{\sum }_{i=1}^{2}{\beta }_{i}{\alpha }_{i}}\right\}}^{2}$$

(18)

In this work, the surface layer of the soil profiles was found to be better described by bimodal hydraulic properties while the subsurface layers were found to have unimodal hydraulic properties.

Study area

The study area is the irrigation district 10 in “Sinistra Ofanto” irrigation system, located at the left bank of “Ofanto” river in Apulia region, Southern Italy. The district is about 2000 ha divided in 19 sectors. The overall objective of the study was to apply the agro-hydrological model to all the 19 sectors to optimize the irrigation scheduling in the whole district. To do that, the sector 6 was identified as a sort of calibration area to analyse preliminary the role of spatial variability of soil HP and VP on the irrigation needs and deep percolation losses, as well as the farmer behaviour related to irrigation, which are essential to constrain the simulations to the real context and to establish the simulation grid size and thus minimum number of soil sites to be characterized in an irrigation sector. Only then, the model was applied to the 18 remaining sectors. The sector 6 is 140 ha, one of the larger sectors, and is divided into 39 fields each provided with an irrigation water delivery hydrant to feed a drip irrigation system. Based on a preliminary pedo-morphological investigation, the sector was found to consist of five soil pedological units (hereafter, PUs).

Ninety well distributed sites were selected to measure HP in the first soil horizon. HP in deeper horizons were measured at 5 sites using 1.5 m depth excavation pits. Those five pits were the same used to classify the soil pedological units in the study area. Based on the 5 excavation pits, the sector 6 was firstly divided into five Thiessen polygons, each assumed to approximately correspond to a PU, by a Geographic Information System (ArcGIS). Thus, these five polygons were further subdivided into 90 Thiessen polygons based on the 90 surface HP measurement sites. HP for deeper horizons were assigned to each of the 90 Thiessen polygon according to the nearest of the five excavation pits. Figure 2 shows the study area, with: (a) the irrigation district with all the 19 sectors and (b) the sector 6: the solid border represents sector 6, the dashed borders represent Thiessen polygons of surface measurement locations, the solid circles represent the 90 surface measurement sites, the red triangles represent the excavation pits and the red borders represent the Thiessen polygons of the 5 PUs.

Fig. 2

The study area. a the irrigation district with all the 19 sectors (labelled by sector numbers) and the excavation pits used to to characterize the pedological soil units; b the sector 6. The solid circles and the red triangles represent, respectively, the ninety surface measurement sites and the five excavation pits. The red solid and the black dashed borders represent, respectively, the five Thiessen polygons based on the five pits and the ninety Thiessen polygons based on the surface measurement locations

Measurement of soil physical and hydraulic properties

Each of the 5 excavation pits was characterized from a pedological point of view. Both disturbed and undisturbed soil samples were collected from each of the pit pedological layer. Samples were also taken at each of the 90 surface measurement sites. All the samples were analyzed to measure their physical parameters in the laboratory, including particle size distribution, aggregate size distribution and soil bulk density, ρ_b. At each of the 90 surface measurement sites, HP were determined using tension infiltrometer (Ankeny et al. 1988; Coppola et al. 2011; Wang et al. 2013). The soil surface was first prepared and levelled, a ring was placed on the surface and then fossil sand was added within the ring to ensure good contact with the infiltrometer’s disc. The Mariotte tower was used to control the pressure head at each site, which was set at 4 pressure head values (– 15, – 10, – 5 and – 1 cm) by progressively raising the bubbling point in the tower. The reservoir’s infiltration depths were measured against time. Soil samples were taken before and after the infiltration experiments to measure the initial and final water contents.

HYDRUS 3D model was used to carry out inverse solutions to Richards 3D infiltration equation using the cumulative infiltration data as inputs (Rassam et al. 2003; Šimůnek et al. 2008; Hassan et al. 2022;). Both the unimodal and the bimodal HP models were used to fit the simulations to the data. The objective function of the inverse solution included 3 parameters for the van Genuchten–Mualem model: α_VG, n, and K₀, and 5 parameters for the Durner–Mualem model: α_VG,1, α_VG,2, n₁, n₂ and K₀. θ₀ was assumed as the total porosity while θ_rwas assumed to be zero, as suggested by most of the fittings of the infiltration curves which resulted in very low or null values. In the bimodal model, the weight β for the macropores was based on both the particle and aggregate size distribution and calculated according to the procedure given in Hassan et al. (2022). The same measurement procedure was applied to the first three layers in the 5 excavation pits. Based on the fitting, the surface HP were better described by a bimodal model (see Hassan et al. 2022), whereas the HP of the deeper layers in the pits were well described by the unimodal van Genuchten–Mualem model.

Obtaining vegetation parameters using GIS and remote sensing applications

Figure 3a shows the cropping pattern of the sector 6 study area. It is relatively uniform, with 58% of the area cultivated with vineyard (wine grape) while the second predominant crop is peach, which represents 19% of the area. As for the root depths and irrigation periods, ArcGIS was used to create an intersection crops, which are the dominant ones in the sector, was considered constant at 80 cm. between the cropping pattern and the Thiessen polygons. The root depth of the tree As for cereals, which occupy a minor part of the sector, the root depth was considered increasing over time based on literature values. If a crop type occupied more than 75% of the Thiessen polygon area, the polygon was considered homogeneous and was assumed to be cultivated only with that crop (Fig. 3b–-i). If else (i.e., no crop > 75%), a fictitious crop was created with a weighted average of the root depth, LAI and Kc and the longer-spanning irrigation period (Fig. 3b–-ii). In any case, it is worth to note here that, in most cases, the conditions of Fig. 3b–-i applied. Besides, one should be aware that this is the condition frequently found in an irrigation sector, which is generally cultivated with a prevailing crop.

Fig. 3

a The cropping pattern of the study area, with the 90 Thiessen polygons. Note that grape and peach are the predominant crops; b Obtaining cropping patterns of Thiessen polygons: (i) an example of a polygon with a predominant crop (occupying more than 75% of the area); (ii) a heterogeneously cultivated Thiessen polygon, with a fictitious crop built by the weighted average of the vegetation parameters of the real crops within the Thiessen polygon

European Space Agency’s (ESA) Sentinel-2 satellite images with 10 m resolutions were utilized to obtain LAI and K_c. Two satellite images were downloaded for each month of the year 2020. Each satellite image provides the reflectance values of 13 different spectral bands. SNAP software, developed by ESA, was first used to resample and subset the satellite images to obtain the reflectance values for different bands. ArcGIS was then used to obtain the normalized difference vegetation index (NDVI) using raster calculator:

$$\mathrm{NDVI}=\frac{\mathrm{NIR}-\mathrm{Red}}{\mathrm{NIR}+\mathrm{Red}}$$

(19)

where NIR and Red are the near-infrared and the red reflectance, respectively.

K_c was then obtained at each pixel using the following formula (Beeri et al. 2020):

$${K}_{\mathrm{c}}=1.1875*\mathrm{NDVI}+0.04$$

(20)

LAI was calculated using the SNAP software’s biophysical processor, which utilises artificial neural networks (ANN) to obtain LAI at each pixel. Zonal statistics tool on ArcGIS was then used to obtain the mean values of K_c and LAI for each Thiessen polygon for each image representing a day of the year. Time series of the evolution of K_c and LAI, along the year 2020, were created for each polygon using linear interpolation. Time series were first checked graphically to remove the outliers. Then, the time series curves were smoothed by calculating 5-day moving averages to reduce the oscillations of the original values.

Agro-hydrological simulations in sector 6

Table 1 provides a list of abbreviations and symbols used in this section.

Table 1 List of abbreviations used in the agrohydrological simulations

Based on the database described above, FLOWS model was used to optimize the irrigation volumes in the study area. Two main scenarios were considered: (1) A reference scenario (DVS), where distributed irrigation volumes and related deep percolation volumes at a depth of 80 cm were calculated on the ninety polygons by accounting for the actual observed variability of soil HP and VP inputs. The irrigation volumes and percolation volumes in the whole sector 6 were obtained by aggregating the distributed outputs obtained in the single polygons; (2) Then, reduced variability scenarios (RVS) were considered where the information on the actual spatial variability of the soil HP and VP was gradually reduced to find the minimum data set needed to have sector scale irrigation volumes and deep percolation volumes comparable to those obtained under the DVS.

In the DVS case, simulations were carried out under two different sub-scenarios:

• (1) First, simulations were carried out under the Irrigation provided by the user (IU) configuration using the registered irrigation volumes as an input. The purpose of these simulations was to understand the farmers’ irrigation management parameters: firstly, the registered volumes allowed to obtain the average time intervals when the farmers are used to irrigate each crop. Also, the output of the simulations using the registered volumes was analyzed graphically by plotting the simulated mean pressure head at different depths to obtain the depth at which the volumes actually provided by the farmers allow keeping the average pressure head above the critical pressure head, h_crit, for the crop. This average depth was then assumed as the irrigation depth, z_c, to be used in all the simulation scenarios to retain the real farmer behaviour in the simulations, making them as much as possible realistic. The critical pressure head was established at h_crit = – 500 cm, which is frequently adopted for both grapes and peach, the prevailing crops in the sector 6. It is an average value between those used for grapes by Phogat et al. (2018) (h_crit =– 600 cm) and Coppola et al. (2019b) (h_crit = – 400 cm). A total of 90 simulations was carried out using this configuration, each representing one of the Thiessen polygons in the whole sector.

• (2) Then, simulations were repeated on the same 90 polygons to optimize the irrigation volumes under the Irrigation computed by the model (IM) configuration, using the irrigation intervals and the z_c deduced from the IU simulation as input. Now, in the IM configuration, the irrigation volumes were not provided as an input but computed by the model according to the irrigation criterion described above in the Sect. “Irrigation criterion used in FLOWS” The results of the two sets of simulations were aggregated at sector scale and then compared in terms of total sector-scale irrigation volumes and related deep percolation volumes during the whole growth season.

In the RVS case, simulations were carried out only under the IM configuration, using sparser data and thus decreasing the number of polygons into which to subdivide the sector. The objective was looking for decreasing the number of simulations in the sector but preserving estimates of sector scale irrigation volumes and deep percolation volumes comparable to those obtained in the DVS scenario under the IM configuration.

Under the RVS configuration, the information on the variability of VP and HP was progressively overlooked according to the following two steps:

• (A) First, the sensitivity of model simulations to VP variability was evaluated by keeping the ninety soil HP and assuming a fictitious crop with unique VP (UVP) for the whole sector. These UVP were obtained as the mean daily values of K_c, LAI and irrigation duration weighted over all the different crop areas in the sector (according to the procedure illustrated in Fig. 3b). This simulation configuration was identified as RVS-UVP and still involved ninety simulations in the sector, one for each soil HP profile, but always with the same VP of the fictitious crop. Again, the sector scale results were obtained by aggregating the outputs in the five polygons. This RVS-UVP configuration allowed quantifying how accounting for the actual vegetation distribution or taking a sort of averaged crop on the whole sector affected the simulation results. Using the fictitious crop approximation allowed to greatly simplify the analyses carried out in the following steps.

• (B) Then, an uncertainty analysis was carried out to further overlook the variability of the input properties by reducing the number of the surface HP inputs to carry out the simulations in IM configurations, anytime comparing the average results to those obtained under the DVS configuration and evaluating the additional estimation uncertainty introduced by progressively ignoring the actual variability. This analysis was implemented according to the following steps:

• An increasing number of HP inputs, N_s, was selected, from 5 to 40 with a step increase of five (N_s = 5, 10, 15, 20, 40.

• Ten combinations of N_s HP inputs were generated from the 90 soil data points in Fig. 1. These combinations were not completely random but were generated so that they were equally distributed among the five PUs (the red borders in Fig. 1). For example, for N_s = 30, 10 combinations of soil sites were generated, each consisting of 6 sites for each of the five PUs; for N_s = 20, 10 combinations were considered with 4 soil sites for each of the five PUs, and so on. This quasi-random sample generation was defined here as controlled random generation, CRG);

• For each of the ten combinations, Thiessen polygons corresponding to the N_s points were delineated and their areas were calculated.

• The optimal irrigation volumes for each of the N_s polygons was obtained as those calculated according to the RVS-UVP scenario. The overall average irrigation volumes of sector 6 for each combination of N_s polygons was obtained as the average of the irrigation volumes of those polygons weighted for their areas. Also, the uncertainty associated to each N_s was evaluated as standard deviation of the volumes coming from the ten groups of simulation in a given N_s value. The reliability of the average estimation for each N_s scenario was evaluated with reference to the average value obtained from the original 90-site simulations (DVS-IM). The average estimation was considered reliable when it was between 0.85 and 1.15 of the DVS-IM volumes. The acceptability of the uncertainty coming from each N_s scenario was evaluated with reference to the case with the largest N_s (N_s = 40). Specifically, the uncertainty was considered acceptable any times it was not larger than 1.25 the uncertainty coming from the N_s = 40 scenario.

• As a final step, five independent simulations were carried out by considering only one soil profile for the whole sector by taking one of the five soil pits HP for each simulation. This step was identified as RVS-UHP where the simulations were carried out using unique HP input for the entire sector.

The authors are aware that using the fictitious crop introduces some baseline uncertainty which would not be present when considering the actual vegetation polygon by polygon. However, one should consider that accounting for the actual vegetation variability would require recalculating each time the VP for each polygon, whose size and shape is changing according to the configuration considered. Just as an example, in the case of N_s = 40, which is simulated 10 times, this would mean having 400 different polygons, each potentially including a different combination of crops and thus with different VP. However, assuming a single fictitious crop at sector scale makes sense as, in general, a quite homogeneous crop distribution is expected at this scale. This is what we actually observed in the whole district 10 under investigation and this is the general situation in the Mediterranean countries, where an irrigation sector is frequently characterized by similar land use and irrigation management practices.

Figure 4 shows an example of the uncertainty analysis with N_s = 5, where 10 different combinations of 5 profiles were generated so that each PU contained one soil profile. The dashed lines represent the Thiessen polygons corresponding to each combination and the red lines represent the five PUs in sector 6.

Fig. 4

An example of the uncertainty analysis carried out to study the effect of the HP spatial variability on the simulations. In this example the number of HP points is 5 and different combinations of 5 profiles were generated so that each PU contained one soil profile. The dashed lines represent the Thiessen polygons corresponding to each combination and the red lines represent the five PUs in sector 6

Figure 5 provides a schematic view of the steps followed to carry out agro-hydrological simulations under the DVS and RVS configurations. Note that the goodness of estimations obtained under the RVS was evaluated by comparing the corresponding average sector volumes to those obtained under the DVS, as well as by quantifying the estimation uncertainty introduced using progressively sparser data. This will be even clearer in the Results section.

Fig. 5

Schematic view of the steps followed to carry out agro-hydrological simulations under the DVS (accounting for the spatial variability of soil and vegetation properties) and RVS (using sparser soil and vegetation data) configurations

Results and discussions

Soil properties in the 90 polygons and in the five soil pits

Table 2 reports the textural data and USDA soil textural classification of the first and second soil horizons of the five pedological soil units (PUs) in sector 6 obtained from the five excavation pits. The table also reports the area of each PU and the number of surface measurement sites within each PU. Table 3 reports the average values and the coefficients of variation for the bimodal Durner–Mualem soil hydraulic parameters in the first layer and the unimodal van Genuchten–Mualem soil hydraulic parameters in the second layer. The table indicates that the lower soil layer is generally more conductive than the upper layer. It also shows that the soil HP variability in the upper layer is higher than that of the lower layer, even if it has to be considered that the variability of the second layer was calculated on only five profiles.

Table 2 Soil textural data and USDA soil texture classification of the first and second horizons of the five pedological soil units in sector 6, as well as, the area of each soil unit and the number of measurement sites within each PU

Table 3 Average values and standard deviations of the bimodal Durner–Mualem soil hydraulic parameters in the first layer and the unimodal van Genuchten–Mualem soil hydraulic parameters in the second layer

Identifying the farmer’s behaviour to be used for the different simulation scenarios

Irrigation periods per crop

Registered irrigation data were used first to identify the irrigation period of each crop during the year. Figure 6 shows the daily irrigation fluxes for ten Thiessen polygons predominantly cropped with wine grapes during 2020. The figure shows that, though with some variability, wine grapes are mostly irrigated between the days 165 and 260 of the year in the study area. This graphical analysis was repeated for each crop type in the sector.

Fig. 6

The registered daily irrigation heights in 2020 for 10 of the Thiessen polygons in sector 6 that were cropped with wine grapes. The labels indicate the Thiessen polygons names

Soil water pressure head along the soil profile induced by the farmer behaviour

The registered irrigation volumes were used as an input to FLOWS simulations using the DVS-IU configuration. As explained, these simulations were run mostly to understand the farmer behaviour, which has to be accounted for realistic modelling. Figure 7 shows a representative example of the output for a Thiessen polygon (S27) that was cropped with grape. The figure shows the inflow by precipitation (solid blue line) and irrigation (solid light blue lines), along with the mean soil water pressure head in the 80-cm root depth (solid black line) and in the first soil horizon (dashed black line), which in the S27 polygon is 46 cm in depth. The solid red line in the figure represents the deep percolation fluxes and the dashed red line represents the soil water fluxes at 47 cm depth below the first soil horizon. The figure indicates that, in the irrigation period (approximately from day 165 to day 260), the farmer tends to keep the water pressure head above the critical value, h_crit, only in the first layer, whereas the average pressure head in the 80 cm root depth may fall well below the h_crit value. This result may be explained by the relatively high-water retention of the first layer, which is able to keep most of the water supplied by the irrigation applications. This high-water retention takes the pressure head in the first layer regularly well above the critical pressure head of – 500 cm, frequently very close to saturation. By contrast, due to the root uptake, the pressure head in the second layer tends to decrease progressively to values far below the critical value, which induces some upward fluxes from the deeper soil not explored by roots and thus with no root uptake (see the positive deep percolation fluxes during the irrigation period). This was evident by looking at Fig. 7b; the water fluxes at 47 cm depth (i.e., from the first to the second soil horizon) during the irrigation season were mostly null or insignificant quite low (2.1 mm/day at most).

Fig. 7

A representative example of the output for a Thiessen polygon (S27) that was cropped with grape. The upper part (a) shows the inflow by precipitation (solid blue line) and irrigation (solid light blue lines), along with the mean soil water pressure head in the 80-cm root depth (solid black line) and in the first soil horizon (dashed black line). The lower part (b) shows the deep percolation fluxes as a solid red line and the soil water fluxes at 47 cm depth as a dashed red line

The farmer probably ignores this actual soil hydraulic behaviour and may perceive the soil root zone well wetted everywhere. In this sense, realistic modelling has to account for this particular behaviour. In fact, a simulation looking for optimising irrigation water volumes, not considering the specific farmer behaviour and looking only for maintaining the pressure head above the h_crit along the whole root zone, would necessarily result in calculated irrigation volumes significantly higher that the registered ones. The results in the next section will show how an appropriate modelling including this consolidated farmer behaviour may still result in comparable root uptake (and thus transpiration) with lower water volumes.

Optimising irrigation volumes for spatially variable HP and VP

Figure 8a shows, as colour-graduated maps, a comparison between the registered irrigation volumes (on the left) and the optimized irrigation volumes obtained using FLOWS model (on the right) for each of the Thiessen polygons for the year 2020. All the volumes are so-called specific volumes in (m³/ha). Figure 8b provides the average values at the sector scale. The mean optimal irrigation specific volumes were lower than the registered ones with a mean value of 1360 m³/ha compared to the registered 1693 m³/ha for the whole irrigation season of the same year. The deep percolation specific volumes at a depth of 80 cm were also calculated for the same irrigation period using FLOWS model for both the registered irrigations and those optimized by the model, with very similar values for the irrigation season (43 and 41 m³/ha, respectively). In both cases, they represent fluxes directed upward. Thus, during the irrigation season, FLOWS predicts water dynamics from the deeper and wetter soil layers to the 0–80 cm root zone. As may be seen in Fig. 9, the root uptake is similar for both the registered and optimized cases. This means that the excess water supplied by the farmers remains stored in the soil profile, which is systematically wetter than the same soil profile under the optimized volumes case (graphs are not shown here for the sake of brevity).

Fig. 8

a: Registered (left) and optimized irrigation specific volumes obtained using FLOWS (right) for the ninety Thiessen polygons under the DVS-IU and DVS-IM respectively, and b: average registered and optimized irrigation fluxes and their corresponding deep percolation fluxes at 80 cm depth for sector 6. All volumes are specific volumes in m³/ha

Fig. 9

The cumulative root water uptake in the first soil layer (black dashed curve) (depth from 0 to 46 cm), the cumulative root water uptake in the second soil layer (black solid curve) (depth from 46 to 80 cm), and the cumulative irrigation depths (blue solid curve) in the Thiessen polygon S27 for two cases: using registered irrigation volumes (left) and using irrigation fluxes optimized by FLOWS model (right), both under the DVS, thus using spatially variable VP and HP

Figure 9 shows the cumulative root water uptake in the first soil layer (black solid curve), the cumulative root water uptake in the second soil layer (black dashed curve) and the cumulative irrigation depths (blue solid curve) in the Thiessen polygon S27 for two cases: using registered irrigation volumes (on the left) and using irrigation fluxes optimized by FLOWS model (on the right), both under the DVS, thus using spatially variable VP and HP. The root water uptake for both cases was calculated by FLOWS model.

Figure 9 indicates that the total root water uptake values along the root zone in both scenarios DVS-IM and DVS-IU are quite close (18.22 and 18.66 cm, respectively) in spite of the total optimized irrigation depth much smaller than the registered one (11.82 and 16.20 cm, respectively). The total root water uptake depth in the second soil layer is 6.62 cm when using registered irrigation volumes while it is slightly smaller (6.41 cm) when using optimized irrigation volumes.

By comparing Figs. 8 and 9, it may be seen that, even if farmers supply an average of 300 m3/ha more, both farmers and optimized irrigation volumes produce practically the same root uptake and similar deep percolation fluxes. This may be explained by looking at the graphs in the Fig. 10a and b. Figure 10a shows the water storage calculated by FLOWS in the first 80 cm of the soil profile (the rooted depth) for both DVS-IU and DVS-IM scenarios. In the first part of the graph, the two storages are perfectly overlapping because the water content evolution is simply induced by rainfall, which is the same for both cases. By contrast, in the period between days 165 and 260 (the irrigation period), the storage coming from farmers’ irrigation volumes is systematically higher than that coming from optimized volumes. Looking at Fig. 10b, showing the average soil water pressure head, h_av, in the first soil layer, it may be observed that the h_av induced by farmer irrigation had regularly high values (close to saturation), whereas those coming from optimized irrigation had much lower values, but always above the critical value, h_crit. This suggests that the farmer tends to keep regularly the first layer at high water content values, which is not really necessary, while the model optimisation aims at keeping h_av just above the h_crit value. It is interesting to note that the increase of pressure head around day 186 in the optimized case (Fig. 10b) corresponds to only a small change in soil water storage (Fig. 10a). This may be explained by the high non-linearity of the water retention curve, which is amplified by the bimodal behaviour in the involved range of pressure heads.

Fig. 10

Time evolution of: a water storage and b average pressure head in the first soil layer for Thiessen polygon S27 for two cases: using registered irrigation volumes (blue lines) and using irrigation fluxes optimized by FLOWS model (orange lines), both under the DVS, thus using spatially variable VP and HP. The difference in water storage and the critical pressure head are also shown in the graphs as black dashed lines. The red dashed lines represent the irrigation period between the days 165 and 260

Overlooking the spatial variability of VP and HP

Figure 11a shows, as colour-graduated maps, a comparison between the optimized irrigation volumes obtained using FLOWS model under the DVS configuration on the left and the optimized volumes under the RVS-UVP configuration (that is, by overlooking the VP variability and using a unique VP input for all the 90 Thiessen polygons) on the right. Figure 11b shows the sector-scale average irrigation and deep percolation specific volumes at 80 cm depth for the two configurations DVS and RVS-UVP in m³/ha. The figure indicates that coalescing the variability of the VP into a unique fictitious VP input for the 140 ha sector did not significantly affect the FLOWS model results. This finding allowed to greatly simplify the analyses carried out in the following steps, where the effect of neglecting the HP spatial variability was evaluated by considering always the same fictitious crop for all the polygons simulated.

Fig. 11

a Optimized irrigation specific volumes obtained using FLOWS for ninety Thiessen polygons under the DVS-IM configuration (left) and under RVS-UVP configuration (right); b the average optimized irrigation specific volumes and their corresponding deep percolation specific volumes at 80 cm depth for sector 6 for the DVS and RVS-UVP configurations, respectively. All specific volumes are in m³/ha

Further reduction in the input variability was considered by analysing the irrigation fluxes results coming from using progressively sparser HP inputs. To do that, the procedure of considering different number of HP input profiles, N_s, was adopted, as described in the Sect. “Agro-hydrological simulations in sector6”. Figure 12 shows the box plots with the weighted averages, the minimum, the maximum, the median, along with first and third quartiles of irrigation specific volumes, corresponding to each selected number of HP input profiles, N_s.

Fig. 12

Box plots with the weighted averages, the minimum, the maximum, the median, along with first and third quartiles of irrigation specific volumes, corresponding to each selected number of HP input profiles, Ns. The additional horizontal axis labels show the average area assigned to each HP profile. The box plot for the 90 sites refers to the deterministic scenario involving simulations carried out using spatially variable soil and vegetation properties (DVS-IM)

Figure 12 is divided into three parts according to the average prediction (the cross symbol) and uncertainty levels (indicated by the Inter-Quartile Range—IQR associated with changing the number of HP inputs: (1) the values in the blue rectangle represent good average predictions and relatively low uncertainty, meaning that HP in 20–40 sites provide more or less similar average irrigation volumes to the DVS configuration and that this is verified practically in all the randomly generated profile combinations. This would support the choice of starting this analysis with no more than 40 sites, which would have been uselessly redundant. Besides, it would be quite impractical having more than 40 real hydro-pedological characterisations in a single irrigation sector. Overall, the distribution of the values more or less symmetrical; (2) the values in the red rectangle still provide good average and median predictions but with additional uncertainty, which comes from using a limited number of HP. In any case, it should be noted that five soil profiles produce approximately the same results coming from using 15 sites, even if the distribution becomes more asymmetrical; and (3) the single value corresponding to simulations using only one soil profile for the entire sector, which produces unacceptable average prediction along with an extremely high uncertainty level. It should be noted that the variability for the 90-site case comes from the real heterogeneity of the polygons within the sector 6, whereas the variability for all the other scenarios (N_s <  = 40) comes from the variability of the mean among the different stochastic configurations within a given N_s. In this sense, a direct comparison of the IQR would be conceptually inappropriate.

In general, the uncertainty resulting from reducing the number of HP inputs stems from overlooking the variability in HP within each PU and/or between different PUs. The low uncertainty levels at N_s ≥ 20 (i.e., N_s ≥ 4 HP inputs per PU) was due to the higher number of HP inputs which was able to capture the HP variability within each PU as well as between different PUs. On the other hand, the moderate uncertainty zone at 5 ≤ N_s ≤ 15 (i.e., at least one HP input per PU) was able to capture the variability of HP between different PUs and only partly that within PUs.

In the remaining of this analysis, using one soil profile for each PU (5 profiles for sector 6 or about 1 profile for about 28 ha) was considered a good compromise between accurate description of the site and costs/feasibility of the samplings, as it allows for drastically reducing the simulation time, with a relatively small increase of uncertainty, comparable to that obtained using 3 profiles per PU (15 profiles for sector 6). Besides, this is also a realistic option, as making more than one soil measurement per PU is normally considered unfeasible.

These results were extended to the remaining 18 sectors of district 10, which were hydraulically characterized by taking one HP profile measurements per PU identified in each sector. Just for completeness of the discussion, Fig. 13 shows the simulation results of the whole district 10 based on this concept.

Fig. 13

Registered (left) and optimized irrigation fluxes (in m.³/ha) obtained using FLOWS (right) for the forty-one Thiessen polygons representing the soil units in all district 10

Conclusions

The main purpose of this study was to investigate the effect of soil and vegetation properties’ spatial variability on optimising irrigation volumes at large scales. Firstly, aggregating VP spatial variability into a single fictitious crop with the average VP parameters did not affect the outputs significantly. This has important implications for irrigation management at large scale by agrohydrological modelling, as this allows focusing only on the spatial variability of soil HP. However, it is worth to note that this result may well be explained by the fact that the model used in this paper is mostly soil-controlled and is, thus, more sensitive to soil HP than VP. Also, the cropping pattern in the study area was relatively homogeneous with the majority of the area cultivated with a single crop type, i.e., wine grapes. In any case, this relatively homogeneous cropping pattern is generally found in irrigation sectors.

Reducing the spatial variability of soil properties by overlooking the soil hydraulic properties (HP) proved to affect the model outputs, especially in terms of uncertainty related to the estimation of optimal irrigation fluxes. This uncertainty partly comes from neglecting the variability of the soil HP within each soil pedological unit (PU) and partly that between different PUs. The results showed that the uncertainty levels are relatively low when accounting for several HP input datasets within each pedological soil unit (PU). Accounting for at least one soil HP input per PU leads to moderate uncertainty levels but still with good average predictions. The simulation results become completely unacceptable when considering only one soil profile for the whole sector.

The results also emphasized the role of farmers’ behaviour in irrigation management. Agrohydrological models should be used first to analyse the farmers’ methods to manage irrigation scheduling, and at least determine: the irrigation period per crop, soil–water pressure head threshold (critical pressure head, h_crit) and the critical depth, zc, at which the farmers maintain h_crit. These parameters are important for optimising irrigation and deep percolation fluxes, just like HP variability. Overall, it is possible to account for one HP and one VP within each PU. Therefore, large sectors and/or districts can be divided into polygons representing the pedological soil units, each with one macrocrop and one soil profile to optimize their irrigation fluxes. In our case, doing that allowed for reducing the number of simulations from 90 to 5 in sector 6, which limited the number of simulations in the whole district by the order of 100.