Introduction

Agriculture is the largest consumer of freshwater and accounts for 70% of current human water use (FAO 2011). The Food and Agriculture Organization of the United Nations (FAO) predicts that if current consumption patterns continue at the present rate, two-thirds of the world’s population could be living in water-stressed countries by 2025. As a result, water will become scarce not only in arid areas but also in regions where precipitation is abundant (Pereira et al. 2002). Within this context, optimal irrigation water strategies are crucial for saving water while guaranteeing maximum crop yields in a near future.

Irrigation scheduling is the process used by irrigation system managers to determine both the correct moment and the required amount of water to irrigate fields while maximizing crop yields with the minimum amount of water applied. Properly defining the moment of irrigation through water content (or pressure head) thresholds has been demonstrated to substantially increase the efficiency of agriculture fields (Dabach et al. 2013). However, irrigation water contains salts and fertilizers that often promote soil salinization, which results in an increase in soil electrical conductivity (\(\mathrm{EC}\)) and a reduction of crop productivity (Machado and Serralheiro 2017). Thus, optimal irrigation strategy should be designed to avoid soil salinization (Pereira et al. 2007).

Because of the importance of irrigation scheduling, a wide range of irrigation methodologies have been recently developed to guarantee crop productivity as well as to avoid soil salinization. We distinguish between methods based on crop water requirements, direct measurements of the plant water status or response to water stress, direct measurements of water content, and numerical modeling of flow and transport of salt concentrations in the vadose zone.

Water requirements determined from crop evapotranspiration calculations (\({\mathrm{ET}}_{c}\)) is the most widely used method for irrigation scheduling (Feng et al. 2007; Orgaz et al. 2006; Salim et al. 1970). \({\mathrm{ET}}_{c}\) is typically estimated as the product of two terms: the reference evapotranspiration (\({\mathrm{ET}}_{0}\)), generally defined for either clipped grass (Allen et al. 1998) or alfalfa (Wright and Jensen 1978), and the specific crop coefficient (\({K}_{c}\)), estimated from different tabulated information such as FAO56 Irrigation and Drainage Paper No. 56 (Allen et al. 1998). In this case, the daily or weekly volume of water evapotranspired is estimated and used to schedule the irrigation in the next days or weeks (Pereira et al. 2009; Sun et al. 2006; Thompson et al. 2007). If water requirements are not combined with soil moisture sensors measurements, this volume of water is arbitrarily applied in the field.

Plant-based methods include direct measurements of the plant water status as well as a number of plant processes that are known to depend on water deficits (Jones 2004). Different types of measurements can be used to determine plant water and salt stress for irrigation scheduling. This includes photosynthesis capacity (Flexas et al. 2004; de Lima et al. 2015; Ribas-carbo et al. 2006), stomatal conductance (Flexas et al. 2004; Jones 1999), leaf water potential (Alberola et al. 2008; Girona et al. 2006; Turner 1990), and crop temperature (Bellvert et al. 2014; DeJonge et al. 2015; Kassie et al. 2018). In general, these methods only provide a direct measurement of the plant status without determining how much water is necessary. Moreover, these methods require sophisticated and calibrated devices (Jones 2004).

Several algorithms for optimizing irrigation can be found in the literature. For instance, Soentoro et al. (2018) optimized irrigation by determining the cropping patterns and planting areas through linear programming. Ortega et al. (2004) and Martínez-Romero et al. (2017) determined irrigation by maximizing the gross margin through genetic algorithms. Noory et al. (2011) proposed a linear and a mixed-integer linear model for optimizing irrigation water allocation and multicrop planning that maximizes the total net benefit. These irrigation scheduling methods are based only on water requirements and do not properly represent the water movement in the vadose zone. Several authors have addressed the problem and made efforts to derive thresholds of irrigation for different soil types with different soil hydraulic properties based on the unsaturated hydraulic conductivity and potential plant water uptake for different root distributions (Jabro et al. 2020; Martínez-Gimeno et al. 2020; Srivastava and Yeh 1991). Even though their work remained theoretical, several authors (Collin et al. 2019; Gendron et al. 2018; Létourneau et al. 2015; Létourneau and Caron 2019; Rekika et al. 2014) have later on derived the theoretical relationship and tested it in the field with success. Campbell (1982) described at which soil water content fringe the crop is under optimal conditions. This soil water content fringe depends on soil hydraulic properties and it is shown to strongly vary in space (Feki et al. 2018). This indicates that it is important to consider soil water movement through the root zone during irrigation scheduling.

Numerical models constitute an efficient tool for assessing irrigation scheduling (Linker et al. 2016; Ma et al. 2015). Among the different models available, HYDRUS (Šimůnek et al. 2008, 2016) is often used to simulate water fluxes, root water uptake, root growth, and solute and heat transport in the vadose zone. Several researchers have used HYDRUS for simulating water content, soil suction and \(\mathrm{EC}\) to improve irrigation scheduling or to provide some information to stakeholders. For instance, Arbat et al. (2008) simulated soil suction with HYDRUS to demonstrate that this model is capable of assessing irrigation scheduling in a given field site. Siyal and Skaggs (2009) and Skaggs et al. (2010) simulated soil moisture distribution patterns during drip irrigation. The water balance components were not analyzed in this study to ultimately define an optimal irrigation strategy. A different point of view was proposed by Twarakavi et al. (2009), who defined the field capacity point by simulating the drainage process with HYDRUS. These authors found that the field capacity point controls the irrigation schedule strategy but important variables such as transpiration were not considered. The role of transpiration was later on studied by Dabach et al. (2013), who analyzed different irrigation scheduling strategies through numerical simulations. In all these works, the assessment of irrigation scheduling disregards the strong link between the water balance components and the crop yield and economical profit.

In this work, we present a combined simulation and optimization framework aimed at obtaining the irrigation scheduling parameters that maximizes crop yield with minimum applied water while guaranteeing maximum net profit without soil salinization. Importantly, in comparison with other methods, the optimal irrigation solution provided here fully couples the water movement in the root zone with irrigation expenses and profits obtained by the soil–water–crop productivity relationship under potential soil salinization conditions. The framework permits to find the optimal control settings of an irrigated field that maximize the net profit obtained in a period of time T, given some forecasted climatic conditions.

