Computers and Electronics in Agriculture 150 (2018) 312–327

Contents lists available at ScienceDirect

Computers and Electronics in Agriculture
journal homepage: www.elsevier.com/locate/compag

Original papers

A hybrid ﬁnite volume-ﬁnite element model for the numerical analysis offurrow irrigation and fertigation

T⁎

Giuseppe Brunettia, , Jirka Šimůnekb, Eduardo BautistacaUniversity of California, Davis, Department of Land, Air and, Water Resources, CA 95616, USA
University of California, Riverside, Department of Environmental Sciences, CA 92521, USAc
USDA, ARS, U.S. Arid Land Agricultural Research Center, Maricopa, AZ 85138, USAb

A R T I C LE I N FO

A B S T R A C T

Keywords:
Surface ﬂow
Subsurface ﬂow
Furrow irrigation
HYDRUSWinSRFRThis study presents a hybrid Finite Volume – Finite Element (FV-FE) model that describes the coupled surfacesubsurface ﬂow processes occurring during furrow irrigation and fertigation. The numerical approach combines
a one-dimensional description of water ﬂow and solute transport in an open channel with a two-dimensionaldescription of water ﬂow and solute transport in a subsurface soil domain, thus reducing the dimensionality ofthe problem and the computational cost. The modeling framework includes the widely used hydrological model,
HYDRUS, which can simulate the movement of water and solutes, as well as root water and nutrient uptake invariably-saturated soils. The robustness of the proposed model was examined and conﬁrmed by mesh and time
step sensitivity analyses. The model was theoretically validated by comparison with simulations conducted withthe well-established model WinSRFR and experimentally validated by comparison with ﬁeld-measured data froma furrow fertigation experiment conducted in the US.

1. Introduction
Agriculture is among the most important human activities due to itsrole in the food supply chain. According to one of the latest FAO reports(FAO, 2017), agricultural production more than tripled between 1960
and 2015. This remarkable expansion has been accompanied with adramatic increase in the use of irrigation and fertilization, and an associated signiﬁcant environmental footprint, thus posing importantsustainability issues. Irrigation for agricultural purposes accounts for
70% of all water withdrawn from aquifers, lakes, and streams (FAO,
2011a). Furthermore, the worldwide demand for fertilizers has grownby 25% in the last decade, and this trend is expected to continue incoming years (FAO, 2015), increasing the risk of nutrient pollution ofwater bodies. Thus, the transition towards more eﬃcient and sustainable irrigation and fertigation strategies is necessary.
Although slowly replaced by pressurized irrigation in developed
regions such as The United States, Europe, and Israel, surface irrigation
systems continue to be a preferred irrigation method in developing
countries. In 2011, surface systems accounted for 96.8% of irrigatedsurfaces in Southern and Eastern Asia (FAO, 2011b). In particular, basin
and furrow irrigation were the most widespread techniques among
farmers. When correctly performed, furrow fertigation can increase theeﬃciency of fertilizer use and crop fertilizer uptake, compared to⁎Corresponding author.
E-mail address: giusep.bru@gmail.com (G. Brunetti).

https://doi.org/10.1016/j.compag.2018.05.013
Received 21 February 2018; Received in revised form 5 May 2018; Accepted 8 May 20180168-1699/ © 2018 Elsevier B.V. All rights reserved.

traditional techniques (Fahong et al., 2004; Horst et al., 2005; Siyal
et al., 2012; Šimůnek et al., 2016a). To be successful, furrow irrigationand fertigation systems should be designed and managed so that theapplication and distribution of water and fertilizer are eﬃcient anduniform, with minimal surface runoﬀ at the lower end of the ﬁeld andwith minimal deep drainage and leaching below the crop root zone(Šimůnek et al., 2016a). To optimize the performance of furrow irrigation systems, numerical models may play an important role since
they represent powerful tools for assessing irrigation and fertigationeﬃciency.
Furrow irrigation and fertigation are coupled surface-subsurface
processes. Water and solute are injected at the soil surface at one side ofan open channel, hence generating a sharp front that moves along a
furrow while water and nutrients inﬁltrate in the underlying soil.
Therefore, furrow irrigation and fertigation are (physically, mathematically, and numerically) described as coupled three-dimensional ﬂow
and transport processes in both surface and subsurface domains.
Nevertheless, solving a fully coupled system of 2D Shallow Water and3D Richards equations would require signiﬁcant computing resources
and would also likely pose substantial stability issues, mainly stemmingfrom the high nonlinearity of the governing equations. This presents a
problem in the modeling of ﬂooding and drying processes over porous
beds.

Computers and Electronics in Agriculture 150 (2018) 312–327

G. Brunetti et al.

Most of the existing furrow irrigation models adopt a Lagrangian
approach, which uses a computational grid that moves along with thewet/dry interface. Although elegant and accurate, this type of approach
can lead to coupling issues between the surface and subsurface models,
mainly because the grid must be regenerated each time the wet/dry
interface moves, and the computational nodes often must be added
during ﬂooding or removed during recession to reduce grid distortion
error. However, while surface processes are discontinuous (surface ﬂow
and transport occur only during irrigation events), subsurface processes
are continuous (subsurface ﬂow and transport continues between irrigation events). Therefore, a ﬁx grid for the subsurface domain is usuallyused, and values of interest (e.g., inﬁltration, soil moisture, pressure
head, etc.) need to be interpolated between surface nodes, thus leading
to a hybrid Lagrangian-Eulerian approach. Nevertheless, results have
proven to be fairly sensitive to the interpolation strategy (Wöhling
et al., 2006). Lazarovitch et al. (2009) applied the moment analysis
techniques to describe the spatial and temporal subsurface wetting
patterns resulting from furrow inﬁltration and redistribution. Furthermore, most of the existing models are based on the ﬁnite diﬀerence
method (e.g., Tabuada et al., 1995; Abbasi et al., 2003; Perea et al.,2010), which can yield spurious oscillations at ﬂow discontinuities
unless ﬁrst-order accurate (upwind) schemes or artiﬁcial dissipation are
employed.
Godunov-type Finite Volume (FV) schemes have been successfully
applied to simulate overland ﬂow over pervious and impervious lands(Bradford and Katopodes, 2001; Bradford and Sanders, 2002; Brufau
et al., 2002; Burguete et al., 2009; Dong et al., 2013). The FV schemes
solve the integral form of the overland ﬂow and solute transport
equations, thus being mass conservative both globally and locally.
Numerical ﬂuxes are evaluated at the cell faces, thus guaranteeing astraightforward and eﬃcient treatment of the dry bed problem by allowing for ﬂooding and drying of ﬁxed computational cells. Furthermore, numerical oscillations near discontinuities can be eliminatedusing ﬂux limiters (Bradford and Katopodes, 2001). However, theirapplication to furrow irrigation and fertigation has been rather limited
and has not involved a coupled mechanistic description of the subsurface domain.
Thus, the main goal of this study is to develop and validate a hybrid
Finite Volume-Finite Element (FV-FE) reduced-order model capable ofdescribing coupled surface-subsurface ﬂow and transport processes involved in furrow fertigation. The proposed numerical approach combines a one-dimensional FV description of coupled water ﬂow and solute transport in the surface domain with a two-dimensional
mechanistic FE description of ﬂow and transport in the variably-saturated zone, thus reducing the dimensionality (3D) of the problem andassociated computational cost. The modeling framework includes the
widely-used FE model, HYDRUS-2D, whose numerical features signiﬁcantly increased the overall modeling ﬂexibility. The proposedmodel is the ﬁrst attempt to include HYDRUS-2D in a coupled surfacesubsurface furrow fertigation model, thus representing a new contribution to this ﬁeld. The problem is addressed in the following way.
First, the hybrid FE-FV model is theoretically validated against the wellestablished model WinSRFR (Bautista et al., 2009) using synthetic validation scenarios. Preliminary mesh and time step sensitivity analyses
are performed to evaluate the accuracy and the robustness of the proposed numerical approach. Next, the model is validated against measured data from an experimental facility in California, US.

Therefore, several models have been proposed in the literature that
reduce the numerical complexity of a fully 3D model. Katopodes andStrelkoﬀ (1977) used a dimensionless formulation of the governing
equations to show that when the Froude number is small, which is typical under irrigation conditions, the inertial terms in the shallow waterequations are negligible. With this in mind, they developed the ﬁrst
zero-inertia model for irrigation (Strelkoﬀ and Katopodes, 1977). In theearly 80s, Walker and Humpherys (1983) proposed a furrow irrigation
model based on a 1D kinematic wave (KW) approximation of open
channel ﬂow coupled with the modiﬁed Kostiakov equation describing
the inﬁltration process. While the model was assessed against experimental data with satisfactory results, the adoption of the KW equation
limited its application to open-ended furrows. Furthermore, the model
did not include any description of the subsurface water dynamics. Lateron, Oweis and Walker (1990) replaced the KW approximation with the
one-dimensional (1D) zero-inertia (ZI) equation, which represents a
more realistic approximation of the Shallow Water equations. A moredetailed fertigation model was ﬁrst proposed by Abbasi et al. (2003)
and later improved by Perea et al. (2010). In these two studies, waterﬂow and solute transport in an open channel were described using the
1D Zero-Inertia and advection-dispersion equations, respectively. The
modiﬁed Kostiakov equation was used to calculate inﬁltration at eachtime step. Although the model was veriﬁed with good results against
experimental data from four experimental sites, subsurface ﬂow processes were again neglected.
To provide a better and more complete description of coupled waterﬂow and solute transport in the soil, several other models have been
proposed in the literature. For example, Zerihun et al. (2005) coupled
the 1D zero-inertia equation with the HYDRUS-1D model (Šimůnek
et al., 2016b), which numerically describes water ﬂow and solutetransport as well as root water and nutrient uptake in 1D variably-saturated porous media. The computational framework, which targetedbasin irrigation, was based on the iterative coupling between the surface and subsurface models and was validated against measured datawith good results. Ebrahimian et al. (2013) extended this concept tofurrows by coupling a 1D furrow fertigation model with the two-dimensional HYDRUS-2D model (Šimůnek et al., 2016b). The coupledmodel satisfactorily reproduced the overland transport as well as thesolute transport in the soil proﬁle. However, the surface and subsurface
components were solved separately, leading to an uncoupled numericalframework. As pointed out by Furman (2008), theoretically, the higher
the level of coupling, the higher the solution accuracy. This is mainly
due to the high nonlinearity of the involved processes, as well as totheir diﬀerent time scales. For instance, overland ﬂow is generally
faster than inﬁltration, thus requiring a diﬀerent temporal resolution.
This temporal misalignment poses signiﬁcant numerical issues since anapproximation is needed to couple surface and subsurface ﬂow. The
accuracy of this approximation strongly inﬂuences the conservativenessof the numerical scheme. Similarly, the overland solute transport needs
to be solved simultaneously with water ﬂow in order to preserve themonotonicity of the solution.
One of the most complete furrow irrigation models was developedby Wöhling et al. (2004). The proposed computational framework
iteratively coupled a 1D analytical zero-inertia equation with HYDRUS2D, thus providing a complete description of surface-subsurface waterﬂow along the furrow. The model was further developed and extended
by Wöhling et al. (2006), Wöhling and Mailhol (2007), and Wöhlingand Schmitz (2007). A similar approach was used by Tabuada et al.
(1995), who coupled a model based on a complete hydrodynamic
equation of overland ﬂow with a two-dimensional Richards equation.
While accurately simulating water ﬂow, neither of the models discussed
above considered solute transport in surface and subsurface domains,
thus restricting their applicability to irrigation. Hence, further development of similar approaches that would include a detailed descriptionof solute transport in the root zone is desirable for both scientists andpractitioners (Ebrahimian et al., 2014).

