Journal of Hydrology 614 (2022) 128539

Contents lists available at ScienceDirect

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

Research papers

A computationally efficient hydrologic modeling framework to simulate
surface-subsurface hydrological processes at the hillslope scale
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 ofPhilip Brunner, Associate Editor

Considering surface and subsurface interactions is imperative to predict water movement, water quantity andquality in the environment. However, substantial execution time and over-parameterization presently limit theapplicability of integrated hydrologic models at larger scales. Herein, a new efficient coupling routine for onedimensional (1D) surface–subsurface modeling was developed by externally coupling two widely used opensource codes. KINEROS2 (K2) solves the 1D kinematic wave equation for overland flow, and HYDRUS-1D(H1D) solves the Richards equation for subsurface flow. A novel approach, combined with water balance andboundary condition switching, is used to account for surface ponding and water exchange between the twomodel domains. A weighting factor related to the surface water depth is used to assign the time series of exchange
fluxes in one H1D profile to multiple surface nodes in the K2 calculation. This novel approach enables us torepresent the entire subsurface below each overland-flow plane by one vertical soil column with different soil
layering. When soil properties vary horizontally, multiple soil columns can represent spatial heterogeneities. The
performance of the coupled H1D-K2 model is examined by comparing simulation results with the HYDRUS-2D(H2D) model for six benchmark problems. The solution’s robustness, stability, and accuracy are assessed for awide variety of cases, including multiple rainstorms, different slopes, heterogeneous subsurface, and different
bottom boundary conditions resulting in 41 cases. The simulated hydrographs, surface water levels, and waterbalance components are all in good agreement. Relative mass balance errors are very small, typically notexceeding 4.0% for the coupled H1D-K2 model. Compared with the H2D model, the coupled model achieves afactor of 1.3–10 speedup by applying one H1D soil profile to an entire overland flow plane rather than to in­dividual surface nodes.

Keywords:
Hydrologic model
Computational efficiency
Surface runoff
Vadose zone flow
Water balance
Surface/subsurface interactions

1. Introduction
Global water demands for industrial, domestic, and agriculturalusers have increased by 600 % over the past century due to population
growth and economic development (UN World Water Development
Report, 2018). Climate change and anthropogenic activities have also
influenced water quantity and quality (Haddeland et al., 2014). For
example, climate change shifts precipitation patterns and increases thefrequency of flooding and droughts (Vorosmarty et al., 2000). Managed
aquifer recharge is a promising adaptation measure to reduce water
vulnerability and maintain its sustainability (Niswonger et al., 2017,
Sasidharan et al., 2018, Marwaha et al., 2021). This water management
approach can maximize water storage in groundwater aquifers and in­crease water supply during droughts (Niswonger et al., 2017, Kourakos

et al., 2019). However, there are limitations in available modeling re­sources for exploring integrated flow processes and impacts associatedwith managed aquifer recharge programs, climate change, etc. It isimperative to understand surface and subsurface water flow interactions
to improve sustainable water resource management (Smith et al., 2004,
Fan et al., 2019).
A large number of models have been developed to predict the hy­
drologic behavior of coupled surface and subsurface systems in response
to human activity and climate variability (Woolhiser et al., 1970, Freeze,
1972, Smith and Hebbert, 1983, Qu and Duffy, 2007). These hydrologic
models can be divided into four broad categories: empirical models,
lumped conceptual models, semi-distributed models (e.g., AGWA/KIN­
EROS2, SWAT), and physically-based distributed models (e.g., HYDRUS(2D/3D), HydroGeoSphere, MIKE SHE, ParFlow-CLM, PIHM,* Corresponding author at: Department of Environmental Sciences, University of California, Riverside, USA.
E-mail address: linc@ucr.edu (L. Chen).
https://doi.org/10.1016/j.jhydrol.2022.128539
Received 8 July 2022; Received in revised form 14 September 2022; Accepted 23 September 2022
Available online 1 November 20220022-1694/© 2022 Elsevier B.V. All rights reserved.

L. Chen et al.

Journal of Hydrology 614 (2022) 128539

WASH123D) (Smith et al., 2004, Goodrich et al., 2012, Lawrence et al.,2019, Liggett et al., 2012, Orth et al., 2015, Šimůnek, 2015, Kollet et al.,2017, Fan et al., 2019). In comparison to empirical and lumped models,
physically-based distributed hydrologic models are valuable tools forconsidering the full complexity of governing physiochemical processes(Singh and Woolhiser, 2002, Smith et al., 2004, Fatichi et al., 2016) andare capable of predicting the spatial and temporal distribution of waterflow and pollutant transport at hillslope to watershed scales (Šimůnek,2015, Kollet et al., 2017, Wen et al., 2021).
Benchmark simulation cases, ranging from simple homogenous caseswith constant rainfall rates to more complex cases with subsurface
heterogeneity and lateral flow, have been utilized to compare the per­formance of different physically-based distributed coupled/integratedhydrologic models (Sulis et al., 2010, Smith et al., 2012, Maxwell et al.,2014, Kollet et al., 2017). Results illustrate that all models can repro­
duce general hydrologic behaviors of benchmark cases (Liggett et al.,2012, Maxwell et al., 2014, Kollet et al., 2017, Wen et al., 2021).
However, fully integrated hydrologic models are computationallyintensive, require extensive parameterization, and need substantial
execution time. These limitations restrict their applicability at large
scales.
The surface–subsurface coupling strategy determines the volume andrate of exchange of water between the two domains (Huang and Yeh,2009), alters the rainfall-runoff behavior, such as peak flow and theonset of water ponding (Kollet et al., 2017), and eventually affects theaccuracy, stability, and robustness of a hydrologic model (Liggett et al.,2012, Coon et al., 2020) because of the faster dynamics of surface flow
compared to subsurface flow. Coupling strategies are broadly classified
into three categories: continuous pressure, first-order exchange (i.e., theexchange flux is determined by a calibrated coefficient and the headdifference between the surface and subsurface nodes), and hybrid ap­
proaches (e.g., boundary condition (BC) switching, and regularization/
penalty methods) (Huang and Yeh, 2009, Coon et al., 2020).
The continuous pressure approach enforces pressure continuity inthe two domains (Kollet and Maxwell, 2006), but rapid changes insurface pressure can cause numerical instabilities at the subsurface
boundary (Huang and Yeh, 2009, Liggett et al., 2012, De Maet et al.,2015). The convergence is especially poor near the transition from dry
surface to surface ponding (Coon et al., 2020). In addition, the same time
step must be used for solving the surface and subsurface equations
simultaneously, which greatly increases the associated computational
cost because smaller time steps are required for surface water flow
(Huang and Yeh, 2009, De Maet et al., 2015).
In contrast, a conceptual interface is considered to represent theconnectivity between the surface and subsurface in the first-order ex­change approach (i.e., Ebel et al., 2009, Liggett et al., 2012). This
coupling approach allows surface and subsurface equations to be solved
separately (Panday and Huyakorn, 2004), and the exchange water flux
through the interface layer is obtained as the product of the pressure
difference in the two domains and an exchange coefficient (Qu andDuffy, 2007, Liggett et al., 2012). A properly calibrated exchange co­efficient can enhance the solver convergence and produce similar results
to the pressure continuity approach (Liggett et al., 2012).
The BC switching approach calculates infiltration rates for stages
with and without surface ponding (Camporese et al., 2010, Camporese
et al., 2014). The subsurface flow equation is first solved for theatmosphere-controlled BC (i.e., specified time-dependent water fluxes),
and then a mass balance calculation is performed to determine theponding condition at each surface node. Once surface ponding occurs,
the land surface BC is switched from an atmospheric-controlled BC to asoil-controlled BC (i.e., specified time-dependent pressure) (Sulis et al.,2010). The BC switching approach causes a time lag between surface
and subsurface flows on an overland plane, influencing the onset of
runoff generation and the movement of the saturation front when inflow
occurs uphill (Sulis et al., 2010).
De Maet et al. (2015) introduced a hybrid coupling approach thatcombines the continuous pressure and first-order exchange approaches.
This hybrid coupling approach can produce similar results to the firstorder exchange approaches and improve the scheme’s robustness.
Although current coupled/integrated hydrologic models reproduce
similar hydrographs (Huang and Yeh, 2009, Kollet et al., 2017, Coonet al., 2020), the impact of coupling approaches has not been extensively
examined under the complex boundary conditions that occur duringmultiple rainstorms.
Furthermore, the application of existing coupled models to largescale problems is limited due to simplified representation of the sub­surface and empirical formulations for infiltration estimation (Giráldez
and Woolhiser, 1996, Singh and Bhallamudi, 1998) or high computa­
tional burden (Morita and Yen, 2002). The computational burden is asignificant obstacle that limits the implementation of hydrologicalmodels at larger scales. The numerical solution of the Richards equation
usually requires fine spatial grids and small time steps, which greatly
increase the computation cost (Niswonger et al., 2006; Šimůnek et al.,2006). Small time steps are required to overcome numerical instabilities
at the subsurface boundary due to rapid changes in surface pressure(Huang and Yeh, 2009, Liggett et al., 2012). Much effort has been
devoted to improving the computational efficiency of integrated hy­
drologic models, including parallel computing (Kollet and Maxwell,2006, Brunner and Simmons, 2012), automatically adjusting time step
size (Panday and Huyakorn, 2004, Simunek et al., 2005, Sulis et al.,2010, Goodrich et al., 2012), nested discretization (Wang et al., 2018),
and dimensionality reduction (Niswonger et al., 2006, Seo et al., 2007,
Brunetti et al., 2018). Different (for overland and subsurface flow) andautomatically adjusted time steps can significantly enhance thecomputational efficiency and robustness of the coupled models (Roberts
2003, Simuunek et al., 2005, Goodrich et al., 2012).
Alternatively, dimensionality reduction of the surface and subsurface
flow domain can provide a comparative advantage. As unsaturated flow
is predominantly vertical (Mantoglou 1992), a single or multiple verticalsoil profiles representing one-dimensional (1D) water flow in the vadosezone can be coupled to a 1D overland flow plane. This simplification anddimensionality reduction could reduce computational time compared tovariably saturated flow simulations in 2D or 3D (Panday and Huyakorn,
2004, Niswonger et al., 2006, Seo et al., 2007). Several conjunctive
surface/subsurface models have been developed to simulate the over­land flow on a plane with 1D or 2D overland flow and 1D, 2D, or 3D flowin the vadose zone (Smith and Woolhiser, 1971, Abbott et al., 1986,
Singh and Bhallamudi, 1998, Morita and Yen 2002, Panday and Huya­
korn, 2004, Thompson et al., 2004, He et al., 2008). Large-scale coupledhydrologic models, such as MIKE SHE (Thompson et al., 2004) and
MODFLOW2000-H1D (Seo et al., 2007, Beegum et al., 2018) are ex­amples of this dimensionality reduction approach.
The general objective of this study is to develop a simple, accurate,and computationally efficient coupled hydrologic model for overland
flow and variably saturated flow through the vadose zone at the hill­slope scale. Two existing open-source 1D codes are coupled for thispurpose: KINEROS2 (referred to below as K2) (Smith et al., 1999,
Goodrich et al., 2012) for overland flow and HYDRUS-1D (referred to
below as H1D) (Simunek et al., 2016) for vertical water flow in thevadose zone. The specific objectives are (i) to develop a novel coupled
H1D-K2 model to reduce the dimensionality of the problem and improvethe computational efficiency; (ii) to evaluate the robustness of the pro­
posed model by comparing its results with the two-dimensional versionof HYDRUS (2D/3D) (referred to below as H2D) on six benchmark
problems; and (iii) to investigate the roles of a rainstorm size and fre­quency, slope, bottom BC, and soil heterogeneity on the surface andsubsurface water flow. We pose the following hypotheses: (i) the
coupled H1D-K2 model accurately represents the surface–subsurface
water flow in response to multiple rainfall storms compared to H2D; (ii)computational time can be reduced by applying one H1D soil profile toan overland flow plane rather than to each individual surface node; and(iii) subsurface heterogeneity and bottom BC affect the runoff2

