water
ArticlePhysics-Informed Data-Driven Models to Predict
Surface Runoff Water Quantity and Quality inAgricultural Fields
Jing Liang 1 , Wenzhe Li 2 , Scott A. Bradford 3 and Jiří Šimůnek 1, *
1
2
3

*

Department of Environmental Sciences, University of California, Riverside, CA 92521, USA;
jing.liang@email.ucr.edu
Department of Computer Sciences, University of Southern California, Los Angeles, CA 90089, USA;wenzheli@usc.edu
USDA, ARS, US Salinity Laboratory, Riverside, CA 92507, USA; Scott.Bradford@ars.usda.gov
Correspondence: jiri.simunek@ucr.edu; Tel.: +1-951-827-7854

Received: 28 December 2018; Accepted: 20 January 2019; Published: 24 January 2019


Abstract: Contaminants can be rapidly transported at the soil surface by runoff to surface water
bodies. Physically-based models (PBMs), which are based on the mathematical description of mainhydrological processes, are key tools for predicting surface water impairment. Along with PBMs,
data-driven models are becoming increasingly popular for describing the behavior of hydrologicaland water resources systems since these models can be used to complement or even replace physically
based-models. Here we propose a new data-driven model as an alternative to a physically-based
overland flow and transport model. First, we have developed a physically-based numerical model tosimulate overland flow and contaminant transport. A large number of numerical simulations wasthen carried out to develop a database containing information about the impact of various relevant
factors on surface runoff quantity and quality, such as different weather patterns, surface topography,
vegetation, soil conditions, contaminants, and best management practices. Finally, the resulting
database was used to train data-driven models. Several Machine Learning techniques were exploredto find input-output functional relations. The results indicate that the Neural Network model with
two hidden layers performed the best among selected data-driven models, accurately predicting
runoff water quantity and quality over a wide range of parameters.
Keywords: machine learning; surface runoff; contaminant transport; physically-based model;
agriculture field; synthetic data; HYDRUS-1D

1. Introduction
A variety of agricultural pollutants resulting from farming and ranching operations (e.g., sediment,
nutrients, pathogens, pesticides, metals, and salts) can lead to impairments of local and far-field water
quality. These diffuse nonpoint sources can directly harm ecosystem and watershed water quality,
and adversely affect our drinking water supply. Agricultural runoff, generated by rainfall or irrigation
events, can transfer pollutants accumulated at the land surface into receiving water bodies. This process
has been identified as one of the major causes of water impairment in agricultural settings [1–3].
The amount of pollutants transported to surface runoff water is dependent on local soil and crop
management, climatic conditions, contaminant properties and environmental factors. The discharge ofpollutants from agricultural activities to surface waters can be minimized by locally implementing
Best Management Practices (BMPs). However, obtaining monitoring data using field investigations isvery time consuming and expensive, and associated with many experimental difficulties. Furthermore,
simple field observations may be difficult to interpret to obtain a complete picture of potential
Water 2019, 11, 200; doi:10.3390/w11020200

www.mdpi.com/journal/water

Water 2019, 11, 200

2 of 21contaminant transport routes and mechanisms, and to extrapolate findings to other environmental
conditions and climates.
Alternatively, physically-based and data-driven models are potential tools for examining andoptimizing the effects of land-use changes and BMPs on surface water quality [4–7]. Indeed,numerical experiments can be cheaply conducted using a physically-based model (PBM) to
predict runoff water quality over a wide range of agricultural fields, weather patterns, and BMPs.
PBMs explicitly account for main hydrologic and contaminant transport processes using mathematicaldescriptions. Complex numerical or analytical techniques are commonly used to solve thesemathematical descriptions, which may require significant computational time as the spatial andtemporal scale of the considered problem increases.
In contrast to PBMs, data-driven models are based on functional relationships between input
and output variables. Artificial intelligence approaches, such as Machine Learning (ML) techniques,
are increasingly being used in data-driven models to quantify input-output functional relations
for complex systems in hydrology [7–19]. Several ML algorithms for modeling rainfall-runoff
processes include Linear Regression (LR), K-Nearest Neighbor Regression (kNN), Feed-forward
Neural Networks (NN), and Support Vector Machine (SVM) models [20–25]. The NN techniqueconsists of a neural network with a single or multiple hidden layer (Deep Neural Network) to describe
complex, nonlinear structures in data. NN with multiple hidden layers have been successfully appliedin many commercial products that are already available to the public, such as computer vision andspeech recognition, language translation, and self-driving cars [26–28]. While there are some studies
investigating the applicability of NN with multiple hidden layers in hydrology/contaminant transport
problems (e.g., [19,29–32]), only limited experience exists with using ML algorithms in environmental
applications involving surface runoff water quantity and quality. It is thus important to investigate thestrengths and limitations of data-driven models and compare their predictions with robust PBMs.
As discussed above, while PBMs are powerful tools for simulating spatial and temporal
distributions of various physical and chemical environmental variables, they require large amountsof data and computational resources. On the other hand, machine learning algorithms trained using
PBMs, or the so-called “metamodels” [33,34], make predictions almost instantly. The overall objectiveof this research thus was to investigate the limitations and strengths of ML-based metamodels incomparison with the HYDRUS overland flow module [35] in predicting contaminant loads in runoff
water from agricultural fields. We synthetically generated an extensive database using the HYDRUS
overland flow module that uses standard descriptions of overland flow and transport processes thathave been extensively verified to properly represent real field processes (e.g., [35–39]). This numericallygenerated database was used to generate metamodels that contained information about the impactof a wide range of physical factors on surface runoff quantity and quality in agricultural fields.
Inputs included information about the agricultural field (slope, length, and area), the soils (type,roughness, initial water content, and initial contaminant concentration), and precipitation events(rainfall intensity and duration). Outputs included the cumulative water volume and the contaminantmass in runoff water. This database was then used in conjunction with ML techniques to develop
correlation relationships between model inputs and outputs using various ML algorithms.
2. Database Preparation
2.1. Physically-Based Model
The one-dimensional diffusion-wave and advection-dispersion equations, which are commonlyused for describing overland flow and solute transport, respectively, have recently been implemented
into the popular HYDRUS-1D code [35]. This code has previously been used to numericallysolve the Richards equation, which simulates water flow, and the convection-dispersion equation,
which simulates contaminant transport in the subsurface [40]. Recent studies indicated that the
governing overland flow equation can be written in a similar mathematical formulation as the Richards

Water 2019, 11, 200

3 of 21equation [35,41–43] so that existing numerical schemes of HYDRUS-1D can be used to solve thefollowing 1D-diffusion wave equation:
∂
∂h
−
∂t
∂x

kh5/3 ∂(zb + h)√∂x
nman S

!

= R−I

(1)

where x is a space coordinate in direction of flow (L; where L denotes units of length), zb is a verticalcoordinate of the soil surface, t is time (T; where T denotes units of time), h is the surface water depth(L), R is rainfall or the irrigation rate (LT−1 ), and I is the infiltration rate (LT−1 ), S is the mean local
slope (defined as dzb /dx) (-), and k is a unit conversion factor (L1/3 T−1 ). The parameter nman is a
Manning’s roughness coefficient for overland flow, and it is dimensionless when k = 1 m1/3 s−1 or hasunits of TL−1/3 when k is not equal to unity. See the graphical abstract for a schematic of some ofthese variables.
In this work, the infiltration rate I is described using the Green-Ampt infiltration model [44].
This model is a simplified physically-based approach that is based on fundamental physics todescribe the infiltration process as a function of the soil suction head, water content, and soil
hydraulic conductivity:ψ∆θI = KS+1(2)
F
where ψ is the wetting front soil suction head (L), ∆θ is the water content difference (L3 L−3 ), KS is thesaturated hydraulic conductivity (LT−1 ), and F is the cumulative depth of infiltration (L). The value of∆θ is defined by using the following equations:∆θ = θe (1 − Se )(3)

θi − ( n − θ e )
θe

(4)

Se =

θ e = n − θr(5)where Se is the effective saturation (-), θ e is the effective water content (L3 L−3 ), n is the porosity (L3 L−3 ),θ i is the initial water content (L3 L−3 ), and θ r is the residual water content (L3 L−3 ).
Solute transport in overland flow is usually described using the advection-dispersion equation(ADE) of the form:∂hc ∂s
∂
∂c
∂Qc
+
=
hD
−
+ hφ + Rcr − Ic
(6)
∂t
∂t
∂x
∂x
∂x
where c is the solute concentration in the aqueous phase (ML−3 ; where M denotes units of mass),
cr is the concentration in rainfall water (ML−3 ); s is the sorbed solute concentration at the soil
surface area (ML−2 ), D is the effective dispersion coefficient accounting for both molecular diffusion
and hydrodynamic dispersion (L2 T−1 ), φ is a sink/source term that accounts for various zero- andfirst-order or other reactions (ML−3 T−1 ), and Q is the runoff flow rate (L2 T−1 ). The parameter Q isgiven as:kh5/3 ∂(h + z)√
Q = hU = −(7)∂xnman S
where U is a depth-averaged velocity (LT−1 ) calculated using the Manning–Strickler uniform
flow formula [41]. The effect of diffusion on the dispersion coefficient can often be ignored andin this case the value of D can be defined as the product of the dispersivity (λ, (L)) and U.
Kinetic sorption/desorption between the solid and aqueous phases can be described using thefollowing equation:∂s= ω (K D c − s)(8)∂tWater 2019, 11, 200