The paper is organized as follows. We first present the formulation of the optimization framework in “Combined simulation–optimization framework”. Then, we apply the method to a specific field site in “Field application”. To do this, we set up a numerical model to simulate flow and transport through the vadose zone at the field site. The model is first shown capable of simulating water content data recorded at different depths over 1 year. After this, we optimize the irrigation problem and we discuss the results in “Simulation–optimization results”. The optimal irrigation solution is compared against traditional methods in this section. Finally, we also use the proposed framework to analyze the impact of soil properties on irrigation scheduling.

Combined simulation–optimization framework

In this section, we introduce the proposed simulation–optimization framework for irrigation scheduling in agriculture fields. In the optimization, the revenue from crop production is to be maximize as a function of a vector of irrigation parameters p. In general, the evaluation of the objective function requires two strongly coupled simulation components: (i) a flow and transport model that predicts the spatiotemporal evolution of water content in the soil and concentrations of chemical components, and (ii) a soil–water–crop productivity model that translates these results into crop yields and economic profit. The optimization problem and the simulation components are presented in general terms here to allow the use of different specific simulation features. Details of an application to a real field site are then provided in “Field application”.

Optimization problem

The irrigation optimization problem permits to find the optimal control settings of an irrigated field that maximize the net profit obtained in a period of time \(T\). Without loss of generality, we consider minimum fluctuations in the water price. That is to say that \(T\) is a short period of time that refers for instance to the duration of a crop growing season. Larger periods of time can be considered by simply incorporating a discount rate in the water price. Let us consider that the schedule of irrigation depends on a vector of control variables \(\mathbf{p}\) that characterize the irrigation rate \({q}_{i}\left(\mathbf{p},t\right)\) [l m−2 h−1]. The components of this vector can include for instance the pressure head threshold \(({h}^{*})\), which indicates the degree to which the soil can dry before irrigation is applied, and the duration of irrigation \((\tau )\) among others. The optimal irrigation scheduling strategy characterized by \(\mathbf{p}\) Eq. (1) is defined here as the one that is most productive and sustainable in terms of minimum water applied. This is mathematically formulated as the problem of finding the set of control variables \(\mathbf{p}\) that maximizes the crop net margin (NM) cost Eq. (2) subject to operational and functional constraints. The NM cost is the revenue from crop production, i.e., the gain margin (GM) Eq. (3) subtracted by the cost of crop production during an operational time \(t=T\). The cost of crop production consists of two main parts: capital costs and operation costs. Infrastructure cost (\(\mathrm{Capex}\)) Eq. (4) are the one-time expenses that usually incur during the purchase of land and equipment, i.e., expenses for bringing the irrigation field to an operable status. This includes the construction and installation of physical facilities such as the watering system, access roads, pipelines, drilling of wells and so on. Operation costs (\(\mathrm{Opex})\) Eq. (5) refers to the cost of specific activities incurred during the crop field lifecycle, which include equipment maintenance, product transport and overheads. From this, the optimization problem is mathematically formulated as

$${\mathbf{p}}_{\mathrm{opt}}=\underset{\mathbf{p}}{\mathrm{max}}\mathrm{NM}\left(\mathbf{p},T\right),$$
(1)

where,

$$\mathrm{NM}\left(\mathbf{p},T\right)=\mathrm{GM}\left(\mathbf{p},T\right)-\mathrm{Opex}\left(\mathbf{p},T\right)-\mathrm{Capex},$$
(2)
$$\mathrm{GM}\left(\mathbf{p},T\right)= {Y}_{a}\left(\mathbf{p},T\right){ C}_{y},$$
(3)
$$\mathrm{Capex}={C}_{m}+{C}_{c },$$
(4)
$$\mathrm{Opex}\left(\mathbf{p},T\right)={C}_{\mathrm{Fix}}+\underset{0}{\overset{T}{\int }}{q}_{i}\left(\mathbf{p},{t}^{{{\prime}}}\right) \left({C}_{\mathrm{wVar}}+{C}_{e}\right) d{t}^{{{\prime}}}.$$
(5)

The \(\mathrm{GM}\) [€ ha−1] is described here as the product between harvest price \({C}_{y}\) [€ t−1] and the actual crop yield \({Y}_{a}\) [t ha−1], \({C}_{\mathrm{Fix}}\) [€ ha−1] is the water fixed cost, \({C}_{\mathrm{wVar}}\) [€ m−3] is the water variable cost, \({C}_{e}\) [€ m−3] is the energy cost, \({C}_{m}\) [€ ha−1] is the irrigation system maintenance cost, and \({C}_{c}\) [€ ha−1] is the capital cost.

In practice, irrigation control settings must satisfy also some operational constraints. The constraints can include limitations on irrigation parameters as well as limitations on leaching water quality (solute concentrations). These features are incorporated into the optimization problem by constraining the solution to practical limitations and requirements, formally written as

$${\mathbf{p}}_{1}\le \mathbf{p}\le {\mathbf{p}}_{2},$$
(6)
$$g\left({\varvec{c}}\right)\le {\varvec{a}}.$$
(7)

The first set of constraints Eq. (6) refers to practical issues such as those determined from the water system capacity installed. Here, we only set an upper and a lower bound of \(\mathbf{p}\) to represent this. The other constraints Eq. (7) refer to solute concentrations, which can be used to limit for instance soil salinization. Note though that this methodology is not affected by the removal or the modification of the constraints included.

Solving this optimization problem results in the selection of optimal irrigation parameters. The solution will identify the maximum profitable solution subject to sustainable and feasible constraints. However, the simulation of crop yields requires two simulation components, one that simulates the distribution of water in the vadose zone as a result of irrigation and environmental conditions, and another that relates soil water content and water availability with crop yields. These two components are discussed in the following sections.

Flow and transport model

The simulation of water availability requires a numerical model capable to predict the flow of water in the unsaturated zone. The unsaturated flow is typically described by Richards’ equation (Richard 1931) Eq. (8), which can be written as

$$\frac{\partial \theta }{\partial t}=\nabla \cdot \left[K\left(h\right)\nabla h+K\left(h\right)\nabla z\right]-S,$$
(8)

