Agricultural Water Management 213 (2019) 792–802

Contents lists available at ScienceDirect

Agricultural Water Management
journal homepage: www.elsevier.com/locate/agwat

A coupled model for simulating water flow and solute transport in furrow
irrigationT

Kun Liua,b, Guanhua Huanga,b, , Xu Xua,b, Yunwu Xionga,b, Quanzhong Huanga,b, Jiří Šimůnekc
⁎aCenter for Agricultural Water Research, China Agricultural University, Beijing, 100083, China
Chinese-Israeli International Center for Research and Training in Agriculture, China Agricultural University, Beijing, 100083, Chinac
Department of Environmental Sciences, University of California Riverside, Riverside, CA 92 521, USAb

ARTICLE INFOABSTRACT

Keywords:
Surface irrigation
Solute transport
Mathematical model
Coupling procedure

For optimal water and fertilizer management under furrow irrigation, it is important to understand the waterand solute dynamics on the land surface and in the subsurface. An efficient mathematical tool is required todescribe these dynamic processes. We propose a coupled model in which surface water flow and solute transport
are described using the zero-inertia equation and the average cross-sectional convection-dispersion equation,
respectively, while the two-dimensional Richards equation and the convection-dispersion equation are used tosimulate water flow and solute transport in soils, respectively. Solutions are computed numerically using finite
differences for surface water flow and finite volumes for solute transports in furrow. Subsurface water flow andsolute transport equations are solved using the CHAIN_2D code. An iterative method is used to couple computations of surface and subsurface processes. Both surface and subsurface water flow and solute transportmodules are coded in program subroutines and functions in the Intel FORTRAN environment. The coupled model
was validated by comparing its simulation results with measured data. Results showed that simulated water front
advances in the furrow and water contents in the soil agreed with the observations reasonably well. Good
simulations can be achieved with a relatively fine temporal resolution. Numerical oscillations can be eliminatedby adopting appropriate time steps. As compared with the traditional furrow irrigation model, the proposedmodel can better quantify soil water and solute dynamics by considering interactions between surface andsubsurface water flow and solute transport processes. The proposed model can be used as a decision tool todesign and manage furrow irrigation.

1. Introduction

required (Crevoisier et al., 2008).
Many field experiments have been conducted to provide management recommendations for increasing the uniformity of fertilizer applications (Mailhol et al., 2001; Halvorson et al., 2002; Abbasi et al.,2003a; Adamsen et al., 2005) and controlling salinity (Evans et al.,
1990; Moreno et al., 1995; Malash et al., 2008; Nagaz et al., 2013; Chen
and Feng, 2013). However, field experiments are time-consuming andexpensive and cannot include all relevant factors. Also, it is quite difficult to obtain general conclusions from field experiments due to thevariability of soil properties, soil infiltration rates, field gradients, soilroughness, etc. (Ebrahimian et al., 2014). Mathematical models canprovide an alternative approach to overcome these limitations.
In recent years, many models have been developed for quantifying
water flow and solute transport both at the soil surface and in thesubsurface under furrow irrigation. These models can be divided into
three categories. The first group of models describes water flow in

Furrow irrigation has been widely adopted for wide-spaced crops
such as corn, cotton, and sunflower all over the world because of its low
investment and low energy requirements (Chen and Feng, 2013). The
application of fertilizers with irrigation water is a common practice infurrow irrigation (Ebrahimian et al., 2014). To evaluate the efficiency
of fertilizer applications, a knowledge of the distribution of the fertilizer
in the soil along the furrow is required. Additionally, poor quality
waters, such as saline waters, have been extensively used for irrigationin arid and semi-arid areas, e.g., in the northwest of China, the Negev
region of Israel, and northwest India (e.g., Oster, 1994; Chen et al.,2015; Mantell et al., 1985; Rajinder, 2004). Improper use of saline
waters for irrigation may cause secondary salinization problems. To
study environmental impacts of irrigation with saline waters, an accurate estimate of water flow and salt fluxes into the soil profile is also⁎

Corresponding author at: Center for Agricultural Water Research, China Agricultural University, Beijing, 100083, China.
E-mail address: ghuang@cau.edu.cn (G. Huang).

https://doi.org/10.1016/j.agwat.2018.11.024
Received 18 April 2018; Received in revised form 13 October 2018; Accepted 23 November 20180378-3774/ © 2018 Elsevier B.V. All rights reserved.

Agricultural Water Management 213 (2019) 792–802

K. Liu et al.