4 of 21

where KD is the linear equilibrium partition coefficient (L) and ω is the first-order desorption rate
coefficient (T−1 ). The product of KD and ω is proportional to the sorption rate coefficient. Details about
the numerical solutions of flow and transport models are provided in the HYDRUS-1D manual [45].
2.2. Numerical Simulations
A large number of numerical simulations was carried out to develop a database containing
information about the impact of various relevant factors on surface runoff quantity and quality.
In agricultural settings, land surface contaminants can be easily picked up and transported by surface
runoff or can seep into the soil by infiltration. Therefore, simulation scenarios considered in this work
mainly focused on the transport of various contaminants in overland flow over different soil and land
surfaces under different rainfall rates.
At first, numerous water flow simulations were run with various input parameters shown in
Table 1. In order to build a realistic training database, a wide range of input parameter values was
selected. In these simulations, precipitation was distributed uniformly over the land surface andlasted one hour. Historical records of precipitation frequency for the San Jacinto, CA watershed(http://hdsc.nws.noaa.gov/hdsc/pfds) were employed to generate precipitation-frequency estimatesthat were employed in simulations. Seven different percent probabilities of exceedance were selected:
1%, 2%, 4%, 10%, 20%, 50%, and 100%. Various agricultural field conditions were obtained byassigning different values of the Manning’s roughness coefficient, slope, and length. As a result,150 different field conditions (6 Manning’s roughness coefficients × 5 slopes × 5 lengths) wereconsidered. Additionally, four different initial water contents were considered: max(θ r , 0.1), 0.2, 0.3,
and the saturated water content (θ s ), for each soil type. Selected Green-Ampt infiltration parameters(soil porosity, residual water content, suction head, and hydraulic conductivity) were taken fromRawls et al. [46] (Table 2). In all water flow simulations, cumulative water volumes at the bottom outlet
at 25 evenly spaced times were outputs. The resulting water flow database thus contained 2,310,000
entries (22 soils × 7 rainfall rates × 150 field conditions × 4 initial water contents × 25 print times).
Table 1. Input parameters used in the numerical simulations.
Input

ValuesRainfall intensity (cm day−1 )Manning’s roughness coefficient
Slope (-)
Length (m)KD (cm)Initial water content, θ i21.6, 28.8, 40.8, 52.8, 69.6, 86.4, 105.6
0.01, 0.02, 0.04, 0.07, 0.15, 0.24
0.01, 0.02, 0.04, 0.08, 0.16
30, 60, 120, 240, 480
0.1, 0.5, 1, 1.5, 2max(θ r , 0.1), 0.2, 0.3, θ s

Table 2. Green-Ampt infiltration parameters for various standard and tilled* soil classes.
Texture

θs

θr

θe

ψ (cm)K (cm s−1 )Sand
Loamy Sand
Sandy Loam
Loam
Silt Loam
Sandy Clay Loam
Clay Loam
Silty Clay Loam
Sandy Clay
Silty Clay
Clay

0.437
0.437
0.453
0.463
0.501
0.398
0.464
0.471
0.430
0.470
0.475

0.020
0.036
0.041
0.029
0.015
0.068
0.155
0.039
0.109
0.047
0.090

0.417
0.401
0.412
0.434
0.486
0.330
0.309
0.432
0.321
0.423
0.385

4.950
6.130
11.0108.890
16.68021.85020.88027.30023.90029.22031.6303.272 × 10−3
8.306 × 10−4
3.028 × 10−4
9.444 × 10−5
1.806 × 10−4
4.167 × 10−5
2.778 × 10−5
2.778 × 10−5
1.667 × 10−5
1.389 × 10−58.333 × 10−6

* Tillage management: all five Green-Ampt infiltration parameters were increased by 20% over values of acorresponding textural class.

Water 2019, 11, 200

5 of 21

Solute transport simulations were subsequently conducted for the previously described waterflow simulations. In these simulations, it was assumed that a unit concentration (ci = 1 g/cm2 ) of
solute was initially distributed along the soil profile. Since the solute transport equation is a linear
equation, the results obtained for a unit initial solute concentrations can be simply multiplied byother initial concentrations to get corresponding results. Runoff water generated by rainfall events“mobilized” solutes from the land surface and transported them to the bottom outlet of the field.
The “mobilization” process was described using different sorption and desorption rates. The exchangeof solute concentrations between solid and liquid phases is controlled by the overland water flow rate
and kinetic sorption/desorption parameters (ω and KD ) in Equation (8). The kinetic sorption model
approaches equilibrium conditions when ω is large, whereas non-equilibrium conditions prevail when
ω is small and solute release is slow. In this study, the database considered transport processes andsolute desorption for near-equilibrium conditions (ω = 8640 day−1 ) to mitigate the complexity ofthese transport processes. Hence, two additional input parameters, variable KD shown in Table 1
and constant ω = 8640 day−1 , were accounted for in these simulations. Solute transport simulations
generated the output that included the cumulative solute mass at the bottom outlet at 25 evenly spaced
times. The inclusion of two additional parameters to describe transport of mobilized solutes (KD andω) resulted in a database that was five times bigger than the water flow database. Furthermore,
the variance in the solute transport database was higher than the water flow database.
3. Data-Driven Models
In this section, we briefly describe the data-driven models used in this paper. They are divided into
two broad classes, linear and nonlinear models. The selected linear data-driven models include Linear
Regression (LR) and Linear Support Vector Machine (SVM-L). Nonlinear models include K-Nearest
Neighbor Regression (kNN), Support Vector Machines with nonlinear kernel (SVM-NL), and Deep
Neural Networks (DNN) with different numbers of hidden layers. We used the previously generated
synthetic databases from the PBM to train the selected data-driven models to predict surface runoff
water quantity and quality. All these models were implemented using the open-source Scikit-learn
and Keras Libraries in Python 3.5 (Python Software Foundation, Wilmington, DE, USA) [47,48].
3.1. Linear Regression
LR is the simplest ML model, which is often treated as the baseline method. It applies linear
mapping to minimize the sum of least square errors between input data (either experimental or derived
from a PBM model) and output data (of the regression model). Suppose we have a set of training dataset (yi , xij ), where i =1, 2, 3, . . . , N and j =1, 2, 3, . . . , k. In here, i is the index for training samples in thedata set, j is the index for the input parameters. A general linear regression problem can be developedby assuming that output variables yi are influenced by input parameters xij as:kyi = β 0 + ∑ β j xij + ei

(9)

j =1

where βj are regression coefficients, β0 is intercept, and ei represent a deviation between actual andpredicted values. The ordinary least squares procedure seeks to minimize the total sum of the residuals.
However, since this approach treats data as a matrix, the process will be very computational intensive
when we have large data sets. Therefore, gradient-based optimization approaches, which optimize thevalues of the coefficients by iteratively minimizing the sum of the squared errors for each pair of input
and output values are used here [49].
3.2. Support Vector Machine
The SVM is another widely used model for both classification and regression tasks. The SVM
model was originally developed for classification problems to find a hyperplane, which can separate aWater 2019, 11, 200

6 of 21

data set into one class from those of another class [50,51]. The SVM models have been successfully
extended to solve regression problems by the introduction of the so-called ε-insensitive loss function.
This loss function is used to penalize errors that are greater than the threshold ε (an insensitivity zone).
A version of the SVM for regression problems is called support vector regression (SVR) and it depends
on a subset of the training data. The ε-SVR aims to find the flattest regression function f (xij ) that has at
most a deviation ε from actual yi for all training points xij :kf ( xij ) = ∑ w j xij + bi

(10)

j =1

where wj are weights and bi is the bias term that should be estimated from the training data. One great
benefit of using ε-SVR is the use of kernels, which inherently map the data into a nonlinear spacedepending on the chosen kernel function. Therefore, the above regression can be formulated as [52]:kf ( xij ) = ∑ w j Φ( xij ) + bi

(11)

j =1

where Φ(xij ) denotes nonlinear transformation (kernel functions). Several nonlinear kernel functions
are available, such as polynomial, radial basis, and sigmoid functions. In this study, we tested both
SVR with linear (SVM-L) and nonlinear (SVM-NL) sigmoid kernel functions.
The optimization of ε-SVR is more complicated than the LR model because there is a linear
constraint. In this case, the penalty factor C and the slack variables ξ i and ξ i∗ are introduced tocope with infeasible constraints [51]. Hence, wj and bi are estimated by minimizing the following
objective function:n1
min |w|2 + C ∑ (ξ i + ξ i ∗ )
(12)
2
i =1
subject to∗ yi − f ( xij ) − bi ≤ ε + ξ if ( xij ) + bi − yi ≤ ε + ξ iξi, ξi ∗ ≥ 0(13)where |w|2 = ∑kj=1 w2j . The penalty factor C determines the trade-off between the flatness of thefunction f (xij ) and the amount up to which deviations larger than ε are tolerated. The slack variables ξ iand ξ i ∗ determine the degree to which data points will be penalized if the error is larger than ε.
3.3. K-Nearest Neighbor Regression
K-Nearest Neighbor (kNN) Regression is another machine learning method for regression
problems. The kNN is a non-parametric regression type of analysis. Predictions are made by computing
the distance between all the data points in the database and searching for the number of neighbors (Knn )which are most similar to a new instance. Once the nearest-neighbor list is obtained, the predictionof output y is achieved by assigning weights to the contributions of the neighbors, so that the nearer
neighbors contribute more to the average than the more distant ones. For example, the Euclidean
distance d between lth training sample xlj and ith test sample xij is defined as [53]:vu kud( xij, xlj ) = t ∑ ( xij − xlj )2
j =1

(14)

Water 2019, 11, 200

7 of 21

The weight wd can be defined by a selected Knn value and a distance between lth and ith samples:
wd ( xij , xlj ) =exp(−d( xij , xlj ))

(15)

Knn

∑ exp(−d( xij , xlj ))l =1

7 of 23

yi = Kwnnd ( xij , xlj ) yl

(16)

Water 2019, 11, x FOR PEER REVIEW

Thus, the output yi can be predicted byKnn=1 ∑ wd ( xij , x lj ) y l
yi l=

(16)

l =1

3.4. Deep Feed-Forward Network
3.4. Deep Feed-Forward Network
Neural networks (NN) are the most popular ML models used nowadays. The classical feed(NN)are thepopularmodelsused nowadays.
Thequality
classical
forwardNeuralnetworknetworks
is a neuralnetworkthat mosthas been
widelyMLused
in the hydrologic
and water
feed-forward
networkisaneuralnetwork
that
has
been
widelyusedinthehydrologic
and
water
problems [22,23,54–56]. The fundamental components of these algorithms include the input layer,
quality problems[22,23,54–56].
Theneurons.
fundamental
componentsof thesealgorithms
includethethat
input
intermediate
layers, output
layer, andAdjustableparameterscharacterizetheneurons
layer,
intermediate
layers,outputlayer,
andneurons.
Adjustableparameterscharacterizetheneurons
are connected to the input and transmit the information from one neuron to another in one direction:
thatconnected
to the input
transmit
information
fromlayer.
one neuron
another in
from
anare
input
layer, throughone orandmorehiddenthelayers,to an output
Figure 1toillustrates
theonedirection:from an input
through
one orhiddenlayers,an output
layer.
Figure 1
weighted
adjustableparameters
for anlayer,
individual
neuron.
In morea network,theinput tovariables
xij areillustrates
theadjustableparametersforanindividual
neuron.
Inanetwork,the
input
variables
xij areand summed up to produce the hidden neurons Sp (p = 1, 2, …, P):
weighted and summed up to produce the hidden neurons Sp (p = 1, 2, . . . , P):k

Sp = w
k jp xij

j =1 ∑ w jp xij
Sp =

(17)
(17)

j =1

where p is the number of hidden neurons. A neuron computes an output, based on the weighted sumof allits inputsp), according
to anneurons.
activation(f(Sp)). Finally,an activation
which
where
p is the(Snumberof hidden
A function
neuron computes
an output,
based onfunction,
the weighted
sum
is applied
toSp,providesthe
final
outputfrom
thislogistic
sigmoid[10](see
Figure1).
Thesigmoidof all its inputs (Sp ), according to an activation function (f (Sp )). Finally, an activation function, which isactivation
is giventhe
as: final output from this logistic sigmoid [10] (see Figure 1). The sigmoid
appliedfunction
to Sp , provides
activation function is given as:
1
f (S p ) =
1
(18)
−S
(18)
f (S p ) =1 + e p −S p
1+eFigure 1. Activation of a single neuron.

Figure 1. Activation of a single neuron.

The so-called Deep Neural Networks (DNN) model is a direct extension of the feed-forward
networkthat contains
multiple
hidden layers
thatare stacked
as shown
in Figure[57]. The DNN
The so-called
Deep Neural
Networks
(DNN)
model
is a direct
extensionof the2feed-forwardmodel
becomes“deeper”
asthenumberofhiddenlayersincreases.
DeeperNNmodelsnetwork that contains multiple hidden layers that are stacked as shown in Figure 2 [57]. Theallow
DNNforsimulating
more
complicated,mapping
betweenDeeper
input anddata.
However,model
becomes“deeper”
as the nonlinear
number of
hidden functions
layers increases.
NNoutputmodels
allow
for
DNN
modelsoften
suffer
from
overfitting
when
using
too
many
hiddenlayers.
This
problemcansimulating more complicated, nonlinear mapping functions between input and output data. be
mitigatedby dropping
layers.suffer
We havetestedthe DNNmodel
withtoo
one,many
two, andthreelayers.
hiddenThis
layersHowever,DNN
models often
from
overfitting
when
using
hiddento determine
optimum
oflayers.
hiddenWelayersandto minimize
overfitting.
problemcan be themitigatedby number
dropping
have
testedthe DNNmodel with one, two, andthree hidden layers to determine the optimum number of hidden layers and to minimize overfitting.

Water 2019, 11, 200
Water 2019, 11, x FOR PEER REVIEW

8 of 21
8 of 23

Figure 2. A feed-forward neural network (a single hidden layer, left) and a deep feed-forward neural
Figure 2. A feed-forward neural network (a single hidden layer, left) and a deep feed-forward neuralnetwork (multiple hidden layers, right).
network (multiple hidden layers, right).

4. Model Training and Evaluation
4. Model Training and Evaluation
The runoff water database generated using the PBM was divided into training (80% of the entire
water
was
divided intotraining(80% ofinput
the entire
database)The
andrunoff
test (20%of database
the entiregenerated
database)using
data thesetsPBMfor thedata-driven
models.
Variable
database)
and
test
(20%
of
the
entiredatabase)
data
sets
for
the
data-driven
models.
Variableparameters and output values were used to train the data-driven models by randomly selecting 18input
parametersandEachoutputvalues wereused
to tunedtrain thedata-driven
modelsandby randomly
selecting 18
out of
22 soil types.
data-driven
model
was
using
cross validation
further evaluated
out
of
22
soil
types.
Eachdata-driven
model
was
tuned
using
cross
validationandfurtherevaluatedusing the left-out test data set. After the training was completed, data-driven models were testedonusing the left-out
dataoverset. many
After thetraining
was completed,data-driven
models were
tested on
the remaining
four soiltest
types
different
scenarios.
The accuracy
and predictive
capability
four
soil types
over
many
different scenarios.
The accuracy
capability
werethefirstremaining
assessed bymakingscatterplots
of cumulative
runoff water
volumesand(Qc )predictive
for data-driven
2
were
first
assessedbymakingscatterplots
ofcumulative
runoffwater
volumes(Qc)
fordata-driven
models and the PBM. The calculated R and the best-fitted regression line were shown in these scatterand the
The calculated
R2 and thebest-fitted
regression line were shown in these scatterplotsmodels
to quantify
the PBM.
goodnessof the data-driven
model
prediction.
plots
toquantify
thegoodnessofthedata-driven
model
prediction.
In addition, two examples were chosen to show the accuracy between observed and predicted
Inrunoff
addition,
two examples
wereto show
the accuracy
observed
predicted
cumulativehydrographs
betweenthechosen
PBM anddata-driven
models,between
respectively,
over aandrange
ofcumulative
runoffhydrographs
betweenthe PBM considered
and data-driven
models, initialrespectively,
over a range
initialwater contents
andsoil types. The
first example
four different
water contents
of initial
water
contents theandrainfall
soil types.
exampleconsidereddifferentwater
and thefollowingconditions:rate ofThe69.6first
cm/dayfor 1 h,the surfacefour
lengthof 240 initialm, a Silt
the following
conditions:
theofrainfall
rate grass),of 69.6andcm/dayforslope
1 h, ofthe16%.
surface
Loamcontents
soil, theand
Manning’s
roughnesscoefficient
0.24 (densea steep
The length
total of
240 m, timea Siltwas
Loamof 0.24 (dense
grass),
a steepinslope
simulation
1 hsoil,andthe10 Manning’s
min so thatroughness
the water coefficientrecession process
could befullyand
captured
16%. The total
time consideredwas 1 h andsurface
10 minso thattheSilt
water
recession
process
could betheseofsimulations.
The simulation
second example
runofffrom
Loam,Clay, tilled
Silt Clay,these thesimulations.
second example
consideredsurfacerunofffromrate
Silt ofLoam,
and fully
tilled captured
Clay soilsinunderfollowingThe
conditions:an initiallysaturatedsoil, therainfallClay,tilled
Silt
Clay,
and
tilled
Clay
soils
under
the
followingconditions:aninitiallysaturated
105.4 cm/day for 1 h, a surface length of 48 m, a slope of 16%, and the Manning’s roughness coefficient soil,
the rainfall rate of 105.4 cm/day for 1 h, a surface length of 48 m, a slope of 16%, and the Manning’sof 0.24.
roughnesscoefficientwereof 0.24.
Statistical measuresused to quantify the ability of various ML models to predict an output
Statisticalmeasureswereto quantify
the abilityof meanvarious
ML modelsto predictan output
from the PBMs. Calculated statisticalused
parameters
included
the rootsquareerror (RMSE),
the mean
fromerror
the PBMs.
statisticalparameters
included
the rootmeanerror Error(RMSE),
absolute(MAE),Calculated
the mean bias
error (MBE),
the model
efficiency(EF),andsquare
the Relative∆i . themean
absolute
error (MAE),
the mean bias
These
statisticalparameters
were determined
as: error (MBE), the model efficiency (EF), and the Relative

Error Δ i . These statistical parameters weres determined as:
2

RMSE =∑iN=1 (Oi N− Pi ) 2
 (Oi − Pi )RMSE = N i =1

(19)

Pi |
∑iN=1 |Oi −N|
Oi − Pi |
MAE = N i =1

(20)

N

MAE =
MBE =

N
∑iN=1 (Oi − Pi )N N

 (O −2 P )
MBE
N =
i

i =1

EF = 1 −Oi )∑iN=1 (Oi −N
EF = 1 −

(20)

(21)

i

∑i=1 (Oi − P
Ni )

(19)

(21)
(22)

2

 (O − P ) (O − O )

2

i =1
N

i

i

i =1

i

i

2

(22)

Water 2019, 11, 200

9 of 21

Relative Error ∆i =( Pi − Oi )( Pi +Oi )2(23)

where Pi and Oi are the output values predicted by data-driven models and observed data obtained
from the PBM, respectively, Oi is the average of observed data, and N is the number of observations.
It should be noted that data with zero values of surface runoff were removed from this analysis formodel inter-comparison.
A database for transport of mobilized solute with runoff water (e.g., the cumulative solute mass atthe outlet) was developed using the PBM that encompassed a wide range in rainfall rates, soil physicalproperties, initial water contents, slopes, field size, Manning’s roughness coefficients, and kinetic
sorption/desorption parameters. As we mentioned before, the solute transport database in this study
considered transport processes and solute desorption for near-equilibrium conditions with a constantdesorption rate coefficient (ω = 8640 day−1 ) and different KD values. The solute transport database
had a similar structure as the water flow database but was five times bigger in size.
The contaminant transport is highly related to water flow. The data-driven model with thebest performance in predicting the runoff water volume was therefore selected for training withthe solute transport databases. Data-driven models were trained against 80% of the entire solutetransport database using Qc , KD , and other input parameters associated with the water flow
database. The remaining 20% of the solute transport database were employed for testing the model.
Two simulation examples were chosen to further investigate the predictive ability of the trained
data-driven models. These simulation scenarios considered the same input variables as for thewater flow model (e.g., differences in initial water contents or soil texture) with KD = 2 cm andω = 8640 day−1 .
5. Results and Discussion
5.1. Surface Runoff Quantity
5.1.1. Linear Regression
A scatter plot of estimated Qc using the LR model versus observed Qc for the PBM test data set
are shown in Figure 3a. We can see that R2 is equal to 0.439, which indicates that the LR model failed
to accurately describe the PBM data. The LR algorithm usually works well with high-bias and low
variance data. However, the PBM describes real-world scenarios with complex, nonlinear relationships
between input and output variables that generate high variance data sets.
The observed (PBM) runoff water volume in Figure 3b decreased and was increasingly delayedwhen the initial soil water content was lower. The observed (PBM) simulations in Figure 3c demonstratethat the Clay and tilled Clay textured soils produced the most runoff water because of their lower
saturated hydraulic conductivity and infiltration. After rainfall stops, the cumulative runoff water
flux continues to increase because water flows slowly over the rough surface before reaching thebottom boundary. The LR model failed to accurately describe runoff water dynamic for these
two selected examples, especially when the initial water content was higher (Figure 3b) and thehydraulic conductivity was lower (Figure 3c). In particular, the LR model tended to overestimate
Qc for early times and underestimate Qc for later times. Other authors (e.g., [18]) similarly foundthat the linear regression model does not provide enough flexibility to represent complex nonlinearhydrological phenomena.

Water 2019,
Water
2019, 11,
11, 200
x FOR PEER REVIEW

10of
of 23
21
10

(a)

(b)

(c)
Figure3. Comparison
cumulative water(Qcc incalculated by
the physically-based
physically-based
Figure 3.
Comparison ofof cumulative
water fluxes
fluxes (Qin liters,
liters, L)L) calculatedby themodels((a)
overall,(b)(b)
different
initialwater
contents,(c)models PBM
PBMandandthetheLinear
LinearRegression
Regression(LR)(LR)model
model
((a)
overall,different
initialwater
contents,
different
textures)).
(c) different textures)).