where \(\theta\) [–] is the volumetric water content, \(h\) [hPa] is the water pressure head, \(t\) is time, \(S\) [T−1] is the source-sink term (includes root water uptake), and \(K\) [LT−1] is the unsaturated hydraulic conductivity. Solving this partial differential equation requires the knowledge of the soil water retention curve (SWRC), which is the relationship between the water content and the pressure head, and the hydraulic conductivity curve (HCC), which is the relationship between the hydraulic conductivity and the pressure head. These curves are characteristic for different types of soils.

The solution of Richards’ equation requires also the knowledge of initial and boundary conditions. For our purposes, the boundary condition specified at the soil surface plays an important role as it defines the amount of water infiltrated into the soil during irrigation. The water flux across the soil surface depends on both external conditions such as irrigation water rates \({q}_{i}\left(\mathbf{p},t\right)\) and the water content conditions in the soil. The corresponding boundary condition should represent for instance that run-off or ponding effect occurs when an irrigation rate exceeds the infiltration capacity of the soil or the fact that evaporation cannot exceed the capacity of the soil to deliver enough water to the soil surface. Solving Eq. (8) subject to a system-dependent boundary condition that limits the surface flux by the following two conditions Eqs. (9, 10) is often used to incorporate these features (Feddes et al. 1974; Neuman et al. 1974),

$$\left|q({z}_{\mathrm{top}},t)\right|\le \left|X(t)\right|,$$
(9)
$${h}_{d}\le h\le {h}_{s},$$
(10)

where \(q({z}_{\mathrm{top}},t)\) is the water flux at the soil top surface, \({z}_{\mathrm{top}}\) is the z-coordinate of the soil top surface, \(X(t)\) is the prescribed maximum potential rate of infiltration or evaporation given by meteorological conditions [LT−1], and \({h}_{d}\) and \({h}_{s}\) refer to dry and saturated water pressure heads, respectively. During evaporation, \(X\left(t\right)>0\) and this value represents the maximum evaporation rate \({E}_{p}\).

The numerical model should also describe the water extraction by plant roots. In this context, the root water uptake \(S\left(h,{h}_{\Phi },z,t\right)\) Eq. (11) is typically determined by the product of the water and salinity stress function \(\alpha \left(h,{h}_{\Phi }\right)\), the root density distribution function \(\beta \left(z,t\right)\) in the vertical direction, and the potential transpiration \({T}_{p}(t)\),

$$S\left(h,{h}_{\Phi },z,t\right)=\alpha \left(h,{h}_{\Phi , }\right){ \beta \left(z,t\right) T}_{p}\left(t\right),$$
(11)

where \({h}_{\Phi }\) is the osmotic pressure head. Several stress models can be found in the literature. Among them, the models presented by Feddes et al. (1978) and van Genuchten (1987) are the most widely used. Essentially, the latter model considers a smooth monotonic function with maximum root uptake at saturated conditions, while the other represents a piecewise linear function with maximum uptake within a pressure (saturation) interval. Importantly, only the formulation presented by Feddes considers a transpiration reduction near saturation. Assuming an isothermal system, the osmotic pressure head depends only on the solute concentrations of the chemical compounds in water, i.e., \({h}_{\Phi }=F({c}_{i})\). This relationship is often determined based on empirical relationships. Notice that the presented methodology can be used with any root water uptake model of the general form given by Eq. (11). For completeness, an example of application is presented later on in “Field application”, where we further describe in detail the root water uptake model for the Foradada field site and its characterization.

The simulation of solute concentrations is needed for two reasons: to evaluate constraints in concentrations (prevent soil salinization) and to estimate osmotic pressures. Solute transport in the unsaturated zone is typically described by the advection–dispersion equation Eq. (12), which is written as,

$$\frac{\partial \left(R\theta {c}_{i}\right)}{\partial t}=-\nabla \cdot \left({\varvec{q}}{c}_{i}\right)+\nabla \cdot \left(\theta \mathbf{D}\nabla {c}_{i}\right)+{f}_{i}, i=1,\dots ,{N}_{s},$$
(12)

where \({c}_{i}\) is the solute concentration of the ith chemical component, \({\varvec{q}}\) is the Darcy flux, \(\mathbf{D}\) is the hydrodynamic dispersion tensor, \({N}_{s}\) is the number of chemical components considered, R is the retardation factor, and \({f}_{i}\) is the concentration source–sink term.

Soil–water–crop productivity relationship

An estimation of the crop yield is required to evaluate the revenue from crop production during optimization. Crop productivity models can be complex as they include the interaction between genetics, physiology and environmental conditions such as water content. To facilitate the estimation procedure, one often considers an empirical model that relates crop water needs to crop yields. The model presented by Stewart et al. (1977) Eq. (13) has been widely accepted and recommended by the Food and Agriculture Organization of the United Nations (FAO) (Doorenbos and Pruitt 1975, hereafter FAO24). In addition, it has been recently used by several authors (Domínguez et al. 2012; Irmak et al. 2016; Martínez-Romero et al. 2017; Saadi et al. 2015). Note though that the methodology proposed is not limited by this model and it can also be used with other soil–water–crop productivity models if needed. This simplified crop yield model represents the seasonal pattern of crop water needs by different growing stages (four in the case of maize). The potential crop yield is penalized at each growth stage depending on the deficit of water, which is estimated by the relative discrepancy between potential (maximum) and actual evapotranspiration. The model determines that

$${Y}_{a}= {Y}_{p}\prod_{{k}=1}^{{{N}}_{{y}}}\left(1-{K}_{{y}_{k}}\left(1-{\left(\frac{{\mathrm{ET}}_{a}}{{\mathrm{ET}}_{c}}\right)}_{{k}}\right)\right),$$
(13)

where \({Y}_{p}\) [t ha−1] is the potential crop yield for the total growing season, \(k\) is the growing stage index, \({N}_{y}\) is the number of growing stages, and \({K}_{{y}_{k}}\) [–] is the crop yield response factor associated with the kth growing stage. \({\mathrm{ET}}_{a}\) and \({\mathrm{ET}}_{c}\) are the actual and potential accumulated crop evapotranspiration in each growing stage \(k\) [mm]. Note that the crop is under stress conditions when \({\mathrm{ET}}_{a}/{\mathrm{ET}}_{{c}} < 1\).