L. Chen et al.

Journal of Hydrology 614 (2022) 128539

generation mechanism, pattern, and magnitude, while the slope
magnitude impacts the timing of the hydrograph and volume of cumu­
lative runoff.

widely used and open-source watershed model representing a watershedby a cascade of overland flow planes and channels (Smith et al., 1999,
Goodrich et al., 2012). The overland flow is described using the kine­
matic wave approximation of the Saint-Venant equations (assuming
uniform flow and a friction slope approximately equal to the planeslope) (Woolhiser et al., 1970, Giráldez and Woolhiser, 1996). The ki­
nematic wave solution is accurate for most overland flow conditions(Woolhiser and Liggett, 1967, Goodrich et al., 2012). The kinematicwave equation can be expressed as follows:

2. Hydrologic models and setups
To achieve these objectives, a novel approach was developed toexternally couple two well-established and widely used open-source
codes to simulate hydrologic processes at the surface and in the sub­surface using a MATLAB (R2021b) script. The K2 model simulates 1D
overland flow via the kinematic wave equation (Goodrich et al., 2006),
and H1D simulates 1D variably saturated flow in the subsurface(Simunek et al., 2005, Simunek et al., 2016). To evaluate the model
performance, we compared the simulation results of the H1D-K2 modelwith the H2D model that includes an overland flow routine (Šimunek,
2003, Simunek, 2015). The H2D model outputs are used as a reference
for model comparison because the model is physically based (solves theRichards and kinematic wave equations simultaneously) and does not
reduce the problem’s dimensionality (contrary to the coupled H1D-K2
model). The conceptual models for H2D and coupled H1D-K2 are pre­sented in Fig. 1c and 1d, respectively.

∂h ∂Q+= q(x, t)
∂t ∂x

(1)

where h is the surface water depth [m], t is time [s], x is the distancecoordinate over the plane [m], Q is the discharge per unit width [m2/s],
and q(x,t) is the change of the lateral inflow rate [m/s], which is equal tolocal precipitation rate minus local infiltration and evapotranspiration
rates.
Discharge per unit width (Q) is related to the surface water depth (h)by a simple power relation (Giráldez and Woolhiser, 1996):
Q = αhm(2)where parameters α [m1/3s− 1] and m [-] are related to the slope, surface
roughness, and flow conditions (laminar or turbulent flow). The
Manning hydraulic resistance law is commonly used to obtain the

2.1. Coupled H1D-K2 model
The KINEmatic runoff and EROSion model (KINEROS2, K2) is a

Fig. 1. (a) The conceptual model of surface/subsurface hydrological processes; (b) a cascade of 1D soil columns connected with an overland flow plane, where eachsoil column is linked to one surface node (j = 1, 2, …, nk); (c) a conceptual model of a hillslope simulated by the H2D model; and (d) a conceptual model of a hillslope
simulated by the coupled H1D-K2 model, where one or multiple H1D columns are linked to multiple surface nodes (one vertical soil column for a homogeneous
system and three columns for a heterogeneous system in this study). Note: the red circle represents the discretization nodes for solving the surface flow equation. The
red and yellow lines represent the applied boundary conditions (BCs) at the land surface and the bottom of a soil profile, respectively. No flow BCs were applied on
the vertical sides of the 2D model. Different color zones represent different soil types at the soil surface (0 ~ 0.05 m), while the subsurface is uniform. Ks1, Ks2, and Ks3
are the near-surface hydraulic conductivities of the uphill (purple zone), middle (red zone), and downhill (blue zone) sections, respectively. (For interpretation of thereferences to color in this figure legend, the reader is referred to the web version of this article.)3

L. Chen et al.

Journal of Hydrology 614 (2022) 128539

parameter α:knα = S1/2 andm =5
3on the runoff generation is relatively small for the considered timescales. The total simulation time is divided into multiple periods (Δtp),
which are usually larger than the individual time steps (Δti) for solving
the surface and subsurface equations. Furthermore, convergence con­ditions are employed to automatically adjust time steps in both the H1D(Simunek et al., 2005) and K2 (Goodrich et al., 2012) models. To
improve computational efficiency, one soil profile is used to represent
the entire subsurface domain when soil characteristics do not varyhorizontally but may change in the vertical direction. When soil prop­
erties vary horizontally, the overland flow plane is divided into multiple
sections, each connecting to one soil type (Fig. 1d).
A MATLAB script executes the H1D and K2 codes sequentially. A
series of MATLAB functions read and write the input and output files,
initialize local and global variables, and update exchange information at
a selected Δtp. Fig. 2 shows a flowchart depicting the coupling approachin the H1D-K2 model. This coupling sequence is summarized below:(1) The H1D code is executed with the applied land surface BCs (i.e.,
rainfall, evaporation) to obtain the actual infiltration and evaporation
rates. The land surface BC of each H1D soil profile is determined using
the surface water depth calculated in the previous Δtp. If there is no
ponding, an atmospheric BC is used. When water starts ponding at thesoil surface, the land surface BC is switched from an atmospheric BC to a
constant pressure head BC. The corresponding pressure head value isequal to the average water depth of surface nodes in an overland flow
plane. If multiple H1D soil profiles are set, a loop is employed to estimate
infiltration for each soil profile (Fig. 2).
(2) A water balance considering infiltration, rainfall, evaporation,and ponded water is performed to determine if the soil infiltration ca­
pacity (qsoil) exceeds the water supply (qsupply, the sum of rainfall, inflow
from a neighboring plane, and ponded water minus evaporation). Whenthe pressure head BC is imposed, the infiltration volume over a Δtp is(3)

where S is the slope [–], and n is a Manning’s roughness coefficient foroverland flow [m1/3s− 1], k is a unit conversion factor between Inter­national System (SI) and English units. k = 1 (m1/3 s− 1) for SI units and1.49 (ft1/3 s− 1) for English units.
The K2 model has been applied to investigate the processes ofinterception, infiltration, runoff generation, erosion, and sedimenttransport (up to five particle sizes) for individual rainfall-runoff events(Giráldez and Woolhiser, 1996, Goodrich et al., 2006, Guber et al., 2014,
Meles Bitew et al., 2020). Infiltration is described using an equation for
cumulative vertical infiltration into a homogeneous (infinite) soil profile
with uniform initial moisture content (Parlange et al., 1982). The K2
model can also consider a two-layer soil profile without a bottom
boundary condition to represent differences in soil texture in the sub­surface. K2 cannot simulate saturation excess overland flow (Fig. S5)
because it assumes an infinitely deep soil profile. Furthermore, the K2
model does not provide information on the spatial and temporal vari­ability of water content and neglects evaporation. These limitations
restrict the ability of the K2 model to investigate surface–subsurface
connectivity.
To overcome the above-discussed limitations, we have developed anovel coupling approach to replace the empirical infiltration equation in
the K2 model with the physically-based H1D model. H1D simulates 1Dvariably saturated water flow in the vadose zone by numerically solving
the Richards equation (Fig. 1d) and using realistic BCs at the soil profile
bottom (Simunek et al., 2005). Evapotranspiration processes areaccounted in H1D but are not considered in this paper since their impact

