Journal of Hydrology 628 (2024) 130495

Contents lists available at ScienceDirect

Journal of Hydrology
journal homepage: www.elsevier.com/locate/jhydrol

Research papers

Coupling water, solute, and sediment transport into a new computationallyefficient hydrologic model
Lin Chen a, *, Jiří Šimůnek a, Scott A. Bradford b, Hoori Ajami a, Menberu B. Meles bab

Department of Environmental Sciences, University of California, Riverside, CA 92521, USA
USDA, ARS, Sustainable Agricultural Water Systems Unit, Davis, CA 95616, USA

A R T I C L E I N F O

A B S T R A C T

This manuscript was handled by Corrado Cor­
radini, Editor-in-Chief, with the assistance ofCarla Saltalippi, Associate Editor

While surface–subsurface interactions strongly impact water quality and quantity in ecosystems, only relatively
little attention has been given to solute mass transfers at this boundary. The KINEROS2 (K2) model solves the 1D
kinematic wave equation for water flow and the advection–dispersion equation for non-absorbing solute (orsediment) transport at the surface, whereas the HYDRUS-1D (H1D) model solves the Richards equation for flow
and the advection–dispersion equation for solute transport in the subsurface. We previously coupled water flowin the surface and subsurface environments using the K2 and H1D models (H1D-K2). Herein, an innovative
framework was developed to additionally couple solute transport at the surface–subsurface boundary in thesimple, accurate, and computationally efficient H1D-K2 model. To test the H1D-K2 model’s accuracy androbustness, simulation results are compared to analytical solutions, experimental data, and results from the twodimensional HYDRUS-2D (H2D) code. The relative differences in water and solute balance components between
the H2D and H1D-K2 models are less than 1% for homogeneous cases and 4% for heterogeneous cases. The
simulated hydrograph and chemograph are almost identical for the two models. Furthermore, simulated Bromide(Br) concentration in runoff is in good agreement with observations. The proposed model also accurately re­produces soil erosion under different rainfall patterns and intensities.

Keywords:
Solute transport
Overland flow
Sediment transport
Vadose flow
Computational efficiency
Surface/subsurface interaction

1. Introduction
Solute transport between the surface and subsurface greatly affects
water quality, nutrient cycling, biodiversity, and ecosystem dynamics(Weill et al., 2011; Kayhanian et al., 2012). Overland flow may carryvarious pollutants, such as phosphorous, nitrogen, pesticides, heavy
metals, sediment, or bacteria, directly to streams, rivers, lakes, and othersurface water bodies (Gao et al., 2011; Kayhanian et al., 2012; Miranda
et al., 2021; Brookfield et al., 2023). Additionally, polluted surface
water can infiltrate into the vadose zone and eventually contaminate
groundwater aquifers (Liggett et al., 2014; Sasidharan et al., 2021).
These processes pose a potential threat to aquatic ecosystems andhuman health. It is, therefore, imperative to understand and predict the
movement, fate, and transport of pollutants at the surface and in thesubsurface.
Hydrological models have been developed to investigate the in­teractions of water and solutes at the surface–subsurface boundary overa wide range of spatial and temporal scales (Havis et al., 1992; Shi et al.,2011; Brookfield et al., 2023). These models can be classified into four