At each growing stage, according to Allen et al. (1998) the potential evapotranspiration \({\mathrm{ET}}_{c}\) Eq. (14) is estimated by the crop coefficient \({K}_{c}\) and the reference evapotranspiration \({\mathrm{ET}}_{0}\). The latter can be calculated by the Penman–Monteith equation (Allen et al. 1998), which is a function of the input daily mean temperature, the wind speed, the relative humidity and the solar radiation.

$${{\mathrm{ET}}_{c}={K}_{c}\mathrm{ ET}}_{0},$$
(14)

The \({\mathrm{ET}}_{a}\) Eq. (15) is estimated from the simulation of flow in the vadose zone, which gives the evaporation at the soil surface and the water uptake by plant roots \(S\left(h,{h}_{\Phi },z,t\right)\). From this, the \({\mathrm{ET}}_{a}\) of the kth growing stage taking place in the time interval \(({t}_{k},{t}_{k+1})\) can be determined by

$${\mathrm{ET}}_{a}=\underset{{t}_{k}}{\overset{{t}_{k+1}}{\int }}q\left({z}_{\mathrm{top}},{t}^{{{\prime}}}\right)I(q\left({z}_{\mathrm{top}},{t}^{{{\prime}}}\right)>0)d{t}^{{{\prime}}}+\underset{{t}_{k}}{\overset{{t}_{k+1}}{\int }}\underset{{z}_{\mathrm{top}}-{L}_{R}}{\overset{{z}_{\mathrm{top}}}{\int }}S\left(h,{h}_{\Phi },z,{t}^{{{\prime}}}\right)dz d{t}^{{{\prime}}},$$
(15)

where \(I(q>0)\) is an indicator function that is equal to one during evaporation (\(q>0\)) and zero otherwise, \({z}_{\mathrm{top}}\) is the z-coordinate of the soil top surface, and \({L}_{R}\) is the vertical length of the root zone.

Field application

Site description

In this section, we illustrate the applicability of the method in a real field setting. The study area considers a 25 ha commercial farm located in Foradada, Spain (Fig. 1). The field is divided into 23 irrigation sectors, which are fully covered by sprinkles. The irrigation rate is 6.5 l m−2 h−1. Two sectors can be irrigated at the same time in 1 day. The soil texture can be classified as a Silty Clay Loam (USDA soil taxonomy) with 28% Clay, 58.4% Silt and 13.6% Sand. Every year two different crops are grown: the first crop is usually canola and the growing season extends between winter and spring, during which the soil is wetted by rainfall events. The second crop is maize, which is grown during the summer and autumn season, when the soil is irrigated. Irrigation scheduling is needed to maximize productivity.

Fig. 1
figure 1

The Foradada field is located within the Segarra Garrigues system (ASG) canal. Solid irrigation sprinkler systems installed and distributed in 23 sectors are represented in different colours. Water content monitoring station and sampling is represented in red

Full size image

A field campaign was conducted aiming to measure soil water dynamics and soil hydraulic properties. For this, soil water content sensors were installed and undisturbed soil core samples were collected. Two EC-5 volumetric water content sensors (METER Group, Pullman, WA, USA) were installed at 10 and 20 cm depth in a representative location (Fig. 1). Sensor data were collected every 5 min with an accuracy of ± 0.03 cm3 cm−3 (Campbell and Devices 1986). One undisturbed soil sample was taken at 10 cm depth using a stainless-steel ring of 250 cm3 capacity. This soil samples were used to measure the SWRC and HCC with high precision (± 1.5 hPa) and over a wide range of pressures. This was achieved by combining the HYPROP, the WP4c, and the KSat devices (METER Group, Pullman, WA, USA). The HYPROP device is capable to measure SWRC and HCC. The WP4c device complements the SWRC data in the dry region. The KSat system does the same for HCC. A comparison of approaches has been reported by Schelle et al. (2013). These authors demonstrated that this combined method shows less noise than other traditional methods. The van Genuchten–Mualem SWRC model was used to fit experimental data using the HYPROP Fit software (METER Group, Pullman, WA, USA). Figure 2 shows SWRC and HCC measured by HYPROP, WP4c and KSat. Note that these devices provide SWRC and HCC with high resolution. SWRC describes a soil with a quite high water content retention capacity and a low air enter potential and slope. This is shown by the shape parameters α and n with values of 0.0678 and 1.186, respectively. HCC is characterized by a low \({K}_{s}\) value with 12 cm d−1, indicating slow wetting front movement during irrigation. To evaluate the vertical distribution of water uptake by plants we measured the root depth by pulling a plant off twice a month during the field campaigns. The maximum root depth registered was 55 cm after 78 days from sowing.

Fig. 2
figure 2

Soil water retention curve (SWRC) and hydraulic conductivity curve (HCC) measured by HYPROP, WP4c and KSat systems and fitted by van Genuchten–Mualem model

Full size image

Model setup

We use the HYDRUS-1D software package (Šimůnek et al. 2008, 2016) for simulating the one-dimensional movement of water and solute transport in variably saturated porous media at the field site. This code solves Richards’ equation to simulate water flow in the unsaturated zone and the advection–dispersion equation to simulate solute transport using numerical methods based on the Galerkin finite element method. We consider a 60 cm vertical soil profile representative of the Foradada field site. The domain was discretized into 101 segment elements. The column represents the movement of water through the soil profile associated with the water content sensors and core samples. A system-dependent boundary condition was imposed at the soil top surface according to Eqs. (9, 10). Since the water table is far below, a free drainage boundary condition was imposed at the column bottom, i.e., \(q=-K(h)\). Considering the soil type and the soil water content measurements, initial conditions were set to \(\theta =0.25\). A multiplicative model was used to represent the water and salinity stress functions, i.e., \(\alpha \left(h,{h}_{\Phi }\right)=\alpha (h)\alpha ({h}_{\Phi })\), where \(\alpha (h)\) and \(\alpha ({h}_{\Phi })\) are the water and salinity stress functions, respectively. In this case, we have considered the Feddes et al. (1978) Eq. (16) model to represent the water stress function,

