The Furrow Modulefor HYDRUS (2D/3D)
A Hybrid Finite Volume-Finite Element Model for the
Numerical Analysis of Furrow Irrigation and Fertigation
Version 1.0

Giuseppe Brunetti1, Jirka Šimůnek2, and Miroslav Šejna3
April 20191

Institute of Hydraulics and Rural Water Management,
University of Natural Resources and Life Sciences, Vienna, Austria
2

Department of Environmental Sciences
University of California Riverside, Riverside, CA, 92521, USA
3

PC Progress, Ltd., Prague, Czech Republic© 2019 J. Šimůnek and M. Šejna. All rights reserved.

2

Table of Contents
Table of Contents................................................................................................................................. 3
List of Figures ...................................................................................................................................... 5
List of Tables ........................................................................................................................................ 6
Abstract ................................................................................................................................................ 7
1. Introduction ..............................................................................................................................9
2. Conceptual Model....................................................................................................................... 13
3. One-dimensional Overland Flow and Solute Transport ....................................................15
3.1.

Governing Equations ...................................................................................................15

3.2.

Spatial Discretization...................................................................................................16

3.3.

Temporal Discretization ..............................................................................................19

3.3.

Boundary Conditions ...................................................................................................21

4. Two-dimensional Subsurface Flow and Solute Transport .................................................23
4.1.

Governing Equations ...................................................................................................23

4.2.

Soil Hydraulic Properties ............................................................................................23

4.3.

Numerical Domain and Boundary Conditions ............................................................24

5. Models Coupling Strategy .....................................................................................................23
5.1.

Wet/Dry Boundary Tracking ........................................................................................25

5.2.

Model Implementation .................................................................................................25

6. Preprocessing and Input Data ..............................................................................................27
6.1.

Main Processes ............................................................................................................27

6.2.

Furrow Parameters ......................................................................................................27

7. Post-Processing.......................................................................................................................31
7.1.

Results - Profile Information .......................................................................................31

7.2.

Results - Other Information .........................................................................................32

8. Example Problems .................................................................................................................33
8.1.

Problem Description ....................................................................................................33

8.2.

Advance and Recession Times .....................................................................................33

8.3.

Water Depths at Different Stations................................................................................34

8.4.

Actual and Cumulative Infiltrations.............................................................................35

8.5.

Actual and Cumulative Water and Solute Outflow ......................................................36

8.6.

Water Content Profile Distributions ............................................................................36

3

8.7.

Concentrations at Different Stations .............................................................................37

8.8.

Actual and Cumulative Solute Infiltrations..................................................................38

8.9.

Mass Balance Calculations..........................................................................................39

References .....................................................................................................................................41

4

List of Figures
Figure 1.

A schematic of the modeled transport domain consisting of the trapezoidal open
channel and soil (Qinj is the solute injection rate [L3/T], Cinj is the soluteconcentration [M/L3], Qir is the irrigation inflow rate [L3/T], ET isevapotranspiration, P is precipitation, h is the water depth in the furrow, hmax isthe 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 sideslope, L is the length of the furrow, and dx is the discretization length). The two
figures at the bottom show boundary conditions for water flow and solutetransport. .....................................................................................................................13

Figure 2.

A schematic of the temporal discretization. ...............................................................20

Figure 3.

A schematic showing the linear interpolation of the water depth in the ghost cellwhen the free overfall BC is used. .............................................................................22

Figure 4.

The Main Processes dialog window where the Furrow module is selected. ..............27

Figure 5.

The Furrow Parameters dialog window. ....................................................................28

Figure 6.

The part of the Edit Bar where one selects which Station and which Time Layer(Print Time) is displayed. ...........................................................................................32

Figure 7.

The Furrow Parameters dialog window for the example problem. ............................33

Figure 8.

Advanced and recession times for a furrow with open (left) and blocked (right)
ends. ............................................................................................................................34

Figure 9.

Water depth at different stations for a furrow with open (left) and blocked (right)
ends. ............................................................................................................................34

Figure 10. Actual (left) and cumulative (right) infiltrations at different stations for a furrow with
open (top) and blocked (bottom) ends. .......................................................................35
Figure 11. Actual (left) and cumulative (right) water (top) and solute (bottom) outflow at the ofa furrow with an open end. .........................................................................................36
Figure 12. Water content distribution after 7,200 s at the first (left) and last (right) station for a
furrow with an open (top) and blocked (bottom) end. ...............................................37
Figure 13. Concentrations at different stations for a furrow with an open (left) and blocked(right) end. ..................................................................................................................38
Figure 14. Actual (left) and cumulative (right) solute infiltration at different stations for a
furrow with an open (top) and blocked (bottom) end. ...............................................395

List of Tables
Table 1.

The list of output files created by the Furrow module. ................................................31

Table 2.

Mass balance calculations for the furrow systems with open and blocked ends .........40

6

Abstract
Brunetti, G., J. Šimůnek, and M. Šejna, The Furrow Module for HYDRUS (2D/3D) Simulating
Water Flow and Solute Transport in a Pseudo-Three-Dimensional Furrow Irrigation System,
Version 1.0, PC Progress, Prague, Czech Republic, 43 pp., 2019.
This manual presents a hybrid Finite Volume - Finite Element (FV-FE) model that describes thecoupled surface-subsurface flow processes occurring during furrow irrigation and fertigation.
The numerical approach combines a one-dimensional water flow and solute transport in an open
channel with a two-dimensional water flow and solute transport in a subsurface soil domain, thus
reducing the dimensionality of the problem and the computational cost. The modelingframework includes the widely used hydrological model, HYDRUS, which can simulate themovement of water and solutes, as well as root water and nutrient uptake in variably-saturated
soils. The robustness of the proposed model was examined and confirmed by mesh and time stepsensitivity analyses reported in Brunetti et al. (2018). The model was theoretically validated bycomparison with simulations conducted with the well-established model WinSRFR andexperimentally validated by comparison with field-measured data from a furrow fertigation
experiment conducted in the US (Brunetti et al., 2018).
The purpose of this report is to documents the Furrow module for the HYDRUS (2D/3D)
software package (Šimůnek et al., 2018) simulating water flow and solute transport in a pseudothree-dimensional furrow irrigation system. The module, as well as the text of the description of themodule in this report, is largely based on the theoretical development of Brunetti et al. (2018). The
report serves as both a Technical Manual and a User Manual, as well as a reference document of
the Graphical User Interface of Furrow related parts of the HYDRUS software package.

7