2. Materials and methods
2.1. Modeling approach
The proposed approach combines a one-dimensional description of
coupled water ﬂow and solute transport in an open channel with a twodimensional description of variably-saturated water ﬂow and solutetransport in soil. As in previous studies (e.g., Tabuada et al., 1995;
313

Computers and Electronics in Agriculture 150 (2018) 312–327

G. Brunetti et al.

Fig. 1. A schematic of the modeled transportdomain consisting of the trapezoidal open
channel and soil (Qinj is the solute injection rate
[L3/T], Cinj is the solute concentration [M/L3],
Qir is the irrigation inﬂow rate [L3/T], ET isevapotranspiration, P is precipitation, h is thewater depth in the furrow, hmax is the maximum
allowed water depth in the furrow, bw is the
furrow bottom width, Tw is the top width, Rw is
the ridge width, FS is the furrow spacing, SS is
the side slope, L is the length of the furrow, anddx is the discretization length). The two ﬁgures
at the bottom show boundary conditions forwater ﬂow and solute transport.

Wöhling et al., 2004), a key assumption in the model is that thetransport processes in the soil domain occur only in the vertical plane,
as a function of conditions at a surface ﬂow cross section, and independent of conditions upstream or downstream from that cross section. However, the main novelty of the present study is the development of a complete furrow fertigation model, as well as the inclusion ofa widely used hydrological model, HYDRUS, in the modeling framework.
The furrow of length L is discretized into N elements, each one
characterized by a two-dimensional vertical cross-section (Fig. 1). The
open channel cross-section is assumed trapezoidal, although the numerical approach can be easily adapted to other geometries as well. Theinterface between the surface and subsurface domains acts as a common
boundary.