$$\alpha \left( h \right) = \left\{ \begin{gathered} \frac{{h - h_{4} }}{{h_{3} - h_{4} }}\quad \quad h_{3} > h > h_{4} \hfill \\ 1\quad \quad \quad \quad \;\; h_{2} \ge h \ge h_{3} \hfill \\ \frac{{h - h_{1} }}{{h_{2} - h_{1} }}\quad \quad h_{1} > h > h_{2} \hfill \\ 0\quad \quad \quad \quad \;\; h \le h_{4} \;{\text{or}}\;h \ge h_{1} \hfill \\ \end{gathered} \right.,$$
(16)

which describes that the plant suffers water stress outside the pressure head range \(({h}_{2},{ h}_{3})\). The water stress reduction decreases linearly from those pressure points and gets to a minimum (\(\alpha =0\)) below or above \({h}_{4}\) and \({h}_{1}\). Another function used is the salinity stress function which is defined using the following threshold-slope salinity stress reduction function (Maas and Hoffman 1977) Eq. (17),

$$\alpha \left( {h_{\Phi } } \right) = \left\{ \begin{gathered} 1,\quad \quad \quad \quad \quad a \le h_{\Phi } \le 0 \hfill \\ 1 + b\left( {h_{\phi } - a} \right),\, a > h_{\Phi } > - \frac{1}{b} \hfill \\ 0,\quad \quad \quad \quad \quad h_{\Phi } \le a - \frac{1}{b} \hfill \\ \end{gathered} \right..$$
(17)

Here, the salinity threshold value \(a\) quantifies the minimum osmotic head above which root water uptake occurs without reduction, and the slope b determines the fractional root water uptake decline per unit increase in salinity below the threshold. The parameters adopted to define these stress functions were chosen from the HYDRUS internal database and are summarized in Table 1. The transport model is simplified to simulate only one representative chemical component, i.e., electrical conductivity \(\mathrm{EC}\), which is used as a surrogate for chloride, which is considered to be a conservative species (non-reactive). We neglect salt precipitation and dissolution processes. The osmotic pressure is assumed to be proportional to \(\mathrm{EC}\). Considering that \({h}_{\phi }\) and \(\mathrm{EC}\) are expressed here in the same units, we essentially have that \(\mathrm{EC}={h}_{\phi }\) (Simunek and Sejna 2014). Based on the root depth measurements taken during the field campaigns, we represent the vertical spatial distribution of water uptake by plants through Hoffman and Van Genuchten model (1983) with a root depth \({L}_{R}\) of 55 cm.

Table 1 Water and salinity stress function parameters
Full size table
$$\beta \left( z \right) = \left\{ \begin{gathered} \frac{{1.66667}}{{L_{R} }}\quad \quad \quad \quad \quad \quad z > z_{{{\text{top}}}} - 0.2L_{R} \hfill \\ \frac{{2.0833}}{{L_{R} }}\left( {1 - \frac{{z_{{{\text{top}}}} - z}}{{L_{R} }}} \right)\quad z \in (z_{{{\text{top}}}} - L_{R} ,z_{{{\text{top}}}} - 0.2L_{R} ) \hfill \\ 0\quad \quad \quad \quad \quad \quad \quad \quad \quad z < z_{{{\text{top}}}} - L_{R} \hfill \\ \end{gathered} \right..$$
(18)

An initial estimation of model properties (e.g., SWRC and HCC) was known from core sample measurements. The model parameters were then calibrated manually to reproduce the recorded soil moisture data obtained at the field site. The calibrated parameters were not substantially different from the initial estimation. Table 2 summarizes the measured and calibrated parameters. Parameters modified in the calibration process were \(\alpha\) and \(n\). The simulation considered 240 days, from February 9th 2017 until October 31st 2017. Two different crops growing at different periods of time were accounted for during the simulation; canola from February, 9th 2017 to May, 31st 2017, and maize from June, 1st 2017 to October, 31st 2017. Meteorological parameters were downloaded from the nearest available weather station, located at 15 kms from the site, to compute \({\mathrm{ET}}_{0}\), which was then converted into daily \({\mathrm{ET}}_{c}\) values with the \({K}_{c}\) coefficients shown in Table 3. \({K}_{c}\) coefficients were extracted from the works of Domínguez et al. (2012) and Martínez-Romero et al. (2017), conducted in a maize field located in Castilla la Mancha, Spain. \({K}_{c}\) coefficients were determined from field temperature and an estimation of the growing degree days (GDD). Weather conditions from Foradada field and Castilla la Mancha are similar. During the field campaigns, we visually corroborated the time duration of the different phenological stages proposed by Domínguez et al. (2012) and Martínez-Romero et al. (2017). The potential evaporation and transpiration values needed in the root water uptake model, Eqs. (9, 10, 11), were calculated by partitioning \({\mathrm{ET}}_{c}\) into potential evaporation \({E}_{p}\) and transpiration \({T}_{p}\) based on the Canopy Cover (Raes et al. 2010), which determines that \({\mathrm{ET}}_{c}=\beta {E}_{p}+{(1-\beta )T}_{p}\), being \(\beta\) the soil cover fraction. Figure 3 compares simulation results with soil moisture field measurements obtained at two different depths. Simulations are in good agreement with soil moisture data, except for a relatively small underestimation of the water content measured at depth 20 cm by a factor of about 1.15 during the first 200 days after sowing. Table 4 shows several goodness-of-fit statistics, such as Willmott index (Willmott 1981), calculated for both depths.

Table 2 Foradada and loamy sand soil hydraulic parameters used in simulations
Full size table
Table 3 Canola crop coefficient (Kc1), maize crop coefficients (Kc2) applied for ETc, where Stage I represents the initial period; Stage II is crop development; Stage III mid-season; Stage IV is late season
Full size table
Fig. 3
figure 3

Comparison between daily soil moisture field measurements and the soil moisture output from the validation model at 10 and 20 cm depth

Full size image
Table 4 Statistical index calculated with observed and simulated water content values
Full size table

Optimal irrigation scheduling problem setup

We applied the simulation–optimization framework presented in “Combined simulation–optimization framework” to the Foradada field site to estimate optimal irrigation scheduling parameters during a growing season under highest \({\mathrm{ET}}_{c}\) demand, defined by the dry conditions without rainfall events taking place between 2008 and 2017. Thus, we considered the most unfavourable weather conditions, even though the methodology allows for other weather conditions as well. The parameters used here were directly adopted from “Site description” and “Model setup”. The corresponding \({K}_{c}\) coefficients are presented in Table 3. The year with more water demand was 2016 with an atmosphere demand of 478 mm and a total rainfall of 80 mm. During this period of time, the maximum and minimum temperatures were 39 °C and 21 °C, respectively. Initial conditions were \(\theta =0.25\) and \(\mathrm{EC}\) = 0.6 dS/m, which represents an average water content of the field site without salinity problems. Simulations considered the entire growing season of maize, from June 15th (sowing) to November 11th (harvesting). This crop is cultivated during the dry season; thus, irrigation should be applied to ensure crop development. The rest of the year, rainfall events maintains the soil under wet conditions and irrigation is not necessary. Consequently, it is not necessary to irrigate during this period of time. Data necessary to evaluate the crop net margin cost \(\mathrm{NM}\) are summarized in Table 5, mostly provided by the Aigües Segarra Garrigues (ASG) company. Two control irrigation parameters were used to characterize water irrigation rates, i.e., the pressure head threshold \({h}^{*}\) observed at a control point and the duration of irrigation \(\tau\). The control point is located at the vertical midpoint of the maximum maize root length (20 cm below soil surface). We considered the depth where the maize main roots are located. The pressure head at the control point triggers irrigation when \({h<h}^{*}\) (when the soil is dry). At this moment, the model applies water at a constant irrigation rate \({q}_{i}\) for an irrigation time \(\tau\). Since the watering capacity is fixed by the type of irrigation equipment installed, the irrigation rate is not considered to be a control parameter in this case (\({q}_{i}\) is set constant to 6.5 l m−2 h−1). Figure 4 presents an illustrative sketch of the irrigation. The salt concentration in the irrigated water was measured with a portable water electrical conductivity meter (Hanna Instruments). The EC of irrigated water was 0.4 dS m−1. The optimization is constrained to fulfil that \(\mathrm{EC}<3.4\) dS m−1 (below this threshold maize is not under salinity stress).

Table 5 Parameters necessary to apply Stewart and net margin (\(\mathrm{NM}\)) equations
Full size table
Fig. 4
figure 4

On the right, boundary conditions imposed in the model where triggered irrigation is a function of \({h}^{*}\) and \(\tau\). On the left, a synthetic case about how the model triggers the irrigation, when \({h}^{*}\) and \(\tau\) are defined

Full size image

A large number of algorithms can be used to maximize the crop net margin cost function \(\mathrm{NM}\) with constraints. Here, we chose to maximize \(\mathrm{NM}\) over a given range by brute force, which simply consists in exploring the parameter space to find the global maximum. This can be inefficient in practical applications but provides detail insights about irrigation scheduling as well as the full shape of the \(\mathrm{NM}\) cost function, which is the objective here. To do this, the parameter space \((\tau , {h}^{*})\) was discretized into a \(4\times 10\) regular mesh, where \(\tau\) ranges between 1 and 4 h and the threshold pressure head \({h}^{*}\) varies between − 100 and − 10 kPa. Results shown in the next section will demonstrate that the optimal set of parameters for soil irrigation lies within this region because pressure heads \({h}^{*}\) below − 100 kPa result in crops under strong water stress conditions. This agrees with the works of Bianchi et al. (2017) and Thompson et al. (2007), who found that the optimal pressure head for maize irrigation is about − 40 kPa.

To analyze the performance of the method, we compared the optimal irrigation results obtained with our proposed framework with those given by a traditional irrigation method. The traditional irrigation scheduling method is based on water requirements and consists in irrigating as much water as that evapotranspirated in the previous week. To this end, the farmer must devise an irrigation calendar. Thus, to simulate the traditional method, the weekly \({\mathrm{ET}}_{c}\) value is calculated and this volume of water is applied over the next week. This amount of water is uniformly distributed over the next week. All other parameters are kept the same. The relative difference between the Net Margin obtained with the traditional method (\({\mathrm{NM}}_{\mathrm{trad}})\) and the optimal one is computed as follows:

$${\Delta }_{r}\mathrm{NM}=\frac{\mathrm{NM}-{\mathrm{NM}}_{\mathrm{trad}}}{\mathrm{NM}}.$$
(19)

Simulation–optimization results

In this section, we present the simulation–optimization results of the Foradada irrigation scheduling problem. Figure 5a shows a map of the net margin \(\mathrm{NM}\) function obtained from all the irrigation strategies simulated as a function of \({h}^{*}\) and \(\tau\). A clear \(\mathrm{NM}\) maximum value of 2791 € ha−1 can be seen in this figure for \(\tau\) = 1 h and \({h}^{*}\) = − 40 kPa. This optimal irrigation strategy represents a short but moderately frequent irrigation, i.e., moderate \({h}^{*}\) value, which results from balancing the gain margin \(\mathrm{GM}\) with the operational expenses \(\mathrm{Opex}\). To better appreciate this, the \(\mathrm{NM}\) map is decomposed into the corresponding \(\mathrm{GM}\) and \(\mathrm{Opex}\) contributions, respectively, shown in Fig. 5b and c. The maximum gain \(\mathrm{GM}\) requires a more frequent irrigation with \(\tau\) = 1 h and \({h}^{*}\) = − 25 kPa (more irrigated water) but the operational expenses \(\mathrm{Opex}\) significantly increase in this region, economically penalizing the gain margin \(\mathrm{GM}\). Thus, even though the irrigation strategy \(\tau\) = 1 h and \({h}^{*}\) = − 25 kPa is the most productive in terms of crop yield, the maximum net margin \(\mathrm{NM}\) is obtained at higher suction values \({h}^{*}\) = − 40 kPa (less irrigation) with the same irrigation frequency \(\tau\) = 1 h. Thus, the expenses related with the applied volume of water makes the optimal irrigation strategy to substantially depart from the most productive strategy. On the contrary, \(\mathrm{Opex}\) is minimum when the frequency of irrigation is very small, implying that less water is used for irrigation. Applying the optimal irrigation strategy, the volume of water applied is 470 mm. Figure 5d shows the relative difference between the traditional method and the optimal solution obtained with the proposed methodology. Results show that the proposed method can increase the net margin by 7%, decreasing by 6% the total amount of water applied at the end of the campaign, and reducing by 5% the costs associated by irrigation. We also note that the \(\mathrm{NM}\) function can give smaller values than the traditional one when either irrigations strategies apply a large volume of water (increasing \(\mathrm{Opex}\)) or \({h}^{*}\) is too small (decrease in \(\mathrm{GM}\)). This situation is far away from the optimal irrigation strategy and thereby the simulation–optimization method seems mandatory in routine field applications.

Fig. 5
figure 5

a Net margin (NM) Foradada soil results and objective functions elements, being b GM, gain margin; c \(\mathrm{Opex}\), operational costs and d \({\Delta }_{r}\mathrm{NM}\), fractional difference between \(\mathrm{NM}\) and \({\mathrm{NM}}_{\mathrm{trad}}\). Dash line in \({\Delta }_{r}\mathrm{NM}\) map represents when the relative increment is zero

Full size image

A more profound understanding of the difference between the traditional irrigation scheme and the optimal irrigation method (\(\tau\) = 1 h and \({h}^{*}\) = − 40 kPa) can be seen in Fig. 6, which compares the temporal evolution of the water content resulting from both methods at four different soil depths during the growing season. The optimal irrigation strategy applies water for 1 h during several days and stops when \({h}^{*}\) approaches − 40 kPa. The traditional method applies water every day (the total of volume of water applied is 478 mm). The main difference is that the optimal solution provides a more uniform variation of the water content over the entire season, i.e., from 10 to 40 cm depth water content oscillates between 0.20 and 0.32 cm3 cm−3, compared with the traditional method, which fluctuates between 0.16 and 0.34. Moreover, the optimal method is capable to increasing the water content globally in the root zone, in particular at 10 cm depth, whereas the traditional method stresses the system during periods with higher \({\mathrm{ET}}_{c}\) demand by letting the water content to decrease up to 0.16 cm3 cm−3. This water stress affects crop productivity because water content conditions are not optimal. We also note that when the volume of water applied is less than 5 mm d−1, the wetting front does not arrive to the deepest soil parts, compromising root water uptake. This way, even though the traditional method applies the total volume of water evapotranspirated during the growing season, this is seen not enough to maintain constant the water content during the growing season. In contrast, optimal irrigation method is based on soil water status and how water moves through the root zone. Results show that the optimal irrigation can better maintain the crop under soil water optimal conditions during the entire growing season.

Fig. 6
figure 6

Irrigation scheduling from both strategies, where irrigation, volumetric water content dynamics at different depths and the potential evapotranspiration demand (\({\mathrm{ET}}_{c}\)) are plotted

Full size image

To better understand the results, Fig. 7 describes the impact that each irrigation strategy produces on \({T}_{a}\), \({E}_{a}\) and soil \(\mathrm{EC}\). The optimal strategy is highlighted in these figures with a black circle. Figure 7a presents the simulated transpiration values obtained as a function of \({h}^{*}\) and \(\tau\). Remarkably, the optimal irrigation strategy is not the one that produces more transpiration but lies within the plant water stress region, i.e., \(\alpha \left({h}^{*}\right)\) is smaller but close to 1. This means that the gain in crop productivity obtained for \(\alpha =1\) does not compensate the expenses associated with the increase in irrigated water. Results also show a decrease in transpiration when irrigation applies water during several hours. This is caused by the saturation of the top horizon of the soil (see Fig. 8).

Fig. 7
figure 7

a Actual transpiration (\({T}_{a}\)), b actual evapotranspiration (\({E}_{a}\)) and c electrical conductivity at the root zone (\(\mathrm{EC}\)) resulting from Foradada soil showing all the irrigation strategies simulated. Circles show the strategy who provides maximum net margin (\(\mathrm{NM}\))

Full size image
Fig. 8
figure 8

Wetting patterns when irrigation strategy is fixed at \(\tau\) = 1, 2, 3, 4 h and \({h}^{*}\) = 10 kPa. Some water stress function parameters (\({h}_{1}\) and \({h}_{2}\)) are plotted indicating when transpiration decreases as a consequence of water logging

Full size image

Figure 7b plots the corresponding evaporation values. According to Philip (1956) and Ritchie (1972), we distinguish between two evaporation regions: Range I represents an energy-limited evaporation process where the soil surface is wetted by irrigation and water evaporates from a thin soil surface layer; Range II represents a falling-rate evaporation process that occurs when water content flows from the soil layer below. Results show that the optimal irrigation strategy is energy-limited because the Foradada soil has a high-water content retention capacity.

Figure 7c presents the simulated \(\mathrm{EC}\) values at the end of the season as a function of \({h}^{*}\) and \(\tau\). These values are estimated as the average from the four observation points inserted at different depths. Even though none of the irrigation strategies exhibit salinity stress (above threshold) an increase in \(\mathrm{EC}\) is seen in all simulations. Consequently, salts are accumulated through the root zone presenting a maximum at the end of the season. This is more pronounced for small \({h}^{*}\) values. This increase could significantly affect the crop productivity in the following seasons and reflects the importance of rainfall events in wet periods (winter and spring). This volume of water should be responsible for flushing out the salts accumulated during dry periods. Our simulated scenarios considered only dry conditions without rainfall events and thereby this flushing mechanism could not be seen.

Figure 8 displays the wetting front resulting from different irrigation events for four different irrigation strategies. One can observe that when water is applied during several hours (i.e., more than indicated by the optimal irrigation strategy), the top horizon of the soil becomes almost saturated and thereby water uptake is negligible in this region (pressure head overpasses the \({h}_{2}\) threshold). This prolonged saturation of the top soil horizon due to poor drainage inhibits crop productivity. This is in turn reflected by smaller \({T}_{a}\) values in Fig. 7a.

To evaluate the impact that soil hydraulic properties have on the optimal irrigation strategy, we also solved the optimization–simulation problem considering a different soil type. The chosen soil hydraulic properties have been downloaded from the Rosetta database (Schaap et al. 2001) available from HYDRUS 1D. Table 2 shows the soil hydraulic parameters of a loamy sand soil. In this specific case, \(\alpha\) and \(n\) parameters represent a soil with less soil moisture retention capacity than the Foradada soil, and \({K}_{s}\) is substantially higher than the soil previously studied (\({K}_{s}=350\) cm d−1). Taking these observations into account, a lower soil moisture retention with a faster wetting front movement is expected. Figure 9 shows the map of the \(\mathrm{NM}\), \(\mathrm{GM}\), and \(\mathrm{Opex}\) function depending on \({h}^{*}\) and \(\tau\). The linear-like shape of these functions is strikingly different from those of Foradada, which exhibited a strong nonlinear behavior close to the optimal value. Consequently, results show that in this case \({h}^{*}\) is mostly controlling \(\mathrm{NM}\) and \(\tau\) has little effect. An optimal irrigation strategy is found for \(\tau\) = 1 h and \({h}^{*}\) = − 10 kPa. The maximum \(\mathrm{NM}\) value is 2719 € ha−1, which is slightly smaller than that of the Foradada soil of 2791 € ha−1. This irrigation strategy represents a short but very frequent irrigation. We conclude then that the optimal irrigation strategy can drastically change from one soil type to another.

Fig. 9
figure 9

a Net margin (\(\mathrm{NM}\)) loamy sand soil results and objective functions elements, being b \(\mathrm{GM}\), gain margin and c \(\mathrm{Opex}\), operational costs, d actual transpiration (\({T}_{a})\), e actual evapotranspiration (\({E}_{a})\) and f electrical conductivity at the root zone (\(\mathrm{EC}\)). Circles show the strategy who provides maximum net margin (\(\mathrm{NM}\))

Full size image

The physical process by which the \(\mathrm{NM}\) function associated with a loamy sand soil is now mostly controlled by the pressure head threshold \({h}^{*}\) and not \(\tau\) can be seen from the dependence of \({T}_{a}\), \({E}_{a}\), and \(\mathrm{EC}\) with \({h}^{*}\) shown in Fig. 10. Note that in all cases these results are substantially different than those of the Foradada soil (Fig. 7). In this case, results indicate that the irrigation water can easily infiltrate and redistribute through the entire root zone, following a falling-rate evaporation process. In contrast to the Foradada soil, a more permeable soil leads to an optimal irrigation strategy (black circle) that favors maximum transpiration rates with zero plant water stress. The optimal strategy also inhibits waterlogging and provides minimum salinization compared to other irrigation strategies obtained with different pressure head thresholds (Fig. 10). Again, this can be explained by an effective percolation of water through the root zone, which gives good internal drainage.

Fig. 10
figure 10

a Actual transpiration (\({T}_{a})\), b actual evapotranspiration (\({E}_{a})\) and c electrical conductivity at the root zone (\(\mathrm{EC}\)) resulting from loamy sand soil showing all the irrigation strategies simulated. Circles show the strategy who provides maximum net margin (\(\mathrm{NM}\))

Full size image

Conclusions

Irrigation scheduling in agriculture is crucial for saving water while guaranteeing maximum crop yields in arid regions as well as in future areas affected by climate change and water scarcity. However, irrigation scheduling is typically conducted based on simple water requirement calculations without accounting for the strong link between water movement in the root zone, crop yields and irrigation expenses. In this work, we have presented a combined simulation and optimization framework aimed at estimating irrigation parameters that maximize the crop net margin cost subject to operational and functional constraints. The simulation component couples the movement of water in a variably saturated porous media driven by irrigation with plant water uptake and crop yields. The optimization component assures maximum gain with minimum cost of crop production during a growing season.

An application of the method was presented in the Foradada irrigation field test site, where soil hydraulic parameters represented a soil with high water content capacity and slow wetting front percolation. The method was demonstrated to yield optimum irrigation parameters at the site (irrigation duration of 1 h and a frequency determined by a pressure head threshold of − 40 kPa). These parameters were substantially different from those estimated with traditional water requirement methods based on previous evapotranspiration values. Results have shown that even though the volume of water evapotranspired during the growing season is fully replaced by the traditional method, the way to schedule irrigation does not guarantee that water content values vary within an optimal water content fringe. Prolonged saturation of the top soil horizon is often shown to occur with the traditional method. On the contrary, the optimal irrigation scheduling solution prevents waterlogging and provides a more constant value of water content during the entire growing season within the root zone. As a result, the crop net margin cost exhibited increase with respect to the traditional method by 7%. The optimal irrigation solution is also demonstrated to properly balance the crop gain margin and operational expenses. Results have shown that even though some strategies can be more productive, irrigation expenses counterbalance the economic benefit ultimately leading to a compromise between them. At this stage, we highlight the practical advantages of the method proposed compared to traditional water requirement methods: (i) agriculture stakeholders can obtain better crop gain margins to achieve the same crop productivity; and (ii) irrigation scheduling parameters can be known prior to the start of the growing season.

To implement the method some measurements are required. First, it is necessary to measure soil hydraulic properties to provide the model the information necessary to simulate soil moisture through the root zone. Second, it is recommended to have a weather station in the study to calculate the potential evapotranspiration demand. Unfortunately, some stakeholders have not the opportunity to have installed a weather station in the field, in this case, weather data must be downloaded from the nearest station. It is also highly recommended to install pressure head potential sensors to calibrate the model and verify that irrigation is triggered at the correct threshold. The impact of soil hydraulic properties has been also analyzed by assuming another soil type (loamy sand soil). Soil hydraulic parameters described a soil with less water retention capacity than the Foradada soil and easier wetting front percolation. Results have shown that the optimal solution strongly depends on the type of soil. In this case, the frequency of irrigation was much larger given by a pressure head threshold of − 10 kPa.

Finally, our results indicate that the irrigation optimization algorithms that will be developed in the next future should account for the water movement trough the root zone instead of water requirements estimations only. Moreover, irrigation scheduling should be based on soil pressure head status to guarantee optimal crop performance during the growing season. To achieve this, an accurate soil hydraulic characterization is crucial as well as the installation of pressure head sensors to monitor the pressure head status during the season.