Water 2019, 11, 200
Water 2019, 11, x FOR PEER REVIEW

11 of 21
11 of 215.1.2.
5.1.2. Support
Support Vector
VectorMachine
MachineSimilarlinearkernelfunctions(SVM-L)
Similar to
tothe
theLR
LRmodel,model,thetheSupport
SupportVector
VectorMachine
Machinemodel
modelwith
withlinearkernelfunctions(SVMproduced
a linear
input-output
relationship.
TwoTwouser-definedparameters
for themodelmodel
wereL) produced
a linear
input-output
relationship.
user-definedparameters
forSVM-L
the SVM-L
thepenaltyparameter
C on the
term,
which
set to1, andthe1,maximum
number ofnumber
iterations,werethe penalty
parameter
C erroron theerror
term,waswhich
was
set toand the maximumofwhich
was which2000 [58].
4a indicatesthis SVM-L
model
had
an evenR2 lower(R2 = 0.386)
iterations,
was Figure
2000 [58].
Figure 4athat
indicatesthat this
SVM-L
model
hadloweran evenR2 (R2
2 = 0.439).
than
the thanLR modelsimulation
examplesshown shownin Figure4b (different
initial= 0.386)
the LR(Rmodel(R2 = Similarly,
0.439). Similarly,
simulation
examplesin Figure4b (different
water
and 4c (different
textures)indicateindicatethat the that
SVM-L
failed to failed
captureinitialcontents)
water contents)
and 4c (different
textures)the model
SVM-Lcompletely
model completely
tothelargethevaluesQc at later
Linearmodelsand SVM-L)
therefore
not well suitedtocapturelargeofvaluesof Qtimes.
c at latertimes.
Linear(LRmodels(LR andare
SVM-L)
are therefore
not welldescriberunoff problems.
suited tocomplex
describenonlinear
complex nonlinear
runoff problems.

(a)

(b)
Figure 4. Cont.

Water 2019,2019, 11,
11, 200
x FOR PEER REVIEW
Water
Water 2019, 11, x FOR PEER REVIEW

12
21
12 of
of 21
12 of 21(c)(c)

Figure 4. Comparison of cumulative water fluxes (Qc in liters, L) calculated by the PBM and the Linear
Figure 4.
Comparison
ofofcumulative
water
fluxes(Qc(Qin(b)
liters,L) calculatedby thePBM
LinearSupportMachine(SVM-L)model((a)
overall,different
initialwater
contents,(c)thedifferent
Figure4.Vector
Comparison
cumulative
water
fluxesliters,L)
calculatedbythe and
PBM
and
the
c in
Support
VectorMachine(SVM-L)model((a)
overall,(b)different
initialwater
contents,(c)different
textures)).
LinearSupport Vector Machine (SVM-L) model ((a) overall, (b) different initial water contents,
textures)).
(c)different textures)).