furrows using either zero-inertia, kinematic-wave, or fully hydrodynamic equations and characterizes solute transport in furrows using
either advection or advection-dispersion equations, while simple empirical functions are used to represent water and solute infiltration.
Based on this modeling concept, many software packages such as
SIRMOD (Walker, 2003) and WinSRFR (Bautista et al., 2009) have beendeveloped for the design and management of surface irrigation. Furthermore, the software package SIFUM (Boldt et al., 1994) was developed for fertilization management under surge irrigation. However, theabove-discussed models cannot describe distributions of water and solutes in soils in detail, which is important for the evaluation of theuniformity of water and fertilizer applications under furrow irrigationand the quantification of deep percolation and leaching. Also, parameters of empirical infiltration functions are event specific and cannot
account, for example, for the effects of the antecedent soil moisture on
infiltration (Zerihun et al., 2005). Thus, the use of empirical infiltration
functions may result in an inaccurate estimation of infiltration.
The second group of models describes subsurface processes, i.e.,water flow and solute transport in soils, using two-dimensional variably-saturated flow and advection-dispersion equations, respectively,
while surface processes are simplified and treated as an upper boundary
condition. Among these models, the HYDRUS (2D/3D) software
package (mainly its two-dimensional version denoted as HYDRUS-2D(Šimůnek et al., 2008, 2016b) with appropriate boundary conditions
has been widely used to evaluate water and solute dynamics in soils
under furrow irrigation. Šimůnek et al. (2016b) identified some 25
papers in which HYDRUS-2D and its predecessors, such as SWMS_2D
and CHAIN_2D, have been used to simulate water flow and solutetransport in furrow irrigation systems. For instance, Crevoisier et al.
(2008) used HYDRUS-2D to simulate water and nitrogen dynamics
under alternate and conventional furrow irrigation. Siyal et al. (2012)and Šimůnek et al. (2016a) developed a new furrow boundary condition to better account for flow and transport processes in furrows.
Šimůnek et al. (2016a) evaluated the effects of different fertigation
strategies and furrow surface treatments on plant water and nitrogenuse. Wang et al. (2014) coupled HYDRUS-2D with the EPIC crop growthmodule to evaluate soil water dynamics and melon growth under
furrow irrigation. However, these types of models, cannot fully describe
the dynamic processes of water advancement, storage, and recession
during a furrow irrigation event. Additionally, interactions between
surface and subsurface processes have not been fully considered in assessing model performance. These simplifications may produce inaccurate simulation results.
The third group of models, which couples two models fully describing the surface and subsurface processes and their interactions,represent an alternative way to simulate furrow irrigation systems.
Some efforts to couple surface and subsurface models for furrow irrigation have been reported in the literature. For example, Wöhling andSchmitz (2007) coupled the zero-inertia model, HYDRUS-2D, and thecrop model to simulate water flow and crop growth under furrow irrigation. However, they did not consider solute transport (e.g., fertilizers and salts) and their effects on crop growth. Ebrahimian et al.
(2013) loosely coupled the surface model of Abbasi et al. (2003b) with
HYDRUS-2D to simulate water and solute dynamics under furrow irrigation. However, their coupling procedure was not very accurate since
it did not consider the mutual feedback between surface and subsurface
systems. Therefore, more accurate coupling methods, i.e., the iterative
coupling and full coupling methods, are required to develop more reliable furrow irrigation models (e.g., Furman, 2008).
The objectives of this study are (1) to develop a coupled model fornumerically simulating water and solute dynamics under furrow irrigation by an iterative coupling of physically-based surface and subsurface models, (2) to validate the coupled model by comparing its simulation results with available experimental data, and (3) to conduct
numerical tests evaluating the accuracy and efficiency of the coupledmodel.

2. Theory and method
2.1. Model description
2.1.1. Surface flow
Water flow in a furrow can be described using the zero-inertia
equation, which is a simplified form of the full hydrodynamic equation
obtained by neglecting the inertial and acceleration terms (Elliott et al.,1982). The mathematical expression of the zero-inertia equation is:

A
Q
Z
+
+
=0
t
x
t
h
= S0
x

(1)

Sf

(2)
2

3 −1

where A is the cross-sectional flow area (m ), Q is the flow rate (m s ),Z is the infiltration volume per unit length of furrow (m2), h is the flow
depth (m), S0 is the bottom slope of the furrow (-), Sf is the friction slope(-), x is the horizontal distance (m), and t is time (s). The friction slope isdefined using the Manning equation:

Sf =Q 2N 2
cu2 A2 R4 3

(3)
1/6

where N is the Manning’s roughness (m ), R is the hydraulic radius(m), and cu is a dimensional constant equal to 1 when units (m1/2 s−1)
are used.
For a furrow of an arbitrary cross-section, the flow depth can be
expressed as a power function of the flow area (Elliott et al., 1982):(4)

h = bAp
1/2p

where b and p are furrow geometry parameters in units of (m) and(-), respectively. Similarly, the wetted perimeter can be written as a
power function of the flow area. Hence, the denominator of Eq. (3)becomes a function of A:(5)

A2 R 4 3 = dA f

5.33−2f

where d and f are hydraulic section parameters in units of (m)
and (-), respectively. Substituting Eqs. (4) and (5) into Eq. (2), one
obtains(bAp )= S0x

Q 2N 2dAf(6)

Eqs. (1) and (6) are the basis of the surface flow model. The initialand boundary conditions are as follows:L(7)

Q (0, t ) = Qo (t ) t tco
Upstream boundary condition: Q (xr , t ) = 0
0 xr
A (x r , t ) = 0
tro < tL and t > tco
trL

Initial condition: Q (x , 0) = 0; A (x , t ) = 0 0x(8)
Downstream boundary condition:
Q (x a , t ) = 0; A (x a , t ) = 0 0 xa
Q (L , t ) = 0
t > taL

A (L , t ) =Q2 (L, t ) N2
S0 d

1 f

t

Landt

taL

for a blocked-end furrow,taL for a free-draining furrow

(9)
3

−1

where L is the furrow length (m), Q0 is the inflow rate (m s ), tco isthe cutoff time (s), xr is the distance of the receding tip from the upstream end (m), tro is the recession time at the upstream end (s), trL is
the recession time at the downstream end (s), xa is the advance distance(m), and taL is the advance time (s).
2.1.2. Solute transport in the furrow
Solute transport in a furrow can be described using a one-dimensional cross-sectional average dispersion equation (Cunge et al., 1980):
793

Agricultural Water Management 213 (2019) 792–802

K. Liu et al.

AC
QC
ZC
C
+
+
=
ADx
t
x
t
x
x

(10)

where C is the solute concentration (kg m−3) and Dx is the dispersion
coefficient (m2 s-1). The dispersion coefficient accounts for the dispersion due to differential advection and turbulent diffusion (Cunge et al.,1980). Dx can be written as (Abbasi et al., 2003b):

Dx =

(11)

x v + Ddwhere αx is the longitudinal dispersivity (m), v is the mean cross-sectional velocity (m s−1), and Dd is the molecular diffusion in free water(m2 s−1). The initial and boundary conditions are as follows:

Initial condition: C (x , 0) = 0 0

x

(12)

L

Upstream boundary condition: C (0, t ) = Co (t ) tC=0x

Downstream boundary condition:

0

xa

(13)

tco
Landt

trL

(14)

Fig. 1. Discretization of the computational domain for the numerical solution of
the zero-inertia equation.

2.1.3. Subsurface flow
Subsurface flow in any cross-section perpendicular to the directionof surface flow can be described by using the two-dimensional Richards
equation:t=

K K ijAxih+ K izA
xj

S

component of the volumetric flux (m s-1), Dij is the dispersion
coefficient tensor (m2 s-1), and cs is the concentration of the sink term
(kg m−3).
The initial conditions for concentrations can either be constant orvary with the depth and distance along the furrow. The first-typeboundary condition reflecting the solute concentration in the furrow is
considered as the upper boundary condition. The third-type boundary
condition with the spatial derivative of the concentration equal to zerois considered as the bottom boundary condition.

(15)

where θ is the volumetric moisture content (-), h is the pressure head(m), K is the unsaturated hydraulic conductivity (m s−1), K ijA are
components of a dimensionless anisotropy tensor KA, xi and xj are thespatial coordinates (m), and S is the sink term accounting for root water
uptake.
The unsaturated hydraulic conductivity function is determinedusing the van Genuchten-Mualem model (Mualem, 1976; van
Genuchten, 1980):

K (h) = Ks Se0.5 [1
Se =

m=1

r
s

(1

2.2. Numerical solution
2.2.1. Numerical solutions for the surface and subsurface models
The zero-inertia equations, i.e., Eqs. (1) and (6), are solved using theimplicit finite difference method. The discretization of the computational domain for Eqs. (1) and (6) is shown in Fig. 1. The implicit finite
difference equations can be found in Appendix A. Whereas the crosssectional average dispersion Eq. (10) is solved using the finite volume
method, and the resulting linear equations can also be found in
Appendix A.
The subsurface flow and transport module is based on the
CHAIN_2D code (Šimůnek et al., 2008), i.e., an open source version of
HYDRUS-2D. Both the two-dimensional subsurface water flow and solute transport equations subject to appropriate initial and boundary
conditions were solved using the Galerkin finite element method.

(16)

Se1 m ) m]2= [1 + ( h) n] m

(17)

r

(18)

1 n
−1

where Ks is the saturated hydraulic conductivity (m s ), Se is the effective water saturation (-), θr and θs are the residual and saturated
water contents, respectively (-), α is the inverse of the air-entry value(m-1), and n is the pore-size distribution index (-).
The initial conditions are given in terms of pressure heads, which
can be either constant or vary with the depth and distance along the
furrow. The time-variable pressure head representing the flow depth in
the furrow is considered as the upper boundary condition in the furrow,
while an atmospheric boundary condition can be considered on the restof the soil surface, including in the furrow when it is empty. A freedrainage is considered as the bottom boundary condition in cases when
the groundwater table is below the domain of interest. Otherwise, theDirichlet boundary condition with a prescribed pressure head can beconsidered as the bottom boundary condition.

2.2.2. Calculation procedures for the coupled model
The coupled model is coded in program subroutines and functions,which integrate water flow and solute transport in the furrow and
CHAIN_2D (Fig. 2). The model is written in FORTRAN 90. As shown in
Fig. 2, the calculation procedures are as follows:
1.First, the parameters of the surface water flow and solute transport modules are initialized and the initial input files of CHAIN_2D are
read.
2.In each time step, the surface flow module is iteratively coupledwith the CHAIN_2D water flow module. First, the flow water depth ineach discrete node along the furrow is calculated using an estimate of
the infiltration amount. The infiltration amount is assumed to be equal
to that at the preceding time. Second, CHAIN_2D calculates the distribution of pressure heads in the soil profile and the cumulative infiltration amount during one time step for each node using the flow
depth as the upper boundary condition. Third, the simulated cumulative infiltration amount is transferred back to the surface flow module
to update the flow depth. These last two steps are repeated until the

2.1.4. Solute transport in subsurface
Solute transport in any subsurface cross-section perpendicular to thedirection of surface flow can be described by using the two-dimensional
convection-dispersion equation:

c
=
t
xi

Dij

c
xj

qi c
xi

Sc s

where c is the solute concentration in the soil (kg m

(19)
−3

), qi is the i-th794

Agricultural Water Management 213 (2019) 792–802

K. Liu et al.

Fig. 2. Flow chart of the calculation procedure of the coupled model.

convergence criterion is fulfilled. The convergence criterion (0.001 m2)
is given as the difference in the cumulative infiltration volumes along
the furrow between two adjacent iterations.
3.After the convergence criterion is satisfied, the flow modules passselected variables (i.e., flow rates, areas, infiltrations, soil water contents,and soil water fluxes) to the solute modules to calculate solute transport
processes. Unlike in water flow calculations, the interaction between thesurface and subsurface solute modules is unidirectional. The surface solute module first calculates solute concentrations for each discrete node
along the furrow. The CHAIN_2D solute module then calculates the distribution of solute concentrations in the soil profiles using solute concentrations in the furrow as the upper boundary conditions.
4.The results, including flow rates and areas and solute concentrations in the furrow as well as the distribution of pressure heads andconcentrations in the soil profiles, are printed after the completion ofsolute calculations.
5.Steps (2) to (4) are repeated until the entire irrigation process issimulated.

Table 1
Parameters of the field experiment of Abbasi et al. (2003b).
Parameter

Value

Field length (m)Field slope (m m−1)
Furrow geometry parametersd (m5.33−2f)f (-)
b (m1/2p)
p (-)
Infiltration parametersa (-)
k (m3 min−a m-1)
c (m3 m−1)
Manning roughness (m1/6)
Cutoff time (min)Solute application time (min)Dispersivity (cm)115
0.0001
0.326
2.789
0.871
0.635
0.75
0.00209
0.006
0.066
140
140
10

Note: adapted from Abbasi et al. (2003b).

2.3. Model validation

The numerical solution for the solute transport model in the furrow
was also evaluated using the analytical solution. For a steady uniform
flow condition and given the initial condition (C(x,0) = 0, for 0 ≤ x ≤L) and boundary conditions (C(0, t) = C0, C(L, t) = 0, for 0 ≤ t ≤ T),
the analytical solution can be expressed as (Warrick, 2003):

2.3.1. Validation of the surface flow and solute transport models
The numerical solution of the zero-inertia equation was validated
against the field experimental data of Abbasi et al. (2003b) for thewater advance and recession phases using the parameters given in
Table 1. The numerical solution of the solute transport in the furrow
was also validated by comparing simulation results with the field experimental data of Abbasi et al. (2003b). The inflow rate of 1.29 L s−1
and the applied solute concentration of 2.36 g L−1 were considered in
the simulation.

C (x , t ) =

C0
x vt
erfc
2
2 Dx t+ e vx Dx erfc

x + vt2 Dx t(20)

where C0 is the inlet concentration, L is the length of the computation
domain, and T is the total time (Warrick, 2003). A hypothetical795

Agricultural Water Management 213 (2019) 792–802

K. Liu et al.

Table 2
Parameters of three irrigation laboratory experiments of Wöhling et al. (2004).
Parameters
Average inflow rate (L s
Cutoff time (h)Furrow length (m)Furrow slope (m m−1)

−1)

Table 4
Parameters for the time resolution tests.

1st irrigation

2nd irrigation

3rd irrigation

Parameter

Value

1.20
2.55
26.4
0.0025

2.14
5.37
26.4
0.0014

2.58
3.62
26.4
0.0015

Field length (m)Field slope (m m−1)
Geometry parametersdfbpManning roughness (m1/6)
Inflow rate (L s−1)
Applied solute concentration (g L−1)
Cutoff time (min)Solute application time (min)Dispersivity (cm)100
0.0015

Note: adapted from Wöhling et al. (2004).

example with the parameter values C0 = 1, L = 1, T = 1, v = 1, and
Dx = 0.1 was set up for this evaluation.
2.3.2. Validation of the coupled model
The coupled model was validated using the experimental data of
Wöhling et al. (2004), who carried out laboratory experiments to investigate surface and subsurface flow patterns during the advancephase of furrow irrigation. A parabolic furrow, with a top width of
0.35 m and a depth of 0.184 m, was formed in an experimental tank.
Fifty tensiometer probes were placed at five cross sections (x = 1.5, 6.3,
12.3, 18.3, and 24.3 m) to measure pressure heads. The soil hydraulicparameters were optimized by a MATLAB routine using the pressure
head measurements at the cross section x = 12.3 m. The resulting
optimized parameters of the van Genuchten-Mualem model were
α = 1.4 m−1, n = 1.25, Ks = 1.95 × 10-5m s-1, θs = 0.41, andθr = 0.14. More details about the experiments can be found in Wöhling
et al. (2004). Three irrigation experiments were conducted with different inflow rates and bottom slopes. The detailed experimental information is shown in Table 2. The initial conditions for the subsurface
flow model were determined by averaging measured pressure heads atfive cross sections, which are shown in Table 3. The simulated water
advance by the coupled model was compared with the measured data
and with that simulated by the Wöhling’s model (Wöhling et al., 2004),
in which the analytical solution of the zero-inertia equation was used todescribe surface water flow. Simulated pressure heads in different crosssections (i.e., x = 6.3 m and x = 18.3 m) were also compared with themeasurements.
The coefficient of determination (R2) and the root mean squareerror (RMSE) were used as indicators to evaluate the performance ofthe model validation.

0.382
2.767
1.150
0.564
0.03
2.0
6.3
60
60
10

only when they fulfill certain stability criteria, often involving limits on
the spatial and temporal discretization, i.e., spatial and time steps.
Generally, the accuracy of simulation results depends greatly on thetime and space resolution. In our proposed model, a constant time stepis used in the numerical solution of the zero-inertia equation. The space
increment is determined by the difference of the calculated water front
positions at two adjacent times. It should be noted that the space increment is not constant during the simulations because the velocity of
the water advance decreases during the advance phase (see Fig. 1). In
order to analyze the impact of the time resolution, numerical tests wereconducted for time steps Δt = 30 s, 60 s, and 120 s A hypothetical casewas set up for these numerical tests with its parameters given in
Table 4. In this case, the sandy soil was considered and its hydraulicparameters were adopted from Vogel and Hopmans (1992). Initial soilpressure heads and solute concentrations were assumed to be uniformly
distributed with values of -2.5 m and 0 g m−3, respectively. The simulation results for different time steps were compared with each other.
The total volume balance errors, δ, were also computed for different
time steps:

=

|Vin

Vtotal |× 100%
Vin(21)

where Vtotal is the sum of the surface water and infiltration volumes(m3) and Vin is the water volume applied into the furrow (m3).
3. Results and discussion

2.4. Numerical tests

3.1. Simulation results of the surface flow and solute transport models

In the proposed coupled model, all governing equations were solved
using a numerical solution. Numerical solutions can generally handle
more complex initial and boundary conditions, as well as heterogeneities in soil properties, than analytical solutions. However, numerical solutions often encounter numerical problems and are stable

Fig. 3 shows the comparison of the water advance and recession
measured and simulated using the numerical zero-inertia model. It canTable 3
Average initial pressure heads at different locations used in the subsurface flow
model.
Soil depth (m)h (m)1st irrigationh (m)2nd irrigationh (m)3rd irrigation

0.2
0.3
0.35
0.45
0.5
0.6
0.7
0.8
1.05

−4.37
−3.33
−2.83
−2.41
−2.25
−2.02
−1.97
−2.07
−3.46

−6.65
−5.45
−6.02
−3.94
−4.09
−4.22
−3.61
−2.72
−3.43

−1.74
−1.45
−1.58
−1.16
−1.19
−1.06
−0.86
−0.64
−0.68

Fig. 3. Comparison of the measured (symbols) and simulated water advance(solid line) and recession (dashed line) using the zero-inertia model.

Note: adapted from Wöhling et al. (2004).
796

Agricultural Water Management 213 (2019) 792–802

K. Liu et al.

variations in the inflow rate, resulting in unstable solute concentrationsin the irrigation water.
Fig. 5 presents the comparison between analytical and numericalsolutions for the cross-sectional average dispersion equation. The concentrations at different locations simulated using the numerical solution
were in very good agreement with those obtained using the analyticalsolution. This indicates that the numerical solution has relatively highaccuracy.
3.2. Performance of the coupled model
A comparison of the water advance measured and simulated using
the coupled model is presented in Fig. 6. It can be seen that theFig. 4. Comparison of measured and simulated solute concentrations at a distance of 10 (top) and 60 (bottom) m from the furrow inlet.

be seen that the recession time was approximately the same along the
furrow. This may be attributed to a small field slope and spatially relatively uniform infiltration. The simulated water advance and recession agree reasonably well with the measurements, indicating that thenumerical zero-inertia model has good accuracy. The differences between the simulations and measurements may be due to the spatial andtemporal variability of infiltration and roughness, and the nonuniformity of the field slope along the furrow (Abbasi et al., 2003b).
The simulated solute concentrations at the distances of 10 and 60 m
from the inlet were compared with the measurements in Fig. 4. The
solute concentration increased rapidly with the water advance and thenreached a constant value. In general, the simulated solute concentrations were in good agreement with the measurements. Deviations between early simulated and measured values were mainly due toFig. 6. Comparison of the water advance measured and simulated using the
new coupled model and the Wöhling’s model (Wöhling et al., 2004) for threefurrow irrigation experiments.

Fig. 5. Comparison of concentrations obtained using either the analytical or
numerical solutions of the cross-sectional average dispersion equation.
797

Agricultural Water Management 213 (2019) 792–802

K. Liu et al.

Fig. 7. Measured and simulated distributions of the pressure head in the soil profiles at distances of 6.3 (left) and 18.3 (right) m from the furrow inlet after 20 min
irrigation time.

measured and simulated water advance values matched reasonably
well, as reflected by the coefficient of determination R2 of 0.99 and the
RMSE of 1.02–2.73 m, implying that the coupled model has reasonableaccuracy. Differences between simulations and observations may be
due to the spatial variability of infiltration and roughness along the
furrow. While for the first irrigation, the coupled model predicted aslightly slower water advance than measured, for the second irrigation,
it predicted a faster water advance. This may be attributed to a significant increase in the initial infiltration rate caused by soil cracks thatdeveloped between the first and second irrigations (Wöhling et al.,2004).
Water advances simulated by both the coupled and Wöhling’s
models (Wöhling et al., 2004) are also presented in Fig. 6 for comparison. The water advance simulated by the coupled model was slightly
slower than simulated by the Wöhling’s model for the first irrigationand slightly faster for the second and third irrigations. This could be
explained by the fact that the Wöhling’s model used the analytical solution of the zero-inertia equation, whereas the numerical solution wasused in the present work. The analytical solution of Wöhling assumed
that the momentum could be expressed as the momentum at the centerof gravity of the water body. This assumption was not needed in our
numerical solution. In addition, the iterative convergence criteria usedin our model was different from that of Wöhling’s model. This might
result in different simulation results.

Fig. 9. Simulated cumulative infiltration along the furrow during the advancephase for time steps Δt = 30, 60, and 120 s.

Fig. 7 shows the simulated and measured pressure head distributions in soils at the cross-sections 6.3 and 18.3 m from the furrow inlet
after 20 min irrigation time. In general, the simulations were in good
agreement with the measurements. The discrepancy between theFig. 10. Simulated solute concentrations along the furrow after 30 min for time
steps Δt = 30, 60, and 120 s.

Fig. 8. Simulated water advance for time steps Δt = 30, 60, and 120 s.
798

Agricultural Water Management 213 (2019) 792–802

K. Liu et al.

simulated and measured results may be attributed to the approximationof the initial conditions for the simulation. Similarly to Wöhling et al.
(2004), the initial pressure heads over the simulation depth were obtained by interpolating measured pressure heads at nine depths (seeTable 3).

3.3. Impact of the time resolution
The simulated water advances for different time step sizes Δt arepresented in Fig. 8, which shows that the water advance was similar fortime steps Δt = 30 and 60 s, whereas it was faster for Δt = 120 s.

Fig. 11. Simulated pressure heads (top) and solute concentrations (bottom) in the soil profile at a distance of 50 m form the furrow inlet after 30 min irrigation time
for time steps Δt = 30 (up), 60 (middle), and 120 (down) s.

799

Agricultural Water Management 213 (2019) 792–802

K. Liu et al.

Infiltration for Δt = 120 s was underestimated in comparison withthose for Δt = 30 and 60 s during the water advance process (seeFig. 9), resulting in a faster water advance. The total volume balanceerror (Eq. 26) showed a significant increase from 0.8% for Δt = 30 s to
3.2% for Δt = 120 s. This indicates that larger time steps and space sizes(the space increment increases with the time step in the coupled model)
may cause an inaccurate evaluation of infiltration and the water advance time. These results are consistent with those of Wöhling et al.
(2004), who reported that the large space sizes could lead to the biased
estimation of infiltration and advance time. Banti et al. (2011) also
indicated that smaller space increments could result in lower volume
balance errors.
Fig. 10 presents simulated solute concentrations along the furrow
after 30 min irrigation time for different time steps. Significant numerical oscillations could be found in the simulation with Δt = 120 s,
whereas there were no visible numerical oscillations in the simulation
with Δt = 30 s. This may be attributed to the fact that larger time steps
resulted in larger Peclet and Courant numbers. In our numerical tests,
Peclet numbers were less than 8.0 for Δt = 30 s and less than 18.0 for
Δt = 120 s. Corresponding Courant numbers were less than 1.5 for Δt= 30 s and less than 1.7 for Δt = 120 s. Similar results were also foundby Abbasi et al. (2003b), who reported that no significant numericaloscillations could be found for simulations with Peclet numbers smaller
than 10 and Courant numbers not larger than 1. Although the Courant
number is larger than 1 in our simulations, a Courant number smaller
than 2 is acceptable as the average dispersion equation (Eq. 10) is
solved using an implicit scheme. It should be noted that Fig. 10 shows adecrease of solute concentrations below c0 at the outflow end. This is
the result of the numerical dispersion and not a physical phenomenon
since there is no physical process in our hypothetical example that
could lower concentrations below c0.
The simulated pressure heads and solute concentrations in the soil
profile after 30 min irrigation time for different time steps are shown in
Fig. 11. It can be seen that simulated pressure head distributions weresimilar for Δt = 30 and 60 s. The wetted area for Δt = 120 s is slightly
smaller than for Δt = 30 and 60 s. This further demonstrates that using
the time step of 120 s may underestimate infiltration. Simulated soluteconcentrations for Δt = 120 s are smaller than those for Δt = 30 and
60 s. This is also attributed to the fact that the infiltration volumes for
Δt = 120 s are smaller.
The required CPU time for different time steps is presented inFig. 12, which shows that the CPU time increased significantly with thedecrease of the time step. The CPU time for Δt = 30 s was 2.5 times
larger than for Δt = 60 s and 9.3 times larger than for Δt = 120 s. Thiswas mainly due to the numerical solution of the two-dimensional Richards equation and the convection-dispersion equation, i.e., for subsurface processes. Using smaller time steps resulted in having moresubsurface profiles perpendicular to the direction of surface flow, requiring many more computations.

There are many ways to improve the computational efficiency. First,
the numerical code can be parallelized to significantly reduce thecomputational time when the code is run on multiple-core processors or
multiple processors. Second, if one only needs to know the cumulative
infiltration amount rather than detailed water and solute distributions
in soils, a simple infiltration model can be used instead of the sophisticated two-dimensional process model. For example, Bautista et al.
(2016) proposed an approximate furrow infiltration model by treating
the two-dimensional infiltration as a sum of a one-dimensional infiltration and edge effects. This model could obtain comparably accurate estimations of cumulative infiltration as the numerical solution of
the two-dimensional Richards equation. Similarly, an empirical surface
flow function can be used instead of a physically based model.
Ghanbarian et al. (2018) provided theoretical insights into the exponent in the empirical power-law advance-time curve, which is widelyused in surface irrigation. However, improving the computational efficiency may not be a big concern in the current situation when significant computational resources are often readily available.
It should be pointed out that the proposed model is highly computationally efficient. This is because the subsurface water flow and solutetransport modules in the coupled program are called directly from theoverland flow module, rather than passing various variables throughexternal HYDRUS files (as done in Wöhling et al., 2004; Zerihun et al.,2005; Wöhling and Schmitz, 2007). In addition, the modular programming method makes it easy to optimize each part of the model
separately. The proposed model can simulate the water and solute dynamics of all phases of a furrow irrigation event by considering interactions between surface and subsurface water flow and solute transport
processes. It can also provide a more detailed distribution of both waterand solute in soils in comparison with traditional furrow irrigation
models. Hence the proposed model is a useful tool for water and fertilizer management under furrow irrigation.
In this study, crop growth processes were not considered in the
coupled model. However, the effect of irrigation and fertigation on crop
growth using different irrigation techniques is always an importantfactor for water management at different scales. Therefore, future research should focus on incorporating the crop growth module into thecoupled model to allow simultaneous simulations of soil water andsolute dynamics, as well as crop growth during the entire growing
season.
4. Summary and conclusions
A coupled model for simulating water flow and solute transport
under furrow irrigation was developed. In the model, surface waterflow and solute transport were described using the zero-inertia equationand the cross-sectional average dispersion equation, respectively. The
two-dimensional Richards and advection-dispersion equations, respectively, were used to simulate water flow and solute transport in soils.
The surface and subsurface models were coupled using an iterative
coupling method and were coded in program subroutines and functionsin the Intel FORTRAN environment. The proposed model was validatedby comparing its simulation results with measured data from laboratory
and field experiments, as well as with other available models.
Numerical tests were conducted to analyze the accuracy and efficiencyof the coupled model. When smaller time steps were used, more robust
simulations without significant numerical oscillations and dispersion
were obtained, albeit at a higher computational cost. The proposedmodel can simulate all phases of an irrigation event and provide detailed quantification of both water and solute dynamics under furrow
irrigation by considering interactions between surface and subsurface
water flow and solute transport processes. The proposed model canserve as a decision tool for the design and management of furrow irrigation.

Fig. 12. CPU times for time steps Δt = 30, 60, and 120 s.
800

Agricultural Water Management 213 (2019) 792–802

K. Liu et al.

Acknowledgments

Program of China (No. 2017YFC0403301) and National Natural Science
Foundation of China (Nos. 51639009, 51621061, 51125036).

This research was jointly supported by the National Key R&D
Appendix A

The implicit finite difference schemes for Eqs. (1) and (6) can be expressed as (Elliott et al., 1982):

Hj Ajk + 1 + Bj Ajk++11 + Cj Qjk + 1 + Dj Qjk++11 = Gj
Hj

Ajk + 1 + Bj

Ajk++11 + Cj

Qjk + 1 + Dj

(A1)

Q jk++11 = Gj

(A2)

where
xj

Hj =

xj

Gj =

xj 1
,
t

+

t

+ (1

t

)

bp (Ajk ) p 1

Hj =

Sfj

Cj = 2

Bj = (1xj 1tf

Gj = b [(Ajk ) p
A = Ak + 1

Sfj, Cj =Ajk, Dj = ,Z jk ) + Q jk

x j , B’j = bp (Ajk+ 1 ) p 1

Akj

Sfj + 1)

Qkj + 1

Sfj

(1

Q k , and Sfj =Q = Qk+1

Q jk+ 1,)f

(1

Sfj + 1
Akj + 1

xj ,

xj ,

(Ajk+ 1 ) p] + xj [S0Ak ,t(Ajk+ 1 + Z jk+ 1

xj , D’j = 2(1

Q kj

xj

)

) Sfj + 1],(Qkj )2N2
d (Akj ) f

,

(A3)

in which j is the node index, k is the timeline index, and φ and ψ are the space and time weighting factors, respectively. In this study, φ and ψ wereconsidered with values of 0.5 and 0.6, respectively.
When the furrow is discretized into N cells, there are 2N linear algebraic equations and 2N unknowns. The linear Eqs. (A1) and (A2), combinedwith the initial and boundary conditions, are solved using the Gaussian elimination method.
Using the finite volume method, the cross-sectional average dispersion Eq. (10) can be decretized as:aW C jk +11 + aP Cjk + 1 + aE C jk++11 = bP

(A4)

where
Qkj 1 + Qkj

aW =
aE =

2
Q kj + 1 + Q kj

aP = Ajk
t

t
2

2

t
2

(

xj 1 + xj
2

Dkj + 1 + Dkj
2•

D kj 1 + D kj
2
Dkj + 1 + Dkj
2

•

)+Z (
k
j

•

Akj 1 + Akj
2

Akj + 1 + Akj

t
,
xj 1

2

t
,
xj

xj 1 + xj
2)+

Q kj + 1 + Q kj 1

Akj + 1 + Akj 1
D kj 1 + D kj Akj 1 + Akj
• 2
+
2
xj
2

and bP = (Ajk C jk + Z jk Cjk )(xj 1 + xj
2

2

t
,
2

1
xj 1

),

(A5)

in which j is the node index and k is the timeline index. Variables A, Z, and Q are obtained from the numerical solution of the zero-inertia equation.
For N cells, there are N linear algebraic equations and N unknowns. Similarly as for water flow, the linear Eq. (A4), combined with the initial andboundary conditions, are solved using the Gaussian elimination method.

applied nitrogen distribution under surge irrigation. Trans. ASAE 37 (4), 1157–1165.
Chen, L.J., Feng, Q., 2013. Soil water and salt distribution under furrow irrigation ofsaline water with plastic mulch on ridge. J. Arid Land 5 (1), 60–70.
Chen, L.J., Feng, Q., Li, F.R., Li, Ch.Sh, 2015. Simulation of soil water and salt transfer
under mulched furrow irrigation with saline water. Geoderma. 241–242 87–96.
Crevoisier, D., Popova, Z., Mailhol, J.C., Ruelle, P., 2008. Assessment and simulation ofwater and nitrogen transfer under furrow irrigation. Agric. Water Manage. 95,
354–366.
Cunge, J.A., Holly, F.M., Verwey, A., 1980. Practical Aspects of Computational River
Hydraulics. Pitman, London.
Elliott, R.L., Walker, W.R., Skogerboe, G.V., 1982. Zero-inertia modeling of furrow irrigation advance. J. Irrig. Drain. Div. 108 (3), 179–195.
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 (3), 301–316.
Ebrahimian, H., Keshavarz, M.R., Playán, E., 2014. Surface fertigation: a review, gaps andneeds. Span. J. Agric. Res. 12 (3), 820–837.
Evans, R.G., Smith, G.J., Oster, J.D., Myers, B.A., 1990. Saline water application effects on
furrow infiltration of red-brown earths. Trans. ASAE 33 (5), 1563–1572.
Furman, A., 2008. Modeling coupled surface–subsurface flow processes: a review. Vadose

References
Abbasi, F., Adamsen, F.J., Hunsaker, D.J., Feyen, J., Shouse, P., van Genuchten, M.T.,
2003a. Effects of flow depth on water flow and solute transport in furrow irrigation:
field data analysis. J. Irrig. Drain. Eng. 129 (4), 237–245.
Abbasi, F., Šimůnek, J., van Genuchten, M.Th., Feyen, J., Adamsen, F.J., Hunsaker, D.J.,
Strelkoff, T.S., Shouse, P., 2003b. Overland water flow and solute transport: model
development and field-data analysis. J. Irrig. Drain. Eng. 129 (2), 71–81.
Adamsen, F.J., Hunsaker, D.J., Perea, H., 2005. Border strip fertigation: effect of injection
strategies on the distribution of bromide. Trans. ASAE 48 (2), 529–540.
Banti, M., Zissis, T., Anastasiadou-Partheniou, E., 2011. Furrow irrigation advance simulation using a surface-subsurface interaction model. J. Irrig. Drain. Eng. 137 (5),304–314.
Bautista, E., Clemmens, A.J., Strelkoff, T.S., Schlegel, J.L., 2009. Modern analysis ofsurface irrigation systems with WINSRFR. Agric. Water Manage. 96, 1146–1154.
Bautista, E., Warrick, A.W., Schlegel, J.L., Thorp, K.R., Hunsaker, D.J., 2016.
Approximate furrow infiltration model for time-variable ponding depth. J. Irrig.
Drain. Eng. 142 (11), 04016045.
Boldt, A.L., Watts, D.G., Eisenhauer, D.E., Schepers, J.S., 1994. Simulation of water801

Agricultural Water Management 213 (2019) 792–802

K. Liu et al.
Zone J. 7 (2), 741–756.
Ghanbarian, B., Ebrahimian, H., Hunt, A.G., van Genuchten, M.Th, 2018. Theoreticalbounds for the exponent in the empirical power-law advancetime curve for surface
flow. Agric. Water Manage. 210, 208–216.
Halvorson, A.D., Follett, R.F., Bartolo, M.E., Schweissing, F.C., 2002. Nitrogen fertilizer
use efficiency of furrow-irrigated onion and corn. Agron. J. 94, 442–449.
Mantell, A., Frenkel, H., Meiri, A., 1985. Drip irrigation of cotton with saline-sodic water.
Irrig. Sci. 6, 95–106.
Mailhol, J.C., Ruelle, P., Nemeth, I., 2001. Impact of fertilization practices on nitrogenleaching under irrigation. Irrig. Sci. 20 (3), 139–147.
Malash, N.M., Flowers, T.J., Ragab, R., 2008. Effect of irrigation methods, management
and salinity of irrigation water on tomato yield, soil moisture and salinity distribution. Irrig. Sci. 26, 313–323.
Moreno, F., Cabrera, F., Andrew, L., Vaz, R., Martin-Aranda, J., Vachaud, G., 1995.
Water-movement and salt leaching in drained and irrigated marsh soils of Southwest
Spain. Agric. Water Manage. 27 (1), 25–44.
Mualem, Y., 1976. A new model for predicting the hydraulic conductivity of unsaturated
porous media. Water Resour. Res. 46 (12), 513–522.
Nagaz, K., Toumi, I., Masmoudi, M.M., Mechlia, N.B., 2013. Comparative effects of drip
and furrow irrigation with saline water on the yield and water use efficiency of potato(solanum tuberosum l.) in arid conditions of Tunisia. Agric. J. 6 (4), 501–510.
Oster, J.D., 1994. Irrigation with poor quality water. Agric. Water Manage. 25 (3),271–297.
Rajinder, S., 2004. Simulations on direct and cyclic use of saline waters for sustaining
cotton–wheat in a semi-arid area of north-west India. Agric. Water Manage. 66,
153–162.
Šimůnek, J., van Genuchten, M.Th., Šejna, M., 2008. Development and applications of the
HYDRUS and STANMOD software packages, and related codes. Vadose Zone J. 7 (2),587–600.

Š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 (1), 53–69.
Šimůnek, J., van Genuchten, M.Th., Šejna, M., 2016b. Recent developments and applications of the HYDRUS computer software packages. Vadose Zone J. 15 (7).
Siyal, A.A., Bristow, K.L., Šimůnek, J., 2012. Minimizing nitrogen leaching from furrow
irrigation through novel fertilizer placement and soil management strategies. Agric.
Water Manage. 115, 242–251.
van Genuchten, M.Th., 1980. A closed-form equation for predicting the hydraulic conductivity of unsaturated soils. Soil Sci. Soc. Am. J. 44, 892–898.
Vogel, T., Hopmans, J.W., 1992. Two-dimensional analysis of furrow infiltration. J. Irrig.
Drain Eng. 118 (5), 791–806.
Walker, W.R., 2003. SIRMOD III—surface irrigation simulation, evaluation and design:user’s guide and technical documentation. Dept. of Biological and Irrigation
Engineering, Utah State Univ., Logan, Utah.
Wang, J., Huang, G., Zhan, H., Mohanty, B.P., Zheng, J., Huang, Q., Xu, X., 2014.
Evaluation of soil water dynamics and crop yield under furrow irrigation with a twodimensional flow and crop growth coupled model. Agric. Water Manage. 141, 10–22.
Warrick, A.W., 2003. Soil water dynamics. Oxford Univ Press, NY.
Wöhling, Th., Singh, R., Schmitz, G.H., 2004. Physically based modeling of interacting
surface-subsurface flow during furrow irrigation advance. J. Irrig. Drain. Eng. 130(5), 349–356.
Wöhling, Th., Schmitz, G.H., 2007. Physically based coupled model for simulating 1D
surface-2D subsurface flow and plant water uptake in irrigation furrows.Ⅰ: model
development. J. Irrig. Drain. Eng 130 (5), 349–356.
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 (2),
111–128.

802