Fig. 2. Flowchart depicting the data flow in the coupled H1D-K2 code. Orange color indicates computations in MATLAB. (For interpretation of the references to color
in this figure legend, the reader is referred to the web version of this article.)4

L. Chen et al.

Journal of Hydrology 614 (2022) 128539

controlled by qsoil rather than by water supply, and calculated cumula­
tive infiltration may be larger than qsupply. The pressure head BC will
overestimate infiltration when the soil infiltration capacity exceeds thewater supply. The pressure head BC for the subsurface equation istherefore switched to an atmospheric BC if this exceedance occurs. The
H1D code will maintain water balance when using an atmospheric BC
that considers the actual infiltration rate and the maximum water supply′
volume during Δtp. The potential surface flux (qsurf ) is set equal torainfall rate (qr) plus the ponding water contribution given as:∑nk′
j=2 h(j, t)
q surf (t) = qr (t) +(nk − 1)Δtp

domain is represented using the source term in the surface flow equation(Eq. (1)), while the applied land surface BCs account for the effects ofsurface water conditions on subsurface flow. The soil-hydraulic func­
tions of van Genuchten (1980) are used to obtain unsaturated soil hy­
draulic properties.
The atmospheric BC with surface runoff is applied on the land surface
to reflect the surface conditions, and BC switching between the variable
pressure head and the atmospheric-controlled conditions (i.e., precipi­
tation and evaporation) are implemented to represent ponding or no
ponding conditions, respectively. At any given time step, the subsurface
equation is solved first for prescribed initial conditions and an atmo­
spheric BC. After that, a water balance calculation is performed to
calculate the exchange flux and distribute it to separate surface andsubsurface nodes. Once surface ponding is reached at any surface node
(once qsoil is smaller than ql(j,t)), overland flow is generated. The surface
water is routed to the lower parts of the hillslope via a kinematic wave
equation. When water is ponded at the soil surface, the land surface BC isswitched to a variable pressure head BC. Surface water depths at eachsurface node are updated at each time step. The output information isprinted for each node where the overland flow becomes active. Details
concerning the theory, numerical solution techniques, and example
applications of H2D are given by Šimunek and van Genuchten (2008)
and Kohne et al. (2011).

(4)

where h(j,t) is the time-series of surface water depths at the surface node,and j is the surface node number starting from node 1 and ending with
node nk. Note that the second term on the right side of Eq. (4) gives themean surface water depth for all surface nodes (j = 1, 2, …, nk) divided
by the Δtp. The subsurface equation is then solved with the surface water
depths at the previous time as an initial condition and the new specified
flux BC. The criterion used to stop each iteration cycle is the difference
between cumulative infiltration and cumulative water supply over thetime interval being less than 10-6 m2.
(3) The time-series of net surface flux ql(j,t) (i.e., the difference be­
tween rainfall, evaporation, and actual infiltration) generated by the
H1D is then applied to the corresponding surface nodes in the K2 code.
When rainfall is larger than actual infiltration, positive ql is redistributedat the land surface and routed down over the overland flow plane. In thecase of negative ql, ponded water infiltrates into the subsurface, and its
absolute value is subtracted from the surface storage of each surface
node while solving the surface water routing equation. Note that ql(j,t) inan H1D profile corresponds to vertical inflow from multiple surface
nodes in the K2 calculation (Fig. 1d). Therefore, a new method is used toassign the time series of ql(j,t) in the overland plane. A weighting factor
is calculated as a ratio of the surface water depth at each surface nodeand the average surface water depth of an entire plane at the end of theprevious Δtp. It is reasonable to assume that cumulative infiltration into
the soil with a higher ponding depth is larger than the one with a smaller
ponding depth. When the surface water depth at a surface node is zero,ql(j,t) is zero as well. Additionally, the K2 code was modified to save andinitialize local and global variables with the previously calculated re­sults so that ql(j,t) is applied to the corresponding surface nodes. The net
surface flux is obtained by:⎧⎪⎨ql (j, t) =3. Benchmark cases
The hydrological processes on a hillslope with a length (L) of 400 m
and subsurface thickness of 1 m (Fig. 1a) were simulated and comparedusing the H2D and coupled H1D-K2 models. The subsurface was dis­cretized into 100 finite elements with the mesh size varying from0.0053 m near the soil surface to 0.012 m at the soil bottom. The hor­izontal mesh size was 2 m for the H2D model and 10 m for the K2 model.
The total number of surface nodes was 41 in the K2 model. The firstsurface node was used to receive inflow from uphill (not considered inthis study). Finer horizontal discretization was used in the H2D model
because coarse (>2 m) spatial discretization causes numerical oscilla­tion problems for the heterogenous hillslope with a steep slope. The
subsurface spatial heterogeneity in the horizontal direction determines
the number of H1D profiles. At least three H1D soil profiles were needed
to reflect the spatial heterogeneity in Fig. 3b and 3c. For the model withthree soil profiles, the uphill, middle, and downhill H1D soil profiles
correspond to 15, 10, and 15 surface nodes, respectively. For the coupled
H1D-K2 model with 40 soil profiles, each surface node relates to one
H1D soil profile. Furthermore, an analytical solution of overland flow on
an impermeable plane served as a base case to evaluate the model
performance. The subsurface with the saturated hydraulic conductivityof 10-11 m/min was assumed to be impermeable, and modeled results
were compared with the analytical solution.
The performance of the coupled H1D-K2 code was tested against
H2D for homogeneous and heterogeneous conditions with horizontallydifferent saturated hydraulic conductivities (Ks), initial water content
distributions, bottom BCs (i.e., no flow and free drainage), slope ratios,
and rainfall intensities and durations. Detailed model setups are listed inTable 1. Four slopes (0.0005, 0.001, 0.0025, and 0.01) were selected totest the impact of the surface flow speed on infiltration and overland
flow, while two bottom BC types were employed to test the feedback
between the subsurface and surface flow processes. The bottom BC isused to represent the exchange flux between the vadose zone andgroundwater. The no-flow BC represents an impermeable bottom layer,
while the free drainage BC allows for water to drain from the bottom of
the soil profile by gravity (Simunek et al., 2005).
Time step and iteration criteria settings in the H1D and H2D models
were the same to obtain comparable results. The pressure head andwater content tolerances were equal to 0.01 m and 0.001, respectively.
The initial pressure head distribution was assumed to be hydrostatic
equilibrium with a given top pressure head. The same initial pressureqr (t) − qi (t)if qr (t)⩾qi (t)

⎪
⎩ ∑nk

h(j, t)h(j, t)/(nk − 1)
j=2

[qr (t) − qi (t) ]if qr (t) < qi (t)(5)

where qi is the infiltration rate.
(4) The K2 code is executed to calculate lateral water flow on theland surface while considering the feedback of the subsurface hydro­
logical processes. Surface water depths at each surface node and theoutflow rate from the overland plane are simulated during Δtp.
(5) Calculated pressure heads are set as the initial condition for the
H1D computation at the next Δtp. The information on water depths is
updated for the K2. After that, a temporal loop is executed until the total
simulation time (ti) reaches the final time (tmax).
2.2. HYDRUS-2D (H2D) model with an overland flow routine
The H2D (version 3) model with an overland flow routine was used
to simulate runoff at the land surface and water infiltration and redis­
tribution processes in the vadose zone (Šimunek, 2003, Simunek, 2015).
Water flow in the vadose zone is described using the Richards equation,
while the overland flow is described using the kinematic wave equation(Eq. (1)). The water exchange rate between the surface and subsurface
5

L. Chen et al.

Journal of Hydrology 614 (2022) 128539

Fig. 3. The conceptual model simulated by the coupled H1D-K2 model for (a) a homogeneous case, (b) heterogeneous case 1, and (c) heterogeneous case 2. The
subsurface of an overland plane is uniform, except for the blue zones with a Ks of 6.94 × 10-6 m/min. The blue points h1, h2, and h3 are observation nodes. nh1d is thenumber of the H1D soil profiles. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)
Table 1
Mesh size, initial and boundary conditions, and parameters for the benchmark cases.
Homogeneous subsurface
Mesh size
HYDRUS-2D
Horizontal mesh size, Δx (m)Coupled H1D-K2 model
Horizontal mesh size, Δx (m)Exchange period, Δtp (min)Number of H1D soil profiles
Initial and boundary conditions
Depth of initial water table, IWT (m)Rainfall rate, qr(×10-4 m/min)
Bottom BCs
Slope, (-)
Manning’s roughness coefficient, (×10-4)
Saturated hydraulic conductivity, Ks (×10-5m/min)van Genuchten parametersθsat (-)
θres (-)
n0 (-)
α0 (1/m)Heterogeneous subsurface

Base case

Excess saturation

Excess infiltration

Multiple rainstorms

Case 1Case 2

2.0

2.0

2.0

2.0

2.0

2.0

10.0
10, 100*
1

10.0
10, 100*
1

10.0
10, 100*1

10.0
2**, 10
1

10.0
10
3

10.0
10
3

0.0
3.3

0.5, 0.75, 1.0
3.3

1.0
3.3

1.0
0.33, 3.3

1.0
3.3

1.0
3.3

no flow
0.0005
3.31
10-6

no flow
0.0005
3.31
69.4