Many researchers have indicated that the SVM model with nonlinear kernels (SVM-NL) can
Many researchers
have
the
model
nonlinear
kernels(SVM-NL)
canresearchers
have indicatedindicated that
the SVM
SVMproblems
model with
withnonlinear
kernels previous(SVM-NL)
cansuccessfully
describe surface/subsurfacehydrology[18,59,60].
However,
SVMsuccessfully
describesurface/subsurfacehydrology
problems[18,59,60].
However,previous
SVMsuccessfully
describesurface/subsurfacehydrology
problems[18,59,60].
However,previous
SVM-NLNL applications have considered only a relatively small training and test data set. Training the SVMNL model
applications
considered
only
a relatively
small
training
data
set.
Trainingthe
SVMapplications
haveconsidered
only
a relatively
small
training
andandtest testdata
set.
Training
the
SVM-NL
NL
takes
ahave
tremendous
amountofcomputationaltimebecauseofthecomplexityof
thismodel.
NL
model
takes
a
tremendous
amountofcomputationaltimebecauseofthecomplexityofmodel.
model
takes
a
tremendous
amountofcomputationaltimebecauseofthecomplexityofthisFor example, we estimated that training the SVM-NL model to our complete runoff water database
For example,
example,estimatedtraining
SVM-NLmodel totoour complete
completeonly
runoffwater
databaseFor
we estimatedthat training
SVM-NLmodel
ourrunoffwater
databasewouldtake somewhere
betweena few thedays
and a month.
Consequently,
about
20%
of thewould taketakesomewhere
betweenfew
days andand aa month.
month.
thewould
somewhere
betweenaa few
days
Consequently,
about model.
20% of The
theoriginal
runoffwater training
data
set were
randomly
selected
to train theonly
SVM-NLoriginal
runoffwater
training
data
set
were
randomly
selected
to
train
the
SVM-NLmodel.
The
original
runoffwater
training
data
set
were
randomly
selected
to
train
the
SVM-NLmodel.
The
trainedtrained model was then used for testing against the test data set with C = 1 and ε = 0.2. Figure 5a
trainedwas
model
then
used
foragainst
testingagainstthe settestwith
dataCsetC
1 and= 0.2.
Figure5amodel
thenwasusedfortesting
thetestSVM-NL
data= 1withandobserved
ε == 0.2.
Figure
5athe
shows
shows
the
scatterplotof estimatedQc using
themodelversusQcεfrom
test the
set.
shows
plot
of estimatedQc using
themuch
SVM-NLmodel
versusQc fromthe testset.
scatterplot
estimated
QSVM-NL
SVM-NLmodel
versusobservedQthetestmodels.
set.
TheTheR2
c using thec fromThe
R2 the(R2 scatter=of
0.952)
for themodel
washigherthan
for theobserved
LR
and SVM-L
2 =R0.952)
2 (R2 =forThe
0.952)
forprovided
the SVM-NL
model
was
much
higherthan
for
theSVM-L
LR andmodels.
SVM-Lmodels.
The(R
the
SVM-NL
model
was
much
higherthan
for theLR
andThe
SVM-NL
SVM-NLmodel
also
a much-improved
descriptionof thePBM
resultsforthe twoexamplesSVM-NLalso
a much-improved
the PBMforthe two
model
also
provided
aprovidedmuch-improved
descriptionof thePBMofresults
forresults
the twoexamples
shown
in
shown
in modelFigures
5b and5c. However,
there
weredescription
stillconsiderable
deviations
betweenthe examples
SVM-NL
shown
in
Figures5b
and5c.
However,
there
were
still
considerable
deviations
betweenthe
SVM-NL
Figure5b,c.
However,
there
were
still
considerable
deviations
betweenthe
SVM-NL
model
and
the
model and the PBM, especially for larger print times, higher initial water contents (Figure 5b), and
model
and
the PBM,especiallyfor larger
print
times,
higherinitial(Figure
water 5b),
contents(Figure
5b),
and
PBM,
for
largerprint
higherinitialwater
contents
and forthe Silt
Loamfor
theespeciallySilt Loamsoil(Figure5c).times,
for theSilt Loamsoil(Figure5c). soil (Figure 5c).

(a)(a)
Figure 5. Cont.

Water 2019, 11, 200
Water 2019, 11, x FOR PEER REVIEW

13 of 21
13 of 21(b)(c)
Figure 5.
5. Comparison
and thetheFigureComparison ofof cumulative
cumulative water
water fluxes
fluxes (Q(Qcc inin liters,
liters, L)L) calculated
calculated byby thethe PBMPBM andNonlinear Support
Machine (SVM-NL)(SVM-NL) model
model ((a)
Nonlinear
Support Vector
Vector Machine((a) overall,overall, (b)(b) different
different initial
initial water
water contents,
contents,(c) different
different textures)).
textures)).