DISCLAIMER
The Furrow module was developed as a supplemental module of the HYDRUS (2D/3D software
package to simulate the transport and reactions of colloid and colloid-facilitated transport in soils
and groundwater. The software has been verified against selected test cases. However, no warranty
is given that the program is completely error-free. If you do encounter problems with the code, finderrors, or have suggestions for improvement, please contact the senior author at
Jirka Šimůnek
Tel/Fax: 1-951-827-7854
Email: jiri.simunek@ucr.edu

8

1. Introduction
Agriculture is among the most important human activities due to its role in the food supply chain.
According to one of the latest FAO reports (FAO, 2017), agricultural production more thantripled between 1960 and 2015. This remarkable expansion has been accompanied by a dramaticincrease in the use of irrigation and fertilization, and an associated significant environmental
footprint, thus posing important sustainability issues. Irrigation for agricultural purposes
accounts for 70% of all water withdrawn from aquifers, lakes, and streams (FAO, 2011).
Furthermore, the worldwide demand for fertilizers has grown by 25% in the last decade, and thistrend is expected to continue in coming years (FAO, 2015), increasing the risk of nutrient
pollution of water bodies. Thus, the transition towards more efficient and sustainable irrigationand fertigation strategies is necessary.
Although slowly replaced by pressurized irrigation in developed regions such as The UnitedStates, 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 irrigated surfaces inSouthern and Eastern Asia . In particular, basin and furrow irrigation were the most widespread
techniques among farmers. When correctly performed, furrow fertigation can increase theefficiency of fertilizer use and crop fertilizer uptake, compared to traditional techniques(Šimůnek et al., 2016a; Siyal et al., 2012; Horst et al., 2005; Fahong et al., 2004). To be
successful, furrow irrigation and fertigation systems should be designed and managed so that theapplication and distribution of water and fertilizer are efficient and uniform, with minimal
surface runoff at the lower end of the field and with 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 forassessing irrigation and fertigation efficiency.
Furrow irrigation and fertigation are coupled surface-subsurface processes. Water and solute are
injected at the soil surface at one side of an open channel, hence generating a sharp front thatmoves along a furrow while water and nutrients infiltrate in the underlying soil. Therefore,
furrow irrigation and fertigation are (physically, mathematically, and numerically) described ascoupled three-dimensional flow and transport processes in both surface and subsurface domains.
Nevertheless, solving a fully coupled system of 2D Shallow Water and 3D Richards equations
would require significant computing resources and would also likely pose substantial stability
issues, mainly stemming from the high nonlinearity of the governing equations. This presents a
problem in the modeling of flooding and drying processes over porous beds.
Therefore, several models have been proposed in the literature that reduce the numericalcomplexity of a fully 3D model. Katopodes and Strelkoff (1977a) used a dimensionless
formulation of the governing equations to show that when the Froude number is small, which istypical under irrigation conditions, the inertial terms in the shallow water equations are
negligible. With this in mind, they developed the first zero-inertia model for irrigation (Strelkoff
and Katopodes, 1977b). In the early 80s, Walker and Humpherys (1983) proposed a furrow
irrigation model based on a 1D kinematic wave (KW) approximation of open channel flow
coupled with the modified Kostiakov equation describing the infiltration process. While themodel was assessed against experimental data with satisfactory results, the adoption of the KW
9equation limited its application to open-ended furrows. Furthermore, the model did not include
any description of the subsurface water dynamics. Later on, 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 more detailed
fertigation model was first proposed by Abbasi et al. (2003) and later improved by Perea et al.
(2010). In these two studies, water flow and solute transport in an open channel were describedusing the 1D Zero-Inertia and advection-dispersion equations, respectively. The modified
Kostiakov equation was used to calculate infiltration at each time step. Although the model was
verified with good results against experimental data from four experimental sites, subsurface
flow processes were again neglected.
To provide a better and more complete description of coupled water flow 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 flow and solute transport as well as root water andnutrient uptake in 1D variably-saturated porous media. The computational framework, whichtargeted basin irrigation, was based on the iterative coupling between the surface and subsurface
models and was validated against measured data with good results. Ebrahimian et al. (2013)
extended this concept to furrows by coupling a 1D furrow fertigation model with the twodimensional HYDRUS-2D model (Šimůnek et al., 2016b). The coupled model satisfactorily
reproduced the overland transport as well as the solute transport in the soil profile. However, thesurface and subsurface components were solved separately, leading to an uncoupled numericalframework. As pointed out by Furman (2008), theoretically, the higher the level of coupling, thehigher the solution accuracy. This is mainly due to the high nonlinearity of the involved
processes, as well as to their different time scales. For instance, overland flow is generally fasterthan infiltration, thus requiring a different temporal resolution. This temporal misalignment
poses significant numerical issues since an approximation is needed to couple surface andsubsurface flow. The accuracy of this approximation strongly influences the conservativeness ofthe numerical scheme. Similarly, the overland solute transport needs to be solved simultaneouslywith water flow in order to preserve the monotonicity of the solution.
One of the most complete furrow irrigation models was developed by Wöhling et al. (2004). The
proposed computational framework iteratively coupled a 1D analytical zero-inertia equation with
HYDRUS-2D, thus providing a complete description of surface-subsurface water flow along the
furrow. The model was further developed and extended by Wöhling et al. (2006), Wöhling and
Mailhol (2007), and Wöhling and Schmitz (2007). A similar approach was used by Tabuada et al.
(1995), who coupled a model based on a complete hydrodynamic equation of overland flow witha two-dimensional Richards equation. While accurately simulating water flow, neither of themodels discussed above-considered solute transport in surface and subsurface domains, thus
restricting their applicability to irrigation. Hence, further development of similar approaches thatwould include a detailed description of solute transport in the root zone is desirable for bothscientists and practitioners (Ebrahimian et al., 2014).
Most of the existing furrow irrigation models adopt a Lagrangian approach, which uses a
computational grid that moves along with the wet/dry interface. Although elegant and accurate,
this type of approach can lead to coupling issues between the surface and subsurface models,10

mainly because the grid must be regenerated each time the wet/dry interface moves, and thecomputational nodes often must be added during flooding or removed during the recession to
reduce grid distortion error. However, while surface processes are discontinuous (surface flow
and transport occur only during irrigation events), subsurface processes are continuous(subsurface flow and transport continue between irrigation events). Therefore, a fixed grid forthe subsurface domain is usually used, and values of interest (e.g., infiltration, 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 theinterpolation 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 infiltration and redistribution. Furthermore, most of the existing models are basedon the finite difference method (e.g., Tabuada et al., 1995; Abbasi et al., 2003; Perea et al., 2010),
which can yield spurious oscillations at flow discontinuities unless first-order accurate (upwind)schemes or artificial dissipation are employed.
Godunov-type Finite Volume (FV) schemes have been successfully applied to simulate overland
flow 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 theintegral form of the overland flow and solute transport equations, thus being mass conservative
both globally and locally. Numerical fluxes are evaluated at the cell faces, thus guaranteeing astraightforward and efficient treatment of the dry bed problem by allowing for flooding anddrying of fixed computational cells. Furthermore, numerical oscillations near discontinuities canbe eliminated using flux limiters (Bradford and Katopodes, 2001). However, their application to
furrow irrigation and fertigation has been rather limited and has not involved a coupledmechanistic 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 of describing coupled surface-subsurface flow
and transport processes involved in furrow fertigation. The proposed numerical approachcombines a one-dimensional FV description of coupled water flow and solute transport in thesurface domain with a two-dimensional mechanistic FE description of flow and transport in thevariably-saturated zone, thus reducing the dimensionality (3D) of the problem and associatedcomputational cost. The modeling framework includes the widely-used FE model, HYDRUS-2D,
whose numerical features significantly increased the overall modeling flexibility. The proposedmodel is the first attempt to include HYDRUS-2D in a coupled surface-subsurface furrow
fertigation model, thus representing a new contribution to this field. 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 andthe robustness of the proposed numerical approach. Next, the model is validated against
measured data from an experimental facility in California, US.
The main purpose of this report is to document a pseudo-three-dimensional numerical module
Furrow that incorporates processes associated with furrow irrigation and fertigation. The
Furrow module was developed specifically for the HYDRUS (2D/3D) software package (Šimůnek
et al., 2008). The general conceptual basis of the module is discussed in the manuscript of11

Brunetti et al. (2018). The documentation in this report focuses mostly only on the surface overland
flow and transport features of the Furrow module. All processes related to variably-saturated waterflow and heat transport are described in detail in the HYDRUS (2D/3D) documentation (Šimůnek etal., 2018), and will not be repeated here. The report serves as both a Technical Manual and a User
Manual, as well as a reference document of the Graphical User Interface of Furrow related partsof the HYDRUS software package.

12

2. Conceptual Model
The proposed approach combines a one-dimensional description of coupled water flow andsolute transport in an open channel with a two-dimensional description of variably-saturated
water flow and solute transport in soil. As in previous studies (e.g., Tabuada et al., 1995;
Wöhling et al., 2004), a key assumption in the model is that the transport processes in the soil
domain occur only in the vertical plane, as a function of conditions at a surface flow crosssection, 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 of a widely used hydrological model, HYDRUS, in themodeling framework.

Figure 1. A schematic of the modeled transport domain 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 inflow rate [L3/T], ET is evapotranspiration, P is precipitation, h is the water
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 sideslope, L is the length of the furrow, and dx is the discretization length). The two figures at thebottom show boundary conditions for water flow and solute transport.

13

The furrow of length L is discretized into N elements, each one characterized by a twodimensional 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.
The interface between the surface and subsurface domains acts as a common boundary.

14

3. One-dimensional Overland Flow and Solute Transport
3.1. Governing Equations
The one-dimensional overland system can be described in the conservative form (Burguete et al.,2009) using the cross-sectional water and solute mass conservation, momentum balance, andinfiltration as follows:

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

(1)

where U is the vector of the conserved variables, F is the flux vector, S is the source term vector,
I is the infiltration vector, and D is dispersion. Vectors can be expressed as:




 Q
0
 −i 


 A 
 0 

2 






Q 
U =  Q  , F =  gI1 +, I =  0  , S =  g I 2 + A( S 0 − S f  , D =  0  (2)


A
∂C 
 − iC  AC 
0




K
A






x
 QC 
∂x 


[

]

where A is the wetted area [L2], Q is the discharge [L3T-1], C is the cross-sectional average soluteconcentration [ML-3], g is the gravitational acceleration [LT-2], S0 is the bed slope [LL-1], Sf is thefriction slope [LL-1], i is the infiltration rate [L3T-1], Kx is the dispersion coefficient [L2T-1], andI1 and I2 are pressure forces. It must be emphasized that the infiltration rate i is estimated using
the HYDRUS-2D model.
Since flow velocities in furrow irrigation are typically small, the resulting Froude number isgenerally less than 0.3. Under such conditions, the inertial terms are negligible and the shallow
water equations can be simplified to the well-established zero-inertia (ZI) approximation(Strelkoff and Katopodes, 1977):

QQm W
∂h
= S0 −∂x
cu2 A10 3
2

43

(3)

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

15




Q = sign(J ) ⋅ k J∂hJ = S0 −
∂x

5

3
 k = A cu2

mW 3


(4)

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

QQm W
∂h
= S0 −∂x
cu2 A10 3
2

43

(5)

was simplified as well by considering only the advective component (eq. 5), as suggested by
Strelkoff et al. (2006). In that study, the results of the pure advection equation were comparedwith the results of the full ADE for different modeling scenarios. The results indicated that theformer formulation was able to provide a sufficiently accurate description of the solute
distribution along the furrow, although the effects of the longitudinal dispersion were appreciable,
especially for blocked-end furrows. Neglecting the dispersion component simplifies thecomputational scheme and the implementation of boundary conditions.
Considering the above-described simplifications, eq. (2) was rewritten as:

 A 
 Q 
 −i 
 , F =  , I = U = 

 AC 
 QC 
 − iC Q = sign( J ) ⋅ k J(6)

3.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 Differential Equations (PDEs)
with an algebraic approximation, which reduces the PDE to an Ordinary Differential 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 a well-established and widely used numerical solution method for a variety ofengineering problems involving PDEs, mainly because of its conservativeness. It has been

16

applied successfully in many studies focused on overland flow and solute transport (e.g., Lal,
1998; He et al., 2008). The integral form of eq. (6) was written as:∂
UdΩ + ∫ ∇ ⋅ FdΩ = ∫ IdΩ
Ω
Ω
∂t ∫Ω

(7)

Application of the divergence theorem to the second term in eq. (7) produces:∂
UdΩ + ∫ F ⋅ ndω = ∫ IdΩ
Ω
ω
∂t ∫Ω

(8)

The contour integral in eq. (8) was approximated using numerical fluxes F at the edge of eachcell j, thus leading to:

dU j
dt

=−


1 
Fj+ 1 − Fj−1  + I j∆x  22

(9)One of the main drawbacks of the FVM is the introduction of artificial diffusion in the solution
when first-order numerical schemes are adopted (e.g., upwind). Indeed, while preserving themonotonicity of the solution (Godunov, 1954), first-order numerical schemes have lower
accuracy compared to the higher-order schemes (e.g., central difference). However, the latter areprone to numerical oscillation (e.g., van Genuchten, 1976, 1978). Therefore, to increase thenumerical accuracy while preserving the stability of the numerical solution, a high-resolution
scheme was adopted in this study.
In this study, the Monotone Upstream-Centered Scheme for Conservation Laws (MUSCL) (vanLeer, 1979) was used to spatially discretize the coupled water-solute equations in theconservation form (eq. 1). The MUSCL methods switch from the first-order scheme in regionswith sharp fronts to the higher-order schemes in areas with a relatively smooth solution, thus
simultaneously increasing the overall accuracy and avoiding non-physical oscillations. A flux
limiter was used to constrain the value of the spatial derivative around sharp discontinuities. In
the present study, the van Leer (1974) flux limiter was used in conjunction with the upwind andcentral differencing schemes. The upwind and central fluxes for overland water flow wereexpressed as: upwind  Q j −1 

= 
F 1

 Q j −1C j −1  j− 2
 1

(Q j −1 + Q j )  F central =  2
 j − 12 1 (Q C + Q C )
j −1 j −1
j j


2


17

(10)

where the 1D cell centered discharge Q was expressed as:

S o j +1 + S o j h j +1 − h j− J j =
∆x
2

k
k
+
+
1
j
j
Q =
J j ⋅ sign(J j )
 j
2

(11)

It is worth noting that channel conveyance and bed slope were both represented with thearithmetic mean of nodal values. Preliminary simulations confirmed that this type ofapproximation combined with the MUSCL scheme leads to more accurate results compared to a
fully upwind discretization of k and S0, which introduces significant numerical diffusion in thesolution. The same discretization is used for the solute transport. Hence, the numerical flux F
obtained using the MUSCL scheme can be written as:

F

1
j −.
2




= F upwind+ φ  F central− F upwind
1
1
1

j−
j−
j−
2
2
 2


(12)

where ϕ is the van Leer flux limiter [LL-1], which was expressed as:
rj + rj
φ j (r ) =
1 + rj lim φ (r ) = 2
∞ j
r

(13)

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

rj =u j − u j −1
u j +1 − u j(14)where u is a dependent variable, which in our case can be the wetted area A for water flow or across-sectional average solute concentration C for solute transport.
Based on eq. (12), the first-order monotonicity preserving upwind was adopted in regions
characterized by steep gradients in the solution, while the higher-order central differencing
scheme was used in regions with a relatively smooth solution, thus increasing the stability andaccuracy of the spatial scheme. To facilitate the handling of the downstream boundary condition,
the upwind scheme was used to compute the numerical flux in the last cell (i.e., j=N). An upwind
discretization was also used for the infiltration term.

18

3.3. Temporal Discretization
Once the spatial discretization is defined, time t remains the only independent variable. Since theinfiltration 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 efficiency of thenumerical scheme is strictly related to the frequency of data exchanges between the surface andsubsurface 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 framework, are only first-order accurate in time.
Thus, due to aforementioned considerations about accuracy and computational efficiency, thefourth-order explicit Runge-Kutta Dormand-Prince (RKDP) (Dormand and Prince, 1980) is
adopted here to solve the coupled water-solute transport in the trapezoidal channel. Morespecifically, 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 flow can be solved separately from the solute
advection equation. However, the solute advection equation depends on the wetted area and flow
velocity, leading to a one-way coupled system of ODEs. More specifically, the temporalvariability of the wetted area strongly influences the solute advection. Thus, a high temporal
resolution scheme for water flow is needed to avoid overshooting as well as inaccuracies in thesolution of solute advection.
The basic method will be first described for the pure conservation law, that is, without source
terms. These terms will be incorporated next and required modifications will be indicated asneeded. Eq. (10) can now be written as:

U ij+1 = U ij −∆tint  i
F 1 − Fi 1 j−∆x  j + 2
2 

(15)

where ∆tint is the internal time step calculated by the adaptive RKDP algorithm. To increase theoverall computational efficiency of the modeling framework, a specific numerical treatment of
the infiltration component is adopted here. As discussed above, the RKDP algorithm uses a fine
temporal discretization near discontinuities, implying that the calculation of infiltration (i.e.,
HYDRUS executions) is carried out with very small time steps, which consequently increases
the computational cost. However, the infiltration vector I of eq. (3) mathematically represents asource/sink term, which is characterized by a moderately low stiffness. Thus, an explicit firstorder approximation of infiltration (Brufau et al., 2002) is used in our model, which reduces thenumber of HYDRUS-2D executions while maintaining good accuracy. This leads to a mixed
explicit time scheme, which uses a 4th order-accurate RKDP time-stepping strategy for the
coupled water-solute transport and an explicit first-order scheme for the infiltration term.
Practically, this is accomplished by dividing the simulation duration into equal time steps ∆t. At
the beginning of each time step, the values of the water depth are passed to HYDRUS-2D, whichcalculates the amount of infiltrated water during ∆t. Next, the vector I is inserted in eq. (1),19

which is solved by the RKDP method by adopting an internal time step ∆tint depending on thesolution. At time t=0, the surface domain is solved assuming zero infiltration. Equation (15) can
now be rewritten by including the infiltration term:

U

i +1
j

 I ij ∆tint  ii =U −F − F 1  − ∆tint   ∆t 
j−∆x  j + 12
2 
 
i
j

(16)

The temporal discretization is depicted in Figure 2. The non-iterative coupling with the external
HYDRUS-2D model significantly decreases the computational cost. However, it must be notedthat an ad hoc sensitivity analysis is recommended for the choice of ∆t.

Figure 2. A schematic of the temporal discretization.
Since eq. (16) is solved explicitly in time, the Courant-Friedrichs-Lewy (CFL) stability condition
must be respected. In particular, the maximum allowable ∆tint calculated by the RKDP algorithm
must be:




∆x


max
∆tint < minii 
 Q + g A 
i
Γi 
A

(17)

where Γ is the top width of the water flow [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:

(

)

ti
i

(Qir + Qinj )dt − ∫0ti FNi + 1 dt − ∫0L Ai dx − ∫0L VWi − VW0 dx − ∫0ti qsoil
dt
∫
0

i
2
⋅100
 RMBEwater (%) =ti(Qir + Qinj )dt
∫
0

ti
ti
L
L
ti
i
i
i i
i
i
0
(
)
−
Q
C
dt
F
C

inj inj
1
1 dt − ∫ A C dx − ∫ VS − VS dx − ∫ q soil c dt
∫
∫
0
0
0
0
0
N+
N+
 RMBE i (%) =22⋅100
solute
ti

(
)
Q
C
dt
ir
inj
∫0


(

20

)

(18)

where ViW and ViS are, respectively, the water volume and solute mass per unit length in the soil
profile 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, fluxes
across the bottom of the soil domain (e.g., deep percolation and solute leaching).
3.4. Boundary Conditions
The problem statement is completed by the definition of the boundary conditions (BCs) at thetwo sides of the numerical domain. In the present study, the ghost cells methodology is used
(LeVeque, 2002). The 1D ZI equation requires the definition of two BCs, namely inflow andoutflow. The specified inflow Q(0,t) determines the numerical flux entering the domain on theleft face of the first discretized volume. Two types of BCs can be used for the outlet dependingon the type of the furrow end. Since the adopted spatial stencil spans two adjacent cells, thevalue of the wetted area A in the right ghost cell has to be adjusted at each time step to match thedesired BC.
If the furrow downstream end is blocked, no outflow is allowed until the water level reaches the
maximum furrow depth. If h>hmax, the model automatically switches to a free overfall BC. Thus,
the ghost cell value has to be set so that the numerical flux in the last cell (i.e., j=N) is zero.
Considering that the upwind scheme is used for the rightmost numerical flux, eq. (12)mathematically reduces to eq. (19):
5
 53

3


A
A
 hghost − h j 
1
j
ghost
 = 0F 1 = 2 + 2  S o − 
j+
2
mc
∆
x
u 



2
3
3
 W j Wghost 

if j = N(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 flow over a brink at theend 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 flow conditions,the upstream flow passed over the brink with a depth hb lower than the critical depth hc. Thisphenomenon was previously investigated by Beirami et al. (2006), who indicated that for
trapezoidal channels hb ≈ 0.75hc. Hence, in the present study, hghost is determined at each time
step by interpolation so that hb = 0.75hc (Fig. 3).
Both the blocked and open end BCs are applied only when the value of h in the last cell is greaterthan 0.001 m. To avoid convergence issues , if h<0.001 m, the water depth in the ghost cell, hghost,is set equal to hj=N.

21

Figure 3. A schematic showing the linear interpolation of the water depth in the ghost cell when
the free overfall BC is used.
The solute advection equation requires the implementation of the boundary condition only on the
inflow side of the domain. In the conservation form, this was accomplished by matching thenumerical flux with the solute injection rate Qinj·Cinj (Fig. 1).

22

4. Two-dimensional Subsurface Flow and Solute Transport
4.1. Governing Equations
The HYDRUS (2D/3D) (Šimůnek et al., 2016b) software simulates water flow and solutetransport in variably-saturated porous media. The software describes water flow using the
Richards equation for multi-dimensional unsaturated flow:∂θ∂  
∂ψ=+ K iz  K  K ij∂t ∂yi  ∂y j
(i, j = 1,2)(20)

where θ is the volumetric water content [L3L-3], t is time [T], z is the vertical coordinate [L], yi
are the spatial coordinate [L], ψ is the soil pressure head [L], K is the hydraulic conductivity [LT1], and Kij are components of the hydraulic conductivity anisotropy tensor. Since we assumed
isotropic porous media with y1 = y and y2 = z being the transverse (horizontal) and verticalcoordinates, respectively, the conductivity tensor is diagonal with the Kyy and Kzz entries equal toone.
Solute transport is described using the advection dispersion equation:

∂θc
∂ 
∂c  ∂qi
=
Dij
θ
∂t
∂yi 
∂y j  ∂yi(i, j = 1,2)

(21)

where c is the solute concentration [ML-3], qi is the water velocity [LT-1], and Dij are componentsof the dispersion tensor (L2T-1). Eq. (21) is valid for non-reactive transport, thus adsorption andprecipitation/dissolution of the solutes are currently ignored. Note that although HYDRUS-2D
considers multiple chemical reactions (e.g., sorption and degradation), 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 Crank-Nicholson scheme for the advection-dispersion equation.
The program also uses an adaptive time-stepping strategy to increase the overall computational
efficiency.
4.2. Soil Hydraulic Properties
The van Genuchten–Mualem (VGM) model (van Genuchten, 1980) was used to describe the soil
hydraulic properties in this study:

23

1

if ψ ≤ 01
1−nΘ =  (1 + (α ψ ) ) n
if ψ > 01θ − θrΘ=θs − θr
1 2

1− 

n
 

 n
 K s Θ l 1 − 1 − Θ n−1   

K =
 

 





K
 s

if ψ < 0

(22)

(23)

(24)

if ψ > 0

where Θ is the effective 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 a pore-size distribution index [-], Ks is the saturated hydraulic conductivity [LT-1],
and l is the tortuosity and pore-connectivity parameter [-].
4.3. Numerical Domain and Boundary Conditions
The soil computational domain was discretized using triangular 2D finite elements. In theexample simulations discussed below, no mesh stretching was used, and the finite 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 infiltration are reported in Figure 1.
The ‘Atmospheric’ boundary condition, which is assigned to the furrow ridge (green lines in Fig.
1), can exist in three different states: (a) precipitation and/or potential evaporation fluxes, (b) azero pressure head (full saturation) during ponding when both infiltration and surface runoff
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 by water fluxes
towards the soil surface. Due to symmetry, the nodes representing the left and right boundary of
the subsurface domain are set as ‘no flux’ boundaries (grey lines in Fig. 1) because no flow orsolute transport occurred across these boundaries.
A hybrid Dirichlet/Neumann (variable pressure/zero flux) 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 specified 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 basedon their elevation relative to the channel 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 the concentration flux along the boundaries (red lines in Fig. 1).

24

5. Models Coupling Strategy
5.1. Wet/Dry Boundary Tracking
In all numerical simulations, the water front is not considered as a moving boundary, andcalculations are carried out in both wet and dry cells. The numerical treatment of wet/dry cells issimilar 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, waterlevels above 0.0001 m indicate wet cells (Bradford and Katopodes, 2001). If hj>ε, the water
depth value is passed to HYDRUS, which calculates infiltration that is then used to solve surface
flow. Two computational procedures can be followed when the cell is dry:
1. Subsurface in steady-state conditions: HYDRUS is completely bypassed and infiltrationis set to 0. For example, during the advance, the last part of the furrow is not simulated,since it is flagged as “dry” and the subsurface flow is negligible. It is worth noting that
this option, while increasing the computational efficiency of the numerical scheme, canbe considered only when the initial condition is in hydrostatic equilibrium and boundary
fluxes are negligible.
2. Subsurface in dynamic conditions: These computations apply during recession. Although
no more water infiltrates, HYDRUS can continue to simulate variably-saturated waterflow 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 final soil condition is stored and used in the next step.
5.2. Model Implementation
The one-dimensional horizontal domain representing the furrow and the two-dimensional crosssectional domains representing the subsurface (Fig. 1) were coupled using a Python script. The
script simultaneously solves overland flow and solute transport in the furrow, interacts with
HYDRUS-2D, and exchanges data between two models. A series of user-defined functions weredeveloped to write and read the input/output files generated by HYDRUS-2D, which was
directly executed from Python, bypassing the HYDRUS (2D/3D) graphical user interface. The
coupling strategy is summarized below:
1. The geometric characteristics of the furrow, soil hydraulic properties, initialconditions and other input data are defined 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 soluteconcentration distributions in the soil are stored in three matrices, which are
overwritten at each time step. Water content is stored only for visualization purposes
and not used directly in the coupling procedure. Similarly, vectors containing wetted
areas, water depths, and solute concentrations in the surface domain are updated at25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 and solute concentrations in
the furrow, respectively, are passed to HYDRUS-2D to setup its boundary conditions.
Soil pressure head and 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 the infiltration vector.
Calculated soil pressure heads and solute concentrations are stored for the next time
step. The coupled ODEs are solved simultaneously using the RKDP algorithm, andnew wetted area, water level, and solute concentration vectors are updated in theirrespective matrices.

26

6. Preprocessing and Input Data
6.1. Main Processes
The Furrow module is activated by checking the "Furrow Irrigation and Fertigation"
checkbox in the "Main Processes and add-on Modules" dialog window (Fig. 4).

Figure 4. The Main Processes dialog window where the Furrow module is selected.

6.2. Furrow Parameters
Input parameters (in addition to those required by the standard HYDRUS (2D/3D) program)required by the Furrow module are specified in the "Furrow Parameters" dialog window (Fig.
5).

27

Figure 5. The Furrow Parameters dialog window.
In the "Furrow Parameters" dialog window, users specify parameters required by the Furrow
Module that are needed to fully define the geometry and parameters of the furrow, the furrow
spatial and temporal discretizations, as well as the boundary conditions defining irrigation andfertigation.
Main Input
Bottom width of the furrow [L]bwSS
Slope of the furrow sides [-] defined as tan(α) where α is the angle of the slope side.
hmax Depth of the furrow [L]L
Length of the furrow [L]Mn
Manning's roughness coefficient [TL-1/3]
Cinj
Solute concentration in the irrigation water [ML3]Spatial and Temporal Discretization
N
Number of Stations, i.e., locations along the furrow with 2D cross-sectional profiles [-].
dt
Time interval at which information is exchanged between the Open Channel and
Subsurface modules [T]Outflow Boundary TypeOpen EndThe end of the furrow is open, allowing water to flow over a bring at the end of
the furrow

28

Blocked End The end of the furrow is blocked, i.e., no outflow is allowed until the water level
reaches the maximum furrow depth hmax
Time Variable Boundary Conditions
Number of Records Number of time intervals with different boundary conditions (when waterof solute inflow changes)
Time Time when the boundary condition changes [T]Irrigation rate [L3T-1]
Qir
Qinj
Fertigation rate [L3T-1]
Spatial Coordinates of Cross-Section Profiles/Stationsx
Distance of the Cross-Section Profile/Station from the beginning of the furrowz
Height of the Cross-Section Profile/Station relative to the end of the furrow

29

30

7. Post-processing
7.1. Results - Profile Information
The Furrow module generates multiple output files that are listed in Table 1.

Table 1. The list of output files created by the Furrow module.
Output File Name

Output File Description

adv.out

Advance and Recession Times, i.e., times when the advanced andrecession fronts arrive at different Stations.

c.out

Concentrations at different Stations. The first column is Time [T], thesecond and further columns are Concentrations at different (N) Stations[L].

c_inf.out

Solute Infiltration Rates [M/L/T] at different Stations versus time.

Conc1-j.out

Subsurface spatial distribution of solute Concentrations [M/L3] at Station
j at different Time Layers (Print Times).

h-j.out

Subsurface spatial distribution of solute Pressure Heads [L] at Station j at
different Time Layers (Print Times).

runoff.out

Water [L3/T] and Solute [M/T] Runoff from the end of the Furrow

th-j.out

Subsurface spatial distribution of Water Contents [L3/L3] at Station j at
different Time Layers (Print Times).

w.out

Water Depths at different Stations. The first column is Time [T], thesecond and further columns are Water Depths at different (N) Stations[L].

w_inf.out

Infiltration Rates [L2/T] at different Stations versus time.

The display of results is divided, similarly as for the standard HYDRUS code into two groups. In
the Results section, users can display spatial distributions of Pressure Heads, Water Contents(Fig. 12), and Concentrations (when solute transport is considered) at different Print Times fordifferent Stations (i.e., locations along the furrow) in the subsurface. This results section thus
display the results for the subsurface flow/transport part of the Furrow module. The selection
which variable is displayed is made at the Navigator bar. The selection which Station andwhich Time Layer is displayed is made on the Edit Bar using either a List Box or a Scroll Bar(Fig. 6).

31

Figure 6. The part of the Edit Bar where one selects which Station and which Time Layer (Print
Time) is displayed.

7.2. Results – Other Information
In the Results - Other Information section, users can display selected results for the overland
flow part of the Furrow Module. The first set of graphs (Stations) displays Water Depths (Fig. 9),
Actual and Cumulative Infiltrations (Fig. 10), Concentrations, and Actual and Cumulative Solute
Infiltration Fluxes (Fig. 14) at different Stations versus time. This set of graphs corresponds withthe standard Observation Points graphs. The second set of graphs (Boundary Fluxes) displays
Runoff, Cumulative Runoff, Solute Runoff, and Cumulative Solute Runoff (Fig. 11) at the end of
the furrow versus time. Finally, the last graph displays Advance and Recession Times (Fig. 8).
Examples of these graphs can be seen in the following Section.

32

8. Example Problems
8.1. Problem Description
Two test examples are used below to demonstrate the use of the Furrow module. Both examples
have the same geometry, the same spatial and temporal discretizations, and the same irrigation
conditions. They differ only in the boundary condition at the end of the furrow and fertigation
conditions. While the 3DOpenFurrow example has an open furrow end and fertigation during
first 20 minutes of the irrigation period of 30 minutes, the 3DClosedFurrow example has the endof the furrow blocked and fertigation during the entire irrigation period. The furrow is 100 mlong and is almost horizontal with a uniform slope of 0.0005 (5 cm/10,000 cm; 0.0286o). The
bottom of the furrow is 16 cm wide, 15 cm deep, and the furrow sides have a slope of 1.2. The
Manning coefficient for the furrow is 0.0004 s/cm1/3. The furrow is discretized into 10 stations
and the time step is 30 s. Water is applied for 1,800 s at a rate of 2,000 cm3/s (2L/s). Solute is
injected at a rate of 14 cm3/s with a concentration of 0.04 g cm3. Figure 7 shows the dialog
window with the additional input parameters required by the Furrow module.

Figure 7. The Furrow Parameters dialog window for the example problem.

8.2. Advance and Recession Times
Figure 8 (data are taken from the adv.out file) shows the advance and recession times for a
furrow with open and blocked ends. While the advance time (movement of the waterfront in the
furrow) is the same in both types of the furrow, the recession front is delayed in the blocked-end
furrow because of the backup of water in the furrow after the waterfront reached the blocked end.
33

Figure 8. Advanced and recession times for a furrow with open (left) and blocked (right) ends.

8.3. Water Depths at Different Stations
Figure 9 (data are taken from the w.out file) shows water depths for a furrow with open andblocked ends. In a furrow with an open end, the water depth is highest at the beginning of the
furrow with the highest water flux and lowest at the end of the furrow with the smallest water
flux. The water flux (and correspondingly the water depth) decreases along the furrow due toinfiltration. A relatively constant water depth is reached at each station after an advanced front
passes through a station. In a furrow with a blocked end, once the waterfront reaches the end of
the furrow, the water depth starts increasing throughout the entire furrow and reaches higher
values than in the open-end furrow. The water depth is the highest at the end of the blocked-end
furrow. The decrease of the water depth is faster in the open-end furrow because of simultaneous
infiltration and outflow, while it is slower in the blocked-end furrow where water is removed
only by infiltration. Also, note that the rate of the water depth decrease at different locations issimilar in the blocked-end furrow due to infiltration.

Figure 9. Water depth at different stations for a furrow with open (left) and blocked (right) ends.

34

8.4. Actual and Cumulative Infiltrations
Figure 10 shows actual and cumulative infiltrations at different stations for a furrow with open
and blocked ends. Infiltration fluxes are data taken from the w_inf.out file, while cumulative
fluxes are obtained by integrating infiltration fluxes. While cumulative infiltration curves are
relatively smooth, the actual infiltration curves are rather erratic. This is a result of the numericalsolution. As the water depth in a furrow increases or decreases, more or fewer boundary nodes
on sides of a furrow become active in the infiltration process. When the water depth in a furrow
is increasing, more boundary nodes are becoming active and when it is decreasing, fewer nodes
are active. This leads to sudden increases and/or decreases in the infiltration rate as additional or
fewer nodes become active, respectively.

Figure 10. Actual (left) and cumulative (right) infiltrations at different stations for a furrow with
open (top) and blocked (bottom) ends.
Figure 10 also shows that in the open-end furrow, significantly more water infiltrates at the
beginning of a furrow than at its end. While water is available at different locations of the furrow
for about the same time period (Fig. 8), a larger cross-sectional area (a higher water depth) is
involved in infiltration at the beginning than the end of the furrow (Fig. 9). The opposite occurs
35

in the blocked-end furrow where more water infiltrates at the end than the beginning of the
furrow. This is not only because a larger cross-sectional area is involved in infiltration at the endof the furrow (Fig. 9), but also because water is available for a longer time period (Figs. 8 and 9).
8.5. Actual and Cumulative Water and Solute Outflow
Figure 11 shows actual and cumulative outflow of water at the of a furrow with an open end.
There is no outflow from a furrow with a block end. The runoff flux is taken from the Runoff.out
file, while the cumulative outflow is obtained by integrating the runoff flux. The solute runoff isobtained by multiplying the runoff flux with the concentration at the outflow station, and thecumulative solute runoff is again obtained by integrating the actual solute runoff.

Figure 11. Actual (left) and cumulative (right) water (top) and solute (bottom) outflow at the of afurrow with an open end.

8.6. Water Content Profile Distributions
Figure 12 (drawn based on the th-j.out files) shows the spatial distribution of water contents at
the first (left) and last (right) station, for a furrow with an open (top) and blocked (bottom) end
after 3,600 and 7,200 s, respectively. Figure 12 shows that more water infiltrates in a furrow witha blocked end than with an open end since in the latter case some water is lost due to runoff at
the end of the furrow.

36

Figure 12. Water content distribution after 7,200 s at the first (left) and last (right) station for a
furrow with an open (top) and blocked (bottom) end.

8.7. Concentrations at Different Stations
Figure 13 (data are taken from the c.out file) shows concentrations at different stations for a
furrow with an open and blocked end. The concentration should go instantaneously from zero to
the max concentration (cmax) with the arrival of the waterfront. This is not possible for anumerical scheme used in the Furrow module, which will always lead to some (minimal)
numerical dispersion. Similarly, the concentration front for the case with the blocked end, when
fertigation is applied till the end of the irrigation event, should be sharp, going directly from themax concentration to zero. On the other hand, for the case with the open end, when fertigation isapplied only during the first half of the irrigation event, the concentration front should be smooth,37

due to diffusion and dispersion at the concentration front at the end of the fertigation event. The
max concentration (cmax) can be calculated from the inflow concentration (cin=0.04 g/cm3) andthe ratio of infiltration (Qir=2000 cm3/s) and fertigation (Qinj=14 cm3/s) fluxes: = cmax = cin
Qinj/(Qir+Qinj) = 0.04*14/(2000+14) = 2.781e-4 g/cm3 = 278.1 g/m3.

Figure 13. Concentrations at different stations for a furrow with an open (left) and blocked(right) end.
8.8. Actual and Cumulative Solute Infiltrations
Figure 14 (data are taken from the c_inf.out file) shows actual and cumulative solute infiltration
fluxes at different stations for a furrow with open and blocked ends. Similarly as in Figure 10,
while cumulative solute infiltration curves are relatively smooth, the actual solute infiltration
curves are rather erratic as a result of the numerical solution (as discussed in text related to
Figure 10).
Figure 14, similarly to Figure 10, also shows that in the open-end furrow, significantly more
solute infiltrates at the beginning of a furrow than at its end. The opposite occurs in the blockedend furrow where more solute infiltrates at the end than the beginning of the furrow.

38

Figure 14. Actual (left) and cumulative (right) solute infiltration at different stations for a furrowwith an open (top) and blocked (bottom) end.
8.9. Mass Balance Calculations
Table 2 shows the mass balance calculations for the furrow systems with open and blocked ends.
The irrigation (I) and fertigation (F) fluxes are taken from input values. Cumulative water
outflow (O) is taken from the Runoff.out file, while cumulative solute outflow (Os) is obtainedby integrating over time the product of water outflow and the concentration at the outflow end.
Cumulative infiltration is calculated in two way. First, simply as a difference between irrigationand outflow (I-O), and second, as a product of average cumulative infiltration (from w_inf.out)
multiplied by the length of the furrow. Similar calculations were also done for solute massbalance (Table 2).

39

Table 2. Mass balance calculations for the furrow systems with open and blocked ends.
Flux
Irrigation, I

Blocked End
Open End
31800 s * (2000+14 cm /s) 1800 s * 2000 cm3/s + 1200 s * 14
= 3.6252 m3
cm3/s= 3.6168 m3
01.381 m3
3.6252 m3
2.236 m3
1
2
3
0.0364 m * 100 m = 3.64 m
0.02231 m3 * 100 m = 2.23 m3

Water outflow, O
Infiltration (I-O)
Infiltration(average
flux1*
furrow length)
Fertigation, F
1800 s * 14 cm3/s * 0.04 g/cm3= 1008 g
Solute outflow, Os
0
Solute Infiltration (I- 1008 g
Os)
Solute Infiltration
9.981 g/m * 100 m = 998 g(average solute flux1*
furrow length)
11200 s * 14 cm3/s * 0.04 g/cm3 = 672
g
268 g
404 g
4.121 g/m * 100 m = 412 g

Calculated as average cumulative (solute) infiltration from values at individual locations.

40

References
Abbasi, F., Simunek, J., van Genuchten, M.T., Feyen, J., Adamsen, F.J., Hunsaker, D.J.,
Strelkoff, T.S., Shouse, P., 2003. Overland water flow and solute transport: Modeldevelopment and field-data analysis. J. Irrig. Drain. Eng. 129, 71–81.
https://doi.org/10.1061/(ASCE)0733-9437(2003)129:2(71)Bautista, E., Clemmens, A.J., Strelkoff, T.S., Schlegel, J., 2009. Modern analysis of surface
irrigation systems with WinSRFR. Agric. Water Manag. 96, 1146–1154.
https://doi.org/10.1016/j.agwat.2009.03.007
Beirami, M.K., Nabavi, S. V, Chamani, M.R., 2006. Free overfall in channels with different
cross sections and sub-critical flow. 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. https://doi.org/10.1061/(ASCE)0733-9437(2001)127:4(216)Bradford, S.F., Sanders, B.F., 2002. Finite-Volume Model for Shallow-Water Flooding of
Arbitrary
Topography.
J.
Hydraul.
Eng.
128,
289–298.
https://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. https://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.
https://doi.org/10.1061/(asce)0733-9437(2002)128:4(212)Brunetti, G., Šimůnek, J., Bautista, E., 2018. A hybrid finite volume-finite element model for thenumerical analysis of furrow irrigation and fertigation. Comput. Electron. Agric. 150, 312–
327. https://doi.org/10.1016/J.COMPAG.2018.05.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 Numerical
Tests. J. Irrig. Drain. Eng. 135, 401–412. https://doi.org/10.1061/(ASCE)IR.19434774.0000097
Dong, Q., Xu, D., Zhang, S., Bai, M., Li, Y., 2013. A hybrid coupled model of surface andsubsurface
flow
for
surfaceirrigation.
J.
Hydrol.
500,
62–74.
https://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. https://doi.org/10.1016/0771-050X(80)90013-3
Ebrahimian, H., Keshavarz, M.R., Playán, E., 2014. Surface fertigation: A review, gaps andneeds. Spanish J. Agric. Res. https://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 flow and nitrate transport in alternate
and conventional furrow fertigation. Irrig. Sci. 31, 301–316. https://doi.org/10.1007/s00271011-0303-3
Fahong, W., Xuqing, W., Sayre, K., 2004. Comparison of conventional, flood irrigated, flat
planting with furrow irrigated, raised bed planting for winter wheat in China. F. Crop. Res.
87, 35–42. https://doi.org/10.1016/j.fcr.2003.09.003
FAO, 2017. The future of food and agriculture – Trends and challenges. Rome.
FAO, 2015. World fertilizer trends and outlook to 2018, Food and Agriculture Organization ofUnited Nations.
FAO, 2011. The State of the World’s land and water resources for Food and Agriculture.

41

Managing systems at risk, Food and Agriculture Organization. https://doi.org/978-1-84971326-9
Furman, A., 2008. Modeling Coupled Surface–Subsurface Flow Processes: A Review. Vadose
Zo. J. 7, 741–756. https://doi.org/10.2136/vzj2007.0065
Godunov, S.K. (Moscow S.U., 1954. Ph. D. Dissertation: Difference Methods for Shock Waves.
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. https://doi.org/10.1061/(ASCE)10840699(2008)13:9(835)Horst, M.G., Shamutalov, S.S., Pereira, L.S., Gonçalves, J.M., 2005. Field assessment of thewater saving potential with furrow irrigation in Fergana, Aral Sea basin. Agric. Water
Manag. 77, 210–231. https://doi.org/10.1016/j.agwat.2004.09.041
Katopodes, N.D., Strelkoff, T.S., 1977. DIMENSIONLESS SOLUTIONS OF BORDERIRRIGATION ADVANCE. J. Irrig. Drain. Eng. 103, 401–417.
Lal, A.M.W., 1998. Weighted Implicit Finite-Volume Model for Overland Flow. J. Hydraul. Eng.
124, 941–950. https://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.
https://doi.org/10.1061/(ASCE)0733-9437(2009)135:1(7)LeVeque, R.J., 2002. Finite Volume Methods for Hyperbolic Problems. Cambridge Univ. Press
54, 258. https://doi.org/10.1017/CBO9780511791253
Oweis, T.Y., Walker, W.R., 1990. Zero-inertia model for surge flow furrow irrigation. Irrig. Sci.
11, 131–136. https://doi.org/10.1007/BF00189449
Perea, H., Strelkoff, T.S., Adamsen, F.J., Hunsaker, D.J., Clemmens, A.J., 2010. Nonuniformand Unsteady Solute Transport in Furrow Irrigation. I: Model Development. J. Irrig. Drain.
Eng. 136, 365–375. https://doi.org/10.1061/(ASCE)IR.1943-4774.0000106
Šimůnek, J., Bristow, K.L., Helalia, S.A., Siyal, A.A., 2016a. The effect of different fertigation
strategies and furrow surface treatments on plant water and nitrogen use. Irrig. Sci. 34, 53–
69. https://doi.org/10.1007/s00271-015-0487-z
Šimůnek, J., Šejna, M., van Genuchten, M.T., 2018. New features of version 3 of the HYDRUS(2D/3D) computer software package. J. Hydrol. Hydromechanics 66, 133–142.
https://doi.org/10.1515/johh-2017-0050
Š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.
https://doi.org/10.2136/vzj2016.04.0033
Šimůnek, J., van Genuchten, M.T., Šejna, M., 2008. Development and Applications of the
HYDRUS and STANMOD Software Packages and Related Codes. Vadose Zo. J. 7, 587.
https://doi.org/10.2136/vzj2007.0077
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 Manag. 115, 242–251. https://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-flow equation. Agric. Water Manag. 116,
160–174. https://doi.org/10.1016/j.agwat.2012.07.008
Strelkoff, T.S., Clemmens, A.J., Perea-Estrada, H., 2006. Calculation of non-reactive chemical
distribution in surface fertigation. Agric. Water Manag. 86, 93–101.

42https://doi.org/10.1016/j.agwat.2006.07.006
Strelkoff, T.S., Katopodes, N.D., 1977. BORDER-IRRIGATION HYDRAULICS WITH ZEROINERTIA. 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 infiltration. Agric. Water Manag. 28, 201–221.
https://doi.org/10.1016/0378-3774(95)01177-K
van Genuchten, M.T., 1980. A Closed-form Equation for Predicting the Hydraulic ConductivityofUnsaturated
Soils1.
Soil
Sci.
Soc.
Am.
J.
44,
892–898.
https://doi.org/10.2136/sssaj1980.03615995004400050002x
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., 1976. On the accuracy and efficiency of several numerical schemes forsolving the convective-dispersive equation, in: Gray, W.. (Ed.), Finite Elements in Water
Resources. Pentech Press, London, p. 1.71-1.90.
van Leer, B., 1979. Towards the ultimate conservative difference scheme. V. A second-ordersequel to Godunov’s method. J. Comput. Phys. 32, 101–136. https://doi.org/10.1016/00219991(79)90145-1
van Leer, B., 1974. Towards the ultimate conservative difference scheme. II. Monotonicity andconservation combined in a second-order scheme. J. Comput. Phys. 14, 361–370.
https://doi.org/10.1016/0021-9991(74)90019-9
Walker, W.R., Humpherys, A.S., 1983. Kinematic‐W ave Furrow Irrigation M odel. J. Irrig.
Drain. Eng. 109, 377–392. https://doi.org/10.1061/(ASCE)0733-9437(1983)109:4(377)
Wöhling, T., Fröhner, A., Schmitz, G.H., Liedl, R., 2006. Efficient Solution of the Coupled OneDimensional Surface—Two-Dimensional Subsurface Flow during Furrow Irrigation
Advance. J. Irrig. Drain. Eng. 132, 380–388. https://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 Flow and Plant Water Uptake in Irrigation Furrows. II: Model Test and
Evaluation. J. Irrig. Drain. Eng. 133, 548–558. https://doi.org/10.1061/(ASCE)07339437(2007)133:6(548)
Wöhling, T., Schmitz, G.H., 2007. Physically based coupled model for simulating 1D surface–
2D subsurface flow and plant water uptake in irrigation furrows. I: Model development. …
Irrig.
Drain.
Eng.
133,
538–547.
https://doi.org/10.1061/(ASCE)07339437(2007)133:6(538)
Wöhling, T., Singh, R., Schmitz, G.H., 2004. Physically Based Modeling of Interacting Surface–
Subsurface Flow During Furrow Irrigation Advance. J. Irrig. Drain. Eng. 130, 349–356.
https://doi.org/10.1061/(ASCE)0733-9437(2004)130:5(349)
Zerihun, D., Furman, a., Warrick, a. W., Sanchez, C. a., 2005. Coupled Surface-Subsurface
Flow Model for Improved Basin Irrigation Management. J. Irrig. Drain. Eng. 131, 111–128.
https://doi.org/10.1061/(ASCE)0733-9437(2005)131:2(111)43