no flow
0.0005
3.31
69.4, 6.94, 0.694

no flow, free drainage
0.0005
3.31
6.94

free drainage
0.01
3.31
0.694, 69.4

no flow
0.00053.31
0.694, 69.4

0.4
0.08
2.0
1.0

0.4
0.08
2.0
1.0

0.4
0.08
2.0
1.0

0.4
0.08
2.0
1.0

0.4
0.08
2.0
1.0

0.4
0.08
2.0
1.0

Notes: Symbol “*” indicates that the exchange period of 100 mins was used to simulate a hillslope with a no-flow BC. Symbol “**” indicates that the exchange period of
2 mins was used to simulate a hillslope with a no-flow BC and a slope of 0.01. The initial water table depth (IWT) is defined as the water table depth relative to theground surface at the beginning of the calculation. θsat and θres are the saturated and residual water contents, respectively. α0 is a parameter corresponding
approximately to the inverse of the air-entry value, cm− 1, and n0 is an empirical shape-defining parameter.

head distribution was used for both the H1D and H2D models. Bothinitial and minimum time steps were set at 10-6 min, and the maximum
time step was set to 1 min. The maximum number of iterations was equal
to 10. All simulations were carried out for 300 mins, except for caseswith multiple rainstorms (1200 mins). The rainfall rate was applied for
200 mins, followed by 100 mins of recession (qr = 0). For the coupled
H1D-K2 model, the duration of water exchange (Δtp) between surface
and subsurface varied from 2 to 100 mins, depending on the slope andheterogeneity of the overland plane. Simulation outputs were reported
at a one-min interval. The van Genuchten parameters are listed inTable 1.
The performance of the coupled H1D-K2 code was evaluated bycomparing its results with the output of the H2D model, including thetime-series of outflow rates, surface water depths at three observation
points (h1, h2, and h3, shown in Fig. 3), and the surface water storage.
Furthermore, the relative mass balance error (MBE) evaluated the errors
caused by the coupling strategy and numerical methods used in the twodomains. The relative MBE was computed for the entire hillslope
domain using the following equation:⃒⃒∫ t⃒ max [q (t) − q (t)− q (t)]dt⋅L− ∫ tmax q (t)dt − Δq − Δq ⃒
r
e
d
out
surf
sub ⃒
⃒
0
∫ tmaxRelativeMBE = ⃒ 0⃒⃒⃒
qr (t)dt
0(6)where Δqsurf (m2) and Δqsub (m2) are the changes in the surface andsubsurface storages over the total simulation time, respectively, qe is theevaporation rate (m/min), qd is the drainage at the soil profile bottom(m/min), and qout is the outflow rate (m2/min).
The ability of the coupled H1D-K2 model to predict spatial andtemporal water content distributions depends on the number of verticalsoil profiles employed. We, therefore, simulated and compared threemodels: H2D (served as a reference), the coupled H1D-K2 model withone (for a homogeneous case) or three (for a heterogeneous case) soil
profiles, and the coupled H1D-K2 model with 40 soil profiles. Kriging
interpolation with a linear semi-variogram model was used to estimate
the water content distribution. Water content variations with depth andtime were also compared to evaluate vadose zone flow (Figs. S1 and S2).
To illustrate the accuracy and robustness of the proposed H1D-K2 modelover a wide range of conditions, we simulated and compared outflow

6L. Chen et al.

Journal of Hydrology 614 (2022) 128539

rates and cumulative runoffs for 32 cases with those obtained by theH2D model. The detailed descriptions of the sensitivity analysis, therelated results, and the discussion are given in the Supporting Material(Fig. S3). In general, simulation results show that the coupled H1D-K2
model can accurately and efficiently simulate the surface and subsur­
face hydrological processes on a hillslope under scenarios with different
spatial heterogeneities (Fig. 3), slopes (i.e., 0.0005, 0.001, 0.0025, and0.01), and bottom BCs (i.e., no flow and free drainage). The overland
flow plane connecting one (homogeneous cases) and three (heteroge­
neous cases) H1D profiles with different soil layering can capture thegeneral hydrological behavior at the surface and in the subsurface(Fig. S3), which is similar to the results of integrated hydrologic models
used in previous studies (e.g., Simunek 2003, Sulis et al., 2010, Kollet
et al., 2017, Coon et al., 2020). Under a uniform rainfall storm, therunoff pattern was sensitive to the spatial heterogeneity, while its
magnitude was related to the drainage condition at the bottom (Fig. S3).

The length of the middle section was 100 m, while it was 150 m forboth uphill and downhill sections. A rainfall rate (qr = 3.3 × 10-4 m/min) was applied for the first 200 mins, followed by the 200 mins ofrecession (qr = 0 m/min). The slopes were 0.01 and 0.0005 for case 1
and case 2, respectively. Free drainage and no flow BCs were applied at
the bottom of the soil profiles in case 1 and case 2, respectively. These
two benchmark cases have been used in the intercomparison of inte­
grated hydrologic models (e.g., Sulis et al., 2010, Simunek 2015, Kollet
et al., 2017, Coon et al., 2020).
3.3. Runoff generation on a homogeneous hillslope with multiple
rainstorms
Four different configurations (slopes of 0.0005 and 0.01, and no flow
and free drainage bottom BCs) with multiple rainstorms were consid­
ered to examine the stability of the coupled H1D-K2 model. The rainfall
rate (qr = 3.3 × 10-4 m/min) was applied for the first 200 mins, thendecreased to 1.0 × 10-4 m/min for the next 200 mins, followed by 200
mins of recession (qr = 0 m/min). The rainfall events were repeated in
the next 600 mins. The total simulation time was 1200 mins. The soil
was assumed to be homogeneous and isotropic with a Ks of 6.94 × 10-5
m/min. The H1D-K2 model only employed one soil profile to represent
vadose zone flow in these simulations (Table 1).

3.1. Runoff generation on a homogeneous hillslope with uniform rainfall
Two homogeneous and isotropic planes with impervious andpervious soils and a slope of 0.0005 were used to evaluate the model
performance in predicting surface runoff in response to uniform rainfall.
For the coupled H1D-K2 model, one soil profile was used to simulate
water flow in the subsurface (Table 1). The rainfall rate (qr = 3.3 × 10-4
m/min) was applied for the first 200 mins, followed by 100 mins ofrecession (qr = 0 m/min). The total simulation time was 300 mins
(Table 1).
An analytical solution of the kinematic wave equation on animpervious rectangular plane (Gottardi and Venutelli, 1993) was used tocompare simulated outflow rates using the coupled H1D-K2 and H2Dmodels. The parameters used in the analytical solution (base case) arelisted in Table 1. The detailed information, including the mathematicalequation, parameters, and BCs can be found in Gottardi and Venutelli(1993) and Jaber and Mohtar (2003).
The mechanism for runoff generation is mainly classified into two
types: saturation excess (Dunne) runoff and infiltration excess (Horton)
runoff (Smith and Woolhiser, 1971, Fatichi et al., 2016). For the Dunnerunoff, surface ponding occurs after the soil saturates and water fills the
subsurface storage. To realize this runoff generation mechanism, thesaturated hydraulic conductivity (Ks = 6.94 × 10-4 m/min) was specified
to be larger than the rainfall rate (qr = 3.3 × 10-4 m/min). Three sim­
ulations with different depths of the initial water table (IWT) were tested(Table 1). Conversely, for the Hortonian overland flow, the surface
ponding occurs when the rainfall rate exceeds the infiltration rate (Ks <qr). Three simulations with the same IWT depth of 1 m and different Ks(6.94 × 10-5 m/min, 6.94 × 10-6 m/min, and 6.94 × 10-11 m/min) wereconsidered.

4. Results
4.1. Surface runoff over a homogeneous subsurface with a uniform
rainfall rate
Simulation results for the H2D and coupled H1D-K2 models haveexcellent agreement with the analytical solution for the lowpermeability subsurface (6.94 × 10-11 m/min) (Fig. 4a) and saturated
condition (Fig. 4b). The average errors were less than 0.5 % in bothcases. Simulated outflow rates for the H2D and H1D-K2 models were ingood agreement for both infiltration excess (Fig. 4a) and saturation
excess (Fig. 4b) conditions. Differences were less than 0.9 % in outflow
rates and less than 0.5 % in the cumulative outflow.
4.2. Surface runoff over a heterogeneous hillslope with a uniform rainfall
rate
4.2.1. Heterogeneous case 1
The hydrograph simulated by the H2D model (Fig. 5a) shows that
outflow began at 3.8 mins, experienced a steep increase to 6.52 × 10-4
m2/s, remained constant until the arrival of the second wave at 145
mins, and then gradually increased and reached peak flow (6.99 × 10-4m2/s). After rainfall stopped (t = 200 mins), the outflow rate decreased
gradually to zero at 229 mins. Similar outflow rates were obtained forcase 1 by the coupled H1D-K2 model when a high-permeability zone wasplaced near the surface in the middle section of the hillslope. The onset
and cessation of surface runoff were almost identical for the three hy­
drological models.
Simulated surface water depths at observation points h1 and h3 were
almost identical for all three models (Fig. 5b). The biggest difference
occurred in predicting the onset of water ponding at h2 node located on
the high permeability section of the hillslope. The corresponding values
for the H2D model, the coupled H1D-K2 model with 40 soil profiles, andthe coupled H1D-K2 model with three soil profiles were 92 mins, 82
mins, and 122 mins, respectively. The maximum differences in thesurface water depth were less than 0.2 mm among the three models.
Table 2 shows that differences in simulated water balance compo­
nents of heterogeneous case 1 were very small among the three models.
Compared to the results produced by H2D, peak flow predicted by the
coupled H1D-K2 model with 40 and three soil profiles was approxi­
mately 8.0 % higher, while the cumulative outflow was quite similar(2.2–3.4 % difference). The shortest computation time was for the