5.1.3.
K-Nearest Neighbor
5.1.3. K-Nearest
Neighbor Regression
Regression
The
The kNNkNN models
models with
with different
different numbers
numbers ofof neighbors
neighbors (the(the KK value)value) were
were tested
tested while
while allall points
pointsineachneighborhood
wereweighted
equally.
Weobservedthattheoptimal
K
value
is
equal
to 5.
5.
in each neighborhood were weighted equally. We observed that the optimal K value is equal to22The
results ofof calculations
calculations with
with thethekNNkNNmodelmodelarearesummarizedsummarizedininFigure
Figure6.6.The
The= 0.957)
The results
R2 R
(R2(R
= 0.957)
for
2 = 0.952). However,
2
for
the
kNN-5
model
was
slightlyhigher
than
that
for
the
SVM-NLmodel(R
the kNN-5 model was slightly higher than that for the SVM-NL model (R = 0.952). However, there
there
were morepredicted
outliers
the observed
Althoughthe overallc was relatively
were morepredicted
outliers
when when
the observed
Qc wasQrelatively
small. small.
Althoughthe overall
runoffrunoffwere correctly
capturedin Figurethe kNN-5
underestimatedthe cumulative
trends trends
were correctly
capturedin Figure6b, the6b,
kNN-5
modelmodel
underestimatedthe cumulative
water
water
flux
at
every
print
time
and
incorrectlypredicted
thatQ
stoppedincreasing
when
rainfallcflux at every print time and incorrectly predicted that Qc stopped increasing when rainfall ended.
ended.
Furthermore,this deviation
to forincrease
forinitialsmaller
initialwater However,contents. However,Furthermore,this deviation
tended totended
increase
smallerwater
contents.
Figure 6c
Figure6c shows
that thekNN-5
didbettera muchjob of predicting
large
QcSVM-NL
than theshows thatthe kNN-5
model
didmodela muchjob betterof predicting
large valuesof values
Qc thanofthe
SVM-NLmodel.
model.

Water 2019,
Water
2019, 11,
11, 200
x FOR PEER REVIEW

14of
of 23
21
16

(a)

(b)

(c)
Figurecumulative
water
fluxes(Qc (Qin cliters,L) calculatedby theby
PBM
theandkFigure 6.
6. Comparison
Comparisonof of
cumulative
water
fluxesin liters,L) calculated
the andPBM
NearestNeighbor
Regression(kNN-5)model((a)
overall,(b)different
initialwater
contents,(c)
the k-Nearest Neighbor Regression (kNN-5) model ((a) overall, (b) different initial water contents,
different
textures)).
(c) different
textures)).

5.1.4. Deep Feed-Forward Neural Network
The DNN models with one, two, and three hidden layers were tested to determine the optimum
number of hidden layers and to minimize overfitting. In this study, the optimal number of neurons
in each hidden layer was determined by trial and error. We tested possible numbers of neuron from

Water 2019, 11, 200