categories: empirical models (e.g., Ahuja, 1986), mixing-layer models(e.g., Ahuja and Lehman, 1983; Havis et al., 1992; Zhang et al., 1999;
Gao et al., 2004; Yang et al., 2016), interfacial diffusion-controlledmodels (e.g., Wallach et al., 1988; Wallach and van Genuchten, 1990),
and physical-based hydrological models (e.g., Weill et al., 2011; Liggett
et al., 2015; Niu and Phanikumar, 2015).
Empirical models can be used to quickly and easily compute soluteconcentrations in runoff at an outlet (Ahuja, 1986; Vadas et al., 2008),
but they do not consider the complexity of physiochemical processes
since the entire field is considered as one unit (Zhang et al., 1997). Also,
most empirical relations are primarily designed to examine the solutetransport process at the surface (Vadas et al., 2008), and the strong in­teractions between the surface and subsurface are neglected. Moreover,this approach commonly requires extensive calibration at each site.
Similar to the soil particle erosion processes (Foster et al., 1977;
Hairsine and Rose, 1991), the mixing or exchange layer concept is
employed to describe solute transfer from the soil to runoff (Shi et al.,2011). The partition of solute between the surface water and infiltrationis largely determined by either the depth of mixing (Ahuja and Lehman,* Corresponding author.
E-mail address: linc@ucr.edu (L. Chen).
https://doi.org/10.1016/j.jhydrol.2023.130495
Received 9 July 2023; Received in revised form 26 September 2023; Accepted 3 November 2023
Available online 20 November 2023
0022-1694/© 2023 The Author(s). Published by Elsevier B.V. This is an open access article under the CC BY-NC-ND license (http://creativecommons.org/licenses/bync-nd/4.0/).

L. Chen et al.

Journal of Hydrology 628 (2024) 130495

1983), the degree of mixing (e.g., Zhang et al., 1999; Walter et al., 2007;
Yang et al., 2016), or the combination of mixing depth and degree(Havis et al., 1992; Shi et al., 2011). Alternatively, several authors haveproposed a mixing model to consider solute release by the ejection of soil
water into runoff due to the raindrop impact and the diffusion of solutein the exchange layer (e.g., Hairsine and Rose, 1991; Gao et al., 2004).
However, these models are mainly tested through laboratory scalemeasurements (less than several m2) (e.g., Ahuja and Lehman, 1983;
Gao et al., 2004; Walter et al., 2007; Yang et al., 2016). Large-scale
simulations using the mixing layer model are hampered by extensive
calibration requirements (e.g., mixing layer thickness and soil detach­ability), spatial distributions of solutes, and failure to adequately
represent the subsurface heterogeneity. Diffusion from the soil surface
into runoff is also neglected in these approaches.
Although interfacial diffusion-controlled models can consider thetransport of solutes due to diffusion and advection with infiltratingwater, the analytical solutions are driven by the assumption of steady
infiltration and an equilibrium state (e.g., Wallach et al., 1988; Wallach
and van Genuchten, 1990). These approaches neglect spatial heteroge­neity at the surface and in the subsurface, causing challenges in simu­
lating complex, nonlinear systems. Thus, hydrological applications onlarge scales are limited.
Physically-based hydrological models (e.g., CATHY, HydroGeo­
Sphere, PAWS) have been developed to capture the complexity ofphysiochemical processes and the coupling of water and solute in­teractions at the surface–subsurface boundary. A series of virtual ex­periments, ranging from hillslope to watershed scales, have been
performed to investigate the transport conditions at the surfa­
ce–subsurface boundary (Liggett et al., 2014; Liggett et al., 2015),
hydrograph separation or transit time distributions (Weiler and
McDonnell, 2004; Weill et al., 2011). To numerically obtain reliable andaccurate solutions, the high nonlinearity in the Richards equation
generally requires fine discretization (usually on the order of centime­ters) and small time steps, which leads to considerable execution times
and requires extensive computational resources and power (Šimůnek
et al., 2016). Additionally, since the surface and subsurface equations
are solved simultaneously, the implementation of integrated/fullycoupled hydrological models at the catchment or continent scales isgreatly restricted due to computation burden and extensive parameter­
ization (Weill et al., 2011; Hwang et al., 2014).
Similar approaches have been developed to understand soil erosion
and sediment yield processes (e.g., Hairsine and Rose, 1991; Gao et al.,2005; Goodrich et al., 2012; Pandey et al., 2016; Tao et al., 2017). Thesephysically-based models exhibit significant variations in terms of input
requirements, process representation, particularly in the subsurface, andrepresenting spatial and temporal variability of various processes(Pandey et al., 2016; Ebel et al., 2023). Many existing models simplify
the infiltration approach by using the Green-Ampt approximation (e.g.,
Kineros2), require extensive calibration (e.g., maximum shield depth),
and neglect subsurface heterogeneity (Merritt et al., 2003; Pandey et al.,2016; Ebel et al., 2023). Additionally, the application of these models tolarge scales remains limited due to the high computational demand andnumerical instability caused by continuous simulations for the entire
system (Brookfield et al., 2023).
Many strategies have been adapted and developed to enhance thecomputational efficiency of hydrological models of surface and subsur­
face water flow at large scales, such as dimensionality reduction (Seoet al., 2007; Chen et al., 2022), domain decomposition (Smith et al.,
1999; Goodrich et al., 2012), high-performance computing (Laloy and
Vrugt, 2012), parallel computing (Kollet and Maxwell, 2006), andmodularization method (Weill et al., 2011; Gatto et al., 2021). The
application of these strategies to coupled solute transport has received
less attention. The K2-O2 code employs domain decomposition, dividing
a large-scale catchment into a cascade of 1D elements (e.g., plane,
channel, etc.) to simulate water flow and nutrient cycling (Massart et al.,2010). Alternatively, complex physically based hydrological models canbe decomposed into well-defined modules or components, which can be
run independently and exchange data through specified interfaces(Gatto et al., 2021). For instance, Weill et al., (2011) developed asequential coupling module with a switching boundary conditions (BCs)algorithm to investigate the partitioning of solute and water between thesurface and subsurface. The four model components (i.e., water flow andsolute transport for the surface and subsurface) are solved sequentially(Sulis et al., 2010; Weill et al., 2011; Gatto et al., 2021). However, we
still lack standard procedures and benchmark test cases to better un­derstand the physiochemical processes at the surface and in thesubsurface.
The overall objective of this research is to couple water flow, solute,and sediment transport processes at the surface–subsurface boundary
into a simple, physically-based, and computationally efficient hydro­
logical model at the hillslope scale. The efficient coupling approach
between KINEROS2 (referred to below as K2) and HYDRUS-1D (referred
to below as H1D) (Chen et al., 2022) will be extended to simulate
advective and diffusive solute transport, soil erosion, and sedimenttransport. The specific objectives are (i) to develop a novel approach toexchange water and solute fluxes between the surface and subsurface;(ii) to evaluate the model performance by comparing its results withanalytical solutions, experiment data, and HYDRUS-2D (referred to
below as H2D) output.
2. Theory and hydrologic models
Water flow and solute transport on a hillslope are complex nonlinear
problems (Fig. 1). Rainfall or water inflow from intermediate hillslopes
either infiltrates into the subsurface or is routed downslope as overland
flow once the soil infiltration capacity is reached. The overland flow
velocity can be one to several orders of magnitude larger than theinfiltration rate (Havis et al., 1992). Additionally, the magnitude andpattern of the solute mass transfer from the soil into runoff are affectedby raindrop strikes, runoff scour, soil erosion, and diffusion, which
introduce further complexity in solving the coupled problem (Gatto
et al., 2021; Miranda et al., 2021).
To reduce the computational complexity in space and time, asequential solution procedure was used to solve the surface and sub­surface flow problems with automatic adjustment of the time step (Chenet al., 2022). The dimensionality reduction approach was applied torepresent the subsurface heterogeneity using multiple H1D soil profiles(Fig. 1d). The empirical infiltration relation in the previous K2 code was
replaced with the H1D model solution of the Richards equation using
realistic top and bottom BCs.
Outputs from the H2D model were employed as a reference forcomparing and evaluating the coupled H1D-K2 model’s performance.
The H2D model does not incorporate dimensionality reduction tech­niques to reduce the complexity of the problem (Fig. 1c). In addition,
H2D simultaneously solves the governing equations of water movement
and solute transport across the entire system (the surface andsubsurface).
2.1. Coupled H1D-K2 model
The coupled H1D-K2 model, newly developed by Chen et al., (2022),
simulates the coupling processes of water at the surface and in thevadose zone under variably saturated conditions. It has been tested
against analytical solutions (Gottardi and Venutelli, 1993) and suc­cessfully benchmarked against other integrated hydrologic models
within the Integrated Hydrologic Model Intercomparison Project (Sulis
et al., 2010; Liggett et al., 2012; Maxwell et al., 2014; Kollet et al., 2017;
Coon et al., 2020). The coupled model has successfully simulated
coupled surface–subsurface flows for homogeneous and heterogeneous
hillslopes over a wide range of conditions with different soil properties,
initial conditions, bottom BCs (e.g., free drainage, no flow), slopes, andmultiple rainstorm events. The coupled model also provides the spatial2

L. Chen et al.

Journal of Hydrology 628 (2024) 130495

Fig. 1. Conceptual model of (a) surface–subsurface solute transport processes, (b) soil erosion and sediment transport, (c) conceptual model of HYDRUS-2D (H2D),
and (d) conceptual model of coupled H1D-K2 model. Color zones in the subsurface represent the different soil properties. Red circles represent the surface nodes forcalculating overland flow and solute transport at the surface. (For interpretation of the references to color in this figure legend, the reader is referred to the web
version of this article.)and temporal variability of water content and can potentially consider
the evapotranspiration processes (Chen et al., 2022). In this paper, we
extended the functionality of the H1D-K2 model to simulate advective
and diffusive solute transport, soil erosion, and sediment transport. We
briefly describe the equations and coupling strategies below.

the discharge per unit width can be expressed as:Q=[ (2.1.2. Solute transport
Several erosion models have been modified to simulate non-sorbed
chemical transfer from the soil to overland flow (i.e., Ahuja, 1986;
Gao et al., 2004). A similar approach is employed herein for solutetransport in the K2 model. Transport of a nonreactive chemical at theland surface and in the subsurface is described by the advec­
tion–dispersion equations (van Genuchten et al., 2004; Šimůnek et al.,2016):

(1)

)]

∂θ ∂
∂hsub
=
K
+1
∂t ∂z
∂z

+ ϕw (z, t)(3)

where S is the mean local slope [-], nman is Manning’s roughness coef­
ficient for overland flow [m1/3s− 1], and ξ is a unit conversion factor. ξ =1 (m1/3 s− 1) for International System units and 1.49 (ft1/3 s− 1) for En­
glish units.

2.1.1. Water flow
The kinematic wave equation (Eq. (1) is used to describe overland
water flow (Woolhiser and Liggett, 1967), while the Richards equation(Eq. (2) is used to simulate variably saturated water flow in the sub­surface (van Genuchten et al., 2004; Šimůnek et al., 2016):∂hTop ∂Q+
= qr − qi + qw (x, t)∂t∂x

ξ 1/2 5/3
S hTopnman(2)where t is time [s], x is the distance coordinate over the soil surface [m],hTop is the unit storage of surface water [m], Q is the discharge per unit
width [m2/s], qr is the rainfall intensity [m/s], qi is the infiltration rate[m/s], qw is the rate of the water inflow for the uphill [m/s], θ is thevolumetric water content [-], hsub is the water pressure head in soil [m], z
is the spatial coordinate [m] (positive upward), ϕw is the sink term [s− 1],
and K is the unsaturated hydraulic conductivity function [m/s].
Discharge per unit width (Q) is related to the slope (S), surface
roughness, and flow conditions (laminar or turbulent flow) (Giráldez
and Woolhiser, 1996). Based on the Manning hydraulic resistance law,

∂hTop cs ∂Qcs+
= es + eh + qr cr − qi cs + ϕc (x, t)
∂t
∂x

(4)

∂θc ∂Js+= ϕs (z, t)
∂t
∂z

(5)

where cs is the solute concentration in surface water [mmol/m3], cr is theconcentration in rainfall water [mmol/m3], qi is subsurface water flux[m/s], ϕc is the rate of the solute inflow for the uphill [mmol/m2/s], c isthe concentration in the liquid phase in the subsurface [mmol/m3], D is3

L. Chen et al.

Journal of Hydrology 628 (2024) 130495

the dispersion coefficient [m2/s], ϕs is the term for mass transfer be­
tween the two regions [mmol/m3/s], eh is the rate of solute undergoingfirst-order decay [mmol/m2/s], and es is the rate of solute undergoingzero-order decay [mmol/m2/s]. When sediment transport is activated,
eh and es account for the erosion due to flowing water and raindrop
impact, respectively. The term Js is the solute mass flux density [mmol/m2/s] given by:
Js = qi c − θD∂c∂xparticle size (d) and density (γs). The shear velocity (v*) is defined as√̅̅̅̅̅̅̅̅̅̅̅̅̅
ghTop S(Woolhiser et al., 1990).
2.1.4. Coupling strategy
A MATLAB script was developed to run the H1D and K2 codes
sequentially with dynamic time steps. The total simulation time isdivided into smaller exchange time periods, Δtp, which are determinedby the subsurface heterogeneity, rainfall intensity and duration, andsurface characteristics (e.g., slope, roughness) (Chen et al., 2022). The
surface and subsurface modules exchange data at the end of each Δtp.
Here we extended the coupling approach to include solute exchange.
Fig. 2 presents a flowchart illustrating the coupling approach and massbalance in the coupled H1D-K2 model. More information about the
coupling strategy and time step setup can be found in Chen et al., (2022).
In general, the coupling sequence is summarized below:(6)

The dispersion coefficient is assumed to be the sum of the molecular
diffusion and hydrodynamic dispersion contributions (Wallach and van
Genuchten, 1990):
D = D0 λ + ε|v|(7)where D0 is the ionic or molecular diffusion coefficient in pure water[m2/s], λ is a tortuosity factor [-], ε is the dispersivity [m], and v is theaverage pore water velocity [m/s].
A modified third-type (solute flux) BC was imposed at the soil surface
to account for the convective and dispersive solute transfers between
surface runoff and the subsurface (Wallach and van Genuchten, 1990):
(
)⃒
∂c ⃒⃒
Js (0, t) ≡ qi c − θD
= − θk(c(0, t) − cr )(8)∂x ⃒z=0(1) Calculating actual water and solute fluxes using the H1D code
At the surface of each H1D soil profile, the BC is switched between an
atmospheric flux and a time-variable pressure head to reflect conditions
without or with ponding, respectively (Chen et al., 2022). When a soluteis released into the runoff, a modified third-type BC (Eq. (8) is applied
for solute transport, while a standard third-type BC is applied for otherscenarios. After that, the H1D code is executed to compute the balance
between the potential and actual fluxes and the solute flux with infil­
trating water. The pressure head and concentration distributions in thesubsurface are obtained under the applied water and solute BCs.

where k is a mass transfer coefficient between the soil and the runoff
water [m/s], which can be approximated by the approach provided byWallach et al., (1988).
Note that no reactions (e.g., retardation and degradation) areconsidered in the overland flow domain but can be considered in thesubsurface with the H1D model. This assumption is justified because ofmuch faster flow rates at the surface than in the subsurface.

(2) Performing water and solute mass balances at the interface
The water mass balance is evaluated first to determine whether thewater supply to the soil surface is greater than the cumulative infiltra­tion flux during Δtp (Fig. 2). When the soil infiltration capacity exceeds
the water supply, the upper BC is switched to an atmospheric BC withsurface runoff. The potential water flux (qsurf ) is obtained as the sum of
rainfall intensity and the average ponded water depth (hTop) divided by
Δtp, shown in Eq. (12) as:∑nk j
j=2 hTop
qjsurf = qjr +(12)(nk − 1)⋅Δtp

2.1.3. Soil erosion and sediment transport
The K2 sediment subroutine was modified to continuously simulate
the movement of eroded soil with surface water. The lateral water flux
was routed down the hillslope by the overland flow subroutine with the
net surface water flux input at each surface node provided by the H1Dmodel. Soil erosion on a hillslope was divided into two parts: raindrop
energy erosion (i.e., splash erosion) and hydraulic erosion caused byflowing water (Bennett, 1974; Woolhiser et al., 1990). The splash
erosion rate (es) is estimated as a function of the square of the rainfall
rate, while the hydraulic erosion rate (eh) is linearly approximated as thedifference between the equilibrium concentration (cm) and the current
sediment concentration (cs) (Woolhiser et al., 1990):{)(
af q2r exp − ah hTop , q > 0
(9)
es =
q⩽0
0,
eh = ag (cm − cs )hTop

where the superscript j is the surface node number starting from node 1
and ending with node nk (the total number of surface nodes). Steps (1)and (2) are repeated until achieving water balance.
When only advective flux is considered, the solute flux BC is appliedon the top surface. The actual solute flux is calculated as the product of
the infiltrating water flux and the solute concentration, shown in Eq.
(13). The solute mass that infiltrates into the subsurface is then removed
from the surface. The advective solute flux (qis ) is calculated as:(10)qjis = qji cjTop

where af is a contact related to soil and surface properties [-], ah is thedamping effectiveness of surface water [-], ag is a transfer rate coefficient[s− 1], and cm is the concentration at the equilibrium transport capacity[g/s]. When deposition occurs, ag is theoretically equal to the ratio of theparticle settling velocity (v) to the hydraulic depth, hTop.
The sediment transport capacity concentration (cm) is estimated as afunction of hydraulic conditions using the Engelund-Hansen relation(Engelund and Hansen, 1967):
cm =0.5vv3*

g2 hTop d(γ s − 1)2(13)

The solute concentration (cTop) in the infiltrating water is determinedusing a solute mass balance calculation at the top of each H1D profile toaccount for the contributions from the rainfall and surface ponding as:∑j=nk j j
∑j=nk j j
j=1 qr cr ⋅Δt +j=1 hTop cs
jcTop = ∑j=nk j(14)∑j=nk j
q
⋅Δt
+
r
j=1
j=1 hTopA loop is performed to compute the net surface water flux and solute flux
for each H1D soil profile.

(11)(3) Distributing net surface water and solute fluxes to each individual
surface node

where v* is the shear velocity [m/s], d is particle diameter [m], γ s is thesuspended specific gravity of the particle [g/m3], and g is gravitational
acceleration [m/s2]. Particle settling velocity (v) is calculated fromAs the exchange flux in an H1D profile corresponds to vertical inflow
4L. Chen et al.

Journal of Hydrology 628 (2024) 130495

Fig. 2. Flowchart illustrating the main processes of the coupled H1D-K2 code. The orange color represents computations in the MATLAB scripts handling thecoupling between the KINEROS2 (K2) and HYDRUS-1D (H1D) models. (For interpretation of the references to color in this figure legend, the reader is referred to the
web version of this article.)

at multiple surface nodes, a weighting factor approach is used to
distribute the corresponding values of water and solute fluxes to indi­vidual surface nodes using Eqs. (14) and (15), respectively.
The net surface flux (qlw) is obtained by:

qjlw (t) =

⎧
⎪
⎪
⎨

qjr − qji ,
hjTop

⎪
⎪
⎩ ∑j=nk jh /(nk − 1)
j=2 Topif qjr ⩾qji]
[ jj
qr − qi , if qjr < qjisurface is the sum of the splash erosion (es) and hydraulic erosion (eh)(Eqs. (9) and (10). A loop through spatial nodes is performed to computeerosion and deposition rates for the current particle class. The concen­
tration distribution at the surface is obtained for the current particle
class. Finally, a loop over different particle sizes is performed.
During the calculation period of the K2 code (i.e., Δtp), the compu­
tation time interval is automatically adjusted to adhere to the Courant
condition, which ensures reliable and accurate numerical solutions thatavoid numerical instability (Goodrich et al., 2012).

(15)

The total solute flux (qls) is the sum of contributions from the rainfall andsurface ponding:qjls (t) =

⎧
⎪
⎪
⎨

qjr cjr − qji cjs ,⎪⎪⎩ ∑j=nk (j=2

if qjr cjr ⩾qji cjTop][ j j
hjTop cjs
qr cr − qji cjTop , if qjr cjr < qji cjTop)/(nk − 1)

(5) Updating the upper BCs and initial conditions
The calculated pressure heads and solution concentrations are set as
the initial conditions for the H1D simulation at the next Δtp. The upper
BCs are updated based on the calculation results of water depth (hTop)and concentration (cTop). A temporal loop is executed until the end of thetotal simulation time (Fig. 2).

(16)

hjTop cjs

(4) Executing the K2 code

2.2. HYDRUS-2D with the overland flow and solute transport module

The surface water and transport equations (Eqs. (1) and (3) are

The H2D (version 3) model with the newly developed overland flow
and solute transport routine was used to simulate the coupled processesof water flow and solute transport over a hillslope (van Genuchten et al.,
2004; Šimůnek et al., 2016). At the surface, the 1D kinematic-wave andadvection–dispersion equations were used for describing overland flow
and solute transport, respectively, with an effective source/sink from thesubsurface (van Genuchten et al., 2004; Šimůnek et al., 2016). In thesubsurface, the Richards and advection–dispersion equations were usedjjsolved using qls (t) and qlw (t) as source terms to compute the surface
water depth and concentration distribution at the surface. A loop
through surface spatial nodes is performed to compute the net surface
solute flux (Fig. 2). When sediment transport is activated, the sedimentconcentration is solved with current lateral discharge rates, current
ponding water depths, and previous values of the sediment concentra­tion (cs) for each particle size class (Eq. (4). The erosion rate at the soil
5

L. Chen et al.

Journal of Hydrology 628 (2024) 130495

to describe variably saturated water flow and transport processes,
respectively (Šimůnek et al., 2016).
A switching BC approach between the variable pressure head
(Dirichlet) and atmospheric-controlled (Neumann) BCs is used to reflect
the surface ponding condition. A free drainage BC was implemented at
the bottom to enable the unrestricted drainage of water from the domain
under the influence of gravity. Additionally, a third-type (Cauchy) BC
was implemented at the top and bottom surface to represent a specified
solute flux along the boundary. This type of BC is mass-conservative andmore physically realistic than a first-type BC (Parker and van Gen­uchten, 1984). At the surface, the surface solute flux (Eq. (6) is the sumof diffusion and advection fluxes. At the domain bottom, when the water
flux is zero or directed out of the region, third-type BCs are automati­
cally switched to Neumann BCs.

variable pressure head BC. A free drainage BC was imposed at the bot­
tom to allow soil water to flow freely out of the domain by gravity. There
was no inflow from uphill. The initial pressure head distribution was
generated by hydrostatic equilibrium with a given top pressure head of− 1 m, while the initial concentration was assumed to be uniform andequal to 0.0 mmol/m3 in the entire domain. The water content andconcentration distributions in the subsurface from the H1D-K2 simula­tions were generated by the Kriging interpolation with a linear semivariogram model. The hydrograph and chemograph at the outlet were
presented and compared to evaluate the model performance. The
simulated surface concentrations at three observation points (p1, p2, andp3) were compared. Locations for p1, p2, and p3 are found in the center ofthe uphill, middle, and downhill sections, respectively.
3.2. Solute release from the soil into surface runoff

3. Model validation and benchmark comparisons

The modified H1D-K2 code considers solute exchange at the inter­face between the overland and subsurface domains. At the top surface,
the solute flux (Eq. (6) is the sum of advection (i.e., infiltration in thisstudy) and diffusive (i.e., the product of the mass transfer coefficient andthe solute gradient) fluxes. To assess the performance of the coupled
H1D-K2 model, analytical solutions (Wallach and van Genuchten, 1990)
and experimental data were employed for comparison. Three different
soil types (i.e., Ruston fine sandy loam, Ruston loam, and Parson’s clay)
under three infiltration conditions (i.e., no infiltration, full drainage,
and reduced infiltration) were tested (Table S1). See the supporting
information for details about the methods. The coefficient of determi­
nation (R2) and root-mean-squared error (RMSE) were used to estimate
the temporal dynamics and error, respectively, in predicting the surface
concentration and surface solute flux compared to the analyticalsolutions.

The proposed model was validated by comparing the outputs of the
proposed model with the results of the virtual experiment simulated by
the H2D model, analytical solutions, and experiment data. Three sce­narios were considered in this study: (i) overland flow with solutetransport due to rainfall, (ii) solute release from the soil into surface
runoff, and (iii) soil erosion and sediment transport.
3.1. Overland flow with solute transport due to rainfall
The coupled processes of water flow and solute transport on a hill­slope were simulated and compared using the coupled H1D-K2 and H2Dmodels. Two runoff generation mechanisms, including infiltration
excess (Horton) and saturation excess (Dunne), were investigated underhomogeneous and heterogeneous conditions. Similar to the dimension­ality reduction approach provided by Chen et al., (2022), only one H1D
profile was employed for the homogeneous case, while three H1D pro­files were selected to reflect the spatial heterogeneity for the heteroge­
neous case (Fig. 1d).
For both homogeneous and heterogeneous cases, the hillslope wassubjected to uniform rainfall with different concentrations at different
hillslope locations. The soil transect was 400 m long with a thickness of
1 m and a slope of 0.0005. Manning’s roughness coefficient, nman, wasequal to 0.01. A uniform rainfall of 3.8 × 10− 4 m/min was applied for
200 mins, followed by a 100 mins recession period. The concentrations
in rainfall were 1.0, 0.5, and 0.2 mmol/m3 for the uphill (cr1), middle(cr2), and downhill (cr3) sections, respectively. The solute can move into
the soil profile with infiltrating water. The soil hydraulic and solutetransport parameters are listed in Table 1.
For the coupled H1D-K2 model, the subsurface domain was divided
into 100 elements. The mesh size varied from 0.0005 m near the surface
to 0.02 m at the bottom. The horizontal mesh size for the overland planemodel was 10 m (Fig. 1a). Conversely, a regular finite element mesh wasused to discretize the domain for the H2D model (Fig. 1c). A uniform
mesh size of 1 m was used in the horizontal direction. The same vertical
discretization as in H1D was used in the vertical direction of the H2Dmodel. The simulation output is reported at 1-min intervals.
BCs on the surface were switched between an atmospheric BC and a

3.3. Soil erosion and sediment transport
The coupled H1D-K2 model was modified to simulate soil erosion
and sediment transport, and its outputs were compared with theexperimental data (Tao et al., 2017). The experimental site was locatedin the rain-fed cropping region of the Loess Plateau in China (35◦ 14′ N,
107◦ 41′ E). Artificial rainfall experiments on small plots (1 × 1 m2) wereconducted to obtain the runoff and erosion information. The soil consistsof approximately 3.47 % clay, 92.26 % silt, and 4.27 % sand with a bulk
density of 1.41 g/cm3. The plot was exposed to a spatially uniform
simulated rainfall for 90 mins. Two rainfall scenarios were selected andcompared to reproduce the processes of runoff generation, soil erosion,and sediment transport. In the first scenario, the rainfall intensity was
uniform, with a value of 130 mm/h. In the second scenario, multiple
rainstorms were applied with a value of 100 mm/h during the first 30
min, 160 mm/h during the second 30-minute period, and 130 mm/h
during the last 30 min. Detailed information can be found in Tao et al.,(2017). Furthermore, the Nash-Sutcliffe efficiency (NSE) and RMSE
between the simulated and observed results were used to evaluate theperformance of the coupled H1D-K2 model in simulating discharge andsediment yield processes.
4. Results

Table 1
Soil hydraulic and solute transport parameters for the virtual experiment.

Soil
1
Soil
2

Ks (m/
min)

θsat
(-)

θres
(-)

n
(-)

α (1/

l (-)

DL
(m)

DT
(m)

6.94 ×10− 4
6.94 ×
10− 6

0.4

0.08

2

1

0.5

0.1

40

0.4

0.08

2

1

0.5

0.1

40

m)

4.1. Overland flow with solute transport due to rainfall
4.1.1. Homogeneous subsurface, non-uniform rainfall concentrations
Fig. 3 shows cumulative water and solute fluxes and concentrations
from the H1D-K2 and H2D models for a homogeneous hillslope sub­jected to uniform rainfall with different concentrations (cr1 = 1.0 mmol/
m3, cr2 = 0.5 mmol/m3, cr3 = 0.2 mmol/m3) at different locations.
Fig. 3a shows that surface water concentrations as a function of time atthree observation nodes were almost identical for the H1D-K2 and H2Dmodels. The chemograph simulated by the coupled H1D-K2 model

Note: n, α, and l are the van Genuchten parameters; DL and DT are the longitu­
dinal and transverse dispersivities, respectively. For the H1D calculation, only
DL is required for the calculation.
6

L. Chen et al.

Journal of Hydrology 628 (2024) 130495

Fig. 3. Cumulative water and solute fluxes and concentrations for a homogeneous hillslope subjected to uniform rainfall with different concentrations (cr1 = 1.0mmol/m3, cr2 = 0.5 mmol/m3, cr3 = 0.2 mmol/m3): (a) concentrations at three observation nodes located uphill (p1), downhill (p3), and in the middle (p2), (b)variations of water and solute fluxes at the outlet, (c) components of the water mass balance per unit width, and (d) components of the solute mass balance per unit
width. Solid and dashed lines represent the simulation results of the coupled H1D-K2 and H2D models, respectively. Note: the uniform rainfall occurred for the first
200 min.

shows that the uphill surface concentration, which was initially equal to
0.0 mmol/m3, had a step increase at the early stage and then experi­
enced a gradual rise. After reaching the maximum of 0.55 mmol/m3 at
237 mins, the surface concentration gradually decreased to 0.50 mmol/
m3. Similar trends were observed for the middle and downhill sections.
The peak times were delayed by 25 and 42 mins in the middle anddownhill sections, respectively. When the rainfall stopped, the concen­
trations in the middle and downhill sections experienced a sudden in­crease since the only source of water was from uphill. Differences in thepeak time between the coupled H1D-K2 and H2D models were 5.6 mins,
4.9 mins, and 3.3 mins for the uphill, middle, and downhill sections,
respectively. The differences in concentrations at peak times at the threelocations were less than 0.02 mmol/m3 between the two models.
The simulated hydrograph and chemograph were identical for thetwo hydrological models (Fig. 3b). The water discharge rate showed asteep rise from 0.0 m2/s at 1.5 mins to 2.15 × 10− 3 m2/s at 72.5 mins,
remained constant until the start of the recession period, and gradually
dropped to zero at 279.0 mins. Similarly, the solute flux at the outlet
experienced a gradual increase at the early stage and then steeply
increased from 2.34 × 10− 4 mmol/m/s at 39.5 mins to 1.14 × 10− 3mmol/m/s at 76.0 mins. It remained constant until the rainfall stopped
and then experienced a gradual decrease to zero at 279.0 mins. The steep
increase at 39.5 mins was due to water mixing with inflow from the
middle and uphill sections. As expected, the magnitude of solute flux

was much smaller than the water flux because the average rainfall
concentration on the surface was smaller than 1.0 mmol/m3.
Fig. 3c shows that the temporal behavior of water balance compo­
nents of overland flow simulated by the two models was identical. The
water balance errors were 0.01 × 103 kg for both the H2D and coupled
H1D-K2 models (Table 2). Relative differences in each water balancecomponent were less than 1.0 % between the two models (Table 2).
Similarly, Fig. 3d shows that the coupled H1D-K2 model exactly
reproduced the temporal behaviors of the solute mass balance compo­
nents simulated by the H2D model. By comparing the two models, thesolute balance errors were about 0.01 mmol, and the relative differences
between different solute balance components were less than 1.0 %(Table 2).
4.1.2. Heterogeneous subsurface, non-uniform rainfall concentrations
Fig. 4 shows cumulative water and solute fluxes and concentrations
simulated by the H1D-K2 and H2D models for a heterogeneous hillslope
subjected to uniform rainfall with different concentrations at different
locations (cr1 = 1.0 mmol/m3, cr2 = 0.5 mmol/m3, cr3 = 0.2 mmol/m3).
Compared to the homogeneous case (Fig. 3), one can see distinct dif­ferences in temporal dynamics of surface concentrations between loca­
tions in the middle and downhill sections for both models, while asimilar behavior is observed in the uphill section (Fig. 4a).
Infiltration excess occurred in the uphill section (Fig. 4a), and the7

L. Chen et al.

Journal of Hydrology 628 (2024) 130495

remarkably similar behavior was observed for the coupled H1D-K2
model in the uphill section (Fig. 4a). The surface concentration
reached its peak value of 0.54 mmol/m3 at 237.0 mins and eventually
decreased to 0.50 mmol/m3 at 300 mins. The two models produced
identical behavior because rainfall was the only water source available
for infiltration in this section (no inflow was considered in this study).
In the middle section, where excess saturation occurred (Fig. 4a), thesurface concentration simulated by the H2D model can be categorized
into four stages. No ponding existed during the first stage, and the sur­face concentration at p2 gradually increased from zero to 0.26 mmol/m3
at 43.0 mins. The soil surface was ponded during the second stage,
which triggered a steep increase in the infiltration rate at 43.0 mins. As aresult, the soil surface concentration started suddenly rising. Rainfall
stopped, and thus overland flow was the only source of infiltration
during the third stage, i.e., in the period between 200 and 300 mins.
Surface water from the uphill section had a high concentration (cr1 =1.0 mmol/m3), and this caused a steep rise of the surface concentration
by 0.07 mmol/m3 from 200.0 mins to 213.9 mins. Overland flow left thispart of the hill quickly, and no water ponded at the top surface after
256.3 mins during the fourth stage.
The surface concentration (Fig. 4a) experienced a slight decrease due
to free drainage at the bottom of the domain, which allows water andsolute to flow freely downward. Conversely, the surface concentration
simulated by the coupled H1D-K2 model experienced a steep increase
and reached 0.75 mmol/m3 at 200.0 mins. During the recession period,
the surface concentration had a sudden increase in a short period,
reached a peak value of 0.83 mmol/m3 at 273.0 mins, then followed a

Table 2
Water and solute balance components for homogeneous hillslope simulated by
two models.

Cum. rain (103 kg)
Surface volume (103 kg)
Cum. runoff (103 kg)
Cum. infiltration (103 kg)
Water balance error* (103kg)
Cum. solute rain (mmol)
Surf. solute volume(mmol)
Cum. solute runoff(mmol)
Cum. solute infiltration(mmol)
Solute balance error*(mmol)

H2D
model

Coupled H1D-K2
model

Relative
Difference

30.40
0.00
23.82
6.59
− 0.01

30.40
0.00
23.74
6.65
0.01

0.00 %
0.00 %
0.35 %
0.97 %
–

17.48
0.00

17.48
0.00

0.00 %
0.00 %

12.29

12.26

0.28 %

5.18

5.23

1.01 %

0.01

0.00

–

Note: water and solute balance errors are calculated for the net cumulativevalues of fluxes across the top boundary. The negative sign illustrates the total
inflow (i.e., cumulative rain) at the surface is smaller than the total outflow (e.g.,
the sum of surface volume, cumulative runoff, and cumulative infiltration).

simulated surface concentration by the H2D model, therefore, experi­
enced a steep increase at the early stage, followed by a gradual increase.
After the disappearance of surface water, the surface concentration
decreased slightly to 0.49 mmol/m3 at the end of the simulation time. A

Fig. 4. Cumulative water and solute fluxes and concentrations for a heterogeneous hillslope subjected to uniform rainfall with different concentrations (cr1 = 1.0mmol/m3, cr2 = 0.5 mmol/m3, cr3 = 0.2 mmol/m3): (a) concentrations at three observation points located uphill (p1), downhill (p3), and in the middle (p2), (b)variations of water and solute flux at the outlet, (c) components of the water mass balance per unit width, and (d) components of the solute mass balance per unit
width. Solid and dashed lines represent the simulation results of the coupled H1D-K2 and H2D models, respectively. Note: the uniform rainfall occurred for the first
200 min.
8

L. Chen et al.

Journal of Hydrology 628 (2024) 130495

slight decrease, and eventually dropped to 0.82 mmol/m3 at 300.0 mins.
Both models exhibited similar behavior in the downhill section(Fig. 4a). The temporal behavior of the surface concentration can bedivided into four distinct stages: (1) the surface concentration experi­
enced a slightly increasing trend within the first 110.0 mins; (2) a steep
increase was observed at 110.0 mins due to the arrival of the second
overland flow wave that was generated by saturation excess in the
middle section; (3) another steep increase occurred at 200.0 mins due tohigh concentration runoff from the mid-section, which was the only
source of infiltrating water after the cessation of rainfall; and finally (4)a slight decrease in concentration was observed after runoff disappeared
because of the free drainage condition at the bottom.
Fig. 4b shows that the hydrograph simulated by the two models
agreed closely. The simulated hydrograph by the H2D model wascharacterized by two “plateaus” of quasi-steady flow: (1) the firstplateau is due to infiltration excess in the downhill section and no inflow
from the upper part. Initially water ponding occurred quickly and wasrouted to the outlet, resulting in a quick increase in the hydrograph(Fig. 4b). After that, outflow at the outlet reached a quasi-equilibrium
state. During this period, no runoff was routed to the downhill section
due to the high soil infiltration capacity in the mid-section; (2) satura­tion was fully reached in the mid-section at 105.8 mins. The runoff wasrouted quickly to the downhill section and caused a step increase in theoutlet flux. After a short rising period, a second quasi-equilibrium state
was achieved. Eventually, the water flux experienced a quick decreaseduring the recession period. The coupled H1D-K2 model reproduced
similar behavior. The arrival time of the second runoff wave caused by
the mid-section was 118.5 mins. The relative difference in peak flows
simulated by the two models was less than 1.87 %.
Chemographs simulated by the two models closely agreed (Fig. 4b).
The solute flux simulated by the H2D model arrived at 2.8 mins, grad­ually increased to 1.54 × 10− 4 m2/min at 37.1 min, and then remained
relatively constant before the appearance of the second runoff wave,
which was generated by excess saturation in the middle section. The
solute flux experienced a gradual increase at 106.7 min, reached its
peak, exhibited a second plateau, and then a slow recession once rainfall
stopped. Only a small difference in the arrival time of the second wave
was observed when comparing the two models. The relative differencein the peak solute flux simulated by the two models was 6.3 %.
Almost identical temporal behaviors of water (Fig. 4c) and solute(Fig. 4d) balance components were obtained for the two models. The
relative difference of each water balance component was less than 3.0%, while the corresponding value for the solute balance components wasless than 4.0 % (Table 3).
Fig. 5 shows water content and solute concentration distributions at
different times simulated by the two models for the heterogeneous
profile scenario. The H2D results show that infiltration excess induced
runoff first in the uphill and downhill sections. At the early stage, rainfall
was the only water source available for infiltration at these two sections(Fig. 5a). Rainfall reached the soil surface and infiltrated into the sub­surface, which increased soil water contents and solute concentrations.
The subsurface concentration in the uphill section was uniformly
distributed in the horizontal direction during the entire simulation time
because there was no inflow from the upslope (Fig. 5b). This section was
not fully saturated during the entire simulation time (Fig. 5a). Water
content and concentration distributions simulated by the H1D-K2 model
were nearly identical to the H2D outputs in the uphill section.
The infiltration capacity is higher than the rainfall intensity in the
middle section. In this case, rainfall rapidly infiltrated into the soil
instead of flowing over the surface (Fig. 5a). The runoff generated in the
uphill section was routed to the middle section and infiltrated into its
subsurface. Additional water from the uphill section contributed to thedevelopment of a saturated zone at the upper part of this area, which
enlarged and extended downhill and towards the bottom over time.
Saturation-excess runoff occurred at 83.7 mins in the middle section.
After the cessation of rainfall (200 mins), water contents experienced aTable 3
Simulated water and solute balance components for heterogenous hillslope by
two models.

Cum. rain (103 kg)
Surface volume (103 kg)
Cum. runoff (103 kg)
Cum. infiltration (103 kg)
Water balance error* (103kg)
Cum. solute rain (mmol)
Surf. solute volume(mmol)
Cum. solute runoff(mmol)
Cum. solute infiltration(mmol)
Solute balance error*(mmol)

H2D
model

Coupled H1D-K2
model

Relative
Difference

30.40
0.00
10.73
19.83
− 0.16

30.40
0.00
10.75
20.36
− 0.700.00 %
0.00 %
0.15 %
2.65 %
–

17.48
0.00

17.48
0.00

0.00 %
0.00 %

2.65

2.68

1.20 %

14.99

15.50

3.37 %− 0.16− 0.70–Note: water and solute balance errors are calculated by the net cumulative valuesof fluxes across top boundary. The negative sign illustrates the total inflow (i.e.,
cumulative rain) at the top surface is smaller than the total outflow (e.g., the sumof surface volume, cumulative runoff, and cumulative infiltration).

slight decrease due to the drainage condition applied at the bottom of
the soil profile. Solute concentrations in the subsurface of the middle
section increased quickly due to high infiltration rates, a mixture of
rainfall and overland inflow with a high concentration (Fig. 5b). The
movement of the saturation front (Fig. 5a) and the high-concentration
zone progressed from the upper right (upslope surface) to lower left
(downslope bottom) portions of the middle section (Fig. 5b). Water
drainage at the bottom caused the downward movement of the con­centration front during the recession period (200 to 300 mins). The
coupled H1D-K2 model accurately captured the horizontally averaged
characteristics of the saturation (Fig. 5c) and concentration (Fig. 5d)
fronts.
Similar water content (Fig. 5c) and concentration (Fig. 5d) distri­butions were observed in the downhill section of both models. In theearly stage, water contents were uniformly distributed in the horizontal
direction, and the saturation front moved gradually toward the bottom(Fig. 5c). A similar pattern was observed for the concentration distri­bution (Fig. 5d). After saturation excess occurred in the middle section
at 83.7 mins, excess water was routed downward to the downhill section
and infiltrated into its subsurface, which triggered non-uniform distri­butions of water contents and solute concentrations (Fig. 5c). A highconcentration zone developed in the upper part of this section andgradually extended to the bottom and downslope (Fig. 5d).
4.2. Solute release from the soil into surface runoff
4.2.1. Comparison with the analytical solution
The simulated solute concentrations in the soil by the coupled H1DK2 model have an excellent agreement (R2 > 0.99, RMSE < 2.7 × 10− 4)with the analytical solutions at different time steps when the infiltration
rate was 0.28 cm/h (Fig. 6). Similarly, solute fluxes from the soil into
runoff water simulated by the coupled H1D-K2 model were almost
identical (R2 > 0.94, RMSE < 6.8 × 10− 3) to the analytical solutions fordifferent infiltration conditions (Fig. 7).
4.2.2. Comparison with the experimental data
Fig. 8 shows observed (Ahuja and Lehman, 1983) and H1D-K2
simulated Bromide (Br) concentrations in the runoff as a function oftime. A generally good agreement was obtained for the no infiltration
and reduced infiltration cases. For the free drainage conditions, thevalue of Br concentration (c/c0) experienced a sharp decrease, dropping
to less than 0.05 within the first 2 mins and eventually reaching a valueof 0.0 after 15 mins. The largest discrepancies between the simulation
9

L. Chen et al.

Journal of Hydrology 628 (2024) 130495

Fig. 5. Simulation results for the heterogeneous subsurface subject to rainfall with different concentrations (cr1 = 1.0 mmol/m3, cr2 = 0.5 mmol/m3, cr3 = 0.2 mmol/m3) at different locations and with a free drainage BC at the bottom, (a) water content and (b) concentration distributions simulated by the H2D model, and (c) watercontent and (d) concentration distributions simulated by the coupled H1D-K2 model for different times (increasing from top to bottom).

results and observations were observed at the early stage (t < 5 mins) of
the runoff process.

5. Discussions
The coupled H1D-K2 model produces similar hydrochemical be­
haviors at the surface and in the subsurface as the H2D model, analyticalsolutions, and experimental data. In the first set of scenarios involving
virtual experiments, the hydrographs and chemographs simulated by thetwo models were almost identical for both homogeneous (Fig. 3) andheterogeneous (Fig. 4) hillslopes subjected to uniform rainfall withspatially different concentrations. Small differences between water andsolute components simulated by the two models (Fig. 3 and Fig. 4)
indicate that the H1D-K2 model can accurately capture the strong in­teractions and movement of water and solutes between the surface andsubsurface in space and time. Additionally, a modified third-type (soluteflux) BC applied at the surface–subsurface boundary ensures mass con­servation (Fig.s 3b, 3c, 4b and 4c, and Tables 2 and 3), as the total mass
exchange between the surface and subsurface is full controlled at a
specific time step. This mass-conservative BC (Eq. (6) can accurately
represent the convective and dispersive solute transfers between surface
runoff and subsurface flow (Fig.s 7 and 8). Therefore, it is more physi­
cally realistic than a first type (concentration) BC, which may cause

4.3. Soil erosion and sediment transport
The coupled H1D-K2 model accurately reproduced the observed
discharge rate and sediment yield rate for uniform rainfall with NSE
values of 0.85 and 0.60, respectively (Fig. 9a). The corresponding valuesof RMSE were 1.44 cm2/min and 6.56 g/min. A sharp increase in therunoff and sediment fluxes was observed at an early stage, followed by aquasi-equilibrium stage (Fig. 9a). Similarly, under multiple rainstorms(Fig. 9b), the simulation results showed good agreement with the ob­servations, with NSE values of 0.87 and 0.53 for the outflow and sedi­ment discharge rates, respectively. The corresponding values of RMSE
were 1.98 cm2/min and 17.52 g/min. Moreover, the differences in
estimated total sediment yields between the model and observations
were 1.15 kg and 0.31 kg for uniform rainfall and multiple rainstorm
scenarios, respectively.

10

L. Chen et al.

Journal of Hydrology 628 (2024) 130495

Fig. 6. Calculated solute concentration distributions versus depth during rainfall-induced runoff at an infiltration rate I of 0.28 cm/h (D = 0.09 cm2/h, k = 0.83 cm/h). Solid lines represent the analytical solutions given by Wallach & van Genuchten (1990), and the dashed lines represent the output from the H1D model.

spatial arrangement of net surface water and solute fluxes within thesurface nodes. Also, the same water and solute volumes enter the sub­surface (into one H1D soil profile) due to the averaging of water andsolute concentrations at multiple surface nodes. Previous studies also
demonstrate similar conclusions regarding the limitation of this timelagged approach (e.g., Huang and Yeh, 2009; Weill et al., 2011; Chenet al., 2022). Therefore, the H1D-K2 model with three soil profiles only
broadly approximates a heterogeneous hillslope’s surface water andsolute front movement.
When a solute release into the runoff is considered, the large
discrepancy for the free drainage conditions at the early stage (Fig. 8) isdue to the uncertainty in the measured hydraulic conductivity caused by
the soil swelling due to wetting. As mentioned in Ahuja and Lehman(1983), soil swelled by nearly 0.5 ~ 2.0 cm upon wetting, which caused
higher values of hydraulic conductivity compared to those measured
before the experiment. Moreover, this discrepancy is also related toaccelerated diffusion due to the mechanical raindrop impact that isneglected in the model (Zhang et al., 1999). The quick decrease insimulated Br concentrations at the early stage in this scenario is attrib­
uted to the dominance of the downward advection flux. Thus, the low
value of the Br concentration in the subsurface soil does not significantly
contribute to the pollution of runoff water (Fig. 8). A similar conclusion
has been reported in previous studies (Wallach and van Genuchten,1990).
While the proposed model captures the horizontally averaged dis­tributions of water contents and solute concentrations (Fig. 5), consid­
ering additional vertical soil columns would yield more precise watercontents and solute concentrations distributions in the subsurface.
Similarly, compared to the H2D results, no steep increase occurred in thesecond period (Fig. 4a) due to the averaging and redistributing of pon­ded water at multiple surface nodes in the mid-section of the H1D soil
profile. In our previously published paper (Chen et al., 2022), we
demonstrated that the coupled model, wherein one soil profile is applied
to each surface node, accurately reproduces the water movement
behavior at the surface and in the subsurface. However, when the spatial
distributions of water contents and solute concentrations in the sub­surface are not a primary concern, the coupled H1D-K2 model, with one
soil profile for the entire overland flow plane, is appropriate formodeling surface–subsurface water and solute transport due to its higher
computational efficiency. This coupled model can sufficiently simulate

Fig. 7. Calculated solute fluxes from the soil into runoff water at infiltration
rates I of 0.0, 0.28, 2.75 cm/h (k = 0.83 cm/h). Solid lines represent theanalytical solutions given by Wallach & van Genuchten (1990), and the dashedlines represent the output from the H1D model.

unrealistic mass movements due to concentration gradients (Liggett
et al., 2014). Previous studies presented similar conclusions (e.g., Parker
and van Genuchten, 1984; Liggett et al., 2014).
Discrepancies in surface (Fig. 4a) and subsurface (Fig. 5) concen­
tration distributions in the middle section are linked to the sequential
coupling approach and dimensionality reduction (i.e., the hillslope
subsurface is represented by one or multiple vertical H1D profiles). First,
the sequential coupling approach causes the time lag (Δtp) in updating
hydraulic and solute transport BCs of the H1D profiles. The subsurface
solutions (Eqs. (2) and (5) are obtained with previous values of hTop andcTop. This also contributes to the time delay in the arrival of the second
runoff wave and solute flux between the two models (Fig. 4). Second,
when dimensionality reduction is applied, one H1D profile corresponds
to multiple surface nodes, which may result in uncertainties in the11

L. Chen et al.

Journal of Hydrology 628 (2024) 130495

Fig. 8. Calculated and measured Br concentrations in the runoff for three different soils at different times, (a) Ruston fine sandy loam (k = 0.83 cm/h), (b) Ruston
loam (k = 0.90 cm/h), (c) Parson’s clay (k = 0.79 cm/h). Dashed lines represent the simulation results by the coupled H1D-K2 model, and symbols represent theexperimental data from Ajuha & Lehman (1983). The infiltration rates (I) for each scenario are listed in Table S1.

Fig. 9. Calculated and measured runoff rates and sediment fluxes under different rainfall events, (a) uniform rainfall intensity (qr = 130 mm/h), and (b) multiple
rainstorms (qr1 = 100 mm/h, qr2 = 160 mm/h, qr3 = 130 mm/h). Solid lines represent the simulation results by the coupled H1D-K2 model, and symbols represent
the experimental data from Tao et al. (2017).

the overall hydrochemical dynamics of the entire system for mostpractical problems.
The capabilities of the H1D-K2 model have been extended to simu­
late soil erosion and sediment transport, in addition to both advective
and diffusive solute transport. The coupled H1D-K2 model can accu­
rately capture the sediment yield process under different rainfall pat­terns and intensities (Fig. 9). The impact of raindrop strikes can be
identified by the sharp increase in the sediment flux at the early stage,
followed by an equilibrium stage when the transport capacity thresholdis reached. Similar behaviors of erosion yield rates were observed else­
where (e.g., Fox and Bryan, 2000; Polyakov et al., 2018). Furthermore,
the soil surface is not eroded uniformly due to the surface and subsurface
heterogeneities, which can explain the deviations between the obser­
vation and simulation results during the quasi-equilibrium stage (Fig. 9).
The combined techniques of the sequential coupling method, dy­
namic time-stepping, and dimensionality reduction lead to an increasein computational efficiency and a more robust model structure.
Compared to fully coupled hydrological models, where all governing
equations are solved simultaneously, the coupled H1D-K2 model re­solves surface and subsurface equations independently and sequentiallywith automatically adjusting time steps, and the results of eachmodeling component are exchanged at the end of Δtp through a MATLAB

script. Additionally, the 2D subsurface domain is represented by mul­
tiple individual H1D profiles (Fig. 1), which facilitates the imple­
mentation of parallel processing and greatly reduces the number ofvariables or parameters. Computationally, the coupled H1D-K2 model
demonstrates significantly higher efficiency compared to the H2Dmodel. The computation time has been significantly reduced in our test
cases, ranging from 30 % to 90 % compared to the H2D model. Similarly,
our previously published paper showed that the coupled H1D-K2 model
exhibits comparable computational efficiency to the H2D model(without dimensionality reduction techniques, i.e., a few soil profiles per
hillslope) across various benchmark cases with different slopes, bottomBCs, soil properties, and rainfall patterns and intensities (Chen et al.,2022). Furthermore, the proposed model has externally coupled thelatest version of the K2 and H1D models, which include a more robustmodel structure and additional enhancements (e.g., geometry of
microtopography at the surface, automatically adjusting time step size,and dispersive solute transfers at the interface) (Wallach and van Gen­uchten, 1990; Goodrich et al., 2012; Ebel et al., 2023). By contrast, thenumerical schemes, parameters, and detailed microtopography consid­
erations could potentially induce numerical oscillations and conver­gence issues and subsequently affect the computational accuracy insolving the 2D kinematic-wave overland flow and advection–dispersion
12

L. Chen et al.

Journal of Hydrology 628 (2024) 130495

equations (e.g., Glaister, 1993; Singh and Bhallamudi, 1998; Jaber andMohtar, 2003; Heppner et al., 2006; Gatto et al., 2021). For example,
Jaber and Mohtar (2003) demonstrated the necessity of using small-time
steps and proper mesh sizes (e.g., a small ratio between horizontal andvertical mesh sizes) to achieve stable and accurate solutions. As the H2Dmodel does not have an erosion module, it is impossible to assess thecomputational efficiency of the coupled H1D-K2 for this process.
The proposed model can be applied to predict the hydrological andchemical processes at the hillslope scale, estimate groundwater
recharge, and assess groundwater pollution resulting from infiltration ofpolluted surface water. The current version of the H1D-K2 model will be
extended to the watershed or regional scales by coupling with the
MODFLOW program to simulate groundwater water flow and solutetransport (e.g., Seo et al., 2007; Beegum et al., 2018). In this approach,
one can divide the watershed into a series of overland plane and channelelements, as is typically done in the K2 model (e.g., Woolhiser et al.,
1990; Goodrich et al., 2012). The coupled H1D-K2 model has important
advantages over the previous K2 model, including the ability to simulate
the saturation excess runoff, to consider a nonreactive chemical transfer,and to provide spatial and temporal distributions of water and solutes.
Also, future developments in the model should take into account the
coupled processes of surface–subsurface contaminant transport model
considering erosion and sorption. Also, the consideration of the colloids(e.g., bacteria and clay) and virus transports into the soil could be ad­vantageous for implementing sustainable water resource management.

CRediT authorship contribution statement
Lin Chen: Conceptualization, Formal analysis, Methodology, Soft­ware, Writing – original draft, Visualization. Jiří Šimůnek: Conceptu­
alization, Funding acquisition, Supervision, Project administration,
Writing – review & editing. Scott A. Bradford: Conceptualization,
Funding acquisition, Supervision, Project administration, Writing – re­view & editing. Hoori Ajami: Conceptualization, Funding acquisition,
Supervision, Project administration, Writing – review & editing. Men­beru B. Meles: Conceptualization, Writing – review & editing.
Declaration of competing interest
The authors declare that they have no known competing financial
interests or personal relationships that could have appeared to influence
the work reported in this paper.
Data availability
Data will be made available on request.
Acknowledgment
This research was supported by the United States Department ofAgriculture, Agricultural Research Service, National Program 211.

6. Conclusions

Appendix A. Supplementary data

Water flow and solute transport are complex nonlinear processes thatinvolve strong interactions between the surface and subsurface domains.
We developed an innovative framework to couple the surface–subsur­
face water-solute flow processes into a simple, accurate, and computa­
tionally efficient hydrological model at the hillslope scale. The coupled
H1D-K2 model can accurately capture the spatial and temporal varia­
tions of solute concentration in runoff and soil water under a range of
infiltration conditions.
Three benchmark problems were conducted, encompassing a broad
range of conditions, to verify, evaluate, and analyze the performance ofthe proposed models. The hydrographs, chemographs, water balancecomponents, and solute balance components simulated by the coupled
H1D-K2 model are in excellent agreement with the outputs of the H2Dmodel, analytical solutions and experimental data. The coupled H1D-K2
model can simultaneously simulate the diffusive mass transfer from thecontaminated soil to runoff due to the concentration gradient andadvection flux triggered by infiltrating water. Furthermore, by replacing
the empirical infiltration equation in the Kineros2 model with theRichards equation in the H1D model, the proposed coupled model
accurately captures the complex interactions between rainfall-induced
runoff, soil detachment, sediment transport, and deposition. Its
computational efficiency is greatly improved by applying the adaptive
time stepping scheme, dimensionality reduction, and sequential
coupling strategy. Furthermore, the generally good agreement between
the proposed model with the H2D model strengthens confidence in theextension and implementation of the coupled H1D-K2 model to large
scales (e.g., catchment) and other conditions.
When detailed spatial variations in soil moisture and concentration
are not the main focus, the coupled H1D-K2 model is more appropriate
for investigating the coupled processes of surface–subsurface water andsolute transport and the estimation of groundwater recharge due to its
higher computational efficiency and less data requirements compared to
the fully integrated hydrological models. The proposed model effec­
tively captures the overall hydrochemical dynamics of the entire system
for most practical problems. Future research focusing on the transport of
the finest sediment fractions (e.g., bacteria and clay) into the soil and theadsorption–desorption behavior of solutes would be beneficial for the
implement of sustainable water management.

Supplementary data to this article can be found online at https://doi.
org/10.1016/j.jhydrol.2023.130495.
References
Ahuja, L., Lehman, O., 1983. The extent and nature of rainfall-soil interaction in therelease of soluble chemicals to runoff. J. Environ. Qual. 12, 34–40. https://doi.org/
10.2134/jeq1983.00472425001200010005x.
Ahuja, L. (1986). Characterization and modeling of chemical transfer to runoff. Advancesin Soil Science: Volume 4, Springer. DOI: 10.1007/978-1-4613-8612-4_3.
Beegum, S., Šimůnek, J., Szymkiewicz, A., Sudheer, K., Nambi, I.M., 2018. Updating the
Coupling Algorithm between HYDRUS and MODFLOW in the HYDRUS Package for
MODFLOW. Vadose Zone J. 17 (1), 180034. https://doi.org/10.2136/
vzj2018.02.0034.
Bennett, J.P., 1974. Concepts of mathematical modeling of sediment yield. Water Resour.
Res. 10 (3), 485–492. https://doi.org/10.1029/WR010i003p00485.
Brookfield, A.E., Ajami, H., Carroll, R.W.H., Tague, C., Sullivan, P.L., Condon, L.E., 2023.
Recent advances in integrated hydrologic models: Integration of new domains.
J. Hydrol. 620, 129515 https://doi.org/10.1016/j.jhydrol.2023.129515.
Chen, L., Šimůnek, J., Scott, B., Ajami, H., Meles, M.B., 2022. A computationally efficienthydrologic modeling framework to simulate surface-subsurface hydrologicalprocesses at the hillslope scale. J. Hydrol. 614, 128539 https://doi.org/10.1016/j.
jhydrol.2022.128539.
Coon, E.T., Moulton, J.D., Kikinzon, E., Berndt, M., Manzini, G., Garimella, R.,
Lipnikov, K., Painter, S.L., 2020. Coupling surface flow and subsurface flow incomplex soil structures using mimetic finite differences. Adv. Water Resour. 144,
103701 https://doi.org/10.1016/j.advwatres.2020.103701.
Ebel, B.A., Shephard, Z.M., Walvoord, M.A., Murphy, S.F., Partridge, T.F., Perkins, K.S.,
2023. Modeling post-wildfire hydrologic response: review and future directions forapplications of physically based distributed simulation. Earth’s Future. 11 (2),e2022EF003038 https://doi.org/10.1029/2022EF003038.
Engelund, F. and E. Hansen (1967). A monograph on sediment transport in alluvial
streams. Technical University of Denmark 0stervoldgade 10, Copenhagen K.
Foster, G.R., Meyer, L., Onstad, C., 1977. An erosion equation derived from basic erosion
principles. Transactions of the ASAE. 20 (4), 678–0682. https://doi.org/10.13031/
2013.35627.
Fox, D.M., Bryan, R.B., 2000. The relationship of soil loss by interrill erosion to slope
gradient. Catena 38 (3), 211–222. https://doi.org/10.1016/S0341-8162(99)000727.
Gao, G., Falconer, R.A., Lin, B., 2011. Numerical modelling of sediment–bacteriainteraction processes in surface waters. Water Res. 45 (5), 1951–1960. https://doi.
org/10.1016/j.watres.2010.12.030.
Gao, B., Todd Walter, M., Steenhuis, T.S., Hogarth, W.L., Parlange, J.Y., 2004. Rainfall
induced chemical transport from soil to runoff: theory and experiments. J. Hydrol.
295 (1), 291–304. https://doi.org/10.1016/j.jhydrol.2004.03.026.
Gao, B., Walter, M.T., Steenhuis, T.S., Parlange, J.Y., Richards, B.K., Hogarth, W.L.,
Rose, C.W., 2005. Investigating raindrop effects on transport of sediment and non-13

L. Chen et al.

Journal of Hydrology 628 (2024) 130495
Niu, J., Phanikumar, M.S., 2015. Modeling watershed-scale solute transport using an
integrated, process-based hydrologic model with applications to bacterial fate andtransport. J. Hydrol. 529, 35–48. https://doi.org/10.1016/j.jhydrol.2015.07.013.
Pandey, A., Himanshu, S.K., Mishra, S.K., Singh, V.P., 2016. Physically based soil erosion
and sediment yield models revisited. Catena 147, 595–620. https://doi.org/
10.1016/j.catena.2016.08.002.
Parker, J.C., van Genuchten, M.T., 1984. Flux-Averaged and Volume-Averaged
Concentrations in Continuum Approaches to Solute Transport. Water Resour. Res. 20(7), 866–872. https://doi.org/10.1029/WR020i007p00866.
Polyakov, V., Stone, J., Holifield Collins, C., Nearing, M.A., Paige, G., Buono, J., GomezPond, R.L., 2018. Rainfall simulation experiments in the southwestern USA using the
Walnut Gulch Rainfall Simulator. Earth Syst. Sci. Data. 10 (1), 19–26. https://doi.
org/10.5194/essd-10-19-2018.
Sasidharan, S., Bradford, S.A., Šimůnek, J., Kraemer, S.R., 2021. Virus transport fromdrywells under constant head conditions: A modeling study. Water Res. 197, 117040
https://doi.org/10.1016/j.watres.2021.117040.
Seo, H.S., Šimůnek, J., Poeter, E.P., 2007. Documentation of the hydrus package formodflow-2000, the us geological survey modular ground-water model. IGWMCInternational Ground Water Modeling Center.
Shi, X., Wu, L., Chen, W., Wang, Q., 2011. Solute transfer from the soil surface tooverland flow: a review. Soil Sci. Soc. Am. J. 75 (4), 1214–1225. https://doi.org/
10.2136/sssaj2010.0433.
Šimůnek, J., M. T. Van Genuchten and M. J. V. Z. J. Šejna (2016). Recent developments
and applications of the HYDRUS computer software packages. Vadose Zone Journal.
15(7): vzj2016. 2004.0033. DOI: 10.2136/vzj2016.04.0033.
Singh, V., Bhallamudi, S.M., 1998. Conjunctive surface-subsurface modeling of overland
flow. Adv. Water Resour. 21 (7), 567–579. https://doi.org/10.1016/S0309-1708(97)00020-1.
Smith, R.E., Goodrich, D.C., Unkrich, C.L., 1999. Simulation of selected events on the
Catsop catchment by KINEROS2: A report for the GCTE conference on catchment
scale erosion models. Catena 37 (3), 457–475. https://doi.org/10.1016/S0341-8162(99)00033-8.
Sulis, M., Meyerhoff, S.B., Paniconi, C., Maxwell, R.M., Putti, M., Kollet, S.J., 2010.
A comparison of two physics-based numerical models for simulating surface
water–groundwater interactions. Adv. Water Resour. 33 (4), 456–467. https://doi.
org/10.1016/j.advwatres.2010.01.010.
Tao, W., Wu, J., Wang, Q., 2017. Mathematical model of sediment and solute transport
along slope land in different rainfall pattern conditions. Sci. Rep. 7 (1), 1–11.
https://doi.org/10.1038/srep44082.
Vadas, P., Owens, L., Sharpley, A., 2008. An empirical model for dissolved phosphorus in
runoff from surface-applied fertilizers. Agr Ecosyst Environ 127 (1–2), 59–65.
https://doi.org/10.1016/j.agee.2008.03.001.
van Genuchten, M. T., J. Šimůnek and R. Feddes (2004). Integrated modeling of vadosezone flow and transport processes. Unsaturated Zone Modelling: Progress,
Challenges and Applications, Eds. R. A. Feddes, G. H. de Rooij, and J. C. van Dam,
Wageningen UR Frontis Series, Vol. 6, Chapter 2, pp. 37- 69, x-xi, Kluwer Academic
Publishers, Dordrecht, The Netherlands, 2004.
Wallach, R., Jury, W.A., Spencer, W.F., 1988. Transfer of chemicals from soil solution tosurface runoff: A diffusion-based soil model. Soil Sci. Soc. Am. J. 52 (3), 612–618.
https://doi.org/10.2136/sssaj1988.03615995005200030002x.
Wallach, R., van Genuchten, M.T., 1990. A physically based model for predicting solutetransfer from soil solution to rainfall-induced runoff water. Water Resour. Res. 26(9), 2119–2126. https://doi.org/10.1029/WR026i009p02119.
Walter, M.T., Gao, B., Parlange, J.-Y., 2007. Modeling soil solute release into runoff with
infiltration. J. Hydrol. 347 (3–4), 430–437. https://doi.org/10.1016/j.
jhydrol.2007.09.033.
Weiler, M., McDonnell, J., 2004. Virtual experiments: a new approach for improving
process conceptualization in hillslope hydrology. J. Hydrol. 285 (1), 3–18. https://doi.org/10.1016/S0022-1694(03)00271-3.
Weill, S., Mazzia, A., Putti, M., Paniconi, C., 2011. Coupling water flow and solutetransport into a physically-based surface–subsurface hydrological model. Adv. Water
Resour. 34 (1), 128–136. https://doi.org/10.1016/j.advwatres.2010.10.001.
Woolhiser, D.A., Liggett, J.A., 1967. Unsteady, one-dimensional flow over a plane—The
rising hydrograph. Water Resour. Res. 3 (3), 753–771. https://doi.org/10.1029/
WR003i003p00753.
Woolhiser, D. A., S. R.E. and G. D.C. (1990). KINEROS, A Kinematic Runoff and Erosion
Model: Documentation and User Manual, Manual. U.S. Department of Agriculture,
Agricultural Research Service, ARS-77, 130 pp.
Yang, T., Wang, Q., Wu, L., Zhao, G., Liu, Y., Zhang, P., 2016. A mathematical model forsoil solute transfer into surface runoff as influenced by rainfall detachment. Sci. Total
Environ. 557–558, 590–600. https://doi.org/10.1016/j.scitotenv.2016.03.087.
Zhang, X.C., Norton, D., Nearing, M.A., 1997. Chemical transfer from soil solution tosurface runoff. Water Resour. Res. 33 (4), 809–815. https://doi.org/10.1029/
96WR03908.
Zhang, X., Norton, L., Lei, T., Nearing, M., 1999. Coupling mixing zone concept withconvection-diffusion equation to predict chemical transfer to surface runoff.
Transactions of the ASAE. 42 (4), 987. https://doi.org/10.13031/2013.13280.

sorbed chemicals from soil to surface runoff. J. Hydrol. 308 (1), 313–320. https://doi.org/10.1016/j.jhydrol.2004.11.007.
Gatto, B., Paniconi, C., Salandin, P., Camporese, M., 2021. Numerical dispersion of solutetransport in an integrated surface–subsurface hydrological model. Adv. Water
Resour. 158, 104060 https://doi.org/10.1016/j.advwatres.2021.104060.
Giráldez, J., Woolhiser, D., 1996. Analytical integration of the kinematic equation for
runoff on a plane under constant rainfall rate and Smith and Parlange infiltration.
Water Resour. Res. 32 (11), 3385–3389. https://doi.org/10.1029/96WR02106.
Glaister, P., 1993. Flux difference splitting for open-channel flows. Int. J. Numer. Meth.
Fluids 16 (7), 629–654. https://doi.org/10.1002/fld.1650160706.
Goodrich, D., Burns, I., Unkrich, C., Semmens, D.J., Guertin, D., Hernandez, M.,
Yatheendradas, S., Kennedy, J.R., Levick, L.R., 2012. KINEROS2/AGWA: Model use,
calibration, and validation. Trans. ASABE 55 (4), 1561–1574. https://doi.org/
10.13031/2013.42264.
Gottardi, G., Venutelli, M., 1993. A control-volume finite-element model for twodimensional overland flow. Adv. Water Resour. 16 (5), 277–284. https://doi.org/
10.1016/0309-1708(93)90019-C.
Hairsine, P., Rose, C., 1991. Rainfall detachment and deposition: Sediment transport in
the absence of flow-driven processes. Soil Sci. Soc. Am. J. 55 (2), 320–324. https://doi.org/10.2136/sssaj1991.03615995005500020003x.
Havis, R., Smith, R., Adrian, D., 1992. Partitioning solute transport between infiltrationand overland flow under rainfall. Water Resour. Res. 28 (10), 2569–2580. https://doi.org/10.1029/92WR01366.
Heppner, C.S., Ran, Q., VanderKwaak, J.E., Loague, K., 2006. Adding sediment transport
to the integrated hydrology model (InHM): Development and testing. Adv. Water
Resour. 29 (6), 930–943. https://doi.org/10.1016/j.advwatres.2005.08.003.
Huang, G., Yeh, G.-T., 2009. Comparative study of coupling approaches for surface waterand subsurface interactions. J. Hydrol. Eng. 14 (5), 453–462. https://doi.org/
10.1061/(ASCE)HE.1943-5584.0000017.
Hwang, H.T., Park, Y.J., Sudicky, E.A., Forsyth, P.A., 2014. A parallel computationalframework to solve flow and transport in integrated surface–subsurface hydrologic
systems. Environ. Model. Softw. 61, 39–58. https://doi.org/10.1016/j.
envsoft.2014.06.024.
Jaber, F.H., Mohtar, R.H., 2003. Stability and accuracy of two-dimensional kinematicwave overland flow modeling. Adv. Water Resour. 26 (11), 1189–1198. https://doi.
org/10.1016/S0309-1708(03)00102-7.
Kayhanian, M., Fruchtman, B.D., Gulliver, J.S., Montanaro, C., Ranieri, E., Wuertz, S.,
2012. Review of highway runoff characteristics: Comparative analysis and universal
implications. Water Res. 46 (20), 6609–6624. https://doi.org/10.1016/j.
watres.2012.07.026.
Kollet, S.J., Maxwell, R.M., 2006. Integrated surface–groundwater flow modeling: A freesurface overland flow boundary condition in a parallel groundwater flow model.
Adv. Water Resour. 29 (7), 945–958. https://doi.org/10.1016/j.
advwatres.2005.08.006.
Kollet, S., Sulis, M., Maxwell, R.M., Paniconi, C., Putti, M., Bertoldi, G., Coon, E.T.,
Cordano, E., Endrizzi, S., Kikinzon, E., Mouche, E., Mügler, C., Park, Y.-J.,
Refsgaard, J.C., Stisen, S., Sudicky, E., 2017. The integrated hydrologic modelintercomparison project, IH-MIP2: A second set of benchmark results to diagnose
integrated hydrology and feedbacks. Water Resour. Res. 53 (1), 867–890. https://doi.org/10.1002/2016WR019191.
Laloy, E., Vrugt, J.A., 2012. High-dimensional posterior exploration of hydrologic
models using multiple-try DREAM (ZS) and high-performance computing. Water
Resour. Res. 48 (1) https://doi.org/10.1029/2011WR010608.
Liggett, J.E., Werner, A.D., Simmons, C.T., 2012. Influence of the first-order exchange
coefficient on simulation of coupled surface–subsurface flow. J. Hydrol. 414–415,
503–515. https://doi.org/10.1016/j.jhydrol.2011.11.028.
Liggett, J.E., Werner, A.D., Smerdon, B.D., Partington, D., Simmons, C.T., 2014. Fully
integrated modeling of surface-subsurface solute transport and the effect ofdispersion in tracer hydrograph separation. Water Resour. Res. 50 (10), 7750–7765.
https://doi.org/10.1002/2013WR015040.
Liggett, J.E., Partington, D., Frei, S., Werner, A.D., Simmons, C.T., Fleckenstein, J.H.,
2015. An exploration of coupled surface–subsurface solute transport in a fully
integrated catchment model. J. Hydrol. 529, 969–979. https://doi.org/10.1016/j.
jhydrol.2015.09.006.
Massart, J., Guertin, D., Smith, R., Goodrich, D., Unkrich, C., Levick, L., 2010. K2–O2(KINEROS-Opus) spatially based watershed hydrologic and biogeochemical model.
Joint Federal Interagency Conf. on Sedimentation and Hydrologic Modeling, Proc.
Maxwell, R.M., Putti, M., Meyerhoff, S., Delfs, J.-O., Ferguson, I.M., Ivanov, V., Kim, J.,
Kolditz, O., Kollet, S.J., Kumar, M., Lopez, S., Niu, J., Paniconi, C., Park, Y.-J.,
Phanikumar, M.S., Shen, C., Sudicky, E.A., Sulis, M., 2014. Surface-subsurface modelintercomparison: A first set of benchmark results to diagnose integrated hydrology
and feedbacks. Water Resour. Res. 50 (2), 1531–1549. https://doi.org/10.1002/
2013WR013725.
Merritt, W.S., Letcher, R.A., Jakeman, A.J., 2003. A review of erosion and sedimenttransport models. Environ. Model. Softw. 18 (8–9), 761–799. https://doi.org/
10.1016/S1364-8152(03)00078-1.
Miranda, L.S., Wijesiri, B., Ayoko, G.A., Egodawatta, P., Goonetilleke, A., 2021. Watersediment interactions and mobility of heavy metals in aquatic environments. Water
Res. 202, 117386 https://doi.org/10.1016/j.watres.2021.117386.

14