3.2. Runoff generation on a heterogenous hillslope
A hillslope transect with two soil materials was simulated to consider
spatial heterogeneity of subsurface hydraulic properties on the runoff
generation. To demonstrate the robustness of the coupled H1D-K2
model, we considered two heterogeneous cases in which high or low
permeability zones are near the surface (0–0.05 m) compared to the restof the soil materials. Heterogeneous case 1, similar to a test problem
introduced in Šimunek (2003), uses a high permeable near-surface(0.05 m) middle section (Fig. 1d and 3) with a Ks value of (6.94 × 104
m/min) that is two orders of magnitude higher than the rest of theuniform subsurface with a Ks of 6.94 × 10-6 m/min. Heterogeneous case
2 is similar to Kollet and Maxwell (2006) with a low permeable nearsurface (0.05 m) soil in the middle of the hillslope, and a Ks value twoorders of magnitude smaller (6.94 × 10-6 m/min) than the remaining
soil (Ks = 6.94 × 10-4 m/min). Both heterogeneous cases employed three
H1D soil profiles to describe subsurface water flow in the hillslope’s
uphill, middle, and downhill sections (see Fig. 1b and S1).
7

L. Chen et al.

Journal of Hydrology 614 (2022) 128539

Fig. 4. Outflow rates for (a) infiltration excess and (b) saturation excess, simulated by the H2D (dashed color lines) and coupled H1D-K2 codes (solid color lines). The
hillslope is homogeneous, with a no-flow BC applied at the bottom. Note: IWT represents the initial water table. The solid black lines represent the base case(analytical solution).

Fig. 5. (a) Outflow rates for heterogeneous case 1 with a slope of 0.01 and the free drainage BC applied at the bottom; and (b) surface water depth at the observation
points of h1, h2, and h3, simulated by the H2D (dash-dotted lines) and coupled H1D-K2 codes with three soil profiles (solid lines) and 40 soil profiles (dashed lines).
Table 2Mass balance components for heterogeneous case 1 with a slope of 0.01 and the free drainage BC applied at the bottom.
Cum.
precip.
(m2)

Cum. outflow
(m2)

Cum.
drain.
(m2)

Surfacestorage(m2)

Subsurface
storage(m2)

Relative MBE(%)

Onset ofoutflow(mins)

H2D

26.40

7.57

28.90

0.00

− 10.04

0.01

3.8

H1D-K2(nh1d* =40)
H1D-K2(nh1d* = 3)

26.40

8.07

30.49

0.00

− 11.97

0.72

3.0

26.40

7.97

29.58

0.00

− 11.090.22

3.0

Peak flow(m2/s)

Comp.
time
(s)

6.99 × 10-

32.49

7.60 × 10-

321.747.60 × 10-

21.36

4
4

4

Note: *nh1d – the number of the H1D soil profiles.

coupled H1D-K2 model with three soil profiles (21.36 s), followed by theH2D model (32.49 s), while the coupled H1D-K2 model with 40 soil
profiles had the longest computation time (321.74 s).
Fig. 6 demonstrates similar subsurface water content distributions
simulated by the three models in response to the rainfall-runoff pro­
cesses in case 1. Since the rainfall intensity exceeds the infiltration ca­
pacity of the lower and upper sections of the hillslope, water ponding
occurs quickly in these two sections and is routed downslope, either off
the hillslope or to the midsection (Fig. 6a). As rainfall intensity does notexceed the infiltration capacity of the midsection, no ponding due torainfall occurs. Instead, the ponded water from the upper part of thehillslope was routed downslope and infiltrated into the subsurface in themidsection causing earlier arrival of the saturation front in the middle
section of the hillslope with higher permeability. After the cessation of
the rainfall, water drained out from the bottom of the soil profile,causing a decrease in the water content. The H2D model and the coupled
H1D-K2 model with 40 soil profiles very well reproduced the overlandand subsurface water movement behavior. Conversely, as expected, only
the horizontally averaged saturation front characteristics were captured
in the middle section by the coupled H1D-K2 model with three soil
8

L. Chen et al.

Journal of Hydrology 614 (2022) 128539

Fig. 6. Simulated water content distributions for heterogeneous case 1 with a free drainage BC at the bottom of the soil profile below the overland flow plane with aslope of 0.01 at times of 0, 60, 90, 110, and 300 mins (top to bottom) simulated by H2D (left), coupled H1D-K2 with 40 soil profiles (middle), and coupled H1D-K2
with three soil profiles (right). The vertical distance is not at scale. nh1d is the number of the H1D soil profiles.

profiles, where this section was represented by only one soil profile.
Although this model significantly simplifies the real conditions, it could
produce the overall hydrological behavior of this entire system (e.g.,
outflow rate and peak, average subsurface saturation) sufficiently wellfor most practical problems.

the hillslope). H2D predicts that the outflow occurs at 108.3 mins, ex­periences a step increase to 1.25 × 10-3 m2/s after several mins, thenremains constant for about 50 mins before gradually reaching the peakflow (2.08 × 10-3 m2/s) and a slow recession once rainfall stops (at 200mins). After cessation of the rainfall (t = 200 mins), the outflow rate
decreased gradually, dropping to 1.59 × 10-4 m2/s at the end of thesimulation. The coupled H1D-K2 model shows very similar outflow
behavior as the H2D model, with the only difference in the early onset ofpredicted surface runoff of 4 and 17 mins for the H1D-K2 models with 40