15 of 215.1.4. Deep Feed-Forward Neural Network
The DNN models with one, two, and three hidden layers were tested to determine the optimum
number of hidden layers and to minimize overfitting. In this study, the optimal number of neurons
Water 2019, 11, x FOR PEER REVIEW
17 of 23
in each hidden layer was determined by trial and error. We tested possible numbers of neuron from
3 to 20 to see if the model performance was improved. The R22 and the Root Mean Squared Error
3 to 20 to see if the model performance was improved. The R and the Root Mean Squared Error(RMSE) were not improved when neuron numbers were increased beyond 16. Therefore, the number(RMSE) were not improved when neuron numbers were increased beyond 16. Therefore, the number
ofofneurons
neuronsinin eacheach hidden
hidden layer
layer waswas setset toto 16.
16.
Figure7ashows
thescatterplotofestimated
Qcc usingusingDNN
DNNmodelsmodelswith
withthree
threedifferent
differentlayers
layersFigure 7a shows the scatter plot of estimated Q
2 were equal to 0.978, 0.981, and 0.959 for theversusobserved
Q
for
the
PBM
data
set.
ValuesofRversus observed Qcc for the PBM data set. Values of R2 were equal to 0.978, 0.981, and 0.959 for the
DNN-1,
models, respectively.
respectively. The
The DNN
DNNmodelsmodelsprovided
providedaadescription
descriptionofofrunoff
runoff
DNN-1,DNN-2,
DNN-2, andand DNN-3
DNN-3 models,
2 values than the other ML models. However, the top
thatwas
fairlyaccurateandobtained
higherR
that was fairly accurate and obtained higher R2 values than the other ML models. However, the top
and
DNN-1 andand DNN-3
DNN-3 models
modelstended
tendedtotounderestimate
underestimate
andbottom
bottom rows
rows ofof Figure
Figure 77 indicate
indicate thatthat thethe DNN-1
andoverestimate
observed
valuesofQ,
respectively,
whereas
the
DNN-2
model
provided
thebestbestcand overestimate observed values of Qc, respectively, whereas the DNN-2 model provided theprediction.
This
resultsuggests
that
the
one
hidden
layer
model
was
not
“deep”
enough,whereasprediction. This result suggests that the one hidden layer model was not “deep” enough, whereas a athree-layer
and results
results inin overfitting.
overfitting. AAmore
morequantitative
quantitative
three-layer model
model probably
probably was
was tootoo complicatedcomplicated andcomparisonofvariousMLmodelsisgivenbelow.
comparison of various ML models is given below.

(a)

(b)

DNN-1

(c)

DNN-1

(c)

(b)

(a)

(c)

(b)

DNN-3

DNN-2

DNN-2

DNN-2

(a)

DNN-1

DNN-3

DNN-3

Figurefluxes (Q(Qccininliters,
liters,L)L)calculated
calculatedby
bythe
thePBM
PBMandandthetheDeep
Deep
Figure7.
7. Comparison
Comparison of cumulative water fluxes
Feed-Forward(DNN-1 (top(top row),row),DNN-2
DNN-2(middle(middlerow),row),and
andDNN-3
DNN-3(bottom(bottomrow))row))
Feed-Forward Neural
Neural Network (DNN-1
modelsdifferent textures)).
textures)).
models((a)((a) overall,overall, (b)(b) different
different initial water contents, (c) different

5.1.5.
Models
5.1.5. Comparison
Comparison ofof Data-Driven
Data-Driven Models
Figureboxplot
plotofofthe
therelative
relative
error
distribution
variousdata-driven
models,Figure 88 presents
presents aa box
error
distribution
forforvariousdata-driven
models,
as
as
well
as
other
statisticalparameters
for
their
agreement
with
the
PBM.
The
LR
and
SVM-L
modelswell as other statistical parameters for their agreement with the PBM. The LR and SVM-L models
show
indicates that
thatthese
thesemodels
modelswere
werenot
notwellwellsuited
suitedtotopredict
predict
runoff.
showbig
big error
error distributions,distributions, which
which indicatesrunoff.
The
SVM-NLmodel
performed
muchbetter
than
the
SVM-L
and
LR
models,
but
not
as
well
as
kNN-5
The SVM-NL model performed much better than the SVM-L and LR models, but not as well as kNNand
DNN
models.
It is Itpossiblethat theof the SVM-NL
model could
further
5 andDNN
models.
is possible
thatperformancethe performanceof the SVM-NL
modelbecould
beimproved
further
improved if the training data size was increased, but this would also dramatically increase thecomputational time. The boxplot of relative errors for the kNN-5 model was narrow, but it included
a lot of outliers which indicated that the model performance may be different under other conditions.
Figures 7a and 8 indicate that outliers of the relative error for DNN-1 appeared when Qc was over
100 L, whereas the DNN-3 model slightly underestimated large Qc values. The DNN-2 model showsWater 2019, 11, 200

16 of 21if the training data size was increased, but this would also dramatically increase the computational
time. The boxplot of relative errors for the kNN-5 model was narrow, but it included a lot of outliers
indicatedthat
the model performance may be different under other conditions. Figures
7a 23
and 8
Waterwhich
2019, 11,
x FOR PEER
REVIEW
18 ofindicate that outliers of the relative error for DNN-1 appeared when Qc was over 100 L, whereas thethe narrowest
box slightly
plot of relative
errors andtheQbest
performancestatistics
forthe considered
ML boxDNN-3 model
underestimatedlarge
The DNN-2
model
shows
the narrowest
c values.
models.
This
indicates
that
thetheDNN-2
model provided
the best
predictionof runoffwater This
volumes.
plot ofrelative
errorsandbest performancestatistics
for theconsidered
ML models.
indicatesthat the DNN-2 model provided the best prediction of runoff water volumes.

Figure 8. Statistical parameters associated with the trained data-driven models of runoff quantity
and8.quality.
The
root meanassociatedsquare error(RMSE),
meandata-driven
bias error (MBE),
absolute
errorand(MAE),
FigureStatisticalparameters
withthe trainedmodelsmeanof runoff
quantity
model
(EF),
and a error
box plotof the mean
relativebias
error
distribution
for data-driven
modelsare given.
quality.
Theefficiency
root meansquare(RMSE),
error(MBE), mean
absolute error(MAE),
Theefficiency
units of RMSE,
MBE,
andplotMAE
therunoff error
quantityand quality
liters andmodels
g, respectively.
model(EF), and
a box
of of
the
relative
distribution
for aredata-driven
are
TheTheupper
andboundaries
of MAEthe boxesshow
the 75thand 25th
whiskersof thegiven.
units
oflower
RMSE,
MBE, andof therunoffquantityand percentile,
quality aretheliters
and g,
box
plot
show
the
maximumandminimum
values,
and
the
red
line
within
the
box
is
the
medianvalue.
respectively. The upper and lower boundaries of the boxes show the 75th and 25th percentile, theBlue dot
symbols
indicatedthemaximum
outliers. and minimum values, and the red line within the box iswhiskersof thebox plotshow thethe median value. Blue dot symbols indicated the outliers.

5.2. Surface Runoff Quality

5.2. Surface
Quality
WeRunoff
discoveredthat the DNN-2 model structure with fine adjusted model parameters performs
best
among
all
selectedmodelsfor predicting
the runoff
water
volume. performs
We therefore
We discovered that the data-driven
DNN-2 model
structurewith fine adjusted
model
parameters
trained
the
DNN-2
model
with
the
near
equilibrium
solutetransport
database.
The
DNN-2
model forbest among all selected data-driven models for predicting the runoff water volume. We therefore
water
could
be directly
adjusted
to produce solutea data-driven
modelthat relates
physicalfactors
trainedtheflow
DNN-2
modelwith thenear equilibriumtransport
database.
The DNN-2
model
for andthecumulative
solutemass.
water flow could be directly adjusted to produce a data-driven model that relates physical factors
Figure 9a presentsa scatter plot of the cumulative solute mass predicted by the DNN-2 modeland the cumulative
solute mass.
andobserved
by
the
PBM
under
constructed
Figure 9a presents a scatterplot nearof theequilibrium
cumulativeconditions.
solute massThe
predicted
by theDNN-2
DNN-2relationshipmodel
betweeninput
and
outputvariablesfornearequilibrium
conditions
wasvery
accurate.
The
goodnessand observed by the PBM under near equilibrium conditions. The constructed DNN-2 relationshipof theinput
DNN-2
for solutetransport
near equilibrium
is reflectedbygoodness
the statisticalbetweenandpredictions
output variables
for nearequilibrium
conditionsconditions
was very accurate.
The
parameters
(RMSE
=
4.77
g,
MAE
=1.57
g,
MBE
=
0.18
g,
and
EF
=
0.96
g)
in
Figure8, the valueof the DNN-2 predictions for solute transport near equilibrium conditions is reflected by the statistical
2 = 0.992 in Figure 9a, and the agreement between observed (PBM) and predicted results for theRparameters (RMSE = 4.77 g, MAE =1.57 g, MBE = 0.18 g, and EF = 0.96 g) in Figure 8, the value R2 =examplesimulations
various watercontents(Figure9b) andand predictedsoil types results(Figurefor
9c).the example
0.992
in Figure9a, and thefor
agreement
betweenobserved(PBM)
simulations for various water contents (Figure 9b) and soil types (Figure 9c).

Water 2019, 11, x FOR PEER REVIEW
Water 2019, 11, 200

19 of 23
17 of 21

(a)

(b)

(c)
Figure9.9. Comparison
Comparison ofof thethe cumulative
cumulative solute
solute mass
mass (in(in grams,
grams, M)M) calculated
calculated byby thethe PBMPBM andand thetheFigure
DNN-2model
modelforfornearequilibrium
trainingdata((a) overall,(b) different
initialwater
DNN-2
thethenear
equilibrium
training
data set((a)setoverall,(b) different
initial watercontents,
contents,(c)different
textures)).
(c) different textures)).