⎛ Q ⎞
0
⎞
⎛
⎛ A⎞
2
⎛ −i ⎞
U = ⎜ Q ⎟ , F = ⎜ gI1 + Q ⎟ , I =, S = ⎜ g [I2+A (S0−Sf ]⎟ , D
⎜ 0 ⎟
A
⎜
⎟
0
⎝− iC ⎠
⎝ AC ⎠⎠⎝⎝ QC ⎠⎛ 0 ⎞=⎜ 0 ⎟⎜ Kx A ∂C ⎟
∂x ⎠
(2)
⎝
where A is the wetted area [L2], Q is the discharge [L3T−1], C is the
cross-sectional average solute concentration [ML−3], g is the gravitational acceleration [LT−2], S0 is the bed slope [LL−1], Sf is the friction
slope [LL−1], i is the inﬁltration rate [L3T−1], Kx is the dispersion
coeﬃcient [L2T−1], and I1 and I2 are pressure forces. It must be emphasized that the inﬁltration rate i is estimated using the HYDRUS-2D
model.
Since ﬂow velocities in furrow irrigation are typically small, theresulting Froude number is generally less than 0.3. Under such conditions, the inertial terms are negligible and the shallow water equations
can be simpliﬁed to the well-established zero-inertia (ZI) approximation(Strelkoﬀ and Katopodes, 1977):

2.2. One-dimensional overland ﬂow and solute transport
2.2.1. Governing equations
The one-dimensional overland system can be described in the conservative form (Burguete et al., 2009) using the cross-sectional waterand solute mass conservation, momentum balance, and inﬁltration as
follows:

∂U
∂F
∂D
+
=I+S+
∂t
∂x
∂x

∂h
Q |Q| m2W 4 3= S0−∂x
cu2 A10 3(3)

where h is the water depth [L], m is the Manning roughness coeﬃcient[TL−1/3], W is the wetted perimeter [L], and cu is a units coeﬃcient(=1.0 in SI units). The 1D unidirectional discharge Q can be obtained
from Eq. (3), and is expressed as:(1)

where U is the vector of the conserved variables, F is the ﬂux vector, S is
the source term vector, I is the inﬁltration vector, and D is dispersion.
Vectors can be expressed as:314

Computers and Electronics in Agriculture 150 (2018) 312–327

G. Brunetti et al.

⎧Q = sign (J )·k J⎪⎪ J = S0− ∂h
∂x
5
⎨
A3 cu
⎪
k=
2
⎪
mW 3
⎩

discretize the coupled water-solute equations in the conservation form
(Eq. (1). The MUSCL methods switch from the ﬁrst-order scheme inregions with sharp fronts to the higher-order schemes in areas with arelatively smooth solution, thus simultaneously increasing the overall
accuracy and avoiding non-physical oscillations. A ﬂux limiter was used
to constrain the value of the spatial derivative around sharp discontinuities. In the present study, the van Leer (1974) ﬂux limiter wasused in conjunction with the upwind and central diﬀerencing schemes.
The upwind and central ﬂuxes for overland water ﬂow were expressed
as:

(4)

where J is the hydraulic head loss [LL−1], and k is the channel conveyance [L3T−1].
The solute advection-dispersion equation

∂AC
∂QC
+
= −iC
∂t
∂x

(5)

Qj − 1 ⎞⎧F upwind= ⎜⎛⎟⎪j − 12
Qj − 1 Cj − 1
⎝
⎠
⎪
1
⎨
⎛ 2 (Qj − 1 + Qj ) ⎞= ⎜1⎪ F jcentral
1⎟(Q C + Qj Cj )⎪ −2
⎝ 2 j−1 j−1⎠⎩

was simpliﬁed as well by considering only the advective component(Eq. (5), as suggested by Strelkoﬀ et al. (2006). In that study, the resultsof the pure advection equation were compared with the results of thefull ADE for diﬀerent modeling scenarios. The results indicated that theformer formulation was able to provide a suﬃciently accurate description of the solute distribution along the furrow, although the effects of the longitudinal dispersion were appreciable, especially forblocked-end furrows. Neglecting the dispersion component simpliﬁes
the computational scheme and the implementation of boundary conditions.
Considering the above-described simpliﬁcations, Eq. (2) was rewritten as:

⎧
⎪U =
⎨
⎪
⎩

where the 1D cell centered discharge Q was expressed as:⎧ Jj =⎨Q =⎩ j⎟

Q = sign (J )·k J(6)2.2.2. Spatial discretization
Overland water and solute transport were solved numerically using
the Method of Lines (MOL). The MOL approach replaces the spatial
derivatives in the Partial Diﬀerential Equations (PDEs) with an algebraic approximation, which reduces the PDE to an Ordinary Diﬀerential
Equation (ODE). Central to the MOL is the numerical approximation of
spatial derivatives. In the present study, the Finite Volume Method(FVM) was used to spatially discretize the surface domain. The FVM is awell-established and widely used numerical solution method for avariety of engineering problems involving PDEs, mainly because of its
conservativeness. It has been applied successfully in many studies focused on overland ﬂow and solute transport (e.g., Lal, 1998; He et al.,2008). The integral form of Eq. (6) was written as:∂∂t∫Ω UdΩ + ∫Ω ∇·FdΩ = ∫Ω IdΩ

∫Ω UdΩ + ∮ω F·ndω = ∫Ω IdΩ

dt

=−

1
⎡F 1 −F 1 ⎤ + Ij
Δx ⎣ j + 2 j − 2 ⎦−hj + 1 − hj
Δx|Jj| ·sign (Jj )

2

⎜

(11)

⎟

where ϕ is the van Leer ﬂux limiter [LL

(12)
−1

], which was expressed as:

r + |r |⎧ φj (r ) = 1j + |r j|j⎨ lim φj (r ) = 2⎩ r →∞(13)

where r is a function that monitors the gradients of the solution and thatcan be written as:

rj =

(7)

uj−uj − 1
uj + 1−uj(14)

where u is a dependent variable, which in our case can be the wettedarea A for water ﬂow or a cross-sectional average solute concentration
C for solute transport.
Based on Eq. (12), the ﬁrst-order monotonicity preserving upwind
was adopted in regions characterized by steep gradients in the solution,
while the higher-order central diﬀerencing scheme was used in regionswith a relatively smooth solution, thus increasing the stability and accuracy of the spatial scheme. To facilitate the handling of the downstream boundary condition, the upwind scheme was used to computethe numerical ﬂux in the last cell (i.e., j = N). An upwind discretization
was also used for the inﬁltration term.

(8)

The contour integral in Eq. (8) was approximated using numericalﬂuxes F at the edge of each cell j, thus leading to:

d Uj2
kj + 1 + kj

upwind ⎞Fj − . 1 = F upwind+ φ ⎛F central
1 −F j − 1
j− 1
2
2
2 ⎠
⎝ j− 2

Application of the divergence theorem to the second term in Eq. (7)produces:∂∂t

Soj + 1 + Soj

It is worth noting that channel conveyance and bed slope were both
represented with the arithmetic mean of nodal values. Preliminary simulations conﬁrmed that this type of approximation combined with the
MUSCL scheme leads to more accurate results compared to a fully upwind discretization of k and S0, which introduces signiﬁcant numericaldiﬀusion in the solution. The same discretization is used for the solutetransport. Hence, the numerical ﬂux F obtained using the MUSCL
scheme can be written as:( ACA ) , F = ⎛⎝QCQ ⎞⎠ , I = (−−iCi )
⎜

(10)

(9)

One of the main drawbacks of the FVM is the introduction of artiﬁcial diﬀusion in the solution when ﬁrst-order numerical schemes are
adopted (e.g., upwind). Indeed, while preserving the monotonicity ofthe solution (Godunov, 1954), ﬁrst-order numerical schemes havelower accuracy compared to the higher-order schemes (e.g., central
diﬀerence). However, the latter are prone to numerical oscillation (e.g.,van Genuchten, 1976, 1978). Therefore, to increase the numerical accuracy while preserving the stability of the numerical solution, a highresolution scheme was adopted in this study.
In this study, the Monotone Upstream-Centered Scheme forConservation Laws (MUSCL) (van Leer, 1979) was used to spatially2.2.3. Temporal discretization
Once the spatial discretization is deﬁned, time t remains the only
independent variable. Since the inﬁltration component in Eq. (5) is
calculated externally using HYDRUS-2D at each time step, and HYDRUS-2D has its own time-stepping algorithm, the computational eﬃciency of the numerical scheme is strictly related to the frequency ofdata exchanges between the surface and subsurface models, which
dramatically increases when fully implicit time stepping schemes are
adopted. Nevertheless, it must be further emphasized that fully implicit
schemes, while generally increasing the stability of the computational
315

Computers and Electronics in Agriculture 150 (2018) 312–327

G. Brunetti et al.

maximum allowable Δtint calculated by the RKDP algorithm must be:framework, are only ﬁrst-order accurate in time. Thus, due to aforementioned considerations about accuracy and computational eﬃciency, the fourth-order explicit Runge-Kutta Dormand-Prince (RKDP)(Dormand and Prince, 1980) is adopted here to solve the coupled watersolute transport in the trapezoidal channel. More speciﬁcally, an RKDP
algorithm with an adaptive time-stepping strategy based on the truncation error control is used. This kind of temporal discretization automatically adjusts the time step depending on the solution, thus providing a dense temporal output near discontinuities.
It is worth noting that the ZI equation for water ﬂow can be solved
separately from the solute advection equation. However, the solute
advection equation depends on the wetted area and ﬂow velocity,
leading to a one-way coupled system of ODEs. More speciﬁcally, thetemporal variability of the wetted area strongly inﬂuences the solute
advection. Thus, a high temporal resolution scheme for water ﬂow isneeded to avoid overshooting as well as inaccuracies in the solution ofsolute advection.
The basic method will be ﬁrst described for the pure conservation
law, that is, without source terms. These terms will be incorporated
next and required modiﬁcations will be indicated as needed. Eq. (10)can now be written as:

U ij+ 1 = U ij−

Δtint ⎡ i
F 1 −F i 1 ⎤Δx ⎢⎦⎣ j+ 2 j− 2 ⎥

⎛
⎞
Δx
max
⎟
Δtint
< min ⎜ii
⎜ Q + gA ⎟
Γi ⎠
⎝ Ai

where Γ is the top width of the water ﬂow [L]. Furthermore, the conservativeness of the scheme is checked by monitoring the Relative Mass
Balance Error (RMBE) for both water and solute at the end of each time
step, which is calculated as:i − V 0 ) dx − ∫ti q i dt
∫0ti (Qir + Qinj ) dt − ∫0ti F i 1 dt − ∫0L Ai dx − ∫0L (VW⎧W
0 soilN+⎪
2i
RMBEwater
(%) =
·100
ti (Q + Q ) dt
⎪
∫0 ir
⎪
inj
ti (Q C ) dt − ti F i
L
L
t
i
i
0
i
i
i
i
⎨
C
dt − ∫0 A C dx − ∫0 (VS − V S ) dx − ∫0 qsoil c dt
∫0 inj inj∫0⎪
N+ 1 N+ 1
2
2
i
(%) =
·100
⎪ RMBEsolute⎪∫0ti (Qir Cinj ) dt⎩

(18)
where ViW and ViS are, respectively, the water volume and solute mass
per unit length in the soil proﬁle at time i, with V0W and V0S the corresponding quantities at the beginning of the simulation. The last two
terms at the numerators of Eq. (18) indicate the water, qsoil, and solute,
qsoil c, ﬂuxes across the bottom of the soil domain (e.g., deep percolationand solute leaching).

(15)

where Δtint is the internal time step calculated by the adaptive RKDP
algorithm. To increase the overall computational eﬃciency of themodeling framework, a speciﬁc numerical treatment of the inﬁltration
component is adopted here. As discussed above, the RKDP algorithm
uses a ﬁne temporal discretization near discontinuities, implying thatthe calculation of inﬁltration (i.e., HYDRUS executions) is carried out
with very small time steps, which consequently increases the computational cost. However, the inﬁltration vector I of Eq. (3) mathematically represents a source/sink term, which is characterized by a moderately low stiﬀness. Thus, a ﬁrst-order explicit approximation of
inﬁltration (Brufau et al., 2002) is used in our model, which reduces thenumber of HYDRUS-2D executions while maintaining a good accuracy.
This leads to a mixed explicit time scheme, which uses a 4th orderaccurate RKDP time-stepping strategy for the coupled water-solutetransport and a ﬁrst-order explicit scheme for the inﬁltration term.
Practically, this is accomplished by dividing the simulation duration
into equal time steps Δt. At the beginning of each time step, the valuesof the water depth are passed to HYDRUS-2D, which calculates theamount of inﬁltrated water during Δt. Next, the vector I is inserted in
Eq. (1), which is solved by the RKDP method by adapting an internal
time step Δtint depending on the solution. At time t = 0, the surface
domain is solved assuming zero inﬁltration. Eq. (15) can now be rewritten by including the inﬁltration term:iΔt⎛ Ij ⎞
U ij+ 1 = U ij− int ⎡F i 1 −F i 1 ⎤−Δtint ⎜ ⎟Δx ⎢Δt⎣ j+ 2 j− 2 ⎥
⎦

⎝

⎠

(17)

2.2.4. Boundary conditions
The problem statement is completed by the deﬁnition of theboundary conditions (BCs) at the two sides of the numerical domain. In
the present study, the ghost cells methodology is used (LeVeque, 2002).
The 1D ZI equation requires the deﬁnition of two BCs, namely inﬂow
and outﬂow. The speciﬁed inﬂow Q(0,t) determines the numerical ﬂux
entering the domain on the left face of the ﬁrst discretized volume. Twotypes of BCs can be used for the outlet depending on the type of the
furrow end. Since the adopted spatial stencil spans two adjacent cells,
the value of the wetted area A in the right ghost cell has to be adjusted
at each time step to match the desired BC.
If the furrow downstream end is blocked, no outﬂow is allowed until
the water level reaches the maximum furrow depth. If h > hmax
(Fig. 1), the model automatically switches to a free overfall BC. Thus,
the ghost cell value has to be set so that the numerical ﬂux in the lastcell (i.e., j = N) is zero. Considering that the upwind scheme is used forthe rightmost numerical ﬂux, Eq. (12) mathematically reduces to Eq.
(19):5

Fj + 1 =
2

5

3
3
⎞
Aghosthghost −hj ⎞1 ⎛ Aj= 0 if j = N⎜ 2 + 2 ⎟ So−⎛
2mcu ⎜ W 3
Δx ⎠
⎟
3
⎝
Wghost
j
⎠
⎝
⎜

⎟

(19)

This nonlinear equation was solved for hghost using the Brent algorithm (Brent, 1974) at each time step.
The open end boundary condition was implemented by allowing
water to ﬂow over a brink at the end of the furrow. The procedure,
which is followed in this study, is similar to the one reported in Soroush
et al. (2013). Soroush et al. (2013) showed that, independently fromﬂow conditions, the upstream ﬂow passed over the brink with a depth
hb lower than the critical depth hc. This phenomenon was previously
investigated by Beirami et al. (2006), who indicated that for trapezoidal
channels hb ≈ 0.75hc. Hence, in the present study, hghost is determinedat each time step by interpolation so that hb = 0.75hc (Fig. 3).
Both the blocked and open end BCs are applied only when the valueof h in the last cell is greater than 0.001 m. To avoid convergence issues,
if h < 0.001 m, the water depth in the ghost cell, hghost, is set equal tohj=N.
The solute advection equation requires the implementation of theboundary condition only on the inﬂow side of the domain. In theconservation form, this was accomplished by matching the numerical
ﬂux with the solute injection rate Qinj∙Cinj (Fig. 1).

(16)

The temporal discretization is depicted in Fig. 2. The non-iterative
coupling with the external HYDRUS-2D model signiﬁcantly decreases
the computational cost. However, it must be noted that an ad hoc
sensitivity analysis is recommended for the choice of Δt.
Since Eq. (16) is solved explicitly in time, the Courant-FriedrichsLewy (CFL) stability condition must be respected. In particular, theFig. 2. A schematic of the temporal discretization.
316

Computers and Electronics in Agriculture 150 (2018) 312–327

G. Brunetti et al.

Fig. 3. A schematic showing the linear interpolation of the water depth in the ghost cell when the free overfall BC is used.

2.3. Two-dimensional subsurface ﬂow and solute transport( (2.3.1. Governing equations
The HYDRUS (2D/3D) (Šimůnek et al., 2016b) software simulates
water ﬂow and solute transport in variably-saturated porous media. The
software describes water ﬂow using the Richards equation for multidimensional unsaturated ﬂow (Eq. (7):∂θ∂ ⎡ ⎛ ∂ψ⎞⎤=
K Kij
+ Kiz ⎟ ⎥ (i, j = 1, 2)∂t∂yi ⎢ ⎜ ∂yj
⎠⎦
⎣ ⎝

(20)

(21)

2.3.2. Soil hydraulic properties
The van Genuchten–Mualem (VGM) model (van Genuchten, 1980)
was used to describe the soil hydraulic properties in this study:

⎧
⎨
⎩

1
1

(1 + (α | ψ |)n)1 − nθ−θΘ = θ − θr
s
r

(23)

2.3.3. Numerical domain and boundary conditions
The soil computational domain was discretized using triangular 2D
ﬁnite elements. In the example simulations discussed below, no meshstretching was used, and the ﬁnite element (FE) mesh was assumed
isotropic. The quality of the FE mesh was assessed at each time step bymonitoring mass balance errors, which were always below 1% during
numerical simulations, indicating a good spatial discretization.
The boundary conditions used to simulate solute and water inﬁltration are reported in Fig. 1. The ‘Atmospheric’ boundary condition,
which is assigned to the furrow ridge (green1 lines in Fig. 1), can exist in
three diﬀerent states: (a) precipitation and/or potential evaporation
ﬂuxes, (b) a zero pressure head (full saturation) during ponding whenboth inﬁltration and surface runoﬀ occur, and (c) an equilibrium between the soil surface pressure head and the atmospheric water vapor
pressure head when atmospheric evaporative demand cannot be met bywater ﬂuxes towards the soil surface. Due to symmetry, the nodes representing the left and right boundary of the subsurface domain are set
as ‘no ﬂux’ boundaries (grey lines in Fig. 1) because no ﬂow or solutetransport occurred across these boundaries.
A hybrid Dirichlet/Neumann (variable pressure/zero ﬂux) BC was
assigned to the border of the trapezoidal channel (orange lines in Fig. 1)
to simulate variations of the water depth in the furrow. The speciﬁed
value of the pressure head (i.e., the water level) was assigned to thechannel bed and the pressure heads in other channel nodes are determined by the software based on their elevation relative to thechannel bed. A Dirichlet BC was used at nodes with positive pressure
heads and an atmospheric BC was assigned to nodes with negative
pressure heads (i.e., above the water table). A ‘Third Type’ Cauchy BC
was set on top and bottom of the numerical domain to simulate theconcentration ﬂux along the boundaries (red lines in Fig. 1).

where c is the solute concentration [ML−3], qi is the water velocity[LT−1], and Dij are components of the dispersion tensor (L2T−1). Eq.
(10) is valid for non-reactive transport, thus adsorption and precipitation/dissolution of the solutes are currently ignored. Note that although
HYDRUS-2D considers multiple chemical reactions (e.g., sorption anddegradation), this study is limited to non-reactive solutes.
HYDRUS advances the solution in time using a fully implicit, temporal discretization scheme for the Richards equation and the CrankNicholson scheme for the advection-dispersion equation. The program
also uses an adaptive time-stepping strategy to increase the overall
computational eﬃciency.

Θ=))where Θ is the eﬀective saturation [L3L−3], α is a shape parameter
related to the inverse of the air-entry pressure head [L−1], θs and θr are
the saturated and residual water contents, respectively [L3L−3], n is apore-size distribution index [–], Ks is the saturated hydraulic conductivity [LT−1], and l is the tortuosity and pore-connectivity parameter [–].

where θ is the volumetric water content [L3L−3], t is time [T], z is thevertical coordinate [L], yi are the spatial coordinate [L], ψ is the soilpressure head [L], K is the hydraulic conductivity [LT−1], and Kij arecomponents of the hydraulic conductivity anisotropy tensor. Since weassumed isotropic porous media with y1 = y and y2 = z being thetransverse (horizontal) and vertical coordinates, respectively, the conductivity tensor is diagonal with the Kyy and Kzz entries equal to one.
Solute transport is described using the advection dispersion equation:∂θc∂ ⎛∂c ⎞ ∂qiθDij(i, j = 1, 2)
=
∂t
∂yi ⎜
∂yj ⎟ ∂yi
⎝
⎠

1 2

1− n
⎧
⎤ if ψ < 0⎪ Ks Θl ⎡ 1− 1−Θ n −n 1
⎥
⎢
K=
⎦
⎣
⎨
⎪
Ks if ψ > 0⎩

if ψ ⩽ 0

1 if ψ > 0
1
For interpretation of color in Fig. 1, the reader is referred to the web version of thisarticle.

(22)317

Computers and Electronics in Agriculture 150 (2018) 312–327

G. Brunetti et al.

concentrations are stored for the next time step. The coupled ODEs
are solved simultaneously using the RKDP algorithm, and newwetted area, water level, and solute concentration vectors are updated in their respective matrices.

2.4. Models coupling strategy
2.4.1. Wet/Dry boundary tracking
In all numerical simulations, the water front is not considered as amoving boundary, and calculations are carried out in both wet and drycells. The numerical treatment of wet/dry cells is similar to Bradford
and Sanders (2002), Bradford and Katopodes (2001) and Brufau et al.
(2002). In order to avoid numerical problems associated with extremely
low values of the water depth at the advancing front, a threshold value
ε is used to identify wet cells. In the present study, water levels above
0.0001 m indicate wet cells (Bradford and Katopodes, 2001). If hj > ε,the water depth value is passed to HYDRUS, which calculates inﬁltration that is then used to solve surface ﬂow. Two computational procedures can be followed when the cell is dry:

2.5. Theoretical validation
2.5.1. Mesh sensitivity analysis
Schlesinger et al. (1979) deﬁned model validation as a “substantiation that a computerized model within its domain of applicability possesses asatisfactory range of accuracy consistent with the intended application of themodel.” There exists a variety of validation techniques that can be used
to assess the model accuracy and robustness (Sargent, 1998). Even so, asensitivity analysis on internal model parameters (e.g., mesh size, timestep) represents the ﬁrst step in a validation framework, sometimes
referred to as model veriﬁcation, that is usually accomplished before
using the model itself.
The computational mesh size is an important component of Eulerianmodeling frameworks, especially when the FVM is adopted, since it isdirectly related to false numerical diﬀusion and accuracy. In regions of
the computational domain where the dependent variables exhibit sharp
gradients, a ﬁne mesh is needed to avoid smearing the solution.
However, as the mesh resolution increases, so does the computationalcost. The mesh resolution should, thus, be properly designed to simultaneously minimize the execution time and maximize the numericalaccuracy. To this end, a mesh sensitivity analysis represents an important prognostic tool.
The inﬂuence of the spatial discretization N on both the computational cost and the accuracy of the proposed hybrid FE-FV furrow fertigation model is investigated ﬁrst. The eﬀect of three diﬀerent meshsizes (i.e., dx = 2, 5, 10 m) on the simulated water and solute ﬂow in a
100 m long blocked-end furrow (a blocked end scenario in Table 1) is
examined. For each numerical simulation, the execution time provides
the measure of the computational eﬀort. A time step dt = 20 s was assumed.

1. Subsurface in steady-state conditions: HYDRUS is completely bypassed and inﬁltration is set to 0. For example, during the advance,
the last part of the furrow is not simulated, since it is ﬂagged as“dry” and the subsurface ﬂow is negligible. It is worth noting that
this option, while increasing the computational eﬃciency of thenumerical scheme, can be considered only when the initial condition is in hydrostatic equilibrium and boundary ﬂuxes are negligible.
2. Subsurface in dynamic conditions: These computations apply during
recession. Although no more water inﬁltrates, HYDRUS can continue to simulate variably-saturated water ﬂow and solute transport
in the soil. In such circumstances, a negative pressure head is assigned to the top BC in HYDRUS, which automatically switches to“Atmospheric BC.” HYDRUS is then executed, and the ﬁnal soil
condition is stored and used in the next step.
2.4.2. Model implementation
The one-dimensional horizontal domain representing the furrow
and the two-dimensional cross-sectional domains representing the
subsurface (Fig. 1) were coupled using a Python script. The script simultaneously solves overland ﬂow and solute transport in the furrow,interacts with HYDRUS-2D, and exchanges data between two models. A
series of user-deﬁned functions were developed to write and read the
input/output ﬁles generated by HYDRUS-2D, which was directly executed from Python, bypassing the HYDRUS (2D/3D) graphical userinterface. The coupling strategy is summarized below:

2.5.2. Time step sensitivity analysis
A core assumption of the proposed modeling framework is the explicit linearization of the inﬁltration vector in Eq. (7). As mentioned
above, this assumption eliminates the need for an iterative execution of
HYDRUS-2D, thus reducing the overall computational cost. However,
the accuracy of the forward Euler approximation of the inﬁltration
component is supposed to deteriorate when the time step dt increases.
The loss of accuracy is directly related to the stiﬀness of the inﬁltration
term, which depends on the scenario analyzed.
Therefore, the inﬂuence of the temporal discretization dt was investigated using three diﬀerent time steps (dt = 10, 20, 30 s). The
blocked-end furrow described in Table 1 was used as a test case. For
each numerical simulation, typical hydraulic information was calculated (i.e., front advance, proﬁles, etc.) and the execution time wasanalyzed. A mesh size dx = 5 m was assumed. It must be emphasized
that, in the developed computational framework, the time step dt inﬂuences not only the computation of inﬁltration, but also the exchange
of information between the surface and subsurface models. On the other
hand, the time step does not inﬂuence the numerical solution of thezero-inertia and advection equations, which are solved using the RKDP
with an adaptive time stepping strategy.

1. The geometric characteristics of the furrow, soil hydraulic properties, initial conditions and other input data are deﬁned in both
Python and HYDRUS-2D. The HYDRUS model is then initialized.
2. The horizontal domain is discretized into N homogenous Finite
Volumes. There is a HYDRUS-2D cross section for each FV. Pressure
head, water content, and solute concentration distributions in thesoil are stored in three matrices, which are overwritten at each time
step. Water content is stored only for visualization purposes and notused directly in the coupling procedure. Similarly, vectors containing wetted areas, water depths, and solute concentrations in thesurface domain are updated at each time step in three additional
matrices, in addition to two columns to handle ghost cells.
3. As described above, the MOL is used to reduce the PDEs to a set of
ODEs using the FVM. Boundary conditions are applied using ghost
cells and the time step dt is set.
4. A for loop is used to iterate through time. At each time step, vectors
h and C containing previously calculated values of water levels andsolute concentrations in the furrow, respectively, are passed to
HYDRUS-2D to setup its boundary conditions. Soil pressure headand solute concentration distributions calculated at the previous
time step are used as initial conditions in HYDRUS for the current
time step.
5. Another for loop is used to iterate through space and calculate theinﬁltration vector. Calculated soil pressure heads and solute

2.5.3. Description of validation scenarios
Preliminary sensitivity analyses help identify the most eﬃcientspatial and temporal discretization. The next step is to properly validate
the model theoretically and experimentally, which involves the assessment of the model structure and its capacity to accurately describe
the investigated system (Sargent, 1998). The former is usually accomplished by comparing the results of the proposed model with the resultsof a well-established, previously validated model using diﬀerent
318Computers and Electronics in Agriculture 150 (2018) 312–327

G. Brunetti et al.

Table 1
Geometric characteristics, soil hydraulic properties, and other input parameters
used in theoretical validation.
Blocked end
Furrow length, L (m)Bottom width, bw (cm)Side slope, SS (–)
Maximum depth, hmax (cm)
Bed slope, S0 (–)
Manning’s coeﬃcient, m (s/m1/3)
Residual water content, θr (m3/m3)
Saturated water content, θs (m3/m3)
Retention function shape parameter, α (1/cm)
Retention function shape parameter, n (–)
Saturated hydraulic conductivity, Ks (cm/min)
Tortuosity, l (–)
Initial water content, θ0 (m3/m3)
Wetting front pressure head, ψf (cm) (Warrick et al.,2007)
Inﬂow rate, Qir (l/s)
Solute injection rate, Qinj (l/s)
Inﬂow concentration, Cinj (g/l)
Longitudinal dispersivity (cm), DL
Transverse dispersivity (cm), DT
Irrigation cutoﬀ time, tw (min)Solute cutoﬀ time, ts (min)

100
16
1.2
15
0.0005
0.04
0.065
0.41
0.075
1.89
0.074
0.5
0.085
−5.4
2
0.014
40
0.5
0.1
30
30

Table 2
Geometric characteristics, soil hydraulic properties, and other input parameters
used in experimental validation.

Open end

Holtville
Furrow length, L (m)Bottom width, bw (cm)Top width, Tw (cm)
Ridge width, Rw (cm)Furrow spacing, FS (cm)Side slope, SS (–)
Maximum depth, hmax (cm)
Bed slope, S0 (–)
Manning’s coeﬃcient, m (s/m1/3)
Initial water content, θ0 (–)
Inﬂow rate, Qir (l/s)
Solute injection rate, Qinj (l/s)
Inﬂow concentration, Cinj (g/l)
Longitudinal dispersivity (cm), DL
Transverse dispersivity (cm), DT
Irrigation cutoﬀ time, tw (min)Solute cutoﬀ time, ts (min)0.001

Table 3The VGM hydraulic parameters used in experimental validation.

20

synthetic scenarios, thus avoiding the interference of diﬀerent sourcesof uncertainty typical of experimental data.
In this study, the proposed hybrid FE-FV coupled water and solutetransport model was theoretically validated against the solute transportcomponent of the WinSRFR software package (Bautista et al., 2009).
This software will be released to the general public with WinSRFR,
Version 5.1. The open-channel ﬂow simulator in WinSRFR (Bautista
et al., 2009), SRFR, utilizes a 1D zero-inertia approximation of thehydrodynamic wave model, which is coupled with the solution of the
1D advection-dispersion equation for solute transport (Perea et al.,2010). In the present study, the dispersion component was neglected
and only advection was simulated using WinSRFR. The explicit inﬁltration function of Warrick et al. (2007), modiﬁed by Bautista et al.
(2014) and Bautista et al. (2016) was used to calculate the ﬂow-depth
dependent furrow inﬁltration. This function has been shown to approximate the solution to the two-dimensional Richards equation reasonably well. The pressure head at the wetting front ψf [L] is calculatedusing the unsaturated hydraulic conductivity function and can be expressed as:

ψf =

∫ψ

0 K (ψ)0Ksdψ

182.88
16
52
55
107
1.2
15
−0.0007 to 0.0007 (Variable)
0.04
0.35
3.2
0.014
40
0.5
0.1
44
37.5

θr (cm3/cm3)

θs (cm3/cm3)

α (1/cm)n (–)

Ks (cm/min)

l (–)

0.1

0.5

0.017

1.26

0.017

0.5

FE model under diﬀerent conditions. The ﬁrst scenario consists of a
blocked-end furrow, with a small slope and equal water and solute
cutoﬀ times. The second scenario considers an open-end furrow, characterized by a moderate slope and diﬀerent water and solute cutoﬀ
times. In both cases, the soil is assumed to be a sandy loam. The VGM
parameters used in HYDRUS-2D are taken from Carsel and Parrish(1988). The subsurface vertical domain is discretized into 2043 triangular elements and 1076 nodes. The eﬀect of groundwater is not simulated in this study and thus a ‘Free Drainage’ BC is assigned to thebottom nodes (z = −80 cm) (blue line in Fig. 1). The geometric characteristics, soil hydraulic properties, and other input parameters arereported in Table 1.

2.6. Experimental validation
Theoretical validation provides a ﬁrst assessment of the model accuracy and robustness. However, the main objective of a numericalmodel is to accurately reproduce the behavior of the system under investigation. In this view, experimental validation, which comparesmodel results and measured data, plays a fundamental role.

(24)

Two synthetic scenarios were developed to validate the hybrid FV-

Fig. 4. Geometric characteristics of the furrow located in Holtville (CA) that were used for experimental validation.
319

Computers and Electronics in Agriculture 150 (2018) 312–327

G. Brunetti et al.

Fig. 5. Simulation results for diﬀerent mesh sizes dx: Water depth proﬁles at time t = 0.17 h (top left); surface solute proﬁles at time t = 0.17 h (bottom left);
execution time (top right), and advance time with distance (bottom right).

calcareous, hyperthermic Typic Torriﬂuvents, [Hendricks, 1985]) as
described by the USDA-NRCS Web Soil Survey. This particle size distribution was provided to the Rosetta pedotransfer module (Schaap
et al., 2001) of HYDRUS-2D to estimate a set of soil hydraulic parameters. However, these parameters produced less than satisfactory irrigation simulation results. A better set of parameters was found using a
diﬀerent soil found in the area, Holtville clay (clayey over loamy,
smectitic over mixed, superactive, calcareous, hyperthermic Typic
Torriﬂuvents [Hendricks, 1985]). For this soil (sand = 12%,silt = 32%, clay = 56%), Rosetta estimated a higher value of the saturated water content, which was more consistent with the measured
pre-irrigation water contents, and also a higher value of the saturated
hydraulic conductivity, consistent with the value reported by the Web
Soil Survey. These parameters are shown in Table 3. The residual watercontent, α, and n computed for the Kofa clay were not very diﬀerent
from those shown in the table.
The 182.88 m long furrow had a blocked downstream end and on
average zero slope, which increases the numerical complexity compared to a constant slope scenario. In particular, the slope varies between −0.0007 and 0.0007 (–) (Fig. 4). This is one of the reasons why
this particular data set was chosen for the experimental model validation.

In the present study, the proposed hybrid FV-FE fertigation model isfurther experimentally validated against ﬁeld data measured in 2001 at
the Yuma Valley Agricultural Center, Yuma AZ, USA (C. Sanchez, 2018,
personal communication).
2.6.1. Site description
A key diﬃculty with this part of the analysis was determining appropriate soil hydraulic parameters. These parameters are diﬃcult tomeasure both in the laboratory and in the ﬁeld. Previous studies haveshown that soil hydraulic parameters under surface irrigation conditions exhibit poor correlation with soil texture and that pedotransfer
functions have wide margins of errors when estimating these parameters (Selle et al., 2011; Zapata and Playán, 2000). Therefore, Zapata
and Playán (2000) recommended ﬁtting irrigation models to irrigation
data as a mechanism for determining these parameters.
The experimental data were collected with the purpose of conducting fertigation simulation studies, but were not intended to be usedin combination with physically based inﬁltration models. As a result,the data set includes only pre-irrigation gravimetric water contents andan average soil texture (33% sand, 32% silt, and 35% clay). The soil
particle size distribution is consistent with one of the soils found in thearea, Kofa clay loam (Clayey over sandy, smectitic over mixed,
320Computers and Electronics in Agriculture 150 (2018) 312–327

G. Brunetti et al.

Fig. 6. Simulation results for diﬀerent time steps dt: Water depth proﬁles at time t = 0.17 h (top left); surface solute proﬁles at time t = 0.17 h (bottom left);
execution time (top right); advance time with distance (bottom right).

The subsurface vertical domain was discretized into 2043 triangularelements and 1076 nodes. The eﬀect of groundwater was not simulated
in this study, and thus a ‘Free Drainage’ BC was again assigned to thebottom nodes (z = −80 cm) (blue line in Fig. 1). The initial volumetricwater content was assumed to be constant in the entire domain and set
equal to 0.30 cm3/cm3. Vertical domains were assumed solute free atthe beginning of the numerical simulation.

A constant inﬂow rate of about 3.2 l/s for 44 min was used to irrigate the furrow, which was simultaneously fertigated with a bromide
tracer. The solute was diluted in a separate tank, resulting in a concentration of approximately 40 g/l, and injected in the furrow with aconstant ﬂow rate of approximately 0.014 l/s. The solute tank emptied
in about 37.5 min. Five measuring stations (grey dots in Fig. 4) wereused to monitor the water level and solute distribution along the furrow
length. Water levels and concentrations were measured with an acquisition frequency of about 5 min. Soil water contents were measured
gravimetrically before the irrigation at ﬁve measuring stations. The
average volumetric water content was 0.35 cm3/cm3. A Manning’s
coeﬃcient of 0.04 s/m1/3 was used as suggested by USDA-NRCS (2012).
Subsequent comparisons of measured with simulated ﬂow depth hydrographs show that this is a reasonable roughness value. The geometric characteristics, soil hydraulic properties, and other input parameters are summarized in Table 2.

3. Results and discussion
3.1. Theoretical validation
3.1.1. Mesh sensitivity analysis
Fig. 5 shows water (top left plot) and solute (bottom left plot)
proﬁles in the furrow at time t = 0.17 h (about 10 min) simulated with
diﬀerent mesh sizes dx. The advance curves (bottom right plot) andassociated computational time (top right plot) are also shown in Fig. 5.
At the ﬁrst inspection, it is evident that the mesh size dx does not
dramatically aﬀect the numerical accuracy of the proposed model,
which exhibits no overshooting and only limited numerical diﬀusion
even for a coarse spatial discretization (i.e., dx = 10 m). Water andsolute proﬁles for diﬀerent dx are relatively close to each other. This is2.6.2. Numerical domain and boundary conditions
Based on the ﬁndings of theoretical validation, the horizontal domain was discretized into 40 Finite Volumes (i.e., dx = 4.57 m) and thetime step dt was set to 20 s. These settings guarantee a good balance
between numerical accuracy and computational cost.
321Computers and Electronics in Agriculture 150 (2018) 312–327

G. Brunetti et al.

Fig. 7. Simulated water (left) and solute (middle) hydrographs, and advance curves (right) calculated using WinSRFR (dashed lines) and hybrid FV-FE (solid lines)
models for the open end (top plots) and blocked end (bottom plots) scenarios.

in the model execution time.

also true for the simulated advance curves, which deviate only at earlytimes. However, as expected, false numerical diﬀusion increases withthe mesh size. For dx = 10 m, the water proﬁle was slightly smeared,
and a faster waterfront advance was predicted for the ﬁrst 60 m of the
furrow length, as indicated by the red line in the bottom right plot ofFig. 5. This eﬀect is more pronounced for the simulated solute proﬁle,
which clearly shows how the waterfront is ahead compared to what is
calculated with dx = 2 or 5 m. This behavior is typical of FV models
when the spatial discretization is coarse. False numerical diﬀusion
tends to vanish when the mesh is reﬁned, as demonstrated by the waterand solute proﬁles as well as the advance curves for dx = 2 and 5 m,which almost overlap. This indicates that further mesh reﬁnement will
not increase the overall accuracy of the numerical solution, only thecomputational eﬀort. Interestingly, the computational time increases
almost linearly with dx. In particular, the execution time for dx = 2 m
was 5.5 times greater than for dx = 10 m, while a mesh size dx = 5 m
reduced the execution time by almost one third in comparisondx = 2 m. As a result of this analysis, a number of conclusions can be
drawn:

3.1.2. Time step sensitivity analysis
The results of the time step sensitivity analysis are shown in Fig. 6.
The graph shows the simulated water (top left plot) and solute (bottom
left plot) proﬁles at time t = 0.17 h (about 10 min) for diﬀerent time
steps dt. The advance curves (bottom right plot) and associated computational cost (top right plot) are also shown in Fig. 6.
It is evident that, compared to the results obtained for diﬀerent
mesh sizes, the model exhibits a lower sensitivity to the time step dt.
Very similar water and solute proﬁles were computed for all tested time
steps (Fig. 6). Likewise, the computed advance curves were in closer
agreement, especially at later times. Furthermore, the execution timeincreased linearly with the time step. It is worth noting that the time
step dt mainly inﬂuences the computation of the inﬁltration vector I,
which is temporally discretized using the forward Euler approximation.
More speciﬁcally, it is assumed that the nodal inﬁltration rate is constant during dt, which is an approximation wherein the accuracy depends on the variability of the inﬁltration rate during a particular time
step. In this regard, the Horton’s inﬁltration theory (Horton, 1933)
states that the inﬁltration rate into an initially dry soil typically decreases with time until reaching a steady-state value. Thus, when thesoil is dry and the water level is small (i.e., the beginning of furrow
irrigation), the inﬁltration rate declines rapidly and a ﬁne time step isneeded to accurately track inﬁltration. Thus, the stiﬀness of the inﬁltration term is appreciable at the beginning of the simulation when
the furrow is dry, and the water level is relatively low. Under suchcircumstances, the vector I makes Eq. (1) moderately stiﬀ. However,this eﬀect is limited only to the ﬁrst part of the simulation and tends tovanish when the water level rises and the soil saturates, as shown in
Fig. 6. This analysis makes clear that the explicit approximation of theinﬁltration term is acceptable, provided that the time step dt is suﬃciently small. Rather than using a constant time step, an adaptive time
stepping discretization of the inﬁltration vector would likely increase
the computational eﬃciency of the model while maintaining a good
accuracy.

- The proposed hybrid FV-FE model maintains a good accuracy evenwhen a coarse mesh is adopted, which demonstrates an overall robustness of the model. This is of practical importance since it indicates that a model with coarse spatial discretization can be conﬁdently used for preliminary design purposes or intensive numericalanalyses, which require running the model multiple times. A classic
example would be the application of a model with coarse discretization as a lower ﬁdelity surrogate model in the Bayesian optimization framework (Razavi et al., 2012).
- The eﬀect of false numerical diﬀusion, although limited, tends tovanish when the mesh is reﬁned. Therefore, when the analysis requires a higher accuracy, a reasonable mesh reﬁnement will increase the overall accuracy of the model.
- A mesh size dx = 5 m oﬀers a good trade-oﬀ between numericalaccuracy and computational cost. Further mesh reﬁnements (i.e.,dx = 2 m) lead to minimal accuracy gains and a signiﬁcant increase322Computers and Electronics in Agriculture 150 (2018) 312–327

G. Brunetti et al.

Fig. 8. Modeled (red lines) and measured (grey dots) hydrographs (plots in the third column) and concentration breakthrough curves (plots in fourth column), andagainst each other (plots in the ﬁrst and second column) at diﬀerent furrow locations (rows from top down for x = 4.57, 45.72, 91.44, 137.16, and 180.74 m) for the
Holtville validation scenario. The dashed black lines indicate the regression lines. (For interpretation of the references to colour in this ﬁgure legend, the reader isreferred to the web version of this article.)(1) would simultaneously increase the numerical accuracy andcomputational eﬃciency. This strategy should consist of small time
steps in regions characterized by high inﬁltration gradients (i.e., drysoils) and larger time steps at later times when inﬁltration stabilizes(i.e., wet soils). Alternatively, future developments could also include the use of a correction factor that initially adjusts the inﬁltration rate and vanishes with time. This could be practically
implemented in the present modeling framework using a scaling
factor for the saturated hydraulic conductivity on the border of the
furrow channel.
- A time step dt = 20 s represents a good trade-oﬀ between the numerical accuracy and computational cost.

The time step sensitivity analysis reveals that the computational
cost is more inﬂuenced by the mesh size dx than by the time step dt.
Indeed, the execution time only triples when dt = 30 s compared towhen dt = 10 s. Overall, a time step dt = 20 s provides a good balance
between the accuracy and computational cost. As a result of this analysis, a number of conclusions can be drawn:- The time step dt has a small inﬂuence on the proposed hybrid FV-FE
model in the analyzed range of time steps, reﬂecting a good robustness of the model. However, a time step sensitivity analysis is
recommended to investigate the stiﬀness of the inﬁltration term,
which increases for dry conditions and highly permeable soils.
- An adaptive temporal discretization of the inﬁltration term in Eq.

323

Computers and Electronics in Agriculture 150 (2018) 312–327

G. Brunetti et al.

Fig. 9. Modeled (solid lines) and measured (dots) water proﬁles (top plot) and advance curves (bottom plot) for the Holtville validation scenario.

make the model open source and available to users in order to receive
continuous feedback and validation.

3.1.3. Validation scenarios
Based on the results of the sensitivity analyses, the mesh size dx andthe time step dt are set to 5 m and 20 s, respectively. Results of thetheoretical validation are reported in Fig. 7. In particular, Fig. 7 shows acomparison between simulated advance curves, and water and solute
proﬁles calculated by WinSRFR (dashed lines) and by the hybrid FV-FE
model (solid lines) for the open end (top plots) and blocked end (bottom
plots) scenarios. At the ﬁrst inspection, it is evident that the results
calculated by the proposed model matched very closely those computedwith WinSRFR for both scenarios. More speciﬁcally, the advance curves
nearly overlapped (right plots in Fig. 7) as did the rising branches of thewater proﬁles (left plots in Fig. 7). However, the hydrographs deviatedslightly from each other during the storage and recession phase. Thiseﬀect is more evident for the blocked-end scenario, for which WinSRFR
predicts a faster recession compared to the hybrid FE-FV model. Thisbehavior is likely explained by the adoption of the explicit inﬁltration
function of Warrick et al. (2007) in WinSRFR. When Warrick et al.
(2007) ﬁrst presented their explicit formulation, they compared thecalculated cumulative inﬁltration for a rectangular channel against
HYDRUS-2D. In particular, two soils were analyzed in that study: loam
and sandy loam. In both cases, the explicit function showed a slightoverestimation of cumulative inﬁltration, which is the same behavior
encountered in the present study. Furthermore, Warrick et al. (2007)
demonstrated that the edge eﬀect increases with the water level andthat it is more pronounced for parabolic and triangular channels thanfor rectangular channels. Thus, the overestimation of inﬁltration is
expected to be substantial for trapezoidal channels, like the one analyzed in the present theoretical validation.
Overall, the theoretical validation suggests that the proposed hybrid
FE-FV model accurately describes the advance, storage, and depletion
phase, as well as the solute transport during furrow irrigation andfertigation. Furthermore, the analysis reveals that the model is oscillation-free in the analyzed cases and stable for a variety of boundary
conditions, thus it can be conﬁdently used for analyzing real systems.
Nevertheless, it must be emphasized that the analyzed scenarios are notexhaustive of all practical modeling situations, therefore we plan to

3.2. Experimental validation
Figs. 8 and 9 compare modeling results (solid lines) with experimental results (circles). The ﬁrst ﬁgure (Fig. 8) depicts hydrographs(left graphs) and concentration breakthrough curves (right graphs),
while Fig. 9 displays water depth proﬁles (upper graphs) and advancecurves (bottom graphs). All ﬁeld measurements were reasonably well
predicted by the model. Diﬀerences between water levels at x = 4.57 m(Fig. 8) can be attributed to furrow survey errors. Nevertheless, thegeneral trend highlights good model performances. This is conﬁrmed by
the model’s ability to accurately reproduce measured concentration
breakthrough curves. Although the solute concentration data are noisy,
both the time arrival of the solute plume and the peak concentrations
are predicted accurately. The solute dilution at x = 4.57 and 45.72 m
induced by solute cutoﬀ is well simulated by the model, thus indicatinga good reliability of the proposed numerical approach.
The advance curve (Fig. 9) is generally reproduced with high accuracy, except at x = 137.16 m where the model predicts a delay in theadvance of the waterfront (Fig. 9). This is conﬁrmed by a slight underestimation (i.e., ≈1 cm) of the water depth at x = 137.16 m (topplot in Fig. 9), which, however, tends to disappear with time. In our
opinion, these small deviations are mainly due to slope uncertainties
and the spatial variability of soil hydraulic properties. Water proﬁles(top plot in Fig. 9) are very well approximated, except for a small underestimation in the ﬁrst part of the furrow. Interestingly, this behavior
is not systematic but limited to particular simulation times (i.e.,t = 42–82 min), likely reﬂecting small measurement errors.
It must be emphasized that soil hydraulic parameters used in thesimulation were generated using the ROSETTA pedotransfer model,
predictions of which for the saturated values can be aﬀected by highuncertainties (Rubio et al., 2008). The bias associated with the estimation of Ks can, therefore, partially explain the model deviations.
Better estimations of soil hydraulic parameters could be obtained by324

Computers and Electronics in Agriculture 150 (2018) 312–327

G. Brunetti et al.

uniform distribution of water in a furrow cross section and instead
assume a uniform water distribution. Clearly, this distribution needs tobe accounted for in order to properly evaluate the irrigation performance (e.g., distribution uniformity, application eﬃciency, requirement eﬃciency). On the other hand, the moisture distribution along the
furrow was homogeneous. Such a uniform longitudinal distribution canbe expected with zero-slope irrigation systems, as diﬀerences in opportunity time along the ﬁeld are small.
On the other hand, a very diﬀerent situation occurs for fertigation.
The solute distribution diﬀers from the water distribution. Solute
mostly accumulates between x = 91.44 m and x = 137.16 m, while it isunevenly distributed in remaining parts of the furrow. This behavior isparticularly exacerbated for x < 91.44 m, where recession times occur
earlier than in the remainder of the furrow. In such circumstances,
plants located in the ﬁrst part of the furrow will receive less fertilizer,
leading to an unbalanced crop production. The analysis reveals that theunevenness of the furrow bottom slope negatively aﬀects the uniformityof solute distribution. However, it must be emphasized that the available dataset did not include measurements of water contents and soluteconcentrations in the soil proﬁle and thus it was not possible to validate
simulated soil water contents and solute concentrations in the soil.
Results demonstrate the usefulness of the proposed modeling approach and its ability to account for the solute and moisture distributions along the entire furrow. In the future, these results could be interestingly used as initial conditions for daily or weekly numericalsimulations with HYDRUS-2D to identify crop stresses and to optimize
the irrigation schedule.
4. Conclusions and summary
The main objective of this study was to develop a numerical model
able to provide an accurate mechanistic description of the coupledsurface-subsurface processes happening during furrow irrigation andfertigation. The proposed modeling framework combines a one-dimensional description of water ﬂow and solute transport in the surface
domain, based on the Zero-Inertia approximation of the hydrodynamic
wave and the advection equation for solute transport, respectively, witha two-dimensional description of water ﬂow and solute transport in thesubsurface domain, based on the Richards and advection-dispersion
equations, respectively. One of the main novelties of this study is theimplementation in the coupled model of the widely used Finite Element
model, HYDRUS, which can describe the simultaneous movement ofwater, heat, and solutes in porous media and which can provide a basis
for analyzing many diﬀerent scenarios and conditions.
The surface domain is solved using the Method of Lines in conjunction with a high-order MUSCL Finite Volume scheme combinedwith the well-established Dormand-Prince (RKDP) temporal discretization method, thus leading to a hybrid FV-FE furrow irrigationand fertigation model. The mesh sensitivity analysis is used to examine
the eﬀects of the coarse spatial discretization on the results of themodel. The analysis reveals that the model guarantees a suﬃcient accuracy even for larger mesh sizes, thus suggesting its application as alower ﬁdelity surrogate in computationally intensive statistical analyses(e.g., Brunetti et al., 2017a, 2017b).
To increase the computational eﬃciency of the proposed numericalframework, the inﬁltration component is explicitly linearized in time,thus avoiding the iterative execution of HYDRUS at each time step. The
validity of this assumption is conﬁrmed by a preliminary time sensitivity analysis, which conﬁrms the robustness of the developed modeland clariﬁes how the stiﬀness of the inﬁltration term decreases withtime. In this light, future work should investigate the use of adaptive
time stepping strategies for the evaluation of inﬁltration.
Theoretical and experimental validations demonstrate the accuracyof the proposed hybrid FV-FV furrow irrigation and fertigation model,
which is able to provide a mechanistic description of the water andsolute distribution in the soil along the furrow. Hence, future

Fig. 10. Simulated distributions of soil water contents (left plots) and soluteconcentrations (right plots) in the vertical 2D domains at diﬀerent measuring
stations.

comparing simulated and measured subsurface quantities (e.g., watercontents or pressure heads). Thus, a complete dataset, which includes
information about both surface and subsurface domains, is recommended for a comprehensive model calibration and an assessmentof its uncertainty. Further plausible sources of uncertainties are related
to the assumption of constant initial conditions and to the adopted
value of the Manning’s coeﬃcient. Indeed, the bias in the simulated
water level at x = 137.16 m could be potentially explained by the underestimation of the near-surface water content at the beginning of theinﬁltration. Similarly, the delay in the simulated advance curve could
be attributed to the assumption of the roughness homogeneity along the
furrow.
One of the key advantages of the proposed hybrid FE-FV model is its
ability to provide a mechanistic description of the transport processes inthe subsurface domain. Fig. 10 shows the simulated post-irrigation(t = 7.7 h) distribution of water contents (left plots) and solute concentrations (right plots) in the soil at diﬀerent measuring stations.
Diﬀerences in the water content are mostly explained by slight variations in ﬁeld elevations – high points in the ﬁeld inﬁltrated less water
than low points. Note that the wetting bulbs of neighboring furrows did
not merge and that a considerable amount of water was stored below
the furrow bed. Conventional furrow irrigation models ignore the non325

Computers and Electronics in Agriculture 150 (2018) 312–327

G. Brunetti et al.

applications and developments should focus on the use and testing of
the proposed model for the numerical investigation of the root soluteand water uptake of typical crops in long-term simulation scenarios.

Horst, M.G., Shamutalov, S.S., Pereira, L.S., Gonçalves, J.M., 2005. Field assessment of
the water saving potential with furrow irrigation in Fergana, Aral Sea basin. Agric.
Water Manage. 77, 210–231. http://dx.doi.org/10.1016/j.agwat.2004.09.041.
Horton, R.E., 1933. The role of inﬁltration in the hydrologic cycle. Trans. Am. Geophys.
Union 445–460. http://dx.doi.org/10.1029/TR014i001p00446.
Katopodes, N.D., Strelkoﬀ, T.S., 1977. Dimensionless solutions of border-irrigation advance. J. Irrig. Drain. Eng. 103, 401–417.
Lal, A.M.W., 1998. Weighted implicit ﬁnite-volume model for overland ﬂow. J. Hydraul.
Eng. 124, 941–950. http://dx.doi.org/10.1061/(ASCE)0733-9429(1998) 124:9(941).
Lazarovitch, N., Warrick, A.W., Furman, A., Zerihun, D., 2009. Subsurface water distribution from furrows described by moment analyses. J. Irrig. Drain. Eng. 135, 7–12.
http://dx.doi.org/10.1061/(ASCE)0733-9437(2009) 135:1(7).
LeVeque, R.J., 2002. Finite Volume Methods for Hyperbolic Problems, vol. 54. Cambridge
Univ. Press, p. 258. http://doi.org/10.1017/CBO9780511791253.
Oweis, T.Y., Walker, W.R., 1990. Zero-inertia model for surge ﬂow furrow irrigation.
Irrig. Sci. 11, 131–136. http://dx.doi.org/10.1007/BF00189449.
Perea, H., Strelkoﬀ, T.S., Adamsen, F.J., Hunsaker, D.J., Clemmens, A.J., 2010.
Nonuniform and unsteady solute transport in furrow irrigation. I: model development. J. Irrig. Drain. Eng. 136, 365–375. http://dx.doi.org/10.1061/(ASCE)IR.19434774.0000106.
Razavi, S., Tolson, B.A., Burn, D.H., 2012. Review of surrogate modeling in water resources. Water Resour. Res. 48, n/a-n/a. http://doi.org/10.1029/2011WR011527.
Rubio, C.M., Llorens, P., Gallart, F., 2008. Uncertainty and eﬃciency of pedotransfer
functions for estimating water retention characteristics of soils. Eur. J. Soil Sci. 59,
339–347. http://dx.doi.org/10.1111/j.1365-2389.2007.01002.x.
Sargent, R.G., 1998. Veriﬁcation and validation of simulation models. Simul. Conf.
Proceedings, 1998. Winter 1, 121–130. http://doi.org/10.1109/WSC.1998.744907.
Schaap, M.G., Leij, F.J., Van Genuchten, M.T., 2001. Rosetta: a computer program forestimating soil hydraulic parameters with hierarchical pedotransfer functions. J.
Hydrol. 251, 163–176. http://dx.doi.org/10.1016/S0022-1694(01)00466-8.
Schlesinger, Crosby, R., Cagne, R., Innis, G., Lalwani, C., Loch, J., Sylvester, R., Wright,
R., Kheir, N., Bartos, D., 1979. Terminology for model credibility. Simulation 32,
103–104. http://doi.org/10.1177/003754977903200304.
Selle, B., Minasny, B., Bethune, M., Thayalakumaran, T., Chandra, S., 2011. Applicabilityof Richards’ equation models to predict deep percolation under surface irrigation.
Geoderma 160, 569–578. http://dx.doi.org/10.1016/J.GEODERMA.2010.11.005.
Šimůnek, J., Bristow, K.L., Helalia, S.A., Siyal, A.A., 2016a. The eﬀect of diﬀerent fertigation strategies and furrow surface treatments on plant water and nitrogen use. Irrig.
Sci. 34, 53–69. http://dx.doi.org/10.1007/s00271-015-0487-z.
Šimůnek, J., van Genuchten, M.T., Šejna, M., 2016b. Recent developments and applications of the HYDRUS computer software packages. Vadose Zo. J. 15, 25. http://dx.
doi.org/10.2136/vzj2016.04.0033.
Siyal, A.A., Bristow, K.L., Šimůnek, J., 2012. Minimizing nitrogen leaching from furrow
irrigation through novel fertilizer placement and soil surface management strategies.
Agric. Water Manage. 115, 242–251. http://dx.doi.org/10.1016/j.agwat.2012.09.
008.
Soroush, F., Fenton, J.D., Mostafazadeh-Fard, B., Mousavi, S.F., Abbasi, F., 2013.
Simulation of furrow irrigation using the Slow-change/slow-ﬂow equation. Agric.
Water Manage. 116, 160–174. http://dx.doi.org/10.1016/j.agwat.2012.07.008.
Strelkoﬀ, T.S., Clemmens, A.J., Perea-Estrada, H., 2006. Calculation of non-reactive
chemical distribution in surface fertigation. Agric. Water Manage. 86, 93–101.
http://dx.doi.org/10.1016/j.agwat.2006.07.006.
Strelkoﬀ, T.S., Katopodes, N.D., 1977. Border-irrigation hydraulics with zero inertia. J.
Irrig. Drain. Eng. 103, 325–342.
Tabuada, M.A., Rego, Z.J.C., Vachaud, G., Pereira, L.S., 1995. Modelling of furrow irrigation. Advance with two-dimensional inﬁltration. Agric. Water Manage. 28,
201–221. http://dx.doi.org/10.1016/0378-3774(95)01177-K.
USDA-NRCS, 2012. Surface irrigation. In: National Engineering Handbook.
van Genuchten, M.T., 1976. On the accuracy and eﬃciency of several numerical schemes
for solving the convective-dispersive equation. In: Gray, W. (Ed.), Finite Elements inWater Resources. Pentech Press, London, pp. 1.71–1.90.
van Genuchten, M.T., 1978. Mass transport in saturated-unsaturated media: one-dimensional solutions, in: Research Rep. No. 78-WR-11. University, Princeton, Princeton.
van Genuchten, M.T., 1980. A closed-form equation for predicting the hydraulic conductivity of unsaturated soils1. Soil Sci. Soc. Am. J. 44, 892–898. http://dx.doi.org/
10.2136/sssaj1980.03615995004400050002x.
van Leer, B., 1974. Towards the ultimate conservative diﬀerence scheme. II. Monotonicity
and conservation combined in a second-order scheme. J. Comput. Phys. 14, 361–370.
http://dx.doi.org/10.1016/0021-9991(74)90019-9.
van Leer, B., 1979. Towards the ultimate conservative diﬀerence scheme. V. A secondorder sequel to Godunov’s method. J. Comput. Phys. 32, 101–136. http://dx.doi.org/
10.1016/0021-9991(79)90145-1.
Walker, W.R., Humpherys, A.S., 1983. Kinematic-wave furrow irrigation model. J. Irrig.
Drain. Eng. 109, 377–392. http://dx.doi.org/10.1061/(ASCE)0733-9437(1983)
109:4(377).
Warrick, A.W., Lazarovitch, N., Furman, A., Zerihun, D., 2007. Explicit inﬁltration
function for furrows. j. Irrig. Drain. Eng. 133, 307–313. http://dx.doi.org/10.1061/(ASCE)0733-9437(2007) 133:4(307).
Wöhling, T., Fröhner, A., Schmitz, G.H., Liedl, R., 2006. Eﬃcient solution of the coupled
one-dimensional surface—two-dimensional subsurface ﬂow during furrow irrigation
advance. J. Irrig. Drain. Eng. 132, 380–388. http://dx.doi.org/10.1061/(ASCE)07339437(2006) 132:4(380).
Wöhling, T., Mailhol, J., 2007. Physically based coupled model for simulating 1D
surface–2D subsurface ﬂow and plant water uptake in irrigation furrows. II: Model
test and evaluation. J. Irrig. Drain. Eng. 133, 548–558. http://dx.doi.org/10.1061/(ASCE)0733-9437(2007) 133:6(548).

Acknowledgements
We thank Dr. Charles Sanchez from University of Arizona (US) for
kindly providing us experimental data used in model validation and Dr.
Honeyeh Kazemi for organizing the experimental data.
References
Abbasi, F., Simunek, J., van Genuchten, M.T., Feyen, J., Adamsen, F.J., Hunsaker, D.J.,
Strelkoﬀ, T.S., Shouse, P., 2003. Overland water ﬂow and solute transport: model
development and ﬁeld-data analysis. J. Irrig. Drain. Eng. 129, 71–81. http://dx.doi.
org/10.1061/(ASCE)0733-9437(2003) 129:2(71).
Bautista, E., Clemmens, A.J., Strelkoﬀ, T.S., Schlegel, J., 2009. Modern analysis of surface
irrigation systems with WinSRFR. Agric. Water Manage. 96, 1146–1154. http://dx.
doi.org/10.1016/j.agwat.2009.03.007.
Bautista, E., Warrick, A.W., Schlegel, J.L., Thorp, K.R., Hunsaker, D.J., 2016.
Approximate furrow inﬁltration model for time-variable ponding depth. J. Irrig.
Drain. Eng. 142, 4016045. http://dx.doi.org/10.1061/(ASCE)IR.1943-4774.
0001057.
Bautista, E., Warrick, A.W., Strelkoﬀ, T.S., 2014. New results for an approximate method
for calculating two-dimensional furrow inﬁltration. J. Irrig. Drain. Eng. 140,
4014032. http://dx.doi.org/10.1061/(ASCE)IR.1943-4774.0000753.
Beirami, M.K., Nabavi, S.V., Chamani, M.R., 2006. Free overfall in channels with diﬀerent
cross sections and sub-critical ﬂow. Iran. J. Sci. Technol. Trans. B, Eng. 30.
Bradford, S.F., Katopodes, N.D., 2001. Finite volume model for nonlevel basin irrigation.
J. Irrig. Drain. Eng. 127, 216–223. http://dx.doi.org/10.1061/(ASCE)07339437(2001) 127:4(216).
Bradford, S.F., Sanders, B.F., 2002. Finite-volume model for shallow-water ﬂooding ofarbitrary topography. J. Hydraul. Eng. 128, 289–298. http://dx.doi.org/10.1061/(ASCE)0733-9429(2002) 128:3(289).
Brent, R.P., 1974. Algorithms for minimization without derivatives. IEEE Trans. Automat.
Contr. 19, 632–633. http://dx.doi.org/10.1109/TAC.1974.1100629.
Brufau, P., Garcia-Navarro, P., Playán, E., Zapata, N., 2002. Numerical modeling of basin
irrigation with an upwind scheme. J. Irrig. Drain. Eng. 128, 212–223. http://dx.doi.
org/10.1061/(asce)0733-9437(2002) 128:4(212).
Brunetti, G., Saito, H., Saito, T., Šimůnek, J., 2017a. A computationally eﬃcient pseudo3D model for the numerical analysis of borehole heat exchangers. Appl. Energy.
http://dx.doi.org/10.1016/j.apenergy.2017.09.042.
Brunetti, G., Šimůnek, J., Turco, M., Piro, P., 2017b. On the use of surrogate-basedmodeling for the numerical analysis of Low Impact Development techniques. J.
Hydrol. 548, 263–277. http://dx.doi.org/10.1016/j.jhydrol.2017.03.013.
Burguete, J., Zapata, N., García-Navarro, P., Maïkaka, M., Playán, E., Murillo, J., 2009.
Fertigation in furrows and level furrow systems. I: model description and numericaltests. J. Irrig. Drain. Eng. 135, 401–412. http://dx.doi.org/10.1061/(ASCE)IR.19434774.0000097.
Carsel, R.F., Parrish, R.S., 1988. Developing joint probability distributions of soil water
retention characteristics. Water Resour. Res. 24, 755–769. http://dx.doi.org/10.
1029/WR024i005p00755.
Dong, Q., Xu, D., Zhang, S., Bai, M., Li, Y., 2013. A hybrid coupled model of surface andsubsurface ﬂow for surface irrigation. J. Hydrol. 500, 62–74. http://dx.doi.org/10.
1016/j.jhydrol.2013.07.018.
Dormand, J.R., Prince, P.J., 1980. A family of embedded Runge-Kutta formulae. J.
Comput. Appl. Math. 6, 19–26. http://dx.doi.org/10.1016/0771-050X(80)90013-3.
Ebrahimian, H., Keshavarz, M.R., Playán, E., 2014. Surface fertigation: A review, gaps andneeds. Span. J. Agric. Res. http://doi.org/10.5424/sjar/2014123-5393.
Ebrahimian, H., Liaghat, A., Parsinejad, M., Playán, E., Abbasi, F., Navabian, M., 2013.
Simulation of 1D surface and 2D subsurface water ﬂow and nitrate transport in alternate and conventional furrow fertigation. Irrig. Sci. 31, 301–316. http://dx.doi.
org/10.1007/s00271-011-0303-3.
Fahong, W., Xuqing, W., Sayre, K., 2004. Comparison of conventional, ﬂood irrigated, ﬂatplanting with furrow irrigated, raised bed planting for winter wheat in China. Field
Crops Res. 87, 35–42. http://dx.doi.org/10.1016/j.fcr.2003.09.003.
FAO, 2011a. The State of the World’s land and water resources for Food and Agriculture.
Managing systems at risk, Food and Agriculture Organization. ISBN: 978-1-84971326-9.
FAO, 2011b. Irrigation in Southern and Eastern Asia in ﬁgures, AQUASTAT Survey – FAO
Water report 37. ISBN: 978-92-5-107282-0.
FAO, 2015. World fertilizer trends and outlook to 2018, Food and Agriculture
Organization of United Nations.
FAO, 2017. The future of food and agriculture – Trends and challenges. Rome.
Furman, A., 2008. Modeling Coupled Surface-Subsurface Flow Processes: A Review.
Vadose Zo. J. 7, 741–756. http://dx.doi.org/10.2136/vzj2007.0065.
Godunov, S.K., Moscow, S.U., 1954. Ph. D. Dissertation: Diﬀerence Methods for ShockWaves. Moscow State University.
He, Z., Wu, W., Wang, S.S., 2008. Coupled Finite-Volume Model for 2D Surface and 3D
Subsurface Flows. J. Hydrol. Eng. 13, 835–845. http://dx.doi.org/10.1061/(ASCE)1084-0699(2008) 13:9(835).
Hendricks, D.M., 1985. Arizona Soils. College of Agriculture, University of Arizona.

326

Computers and Electronics in Agriculture 150 (2018) 312–327

G. Brunetti et al.

349–356. http://dx.doi.org/10.1061/(ASCE)0733-9437(2004) 130:5(349).
Zapata, N., Playán, E., 2000. Elevation and inﬁltration in a level basin. I. Characterizing
variability. Irrig. Sci. 19, 155–164. http://dx.doi.org/10.1007/s002710000017.
Zerihun, D., Furman, A., Warrick, A.W., Sanchez, C.A., 2005. Coupled surface-subsurface
ﬂow model for improved basin irrigation management. J. Irrig. Drain. Eng. 131,
111–128. http://doi.org/10.1061/(ASCE)0733-9437(2005)131:2(111).

Wöhling, T., Schmitz, G.H., 2007. Physically based coupled model for simulating 1D
surface–2D subsurface ﬂow and plant water uptake in irrigation furrows. I: Modeldevelopment. Irrig. Drain. Eng. 133, 538–547. http://dx.doi.org/10.1061/(ASCE)
0733-9437(2007) 133:6(538).
Wöhling, T., Singh, R., Schmitz, G.H., 2004. Physically based modeling of interacting
surface-subsurface ﬂow during furrow irrigation advance. J. Irrig. Drain. Eng. 130,327