4.2.2. Heterogeneous case 2
Fig. 7a shows the hydrograph simulated by the three models forheterogeneous case 2 (low-permeability zone located in the middle of9L. Chen et al.

Journal of Hydrology 614 (2022) 128539

Fig. 7. (a) Time-series of the outflow rate and (b) variations of surface water depths at observed points h1, h2, and h3 for a heterogeneous hillslope (case 2) with aslope of 0.0005 and a no-flow BC applied at the bottom, simulated by the H2D (dash-dotted lines), and coupled H1D-K2 code with three soil profiles (solid lines) and40 soil profiles (dashed lines).

and three soil profiles, respectively. All three models predicted the same
trend of surface water depths at points h1 and h2 (Fig. 7b), with the only
minor differences in the onset of runoff at point h3 (Fig. 6b).
Fig. 8a shows the water content distributions simulated by H2D atselected times. Ponded water from the low-permeability zone (themiddle section of the hillslope) was routed toward the outlet and infil­
trated into the downhill portion of the hillslope. The saturation front
gradually moved downward and toward the outlet. Because of theadditional water source from the midsection, full surface saturation wasachieved faster in the downhill section than in the uphill portion of thehillslope. Similar water content distributions were observed for the H2Dmodel and the H1D-K2 model with 40 soil profiles (Fig. 7a and 7b).
Conversely, the coupled H1D-K2 model with three soil profiles could not
fully capture the same water content distribution dynamics. In partic­
ular, the routed surface water infiltrated and redistributed in the verticaldirection and reached saturation after 90 mins in the downhill portion of
the hillslope.
Table 3 summarizes components of the simulated water balance forthe three models for heterogeneous case 2. Differences in cumulative
outflow and peak flow were less than 0.3 % for the coupled H1D-K2
model with 40 and three soil profiles compared to H2D. The simulated
peak flow and time were the same for all three models. Moreover, thecoupled H1D-K2 model with three soil profiles had the shortest
computation time (16.71 s), followed by the H2D model (24.35 s). The
coupled H1D-K2 model with 40 soil profiles had the longest computa­tion time (145.92 s).

the coupled H1D-K2 and H2D models (Fig. 9c). The biggest difference inponding depths between the two models was observed during the second
rainfall event (about 1 mm) because H2D predicts ponding at the lower
part of the hillslope.
4.3.2. Free drainage bottom BC
Fig. 10 illustrates that the coupled H1D-K2 model reproduced similarhydrographs to the H2D model for homogeneous hillslopes with a freedrainage bottom BC subjected to multiple rainstorms. The free drainage
BC causes dryer soil conditions and facilitates surface water infiltrationcompared to a no-flow bottom BC. The free drainage BC resulted indecreased peak amplitudes of surface runoff for all rainfall events(Fig. 10a) compared to the no-flow BC (Fig. 9a). During the low-intensity
rainfall (0 < qr < Ks), the outflow rate decreased gradually. The second
wave of surface runoff (present in the previous example due to satura­tion excess) did not occur because the subsurface storage was notentirely filled due to the continued draining of the subsurface.
Fig. 10b shows the temporal evolution of the surface water depth at
observation points h1, h2, and h3. While the results simulated by the
H1D-K2 model agree closely with the H2D model’s results for all test
cases, there was a small difference in the timing of surface runoff
cessation. The water depth at the three observation points decreased tozero simultaneously in the H1D-K2 model, whereas they gradually andsuccessively dropped to zero in the H2D model.
5. Discussion
The H2D and coupled H1D-K2 models predict similar hydrologicalbehaviors at the land surface and in the subsurface of homogeneous andheterogeneous hillslopes subject to a uniform rainfall (Figs. 4, 5, and 7)
or multiple rainstorms (Figs. 9 and 10). The simulated components of
the water balance are very close to each other (Tables 2, 3, and S1).
Discrepancies in the onset and cessation of surface runoff and spatial andtemporal water content distributions simulated by the two models arecaused mainly by the coupling strategy, the time interval Δtp, and thenumber of the H1D soil profiles (nh1d).
First, the coupled H1D-K2 model uses the BC switching approach to
execute the two codes sequentially. This results in the subsurface flow
equation being solved before the surface flow equation, with subsurface
flow affecting surface flow at the current time level while surface flow
affecting subsurface flow only at the next time level.
Second, the fixed Δtp also contributes to discrepancies because of atime lag in updating BCs. The H1D model calculates the actual infil­
tration rate with applied surface conditions obtained at the end of theprevious period (Fig. 2). Although water inflows and outflows are

4.3. Surface runoff over a homogeneous subsurface with multiple
rainstorms
4.3.1. No flow bottom BC
Fig. 9a shows hydrographs simulated by the H2D and H1D-K2 (witha single soil profile) models in response to multiple rainstorms when theslope is 0.0005 and 0.01, the subsurface is homogeneous, and a no-flowbottom BC is applied. Even though a single soil profile was employed in
the coupled H1D-K2 model, the hydrograph shapes were very similar forboth models. The only discrepancy was in the cumulative runoff for
relatively flat (0.0005) and steep (0.01) slopes was only 0.14 % and 0.22%, respectively (Table S1). There was also a slight difference (less than
5.0 × 10-5 m2/s) in peak outflow during the second rainfall event (t =200–400 mins). Similarly, the simulated water depths by the coupled
H1D-K2 model were in good agreement with those obtained by the H2Dmodel at the three observation points (Fig. 9b and 9c for slopes of 0.0005
and 0.01, respectively). Furthermore, the timing of the predicted onset
and cessation of surface runoff were the same for all surface points using
10

L. Chen et al.

Journal of Hydrology 614 (2022) 128539

Fig. 8. Simulated water content distributions for heterogeneous case 2 with a no-flow bottom BC at times of 0, 60, 90, 110, and 300 mins (top to bottom) simulatedby H2D (left), coupled H1D-K2 with 40 soil profiles (middle), and coupled H1D-K2 with three soil profiles (right). The vertical distance is not at scale. nh1d is thenumber of the H1D soil profiles.

balanced during particular time steps of the simulation, the time-lagged
approach may cause uncertainty in the distribution of net surface fluxes
in surface nodes (Eq. (5)), which eventually affects the accuracy of thesimulated results. Despite adopting the same sequential solution pro­cedure, the H2D model updates BCs in each individual time step, which
is much smaller than Δtp. This is the main reason for discrepancies in thesimulated outflow rate between two models (Fig. 5a and S3). Similar
conclusions about the limitation of this time-lagged approach have been
discussed in the literature (e.g., Panday and Huyakorn, 2004, Huang andYeh, 2009).
Third, the number of H1D soil profiles causes discrepancies in spatial
and temporal water content distributions (Figs. 6 and 8), onset of thesecond outflow wave, and the surface water flow depths (Fig. 5b and 7).
When only one H1D soil profile is used to represent the entire subsur­
face, the volume of infiltrated water contributed by ponding at multiple
surface nodes is redistributed to one soil profile. This process causes
higher water contents at an earlier stage (Fig. S1c and S2c) and delays
soil saturation (Fig. S1e and S2e). Water fluxes at multiple surface nodes
11

L. Chen et al.

Journal of Hydrology 614 (2022) 128539

Table 3Mass balance components simulated by the coupled H1D-K2 and H2D models for heterogeneous case 2 with a slope of 0.01 and a no-flow BC applied at the bottom.

H2D
H1D-K2(nh1d* = 40)
H1D-K2(nh1d* = 3)

Cum. precip.
(m2)

Cum. outflow(m2)

Surface storage(m2)

Subsurface storge(m2)

Relative MBE(%)

Onset of outflow (mins)

Peak flow(m2/s)

Comp. time
(s)

26.40
26.40

12.65
12.69

0.68
0.70

13.35
13.44

0.54
1.64

108.3
103.0

2.08 × 10-3
2.10 × 10-3

24.35
145.92

26.40

12.65

0.70

13.45

1.50

91.0

2.10 × 10-3

16.71

Note: *nh1d – the number of the H1D soil profiles.

Fig. 9. (a) Comparison of outflow rates for the homogeneous hillslopes with a no-flow BC and a slope of 0.0005 and 0.01 simulated by the H2D (dash lines) andcoupled H1D-K2 (solid lines) models; (b) simulated surface water depths for the hillslope with a slope of 0.0005; and (c) simulated surface water depths for thehillslope with a slope of 0.01. h1, h2, and h3 are three observation points located at 320 m (uphill), 260 m (middle), and 0 m (downhill) away from the outlet,
respectively.

are also applied to only one soil profile (Eq. (5)), causing ponding(Fig. 9) or infiltration (Fig. 10), in all nodes.
Furthermore, differences between H1D-K2 and H2D models for ahomogeneous hillslope (Fig. 4) were much smaller than for the hetero­geneous hillslope (Figs. 5 and 7). The main reason is that the direction of
unsaturated water flow in a homogeneous hillslope is predominantly
vertical, and water starts ponding on the entire plane at the same time.
Also, a relatively small surface flow depth (a maximum of 3.0 cm) does
not significantly impact the infiltration rate (Fig. S4), resulting in anearly uniform infiltration rate along the homogeneous hillslope. This
indicates that one H1D soil profile is sufficient to accurately simulate
water flow in the homogeneous hillslopes. Previous studies have re­ported similar conclusions (e.g., Smith and Woolhiser, 1971, Singh andBhallamudi, 1998). On the other hand, the H1D-K2 model with threesoil profiles only broadly approximates the movement of the surfacewater front at a heterogeneous hillslope (Figs. 6 and 8). The H1D-K2
model can fully capture (similarly to H2D) the general behavior ofheterogeneous hillslopes when more vertical H1D soil profiles are
employed (Fig. 6b and 8b). When users do not need accurate spatial
moisture distribution (Figs. 6 and 8), it is better to use the H1D-K2 modeldue to the higher computational efficiency.
The simulation results of the H1D-K2 model also agree closely withthose presented in an integrated hydrologic model intercomparison
project (Maxwell et al., 2014). In this project, seven broadly used inte­
grated hydrologic models were compared, including ATS (Coon et al.,2020), CATHY (Sulis et al., 2010), HydroGeoSphere (Brunner and
12

L. Chen et al.

Journal of Hydrology 614 (2022) 128539

Fig. 10. (a) Comparison of outflow rates for the homogeneous hillslopes with a free drainage BC and a slope of 0.0005 and 0.01 simulated by the HYDRUS-2D (dashlines) and coupled H1D-K2 (solid lines) models; (b) simulated surface water depths for the hillslope with a slope of 0.0005; and (c) simulated surface water depths forthe hillslope with a slope of 0.01. h1, h2, and h3 are three observation points located at 320 m (uphill), 260 m (middle), and 0 m (downhill) away from the outlet,
respectively.

Simmons, 2012), MIKE SHE (Thompson et al., 2004), and ParFlow(Kollet and Maxwell, 2006). For the homogeneous cases (infiltration
excess and saturation excess, Fig. 4), the peak flow and the onset of flow
are almost identical to the results obtained by the seven integrated hy­
drologic models (Maxwell et al., 2014). For heterogeneous case 2 (Fig. 7
and Table 3), the peak flow (mean relative error of 7.1 % with a devi­
ation of 4.8 %) and the onset of flow (mean value of 17.7 % with adeviation of 13.4 %) fell within the range of results presented in Maxwell
et al. (2014). Similar to our findings, these discrepancies have mainly
been attributed to differences in the coupling strategy, Δtp, and thenumber of soil profiles (Huang and Yeh, 2009, Sulis et al., 2010,
Maxwell et al., 2014).
Compared to H2D, the H1D-K2 code can reduce the computationaltime (Tables 2, 3, and S1) because the dimensionality reduction de­creases the number of elements and nodes for calculation. The BC
switching approach allows surface and subsurface flow equations to be
separately solved using dynamic time steps (Fig. 1), which overcomesnumerical instabilities of subsurface solutions caused by fast surface
flow. Our approach calculates infiltration rates for the entire plane when
ponding occurs using the average surface water depth (Eq. (4)) andrearranges the net surface flux to surface nodes based on individual
surface water depths (Eq. (5)). Similarly, Singh and Bhallamudi (1998)
demonstrated up to 10 times speedup when using a uniform infiltration
rate compared to a non-uniform infiltration rate at surface nodes.
The surface flow on a hillslope can be represented by a cascade of 1D

overland flow planes, while the subsurface flow is assumed to occur onlyin the vertical direction. Our water balance and BC switching approachenables the H1D-K2 model to accurately consider the feedback of sub­surface processes (e.g., evaporation, infiltration, drainage) on surface
water flow. Furthermore, by replacing the empirical infiltration equa­tion in the K2 model with a physically based H1D model, the H1D-K2
model can overcome the limitation of predicting the surface runoff
generated by the saturation excess conditions (Fig. S5) and can provide
spatial and temporal water content distributions (Figs. 6, 8, S1, and S2).
The H1D-K2 model can consider multiple soil layers with different soil
textures, whereas the standard K2 model can only consider two-soil
layers, while the lower layer must be infinite. Spatial variations in thesubsoil characteristics in the vertical direction can be accurately
captured by the H1D-K2 model, which presently does not account forlateral water flow in the vadose zone.
Note that MATLAB only serves as a platform to exchange information
between the H1D and K2 codes and assignes exchange fluxes to separate
surface and subsurface nodes based on the surface water depth. When
access to a comprehensive set of MATLAB and Simulink products is
restricted, one could use other programming languages (e.g., R, Python,
C++, etc.) to create the coupling routine with the same structure and thespatial and temporal loops described in section 2.1 and Fig. 2. Our
external coupling approach keeps two codes independent, allowing
users to use their latest versions in the coupled H1D-K2 code. Anotheradvantage of this approach is the ease of parallel computing
13

L. Chen et al.

Journal of Hydrology 614 (2022) 128539

implementation since hydrological processes in an H1D soil profile are
independent of other soil profiles. Furthermore, the users can also run
the coupled H1D-K2 code on high-performance computing cluster sys­
tems after integrating with the MATLAB Parallel Server.
However, the main advantage of the H1D-K2 model compared to K2
is that it can consider complex meteorological conditions (rainfall
variability), different soil textures, different initial water conditions,different bottom BCs (e.g., free drainage, seepage face, variable fluxes/pressures), and the effects of vegetation, which are the dominant factors
for subsurface flow conditions. Moreover, since the physical nonequi­
librium flow model (e.g., dual-porosity and dual-permeability) has beendeveloped in the computational module of H1D (Simunek et al., 2005),
the same conceptualization can be used in the coupled H1D-K2 model tosimulate the impacts of physical nonequilibrium on surface water flow.
Future research will extend this computationally efficient model to awatershed scale by dividing the domain into a cascade of overland flow
planes and channels.

Data availability
Data will be made available on request.
Acknowledgements
This research was supported by the United States Department ofAgriculture, Agricultural Research Service, National Program 211. The
authors thank Prof. Carl L. Unkrich (USDA-ARS, Tucson, AZ) for
providing insightful comments and suggestions for the modification of
the KINEROS2 code.
Open research
The HYDRUS-1D software is publicly available on the PC-Progress
website; its standard computational module is open source (https://www.pc-progress.com/en/Default.aspx?hydrus-1d). The KINEROS2
model is also open-source and available on https://www.tucson.ars.ag.
gov/kineros/. A HYDRUS (2D/3D) license would be required to run the
HYDRUS-2D model. The MATLAB script is available by contacting Dr.
Lin Chen (linc@ucr.edu, chl_sc@163.com).

6. Conclusions
A novel approach was developed to couple open-source codes forvadose zone flow (HYDRUS-1D) and overland flow (KINEROS2),
resulting in the H1D-K2 model. This new coupling routine sequentiallysolved the Richards equation for subsurface flow and the kinematicwave equation for overland flow. A water balance and BC switchingapproach was used to couple the surface and subsurface hydrologicalprocesses. The water exchange flux in an H1D profile was assigned tomultiple surface nodes in the K2 calculation using a weighting factor
related to the water depth.
Six benchmark problems were performed over a wide range of con­ditions with different saturated hydraulic conductivities, initial condi­
tions, bottom BCs, slopes, and multiple rainstorm events. The simulatedhydrographs, surface water depth, and water components by the
coupled H1D-K2 model closely agree with those obtained by the H2Dmodel. The coupled H1D-K2 model is a simple, accurate, and compu­
tationally efficient hydrologic model. It can simulate surface runoff
generated by infiltration excess and saturation excess in response to
rainfall or inflow. It can also provide spatial and temporal water content
distributions and surface water depths. Its computational efficiency isimproved by adopting a dynamic time step approach and applying one
H1D soil profile to an entire overland flow plane rather than to indi­vidual surface nodes. A heterogeneous hillslope can be divided into acascade of planes with different slopes, spatial variations of soil char­acteristics, and different drainage conditions at the bottom. The novel
coupling approach and the H1D-K2 model can be extended to larger
scales by dividing a watershed into a cascade of overland flow planes
and implementing parallel computing to improve computational
efficiency.

Appendix A. Supplementary data
Supplementary data to this article can be found online at https://doi.
org/10.1016/j.jhydrol.2022.128539.
References
Abbott, M.B., Bathurst, J.C., Cunge, J.A., O’Connell, P.E., Rasmussen, J., 1986. An
introduction to the European Hydrological System — Systeme Hydrologique
European, “SHE”, 2: Structure of a physically-based, distributed modelling system.
J. Hydrol. 87 (1), 61–77. https://doi.org/10.1016/0022-1694(86)90115-0.
Beegum, S., Simunek, J., Szymkiewicz, A., Sudheer, K., Nambi, I.M., 2018. Updating thecoupling 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.
Brunetti, G., Simunek, J., Bautista, E., 2018. A hybrid finite volume-finite element modelfor the numerical analysis of furrow irrigation and fertigation. Comput. Electron.
Agric. 150, 312–327. https://doi.org/10.1016/j.compag.2018.05.013.
Brunner, P., Simmons, C.T., 2012. HydroGeoSphere: a fully integrated, physically basedhydrological model. Groundwater 50 (2), 170–176. https://doi.org/10.1111/j.17456584.2011.00882.x.
Camporese, M., Paniconi, C., Putti, M., Orlandini, S., 2010. Surface-subsurface flow
modeling with path-based runoff routing, boundary condition-based coupling, andassimilation of multisource observation data. Water Resour. Res. 46 (W02512)
https://doi.org/10.1029/2008WR007536.
Camporese, M., Daly, E., Dresel, P.E., Webb, J.A., 2014. Simplified modeling ofcatchment-scale evapotranspiration via boundary condition switching. Adv. Water
Resour. 69, 95–105. https://doi.org/10.1016/j.advwatres.2014.04.008.
Coon, E.T., Moulton, J.D., Kikinzon, E., Berndt, M., Manzini, G., Garimella, R., Painter, S.
L., 2020. Coupling surface flow and subsurface flow in complex soil structures using
mimetic finite differences. Adv. Water Resour. 144, 103701 https://doi.org/
10.1016/j.advwatres.2020.103701.
De Maet, T., Cornaton, F., Hanert, E., 2015. A scalable coupled surface–subsurface flow
model. Comput. Fluids 116, 74–87. https://doi.org/10.1016/j.
compfluid.2015.03.028.
Ebel, B.A., Mirus, B.B., Heppner, C.S., VanderKwaak, J.E., Loague, K., 2009. First-orderexchange coefficient coupling for simulating surface water–groundwaterinteractions: parameter sensitivity and consistency with a physics-based approach.
Hydrol. Process. 23 (13), 1949–1959. https://doi.org/10.1002/hyp.7279.
Fan, Y., Clark, M., Lawrence, D.M., Swenson, S., Band, L., Brantley, S.L., Grant, G., 2019.
Hillslope hydrology in global change research and Earth system modeling. Water
Resour. Res. 55 (2), 1737–1772. https://doi.org/10.1029/2018WR023903.
Fatichi, S., Vivoni, E.R., Ogden, F.L., Ivanov, V.Y., Mirus, B., Gochis, D., Tarboton, D.,
2016. An overview of current applications, challenges, and future trends indistributed process-based models in hydrology. J. Hydrol. 537, 45–60. https://doi.
org/10.1016/j.jhydrol.2016.03.026.
Freeze, R.A., 1972. Role of subsurface flow in generating surface runoff: 2. Upstream
source areas. Water Resour. Res. 8 (5), 1272–1283. https://doi.org/10.1029/
WR008i005p01272.
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.
Goodrich, D., Burns, I., Unkrich, C., Semmens, D.J., Guertin, D., Hernandez, M.,
Levick, L.R., 2012. KINEROS2/AGWA: model use, calibration, and validation. Trans.
ASABE 55 (4), 1561–1574.

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.
14

L. Chen et al.

Journal of Hydrology 614 (2022) 128539

Goodrich, D.C., Unkrich, C.L,. Smith, R.E., Woolhiser, D.A., 2006. KINEROS2-newfeatures and capabilities. Proceedings of the Joint Federal Interagency Conference:
April 2-6, 2006, Reno, NV: proceedings of the 3rd Federal Hydrologic Modeling
Conference: proceedings of the 8th Federal Interagency Sedimentation Conference,[Sl]: Subcommittee on Hydrology, 2006.
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.
Guber, A.K., Pachepsky, Y.A., Yakirevich, A.M., Shelton, D.R., Whelan, G., Goodrich, D.
C., Unkrich, C.L., 2014. Modeling runoff and microbial overland transport with
KINEROS2/STWIR model: Accuracy and uncertainty as affected by source of
infiltration parameters. J. Hydrol. 519, 644–655. https://doi.org/10.1016/j.
jhydrol.2014.08.005.
Haddeland, I., Heinke, J., Biemans, H., Eisner, S., Flörke, M., Hanasaki, N.,
Konzmann, M., Ludwig, F., Masaki, Y., Schewe, J., Stacke, T., Tessler, Z.D., Wada, Y.,
Wisser, D., 2014. Global water resources affected by human interventions andclimate change. PNAS 111 (9), 3251–3256.
He, Z., Wu, W., Wang, S.S., 2008. Coupled finite-volume model for 2D surface and 3D
subsurface flows. J. Hydrol. Eng. 13 (9), 835–845. https://doi.org/10.1061/(ASCE)1084-0699(2008)13:9(835).
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.
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.
Köhne, J.M., Wöhling, T., Pot, V., Benoit, P., Leguédois, S., Le Bissonnais, Y., Simunek, J.,
2011. Coupled simulation of surface runoff and soil water flow using multi-objective
parameter estimation. J. Hydrol. 403 (1–2), 141–156. https://doi.org/10.1016/j.
jhydrol.2011.04.001.
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., Sudicky, E.,
2017. The integrated hydrologic model intercomparison 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.
Kourakos, G., Dahlke, H.E., Harter, T., 2019. Increasing groundwater availability andseasonal base flow through agricultural managed aquifer recharge in an irrigatedbasin. Water Resour. Res. 55 (9), 7464–7492. https://doi.org/10.1029/
2018WR024019.
Lawrence, D.M., Fisher, R.A., Koven, C.D., Oleson, K.W., Swenson, S.C., Bonan, G.,
Zeng, X., 2019. The community land model version 5: description of new features,
benchmarking, and impact of forcing uncertainty. J. Adv. Model. Earth Syst. 11 (12),4245–4287. https://doi.org/10.1029/2018MS001583.
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.
Mantoglou, A., 1992. A theoretical approach for modeling unsaturated flow in spatiallyvariable soils: Effective flow models in finite domains and nonstationarity. Water
Resour. Res. 28 (1), 251–267. https://doi.org/10.1029/91WR02232.
Marwaha, N., G. Kourakos, Levintal, E., Dahlke, H.E., 2021. Identifying agricultural
managed aquifer recharge locations to benefit drinking water supply in ruralcommunities. Water Resour. Res. 57(3): e2020WR028811. 10.1029/
2020WR028811.
Maxwell, R.M., Putti, M., Meyerhoff, S., Delfs, J.-O., Ferguson, I.M., Ivanov, V., Sulis, M.,
2014. Surface-subsurface model intercomparison: A first set of benchmark results todiagnose integrated hydrology and feedbacks. Water Resour. Res. 50 (2),1531–1549. https://doi.org/10.1002/2013WR013725.
Meles Bitew, M., Jackson, C.R., Goodrich, D.C., Younger, S.E., Griffiths, N.A., Vaché, K.
B., Rau, B., 2020. Dynamic domain kinematic modelling for predicting interflow
over leaky impeding layers. Hydrol. Process. 34 (13), 2895–2910. https://doi.org/
10.1002/hyp.13778.
Morita, M., Yen, B.C., 2002. Modeling of conjunctive two-dimensional surface-threedimensional subsurface flows. J. Hydraul. Eng. 128 (2), 184–200. https://doi.org/
10.1061/(ASCE)0733-9429(2002)128:2(184).
Niswonger, R.G., Prudic, D.E., Regan, R.S., 2006. Documentation of the unsaturated-zone
flow (UZF1). Package for Modeling Unsaturated Flow Between the Land Surface andthe Water Table With MODFLOW-2005.
Niswonger, R.G., Morway, E.D., Triana, E., Huntington, J.L., 2017. Managed aquifer
recharge through off-season irrigation in agricultural regions. Water Resour. Res. 53(8), 6970–6992. https://doi.org/10.1002/2017WR020458.
Orth, R., Staudinger, M., Seneviratne, S.I., Seibert, J., Zappa, M., 2015. Does model
performance improve with complexity? A case study with three hydrological models.
J. Hydrol. 523, 147–159. https://doi.org/10.1016/j.jhydrol.2015.01.044.
Panday, S., Huyakorn, P.S., 2004. A fully coupled physically-based spatially-distributedmodel for evaluating surface/subsurface flow. Adv. Water Resour. 27 (4), 361–382.
https://doi.org/10.1016/j.advwatres.2004.02.016.
Parlange, J.-Y., Lisle, I., Braddock, R., Smith, R., 1982. The three-parameter infiltration
equation. Soil Sci. 133 (6), 337–341.

Qu, Y., Duffy, C.J., 2007. A semidiscrete finite volume formulation for multiprocesswatershed simulation. Water Resour. Res. 43 (W08419) https://doi.org/10.1029/
2006WR005752.
Report, U. W. W. D. (2018). 2018 UN World Water Development Report, Nature-based
Solutions for Water, www.unwater.org/publications/world-water-developmentreport-2018/. United Nations Educational, Scientific and Cultural Organization, New
York, United States.
Roberts, A., 2003. A holistic finite difference approach models linear dynamics
consistently. Math. Comput. 72 (241), 247–262. https://doi.org/10.1090/S00255718-02-01448-5.
Sasidharan, S., Bradford, S.A., Simunek, J., DeJong, B., Kraemer, S.R., 2018. Evaluating
drywells for stormwater management and enhanced aquifer recharge. Adv. Water
Resour. 116, 167–177. https://doi.org/10.1016/j.advwatres.2018.04.003.
Seo, H.S., Šimůnek, J., Poeter, E.P., 2007. Documentation of the HYDRUS package for
MODFLOW-2000, the US geological survey modular ground-water model. IGWMCInternational Ground Water Modeling Center.
Šimůnek, J., 2015. Implementation of Overland Flow into HYDRUS (2D/3D), HYDRUS
Software Series 6b. Department of Environmental Sciences, University of California
Riverside, Riverside, CA, 44 pp., May 2015.
Simunek, J., van Genuchten, M.T., Šejna, M., 2016. Recent developments andapplications of the HYDRUS computer software packages. Vadose Zone J. 15 (7), 25.
https://doi.org/10.2136/vzj2016.04.0033.
Šimůnek, J., Van Genuchten, M. T., Šejna, M., 2006. The HYDRUS software package forsimulating two- and three-dimensional movement of water, heat, and multiple
solutes in variably-saturated media. Technical manual, version 1: 241.
Šimůnek, J., 2003. HYDRUS-2D Code Modification: Modeling Overland Flow and
Dynamic Interactions Between Plants and Water Flow in a Hillslope Transect. Idaho
National Environmental Laboratory, Final report for project under contract: 55.
Šimunek, J., van Genuchten, M.T., 2008. Modeling nonequilibrium flow and transport
processes using HYDRUS. Vadose Zone J. 7 (2), 782–797. https://doi.org/10.2136/
vzj2007.0074.
Šimunek, J., Van Genuchten, M.T., Sejna, M., 2005. The HYDRUS-1D software package
for simulating the one-dimensional movement of water, heat, and multiple solutes invariably-saturated media. Univ. California-Riverside Res. Rep. 3, 1–240.
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.
Singh, V.P., Woolhiser, D.A., 2002. Mathematical modeling of watershed hydrology.
J. Hydrol. Eng. 7 (4), 270–292. https://doi.org/10.1061/(ASCE)1084-0699(2002)7:4(270).
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/S03418162(99)00033-8.
Smith, R.E., Hebbert, R.H.B., 1983. Mathematical simulation of interdependent surface
and subsurface hydrologic processes. Water Resour. Res. 19 (4), 987–1001. https://doi.org/10.1029/WR019i004p00987.
Smith, M.B., Seo, D.-J., Koren, V.I., Reed, S.M., Zhang, Z., Duan, Q., Cong, S., 2004. The
distributed model intercomparison project (DMIP): motivation and experiment
design. J. Hydrol. 298 (1), 4–26. https://doi.org/10.1016/j.jhydrol.2004.03.040.
Smith, M.B., Koren, V., Reed, S., Zhang, Z., Zhang, Y., Moreda, F., Cosgrove, B.A., 2012.
The distributed model intercomparison project – Phase 2: Motivation and design of
the Oklahoma experiments. J. Hydrol. 418–419, 3–16. https://doi.org/10.1016/j.
jhydrol.2011.08.055.
Smith, R.E., Woolhiser, D.A., 1971. Overland flow on an infiltrating surface. Water
Resour. Res. 7 (4), 899–913. https://doi.org/10.1029/WR007i004p00899.
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.
Thompson, J., Sørenson, H.R., Gavin, H., Refsgaard, A., 2004. Application of the coupled
MIKE SHE/MIKE 11 modelling system to a lowland wet grassland in southeast
England. J. Hydrol. 293 (1–4), 151–179. https://doi.org/10.1016/j.
jhydrol.2004.01.017.
van Genuchten, M.T., 1980. A closed-form equation for predicting the hydraulicconductivity of unsaturated soils. Soil Sci. Soc. Am. J. 44 (5), 892–898. https://doi.
org/10.2136/sssaj1980.03615995004400050002x.
Vorosmarty, C.J., Green, P., Salisbury, J., Lammers, R.B., 2000. Global water resources:
vulnerability from climate change and population growth. Science 289 (5477),
284–288. https://doi.org/10.1126/science.289.5477.284.
Wang, D., Liu, Y., Kumar, M., 2018. Using nested discretization for a detailed yet
computationally efficient simulation of local hydrology in a distributed hydrologic
model. Sci. Rep. 8, 5785. https://doi.org/10.1038/s41598-018-24122-7.
Wen, H., Brantley, S.L., Davis, K.J., Duncan, J.M., Li, L., 2021. The Limits ofhomogenization: what hydrological dynamics can a simple model represent at thecatchment scale? Water Resour. Res. 57 (6) https://doi.org/10.1029/
2020WR029528 e2020WR029528.
Woolhiser, D., Hanson, C., Kuhlman, A., 1970. Overland flow on rangeland watersheds.
J. Hydrol. 9 (2), 336–356.
Woolhiser, D.A., Liggett, J.A., 1967. Unsteady, one-dimensional flow over a plane—therising hydrograph. Water Resour. Res. 3 (3), 753–771. https://doi.org/10.1029/
WR003i003p00753.

15