Water 2019, 11, 200

18 of 21

6. Conclusions and Outlook
This study tested the ability of data-driven models to mimic PBMs to predict surface runoff
water quantity and quality in agricultural settings. A physically-based overland flow and transportmodel was used to develop a large database containing information about the impact of various
factors on surface runoff quantity and quality. In order to build a realistic training database,a wide range of input parameter values was selected. The input factors included: (a) rainfall
intensities (for seven different percent probabilities of exceedance of 1%, 2%, 4%, 10%, 20%, 50%,
and 100% for the San Jacinto Watershed, Southern California), (b) Manning’s roughness coefficients(reflecting management practices; 6 values), (c) field slopes (5 values), (d) field lengths (5 values),(e) soil properties (affecting infiltration; 11 soil types with and without tillage), (f) initial water contents(4 values), and (g) solute sorption parameters (5 conditions). The resulting water flow database thus
contained 2,310,000 entries, whereas the solute transport database was 5 times larger.
Multiple data-driven models were tested on their ability to relate/correlate model inputs with
outputs. The following data-driven models were used in the analysis: Linear Regression (LR),k-Nearest Neighbor Regression (kNN), Support Vector Machine with linear (SVM-L) and nonlinear(SVM-NL) kernels, and Deep Neural Networks (DNN) (Neural Networks with multiple hidden
layers). The LR and SWM-L models failed to accurately describe runoff water dynamics, having theregression coefficient (R2 ) less than 0.5. In particular, both linear models tended to overestimate therunoff water volume for early times and underestimate it for later times. Both SVM-NL and kNNmodels performed much better than linear models, having R2 of 0.93 and 0.96, respectively. However,there were still considerable deviations between the SVM-NL and kNN models and PBM predictions,especially for larger times, higher initial water contents, and for some textures (such as Silt Loam).
Finally, the DNN models with one, two, and three hidden layers (DNN-1, DNN-2, and DNN-3) were
tested to determine the optimum number of hidden layers and to minimize overfitting of output
variables. The best performance was obtained by the DNN model with two hidden layers (DNN-2)(R2 = 0.98). The DNN-1 and DNN-3 models tended to underestimate, and overestimate observed
values of runoff water volumes, respectively.
In conclusion, DNN techniques exhibited a better capability to reproduce the results of the PBM
compare with traditional ML techniques. In particular, results indicate that the simple linear models
are not well suited to develop correlations for runoff from agricultural lands because runoff is acomplex nonlinear hydrological process. Although nonlinear models like kNN and SVM-NL can fairly
accurately capture the surface water flow dynamic, the DNN models are more promising for dealing
with complex hydrologic problems. Applications of the DNN models with additional hidden layers
provide a possibility to further improve the model performance. Similarly, the trained DNN-2 model
provided an excellent prediction of near equilibrium solute transport.
It should be noted that the training databases from the PBM could also be augmented toinclude real monitoring data and/or more complex model formulations (e.g., multidimensional
simulations, stochastic simulations, other reactive transport processes, and physical non-equilibrium
flow and transport). This approach would preserve the benefits of data-driven models over PBMs
(e.g., faster execution times and ability to incorporate real data), while constraining data-driven models
to the underlying physics of hydrological processes.
Author Contributions: J.L. contributed to the development of the physically based and machine learning models,
the data analysis, and the manuscript writing. W.L. contributed to the machine learning models and the data
analysis. S.A.B. and J.Š. contributed to the physically based model, the data analysis, and the manuscript writing.
Funding: This research was supported by the USDA, ARS, National Program 211 and the Western Riverside
County Agricultural Coalition (WRCAC). Although partial funding for this work was provided by WRCAC,
no site-specific data from the San Jacinto Watershed have been used in this publication.
Conflicts of Interest: The authors declare no conflict of interest.

Water 2019, 11, 200

19 of 21

References
1.
2.
3.
4.
5.
6.
7.
8.
9.
10.
11.
12.

13.
14.
15.
16.
17.
18.
19.
20.
21.
22.
23.
24.
25.

Carpenter, S.R.; Caraco, N.F.; Correll, D.L.; Howarth, R.W.; Sharpley, A.N.; Smith, V.H. Nonpoint pollutionof surface waters with phosphorus and nitrogen. Ecol. Appl. 1998, 8, 559–568. [CrossRef]
United States Geological Survey (USGS). The Quality of Our Nation’s Waters-Nutrients and Pesticides;
US Geological Survey Circular 1225: Washington, DC, USA, 1999; p. 80.
Tyrrel, S.F.; Quinton, J.N. Overland flow transport of pathogens from agricultural land receiving faecal
wastes. J. Appl. Microbiol. 2003, 94, 87–93. [CrossRef]
Kirkby, M.J.; Beven, K.J. A physically based, variable contributing area model of basin hydrology.
Hydrol. Sci. J. 1979, 24, 43–69.
Donigian, A.S.; Love, J.T. The Connecticut Watershed Model—A tool for BMP impact assessment. Proc. Water
Environ. Fed. 2002, 2, 605–624. [CrossRef]
Borah, D.K.; Bera, M. Watershed-scale hydrologic and nonpoint-source pollution models: Review ofmathematical bases. Trans. ASAE 2003, 46, 1553. [CrossRef]
Roz, E.P. Water Quality Modeling and Rainfall Estimation: A Data Driven Approach; The University of Iowa:
Iowa City, IA, USA, 2011.
Hsu, K.L.; Gupta, H.V.; Sorooshian, S. Artificial neural network modelling of the rainfall-runoff process.
Water Resour. Res. 1995, 31, 2517–2530. [CrossRef]
Minns, A.W.; Hall, M.J. Artificial neural network as rainfall-runoff model. Hydrol. Sci. J. 1996, 41, 399–417.
[CrossRef]
Dawson, C.W.; Wilby, R. An artificial neural network approach to rainfall-runoff modelling. Hydrol. Sci. J.
1998, 43, 47–66. [CrossRef]
Dibike, Y.B.; Solomatine, D.; Abbott, M.B. On the encapsulation of numerical-hydraulic models in artificial
neural network. J. Hydraul. Res. 1999, 37, 147–161. [CrossRef]
Abrahart, R.J.; See, L. Comparing neural network and autoregressive moving average techniques for theprovision of continuous river flow forecast in two contrasting catchments. Hydrol. Process. 2000, 14,2157–2172. [CrossRef]
Mjolsness, E.; DeCoste, D. Machine learning for science: State of the art and future prospects. Science 2001,
293, 2051–2055. [CrossRef] [PubMed]
Govindaraju, R.S.; Rao, A.R. Artificial Neural Networks in Hydrology; Springer Science & Business Media:
Berlin, Germany, 2013.
Solomatine, D.P.; Ostfeld, A. Data-driven modelling: Some past experiences and new approaches.
J. Hydroinform. 2008, 10, 3–22. [CrossRef]
Remesan, R.; Mathew, J. Hydrological Data Driven Modelling: A Case Study Approach; Springer: Berlin, Germany,
2014; Volume 1.
Bai, Y.; Chen, Z.; Xie, J.; Li, C. Daily reservoir inflow forecasting using multiscale deep feature learning withhybrid models. J. Hydrol. 2016, 532, 193–206. [CrossRef]
Karandish, F.; Šimůnek, J. A comparison of numerical and machine-learning modeling of soil water content
with limited input data. J. Hydrol. 2016, 543, 892–909. [CrossRef]
Fang, K.; Shen, C.; Kifer, D.; Yang, X. Prolongation of SMAP to Spatiotemporally Seamless Coverage ofContinental US Using a Deep Learning Neural Network. Geophys. Res. Lett. 2017, 44, 11-030. [CrossRef]
Loague, K.M.; Freeze, R.A. A comparison of rainfall-runoff modeling techniques on small upland catchments.
Water Resour. Res. 1985, 21, 229–248. [CrossRef]
Karlsson, M.; Yakowitz, S. Nearest-neighbor methods for nonparametric rainfall-runoff forecasting.
Water Resour. Res. 1987, 23, 1300–1308. [CrossRef]
ASCE Task Committee on Application of Artificial Neural Networks in Hydrology. Artificial neural networks
in hydrology—I: Preliminary concepts. J. Hydrol. Eng. 2000, 5, 115–123. [CrossRef]
ASCE Task Committee on Application of Artificial Neural Networks in Hydrology. Artificial neural networks
in hydrology—II: Hydrologic applications. J. Hydrol. Eng. 2000, 5, 124–137. [CrossRef]
Lin, J.Y.; Cheng, C.T.; Chau, K.W. Using support vector machines for long-term discharge prediction.
Hydrol. Sci. J. 2006, 51, 599–612. [CrossRef]
Nourani, V.; Komasi, M.; Mano, A. A multivariate ANN-wavelet approach for rainfall–runoff modeling.
Water Resour. Manag. 2009, 23, 2877–2894. [CrossRef]

Water 2019, 11, 200

26.

27.
28.
29.
30.

31.
32.
33.
34.

35.
36.
37.

38.
39.

40.
41.

42.
43.
44.
45.

46.
47.

48.

20 of 21

Hinton, G.; Deng, L.; Yu, D.; Dahl, G.E.; Mohamed, A.R.; Jaitly, N.; Senior, A.; Vanhoucke, V.; Nguyen, P.;
Sainath, T.N.; et al. Deep neural networks for acoustic modeling in speech recognition: The shared views offour research groups. IEEE Signal Process. Mag. 2012, 29, 82–97. [CrossRef]
LeCun, Y.; Bengio, Y.; Hinton, G. Deep learning. Nature 2015, 521, 436. [CrossRef] [PubMed]
Goodfellow, I.; Bengio, Y.; Courville, A.; Bengio, Y. Deep Learning; MIT Press: Cambridge, UK, 2016; Volume 1.
Sahoo, G.B.; Ray, C.; Mehnert, E.; Keefer, D.A. Application of artificial neural networks to assess pesticidecontamination in shallow groundwater. Sci. Total Environ. 2006, 367, 234–251. [CrossRef]
Stenemo, F.; Lindahl, A.M.L.; Gärdenäs, A.; Jarvis, N. Meta-modeling of the pesticide fate model MACRO
for groundwater exposure assessments using artificial neural networks. J. Contam. Hydrol. 2007, 93, 270–283.
[CrossRef] [PubMed]
Fienen, M.N.; Nolan, B.T.; Feinstein, D.T. Evaluating the sources of water to wells: Three techniques formetamodeling of a groundwater flow model. Environ. Model. Softw. 2016, 77, 95–107. [CrossRef]
Marçais, J.; de Dreuzy, J.R. Prospective interest of deep learning for hydrological inference. Groundwater
2017, 55, 688–692. [CrossRef] [PubMed]
Blanning, R.W. The construction and implementation of metamodels. Simulation 1975, 24, 177–184. [CrossRef]
Schoumans, O.F.; Mol-Dijkstra, J.; Akkermans, L.M.W.; Roest, C.W.J. SIMPLE: Assessment of non-point
phosphorus pollution from agricultural land to surface waters by means of a new methodology.
Water Sci. Tech. 2002, 45, 177–182. [CrossRef]
Liang, J.; Bradford, S.A.; Šimůnek, J.; Hartmann, A. Adapting HYDRUS-1D to simulate overland flow andreactive transport during sheet flow deviations. Vadose Zone J. 2017, 16, 18. [CrossRef]
Ahuja, L.R.; Sharpley, A.N.; Lehman, O.R. Effect of Soil Slope and Rainfall Characteristics on Phosphorus in
Runoff. J. Environ. Qual. 1982, 11, 9–13. [CrossRef]
Sharda, V.N.; Singh, S.R.; Sastry, G.; Dhruvanarayana, V.V. A finite element model for simulating runoff
and soil erosion from mechanically treated agricultural lands: 2. Field validation and applications.
Water Resour. Res. 1994, 30, 2299–2310. [CrossRef]
Wallach, R.; Grigorin, G.; Rivlin, J. A comprehensive mathematical model for transport of soil-dissolved
chemicals by overland flow. J. Hydrol. 2001, 247, 85–99. [CrossRef]
Goodrich, D.C.; Burns, I.S.; Unkrich, C.L.; Semmens, D.J.; Guertin, D.P.; Hernandez, M.; Yatheendradas, S.;
Kennedy, J.R.; Levick, L.R. KINEROS2/AGWA: Model use, calibration, and validation. Trans. ASABE 2012,
55, 1561–1574. [CrossRef]Šimůnek, J.; Van Genuchten, M.T.; Šejna, M. Recent developments and applications of the HYDRUS computer
software packages. Vadose Zone J. 2016, 15, 25. [CrossRef]
Hromadka, T.V., II; Lai, C. Solving the two-dimensional diffusion flow model. In Hydraulics and Hydrology in
the Small Computer Age, Proceedings of the Specialty Conference on Hydraulics and Hydrology in the Small Computer
Age, Lake Buena Vista, FL, USA, 12–17 August 1985; ASCE: New York, NY, USA, 1985; pp. 555–562.
Panday, S.; Huyakorn, P.S. A fully coupled physically-based spatially-distributed model for evaluating
surface/subsurface flow. Adv. Water Resour. Res. 2004, 27, 361–382. [CrossRef]
Weill, S.; Mouche, E.; Patin, J. A generalized Richards equation for surface/subsurface flow modelling.
J. Hydrol. 2009, 366, 9–20. [CrossRef]
Green, W.H.; Ampt, G.A. Studies in soil physics. I. The flow of air and water through soils. J. Agric. Sci. 1911,4, 1–24.
Šimůnek, J.; Šejna, M.; Saito, H.; Sakai, M.; van Genuchten, M.T. The HYDRUS-1D Software Package forSimulating the Movement of Water, Heat, and Multiple Solutes in Variably Saturated Media, Version 4.0. HYDRUS
Software Series 3, Department of Environmental Sciences, University of California Riverside: Riverside, CA,
USA, 2008; 315.
Rawls, W.J.; Brakensiek, D.L.; Miller, N. Green-Ampt infiltration parameters from soils data. J. Hydrol. Eng.
1983, 109, 62–70. [CrossRef]
Pedregosa, F.; Varoquaux, G.; Gramfort, A.; Michel, V.; Thirion, B.; Grisel, O.; Blondel, M.; Prettenhofer, P.;
Weiss, R.; Dubourg, V.; et al. Scikit-learn: Machine learning in Python. J. Mach. Learn. Res. 2011, 12,2825–2830.
Chollet, F. Keras: Deep Learning Library for Theano and Tensorflow. 2015. Available online: https://keras.io/(accessed on the 21 January 2019).

Water 2019, 11, 200

49.
50.
51.
52.

53.
54.
55.
56.
57.
58.
59.
60.

21 of 21

Ticlavilca, A.M.; McKee, M. Multivariate Bayesian regression approach to forecast releases from a system ofmultiple reservoirs. Water Resour. Manag. 2011, 25, 523–543. [CrossRef]
Vapnik, V.N. Statistical Learning Theory; Wiley: New York, NY, USA, 1998.
Vapnik, V. The Nature of Statistical Learning Theory; Springer Science & Business Media: Berlin, Germany, 2013.
Van Looy, K.; Bouma, J.; Herbst, M.; Koestel, J.; Minasny, B.; Mishra, U.; Montzka, C.; Nemes, A.;
Pachepsky, Y.; Padarian, J.; et al. Pedotransfer functions in Earth system science: Challenges and perspectives.
Rev. Geophys. 2017, 55, 1199–1256. [CrossRef]
Friedman, J.; Hastie, T.; Tibshirani, R. The Elements of Statistical Learning; Springer Series in Statistics: Berlin,
Germany, 2001; Volume 1, pp. 337–387.
Maier, H.; Dandy, G.C. Neural networks for the prediction and forecasting of water resources variables:
A review of modeling issues and applications. J. Environ. Model. Softw. 2000, 15, 101–124. [CrossRef]
Sharma, V.; Negi, S.C.; Rudra, R.P.; Yang, S. Neural networks for predicting nitrate-nitrogen in drainage
water. Agric. Water Manag. 2003, 63, 169–183. [CrossRef]
Khalil, B.; Ouarda, T.B.M.J.; St-Hilaire, A. Estimation of water quality characteristics at ungauged sites using
artificial neural networks and canonical correlation analysis. J. Hydrol. 2011, 405, 277–287. [CrossRef]
Schmidhuber, J. Deep learning in neural networks: An overview. Neural Netw. 2015, 61, 85–117. [CrossRef][PubMed]Joachims, T. Learning to Classify Text Using Support Vector Machines: Methods, Theory and Algorithms;
Kluwer Academic Publishers: Norwell, MA, USA, 2002; Volume 186.
Dibike, Y.B.; Velickov, S.; Solomatine, D.; Abbott, M.B. Model induction with support vector machines:
Introduction and applications. J. Comput. Civ. Eng. 2001, 15, 208–216. [CrossRef]
Bray, M.; Han, D. Identification of support vector machines for runoff modelling. J. Hydroinform. 2004, 6,
265–280. [CrossRef]© 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open accessarticle distributed under the terms and conditions of the Creative Commons Attribution(CC BY) license (http://creativecommons.org/licenses/by/4.0/).