Journal of Hydrology 640 (2024) 131621

Contents lists available at ScienceDirect

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

Research papers

Computationally efficient Watershed-Scale hydrological Modeling:
Integrating HYDRUS-1D and KINEROS2 for coupled
Surface-Subsurface analysis
Menberu B. Meles a,* , Lin Chen b, Carl Unkrich c , Hoori Ajami b, Scott A. Bradford a ,Jiří Šimůnek b, Dave C. Goodrich caUSDA, ARS, Sustainable Agricultural Water Systems Unit, Davis, CA 95616, USA
Department of Environmental Sciences, University of California, Riverside, CA 92521, USAc
USDA, ARS, Southwest Watershed Research Center, Tucson, AZ 85719, USAb

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 Corradini, Editor-in-Chief.

Watershed flow processes consist of partitioning, movement, storage, and redistribution of water fluxes in space
and time. However, the integrated modeling of these processes is challenging due to computational burden,
extensive data requirements, and/or reliance on simplifying assumptions. This study introduces a novel andcomputationally efficient modeling framework that leverages two state-of-the-art process-based models:HYDRUS-1D (H1D) for unsaturated flow and KINEROS2 (K2) for overland flow. The framework extends ahillslope-scale coupled H1D-K2 model to simulate watershed-scale processes, where H1D replaces the threeparameter Parlange’s infiltration equation in the event-based K2 model. Boundary condition switching is
employed to account for surface ponding and water exchange between the two model domains. The structure of
the coupled watershed-scale H1D-K2 model consists of a cascade of connected rectangular planes, channel elements, and 1D soil profiles to simulate 1D overland flow, infiltration, unsaturated zone flow, and recharge.
Computational efficiency relative to HYDRUS-2D is achieved through a dynamic time-stepping approach anddimensionality reduction. The watershed scale model improved the computational time by 14.3% and 47%
compared to the corresponding Hillslope scale H1D-K2 and HYDRUS-2D, respectively. The performance andefficiency of the new watershed model are demonstrated using benchmark watershed and/or hillslope simulations. Calibrated hydrographs and water balance components using Walnut Gulch Experimental Watershed data,
showed excellent agreement with observed data with a Nash-Sutcliffe coefficient (NSC) greator than 0.8 andPearson correlation coefficient (R2) greater than 0.92.

Keywords:
Coupling hydrologic model
Computational efficiency
HYDRUS-1D / KINEROS2
Soil moisture redistribution
Watershed model
Surface/subsurface interactions

1. Introduction
Process-based models that solve partial differential equations forshallow surface water (de Saint-Venant, 1871) and variably saturated
subsurface flow (Richards, 1931) provide invaluable insights into hydrological processes. However, the numerical solution of these governing equations is computationally demanding because non-linear surface
and subsurface processes need to be solved at high spatial and temporal
resolution using iterative time-stepping methods. Moreover, the inclusion of diverse data sources and inherent heterogeneity in the landscape
properties add computational burden to the preparation of the model
input. Computational demands have, therefore, limited the widespread

application and adaptability of process-based models at watershedscales (Singh and Woolhiser, 2002; Wood et al., 2011; Bierkens et al.,2015). Advances in computing resources and evolving numericalmethods aim to solve increasingly complex problems (e.g., Maxwell
et al., 2014; Clark et al., 2011; Brookfiled et al., 2023).
Computational requirements of watershed-scale models may also
extend beyond the mere representation of complex physical processes
required to understand, quantify, and reduce uncertainty (Liu and Gupta,
2007). Considerable efforts have been directed towards coupled model
development (Shi et al., 2013; Shelia et al., 2018; Sophocleous et al.,
1999; Wagner et al., 2016; Zhang et al., 2013) that can simulate interactions between different hydrological processes, increasing the* Corresponding author.
E-mail address: menberu.meles@usda.gov (M.B. Meles).
https://doi.org/10.1016/j.jhydrol.2024.131621
Received 12 April 2024; Received in revised form 10 June 2024; Accepted 20 June 2024
Available online 3 July 20240022-1694/Published by Elsevier B.V. This is an open access article under the CC BY-NC-ND license (http://creativecommons.org/licenses/by-nc-nd/4.0/).

M.B. Meles et al.

Journal of Hydrology 640 (2024) 131621

functionality and utility of the models by leveraging the capabilities ofindividual models. In addition, watershed-scale models need to consider
optimization strategies for parameter identifiability and sensitivity
analysis (Gupta & Razavi, 2017; Harada and Alba, 2020), adaptation ofadvanced data structures and assimilation techniques (Vrugt et al.,2005; Moradkhani et al. 2005), and model calibration (Razavi and Tolson, 2013; Zhang et al., 2009; Kan et al., 2019). For example, large-scalewatershed models (Shen & Phanikumar, 2010) are commonly used tosustainably manage water resources by studying resource allocation,
reservoir operations, and infrastructure planning (e.g., Kurtz et al.,2017) to predict and design early warning systems in susceptible regions(e.g., Chen et al., 2017) for climate change adaptation and environmental protection. High computational demands frequently necessitate
simplifying process representation such as preferential flow pathways(Tetzlaff et al., 2015), subsurface heterogeneity (Hartmann, 2016), etc
in large watersheds and limit our ability to scrutinize behavior at specific locations.
This research advances a novel computationally efficient watershedscale hydrologic modeling framework (e.g., Ajami et al., 2016). In
particular, a hillslope scale surface–subsurface coupling framework
developed by Chen et al. (2022) is extended to the watershed scale. Chenet al. (2022) externally coupled the KINEROS2 (K2) model (Goodrichet al., 2012), which solves the 1D kinematic wave equation for overland
flow (Woolhiser et al., 1990), with the HYDRUS-1D (H1D) model(Šimůnek et al., 2005, 2016b) that solves the Richards equation for thevariably saturated subsurface flow using a boundary switching algorithm. H1D and K2 are well-tested, state-of-the-art models tailored tosimulate specific hydrologic components (processes), often equippedwith highly refined and sophisticated analysis tools to provide invaluable insights for targeted investigation.
The overarching objective of coupling the H1D and K2 models is todevelop a straightforward, accurate, and computationally efficient hydrologic model at the watershed scale, enabling accurate simulation of
surface–subsurface hydrological processes at large scales. The computationally efficient watershed-scale model developed in this work is
referred to as the H1DK2-WS model. The coupled model integrates
specialized modeling and analysis capabilities of both H1D and K2
models while simulating hydrologic processes at the watershed scale. In
this approach, each watershed is subdivided into a cascade of 1D
overland-flow planes to route surface flow from hillslopes into channels
and ultimately to the watershed outlet using K2. Each hillslope andchannel reach are externally coupled to one or multiple vertical H1D soilcolumns in the subsurface, accommodating both homogeneous andheterogeneous systems. We assess the robustness and efficacy of themodeling approach through comparisons with benchmark simulations(Chen et al., 2022; Maxwell et al., 2014; Kollet and Maxwell, 2006) andevaluate the model’s capability to replicate observed hydrographs froma small semi-arid watershed divided into 22 plane and channel elements
within the Walnut Gulch Experimental Watershed WGEW (Goodrichet al. 2022; Moran et al. 2008).
The H1DK2-WS model can provide a comprehensive understandingof water fluxes and soil water dynamics across hydrologic componentsthat require an in-depth and nuanced understanding of specific processes. For example, the model may be used to study the effectiveness ofmanaged aquifer recharge (MAR) or changes in rechange due to factors
like climate and land use. This model also enables adaptation to evolving
conditions that standalone integrated models may struggle to achieve.
Consequently, the H1DK2-WS model provides a valuable tool for researchers, watershed managers, and policymakers striving to quantify
real-world conditions (e.g., soil moisture estimation, stream flow forecasting, and estimation of water fluxes and recharge) and to manage
environmental resources effectively (urban storm management, landuse planning, water management, and conservation planning). However, it is also acknowledged that process-based models are resourceintensive and challenging to work with (Jajarmizad et al., 2012), andthis sometimes makes them less flexible for site-specific and process-specific analyses.
2. Methods
2.1. Summary of K2 and H1D models
The K2 model is an open-source code that solves the 1D kinematicwave equation for surface runoff at various spatial scales (Morris andWoolhiser, 1980). It includes the Parlange three-parameter equation for
estimating infiltration amounts (Smith & Parlange, 1978). The kinematic wave equation in K2 is given as:∂h∂h+ αmhm− 1 = q(x, t) = P − F
∂t
∂x

(1)

where t is time [T], x is the distance along the slope direction [L], h is thedepth of surface water [L], and q is the source and sink term [LT− 1]
represented as a function of precipitation (P) [LT− 1] and infiltration (F)[LT− 1] components. The parameters α and m are flow parameters related
to surface roughness, flow regime, and slope.
K2 has been applied to a wide array of hydrological and environmental problems (Goodrich et al., 2006; Goodrich et al., 2000). The
Automated Geospatial Watershed Assessment (AGWA) tool (Miller et al.,2007; https://www.tucson.ars.ag.gov/agwa/) is a GIS-based hydrologic
modeling tool (ArcGIS 10.x) that uses commonly available GIS data
layers to efficiently processes necessary datasets that delineate watersheds, extract initial parameters from national databases, executes themodel, and supports visualizing simulation results in the GIS. K2 includes well-structured workflow and analysis tools for event-scale runoff
predictions (Goodrich et al., 1997) and uses precipitation data fromvarious sources (Yatheendradas et al., 2005). A continuous version of K2
is now available that also includes the treatment of snowfall and melt.
Hillslope model elements in K2 can be planar or curvilinear (especially
important for erosion and sediment transport) and can become narrower
or wider in the downslope direction. It has been successfully applied toassess the impacts of land use changes, urbanization, and green infrastructure (Hernandez et al., 1999; Kennedy et al., 2013; Korgaonkar
et al., 2020), as well as the effects of wildfires on watershed responses(Sidman et al., 2015; Guertin et al., 2019). The model provides highresolution spatially distributed simulation outputs across the study
watershed at a 1-minute or user-defined time step.
The H1D model is an open-source code that solves the Richards
equation for saturated–unsaturated flow in porous media (Šimůnek
et al., 2005, 2016b). The Richards equation is given as:
[
]
∂θ ∂
∂h
=
K(h) + K(h) − S(h)
(2)
∂t ∂z
∂z
where h is the water pressure head [L], θ is the volumetric water content[L3L-3], z is the spatial coordinate in the vertical direction [L] (positive
upward), S is the sink term [L3L-3T− 1], and K is the unsaturated hydraulic conductivity function [LT− 1].
The HYDRUS models (Šimůnek et al., 2016a, 2024) are among themost widely used numerical models for simulating various transport
processes, such as water flow and heat and solute transport in variably
saturated porous media (e.g., soils or the vadose zone). The HYDRUS
software includes an interactive graphics-based interface for data preprocessing, soil profile discretization, and graphic presentation of theresults. While HYDRUS-1D simulates flow and transport processes in
one dimension and is a public domain software, HYDRUS (2D/3D) (andearlier HYDRUS-2D) extends the simulation capabilities to the second
and third dimensions and is distributed commercially. There are many
thousands of HYDRUS users worldwide, with many applications (alsoseveral thousand) appearing in peer-reviewed literature and manytechnical reports (Šimůnek et al., 2016b, 2024). This model provides acomprehensive grasp of water dynamics in diverse soil and geologicalconditions (Kornelsen & Coulibaly, 2014; Chen et al., 2014), enabling2

M.B. Meles et al.

Journal of Hydrology 640 (2024) 131621

applications such as groundwater recharge estimation (Han et al., 2015;
Tonkul et al., 2019).
2.2. Overview of the H1D-K2 coupling framework
The hillslope H1D-K2 model (Chen et al., 2022) has undergone
rigorous testing and evaluation across a wide range of scenarios,
including various rainstorm events, different slopes, subsurface heterogeneity, and varying bottom boundary conditions (BCs). The workflow
also includes simulation of flow processes, computing the water balance
for each surface node along the single hillslope (plane element) atselected time steps to exchange the information between the H1D andK2 models. Its performance was then compared with the HYDRUS-2D
model (Šimůnek et al., 2006; Šimůnek, 2015), serving as a benchmark
for the evaluation. Results demonstrated that the coupled H1D-K2
model performs exceptionally well in reproducing hydrographs, surface water levels, soil moisture, and water balance component errors
compared to the HYDRUS-2D model. The model also maintains highcomputational efficiency due to a dynamic time-step approach and adimensionality reduction method where a single H1D soil profile represents an entire overland flow plane. The improvement in the computational time ranges from a factor of 1.3 to 10 to that of the HYDRUS-2D
model (Chen et al., 2022).
The watershed-scale coupling involves using a procedure initiallydesigned for the single hillslope scale in Chen et al. (2022). Expanding
the H1D-K2 model to a watershed scale model (H1DK2-WS) involves
dividing the watershed into a series of interconnected rectangular planeand channel elements. Each channel and plane element is associatedwith at least one H1D column to simulate water movement in the subsurface. The H1D column can be homogeneous or heterogeneous. The
planes and channel elements have unique IDs, so their hierarchicalconnections with neighboring elements, upstream, lateral, and downstream, are maintained to describe the flow.
The MATLAB (R2021b) code was derived from the hillslope scale andis designed to read and apply climate inputs to the unique element IDs(plane and channel). As described in the workflow (Fig. 1), the MATLAB
script executes the H1D and K2 codes sequentially using a series offunctions to simulate subsurface and overland flow processes, respectively. The H1D and K2 Fortran codes were modified to write temporary
files to allow the exchange of information between K2 and H1D.
The exchange of information between the two models occurs at
different time scales. The H1D model runs at a coarser time scale, Δthyd ,compared to the K2 model, which runs at a relatively finer resolution(Δtk2 : K2 time step). The H1DK2-WS code first executes the H1D model
for all elements (plane and channel) in the watershed and computes
infiltration at the surface, soil moisture redistribution across the soil
matrix, and recharge (outflow at the bottom of the soil column) at theΔthyd time step. H1DK2-WS then computes precipitation minus infiltration values at every sub-Δthyd time step (i.e., Δtk2 equal to 1 min), at
which K2 computes overland flow rates and surface water levels in eachelement. Δtk2 and Δthyd are multiples of the total simulation time whereΔthyd varies from 2 min to the maximum total simulation time dependingon the model convergence and subsurface soil column resolution.
Δthyd = Total Simtime /nk

Fig. 1. Flowchart depicting the workflow in the coupled H1DK2-WS code. The
MATLAB code operates within a sub-simulation time step (Δthyd : H1D simulation time step); the code loops through individual element IDs where H1D
executes computation of infiltration, surface water level, surface flux, soil
moisture redistribution (SM re-distn), and recharge information. Then the code
executes K2 every time step (i.e., Δtk2 : K2 simulation time step).

water depth (h(t)) at every node in each element. The K2 code runs therouting scheme across the entire watershed at every K2 time step Δtk2
between every simulation pause and start period Δthyd , in which theoutflow rate from the upstream hillslope (plane element) at every time
step is used as an inflow to the downstream plane or channel element inaddition to the precipitation input.
While executing the H1D and K2 codes and implementing the information exchange schemes, the H1DK2-WS model computes the waterbalance on the surface, which dictates the initial condition and infiltration rates to simulate flow for the next time step. In the interconnected K2 elements (plane and/or channel), outflow from upper
hillslopes serves as an inflow for lower elements. Therefore, at every
time step Δtk2 , the inflow rate in the downstream plane or channel
comprises both the outflow rate from upstream planes and qrmi . The key
functions in the exchange are the boundary condition switching functionin the water balance analysis at every time step and the computation of
the water depth across the hillslope (for details of the water balance andwater depth analysis, see Supplementary Material I).

(3)

where nk is an integer number indicating the number of H1D runs foreach watershed element. The Δtk2 time step ranges from the sub-minutescale (usually a small simulation time step in K2 runs) to larger time
steps that are a multiple of Δthyd .
The key component in the H1DK2-WS code was the replacement of
the Parlange infiltration equation (Parlange et al., 1982) by the Richards
equation in the H1D model (Šimůnek et al., 2005). Instead of subtracting
infiltration from surface water depth and precipitation input during flowrouting, the H1DK2-WS model routes the rainfall minus infiltration
amount, qrmi , to compute runoff outflow at every 1-minute interval and2.3. H1DK2-WS model setup
Setting up the H1DK2-WS model involves creating inputs anddomain representations for the K2 and H1D models. The model domains
include the rectangular plane and channel elements, with associatedparameter values describing the geometry, roughness, slope, soil hydraulic conductivity, and soil moisture at saturation and wilting points.
It also consists of the corresponding geometry and properties of the soilcolumns of all the elements representing the subsurface structures.
3

M.B. Meles et al.

Journal of Hydrology 640 (2024) 131621

This process varies depending on the domains used in this study.
These are:graphical user interface. Here, various parameters are defined, encompassing the type of subsurface material, boundary conditions, vertical
discretization, initial conditions, and time steps. H1D input requirements involve surface fluxes (evapotranspiration and rainfall). Soilproperties such as field capacity, wilting point, density, and sand, silt,
and clay fractions can be used to estimate soil hydraulic parameters forthe van Genuchten hydraulic model (van Genuchten, 1980) using aneural network prediction method—specifically, the Rosetta Lite module, version 1.1 (Schaap et al., 2001). This method relies on pedotransfer functions (PTFs) to estimate residual and saturated water contents, saturated hydraulic conductivity, and other soil hydraulic parameters. Evapotranspiration usually has an insignificant influence onrunoff generation because of the short simulation time, and it wastherefore neglected in this work, although the coupled model generally
accounts for it.

1. Hypothetical Domain: This includes simple geometric watersheds
used in previous model development and evaluation studies. Model
geometric domains and associated input files, parameters, and BCs
are defined using information available in the literature.
2. Lucky Hills Domain: Lucky Hills is an experimental watershed within
the WGEW in southeastern Arizona. Its modeling domain development requires using a Geographic Information System (GIS) userinterface − the Automated Geospatial Watershed Assessment(AGWA: Semmens et al., 2007) tool for the K2 setup.
The hypothetical domains used in this study include a single hillslope(e.g., Chen et al., 2022), three interconnected, cascading heterogenous
hillslopes, and a tilted V-catchment (Maxwell et al., 2014). These domains are easy to generate and set up without external pre-processing
tools. For the Lucky Hills domain, the AGWA GIS interface for K2,
developed in 2002 by USDA-ARS and the US Environmental Protection
Agency (Semmens et al., 2007; Miller et al., 2007), provides a preprocessing and an initial parametrization scheme. This interface offers auser-friendly platform that considers landscape information to derive
spatial model inputs and visualizes model results. AGWA relies on adigital elevation model (DEM) to obtain geomorphologic and geometricparameters, surface soil data (e.g., SSURGO data), and land use data
layer from NRCS, USGS, and NASA. This information contributes todeveloping a hierarchical structure and setting up the K2 model byfacilitating the subdivision of watersheds into rectangular overland flow
planes and channel elements, establishing IDs based on the hierarchicalorder for top-down flow routing. It extracts soil, slope, Manning’s coefficient, and land-use parameters and creates initial estimates of model
parameters and the associated file format for K2.
The setup process for the H1D model is conducted individually foreach model element, using the unique ID numbers through the H1D

2.4. Benchmark watersheds
The capabilities and performance of the coupled model were evaluated based on four benchmark watersheds (Figs. 2 and 3), representing adiverse range of physical processes and flux components that define theexpected standard simulation outcomes. In this coupling, we closely
examined the properties and behavior of simulation results, comparing
them with published hypothetical and observed results that allowed for
a comprehensive analysis of key differences and similarities. Water
balance errors and execution times for benchmark simulations were
calculated and compared to literature results for equally complex
modeling scenarios.
2.4.1. Hypothetical benchmark watersheds
Numerical experiments are a useful tool for model intercomparison.
This study explored similarities and differences between H1DK2-WS andother models using hypothetical benchmark simulations, which include
H1D-K2 (Chen et al., 2022), HYDRUS-2D (Šimůnek et al., 2006;

Fig. 2. A) a single hillslope domain for the infiltration excess, and saturation excess flow test cases with a 400 m long slope length, b) a connected heterogeneous
hillslope domain representing varying hillslope properties along the slope length, and c) a tilted V-catchment domain for overland and channel flow tests (after Kollet
and Maxwell, 2006; Sulis et al., 2010) (Sx and Sy are slope inclinations in the x and y directions).
4

M.B. Meles et al.

Journal of Hydrology 640 (2024) 131621

Two types of homogenous hillslopes were examined with a 400meter slope length, a unit width, and a subsurface thickness of 1 m(Fig. 2a). This setup aligns with the descriptions provided by Chen et al.
(2022) and has certain similarities to the hillslope in Maxwell et al.
(2014). The domain incorporates a no-flow BC, representing an impermeable base. The subsurface was discretized into 101 elements, withmesh sizes ranging from 0.0053 m near the surface to 0.012 m at thebottom. In the K2 model, there were a total of 41 surface nodes along thehillslope’s length. To initiate runoff generation, a rainfall event wasapplied at a rate of 3.3 × 10-4 m/s for 200 min, followed by a 100-minuterecession period. For this simulation, we considered two distinct runoff
generation scenarios based on the experimental setup:a. Infiltration Excess (Hortonian) Runoff.
One homogeneous hillslope scenario was designed to ensure pondingon the surface of a 0.0005 slope. To achieve the Horton runoff generation, we specified a smaller Ks (saturated hydraulic conductivity) value
than the rainfall rate (3.3 × 10–4 m/min) to allow runoff generation
before complete saturation of the soil column. In this scenario, we
considered free drainage bottom boundary condition and three different
Ks values of 6.69 × 10-11, 6.69 × 10-6, and 6.69 × 10-5 m/min to change
the amounts of infiltration and runoff (see Table 1 for details). The Ks of
10-11 m/min was assumed to be comparable to an impervious surface.
b. Saturation Excess (Dunn) Runoff.
The second homogeneous hillslope scenario was designed to ensure
ponding on a surface of a 0.0005 slope following complete saturation of
the soil column. To achieve the Dunn runoff generation, we specified alarger Ks (saturated hydraulic conductivity) value than the rainfall rate
(3.3 × 10–4 m/min) to allow complete saturation of the soil column
before surface ponding. In this scenario, we considered a no-flow bottom
boundary condition at the bottom and specified a Ks equal to 6.69 × 10-5
m/min. Runoff occurs after the soil profile becomes fully saturated (seeTable 1 for details).

Fig. 3. The Lucky Hills watershed was subdivided into 22 rectangular planeand channel elements representing a 0.04 km2 K2 watershed area. It is a small
headwater watershed in the Walnut Gulch Experimental Watershed (WGEW) insoutheast Arizona. Rainfall is recorded using two rain gauges (black circle), andrunoff is recorded using three flumes (open triangle).

Šimůnek, 2015), and other models whose results were reported by
Panday and Huyakorn (2004), Maxwell et al., (2014), and Sulis et al.
(2010). The hypothetical benchmark watershed cases intentionallyminimized the complexity in domain geometries and other features,such as topography, hydraulic and hydrogeological properties, and atmospheric forcing. However, they are complex enough to provide athorough basis for comparing model behavior.
The coupled H1DK2-WS code can simulate domains and flow conditions with varying model parameters, BCs, and initial conditions (seeTable 1). This includes different saturated hydraulic conductivities (Ks),
initial water content distributions, initial water tables, bottom BCs (i.e.,
no flow and free drainage), slope ratios, and rainfall intensities anddurations. The bottom BC represents the exchange flux between thevadose zone and groundwater. The no-flow BC represents an impermeable bottom layer, while the free drainage BC allows for water todrain from the bottom of the soil profile by gravity (Šimůnek et al.,2005). Model setups were selected in such a way as to match test scenarios that evaluate the impact of surface flow speed on infiltration andoverland flow and different bottom BC types affecting the subsurface
and surface flow processes. In this study, we used three hypotheticalmodel domains (Fig. 2 and Table 1).

3. Heterogeneous hillslope
Three hillslopes are interconnected in Fig. 2b to illustrate the influence of variability in subsurface materials and slopes. The H1DK2-WS
model domain encompasses three hillslopes and the H1D soil columns
with ID numbers 1 (uphill), 2 (middle), and 3 (downhill). Each ID had
distinct subsurface materials to better represent the spatial heterogeneity across the domain. The uphill and downhill H1D profiles were set
with a Ks of 6.69 × 10-6 m/s, while the middle hillslope had a Ks of 6.69
× 10-4 m/s. The slope lengths were 150 m, 100 m, and 150 m,

1. Homogeneous Hillslope

Table 1Model and watershed properties used in the benchmark watershed simulations. *shows specified calibrated values. **shows the use of observed data for the evaluationof the model.
Properties
H1DK2-WS Model structure
Horizontal mesh, Δx (m)− Exchange period Δthyd min
Watershed properties− Number of elements− Bottom BC− Initial water level (IWL) − m− Manning’s roughness (*10-4)− Rainfall rate (*10-4 m/min)− Slope (− )Soil properties− Number of soil profile− Ks (*10-5 m/min)− van Genuchten parameters− θsat (− )− θres (− )− n0 (− )− α0 (1/m)
Reference simulations

Homogeneous Hillslope
Saturation excess

Infiltration excess

Heterogeneous
Complex

Tilted V-catchment
Impervious layer

Lucy Hills Watershed
Observed catchment

10
10, 100, 300

10
10, 100, 300

10
10, 100, 300

10
10, 100, 300

10
10, 100, 300

1
No flow
0.5, 0.75, 1
3.31
3.30.0005

1
No flow1
3.31
3.3
0.0005

3
No flow1
3.31
3.3
0.0005, 0.01

3
No flow
1
3.31
1.8
0.05, 0.02

22
free drainage−*100–––1.3**Observed events*0.02–0.06

1
69.4

1
0.0694, 6.94, 69.4

1
0.694, 69.4

10.000694

1*0.06–––0.4

0.4
0.08
2
1
Chen et al. (2022)

0.4
0.08
2
1

0.4
0.08
2
1

0.4
0.08
2
1
Maxwell et al. (2014)*0.4 – 0.5*0.12 – 0.02
*2
*1
Meles et al. (2021)

5M.B. Meles et al.

Journal of Hydrology 640 (2024) 131621

respectively, with corresponding water depth calculation nodes of 15,
10, and 15 in K2 simulations. Additional simulations were conducted to
explore the impact of varying slopes. In this case, results from otherconditions were compared to a simulation with a slope of 0.0005 and aroughness coefficient of 0.02.

performed with the results associated with the OpenGeoSys (Delfs et al.,2013, 2012) and HydroGeoSphere (Therrien et al., 2012) hydrologic
simulation models for which the tilted V-catchment simulations areavailable. Both HydroGeoSphere (HGS) and OpenGeoSys (OGS) are oneof the integrated hydrological modeling systems extensively utilized for
simulating watershed processes. Both HGS and OGS employ a 3D
control-volume finite element groundwater model that rigorously conceptualizes interactions within groundwater and surface water systems
capable of solving diffusive-wave overland and Richards equations. In
addition to the flow modeling, the OGS model also simulates various
physics-based processes, including thermal, mechanical, and chemical
processes in porous and fractured media.
The visual graphical inspection was also carried out using various
simulation outputs related to the HSPF (Bicknell et al., 1993), HEC-1(USACE, 1998), ParFlow (Kollet and Maxwell, 2006), and CATHY(Camporese et al., 2010) models. The evaluation metrics included how
the simulated hydrographs and soil moisture values of the coupled
H1K2-WS model tracked the time series of soil moisture redistribution,
the peak flow rate, and the time to peak runoff.
2) Statistical predictive skills: This evaluation used statistical metrics
to evaluate the model’s ability to simulate observed hydrographs.
Goodness-of-fit statistics such as correlation coefficient, Nash-Sutcliffe
Efficiency (NSE), Percent Bias, and Root Mean Square Error (RMSE)
were computed to objectively assess the model using the calibrated
model outputs. In this process, we implemented a manual calibration
utilizing selected parameters that were deemed to influence the twomodels independently. The small size of the watershed facilitated astrategic and efficient calibration process, requiring minimal time andeffort despite using the manual calibration approach. It’s crucial torecognize that manual calibration inherently relies on the modeler’s
expertise and understanding of the model’s fundamental mechanisms
and interactions (Eckhardt & Arnold, 2001). This introduces a degree ofsubjectivity and provides faster calibration time. Ultimately, parameter
sets yielding the best goodness of fit were selected for further analysis.
Moving forward, future versions of the coupled model will include thecoupling with the sensitivity analysis and automated calibration
approaches.
3) Computational efficiency and water balance errors: Computational efficiency was gauged by assessing the computational timerequired for each analysis scenario and comparing it to available data.
We also computed the total water balance errors in each simulation.

3.1. A tilted V-catchment
The tilted V-catchment was described in previous works by Gottardi
and Venutelli (1993), Panday and Huyakorn (2004), Kollet and Maxwell(2006), Sulis et al. (2010), and Maxwell et al. (2014). It is formed byjoining two inclined planar rectangles to a central channel. Each of theserectangles has a width of 800 m and a length of 1000 m, and the channel
has a width of 20 m (Fig. 2c). The geometry of the tilted V-catchment isgiven in Table 1 and Fig. 2c.
This test case is utilized to evaluate the capability of the H1DK2-WS
code to route the surface component. This scenario exclusively focuseson surface processes and excludes subsurface contributions (i.e., minimal infiltration). The saturated hydraulic conductivity value of thewatershed material was set to 6.69 × 10-11 m/s to represent an impervious surface described in the Maxwell et al. (2014) setup. The simulation scenario involves a 90-minute uniform rainfall event of 1.8 × 10-4
m/min, followed by a 90-minute drainage phase (i.e., a total simulation
period of 180 min). Results are compared with other models in theliterature (e.g., Maxwell et al., 2014).
3.1.1. Observed Experimental watershed − Lucky Hills watershed
In this benchmark simulation, we considered a relatively complex
headwater watershed called Lucky Hills within Arizona’s Walnut Gulch
Experimental Watershed (WGEW) (Goodrich et al., 2008). Lucky Hills isa small sub-watershed (Fig. 3) with an area of 0.04 km2. The watershed
has two weighing-type recording rain gauges, RG83 and RG384, with a
precision of 0.25 mm and a 1-minute time step for rainfall recording.
The runoff is measured with an H flume and a supercritical Santa Rita
flume at two internal gauging stations (F106 and F102, respectively) anda Santa Rita flume at the outlet (F104) (Stone et al., 2008; Smith, 1981).
Soils in the study watersheds are dominantly sandy loam (Ustochreptic
Calciorthid) − with a desert shrub-dominated cover condition (e.g.,
Skirvin et al., 2008).
The Lucky Hills watershed was discretized using 22 model elements(planes or channels) and 22 H1D profiles with 1 m depths and freedrainage bottom boundary conditions. We simulated three rainfall
events with runoff observations described by Meles et al. (2021). The
event data were selected from summer rainfall: July 25, 2010, July 26,
2010, and July 14, 2015, with a dry initial condition defined as a minimum of 24 h without rain with observed pre-storm moisture content
values of 6–7 % VMC (volumetric water content) based the soil moisturesensors installed 5 cm below the surface in the watershed. These data are
publicly available in the SWRC-DAP database (https://www.tucson.ars.
ag.gov/dap).

4. Results and discussion
4.1. Homogeneous Hillslope: Infiltration and saturation excess scenarios
The results of the infiltration and saturation excess test cases areshown in Figure SII1, illustrating model parameters and scenarios outlined in Table 1. The surface outflow rate and the soil moisture redistribution with depth were plotted across the simulation period for acomparative analysis between the time series of the hillslope scalemodel (Chen et al., 2022) and HYDRUS-2D (Šimůnek et al., 2006;
Šimůnek, 2015) model simulations. In both infiltration and saturation
excess scenarios, the coupled model H1DK2-WS accurately reproduced
the benchmark simulation (for details, see Supplementary Material II).

3.2. Model evaluation
In this study, we used a comprehensive framework aimed to thoroughly evaluate the coupled model’s performance. This encompassed acomparison with established results, statistical measures of predictive
skills, and considerations of computational efficiency and water balanceerrors.
1) Comparative assessment: This evaluation involved a comparisonof published results from Chen et al. (2022), the HYDRUS-2D model(Šimůnek et al., 2006; Šimůnek, 2015), and an analytical solution(Gottardi and Venutelli, 1993) for both a single homogenous hillslope(Fig. 2a) and the three connected heterogeneous hillslopes (Fig. 2b). For
evaluating the tilted V-catchment watershed (Fig. 2c), the comparison
was conducted through graphical inspection using published simulation
results (e.g., Maxwell et al., 2014). Here, the graphical comparison was

4.2. Heterogeneous hillslope with a uniform rainfall rate
The heterogenous hillslope (Fig. 2b) represented scenarios that
capture variations in topographic and geomorphic features, geometricdimensions, surface and subsurface soil materials, and cover conditions.
In this exercise, the results of the three-plane element-based H1DK2-WS
model configuration accurately reproduced the simulated hydrograph ofthe single hillslope-based HYDRUS-2D model (Šimůnek et al., 2006) andthe H1D-K2 hillslope model (Chen et al., 2022) with heterogeneousmaterials. Compared to both reference simulations, the coupled
6M.B. Meles et al.

Journal of Hydrology 640 (2024) 131621

watershed scale model has exhibited precise tracking of water waves,
except for an early arrival of the second flow wave from the upstream
planes (Fig. 4).
Compared to the HYDRUS 2D model, the second flow wave from the
upstream hills, simulated using the H1DK2-WS and H1K2 models, arrives 7 min and 4 min quicker, respectively. The uphill (Plane 1) anddownhill (Plane 2) planes are the low-permeability plane elements withKs equal to 6.94 × 10-6 m/min, resulting in immediate runoff generation. The middle hillslope, with a highly conductive material (Ks = 6.94
× 10-4 m/min), saturates first before runoff generation starts at about 80
min, initiating the second flow wave at the outlet (at the end of thedownhill plane). The hydrograph simulated by H1K2-WS shows that
outflow began at 3.8 mins, increased to 6.52 × 10-4 m2/min over 50
mins, and then gradually increased until the arrival of the second wave
at about 110 mins (Fig. 4). The time to peak in all three models is equal
at t = 200 min.
Additionally, slight disparities were observed in both the arrival timeof the second flow wave and the total runoff volume between themodels. These differences were 1.2 % and 2.4 % for H1D-K2 and H1DK2WS, respectively, compared to the HYDRUS-2D runoff volume of 18.46
m3 during the simulation period. However, the peak flow, time to peak,runoff initiation time, and runoff rate (over 98 % of the simulation time)remained consistent across the HYDRUS-2D, H1D-K2, and coupled
H1DK2-WS models. The minor discrepancies in runoff between themodels can be attributed to differences in the model domain setup andnumerical solutions. One of the factors leading to the mismatch of thehydrographs plotted in Fig. 4 and subsequent runoff figures could be thedefinition of kinematic wave solution that usually does not include
diffusion. However, an off-centered finite difference scheme that is of
the first order, as described by Ponce (1991), introduces varying degrees
of numerical diffusion depending on the size of Δx and Δt.
We compared the soil water content profiles simulated using theH1K2-WS model throughout the simulation period to the H1D-K2 andHYDRUS-2D codes in Fig. 5. The analysis of water movement during thesimulation period at various points along the hillslope − downhill(Fig. 5a), middle (Fig. 5b), and uphill (Fig. 5c) − demonstrated a similar
distribution pattern between the H1D-K2 and HYDRUS-2D models.
Moreover, the examination across multiple plane elements indicated abroad approximation of soil profile conditions and water movement

across the heterogeneous hillslope (see Figure SII2 and 5). It is important
to note that the H1DK2-WS model effectively captures the generalbehavior of heterogeneous hillslopes similar to H1D-K2 and HYDRUS2D, particularly in the central points of individual hillslopes. In the
middle of the hillslope, the high conductive subsurface material (Ks =6.94 × 10-4 m/min) attained saturation throughout the columns in lessthan 100 min before saturation. In contrast, the uphill and downhill
planes (Ks = 6.94 × 10-6 m/min) registered water content changed only
through the top 20 cm depth over the entire simulation period.
4.3. Tilted V-catchment
Fig. 6 (black and gray hydrographs) shows the runoff outflow rates
computed by various models for the tilted V-catchment with a rainfall
intensity of 1.8 × 10-4 mm/min over 90 min. A very low hydraulicconductivity of Ks = 6.94 × 10-11 m/min was assigned to all three elements to represent an impervious surface. Runoff outflow in m3/s wascompared visually and graphically to the results presented by Maxwell
et al. (2014), Kollet et al. (2006), and Panday & Huyakorn (2004). The
published findings revealed significant variability between model results, underscoring the greater complexity of the domain and the diverse
overland and channel routing schemes applied in each model.
In this analysis, we compared simulation results from two models,
OpenGeoSys (Delfs et al., 2013, 2012) and HydroGeoSphere (Therrien
et al., 2012), alongside our coupled H1DK2-WS model. These simulations were conducted within the same model domain (refer to Fig. 2b),
employing identical parameter values and initial conditions (seeTable 1). The evaluation included the similarities and differences
through a graphical comparison of hydrographs (see Fig. 6), focusing on
runoff initiation, the slope of the rising limb, steady-state flow, cessationof steady-state flow, and the slope of the recession side. The results
showed that the performance of the H1DK2-WS model is well within thebounds of the OpenGeoSys and HydroGeoSphere models.
Furthermore, the visual inspection of published simulated runoff
rates compared to model outputs such as HSPF (Bicknell et al. 1993),
HEC-1 (USACE, 1998), ParFlow (Kollet and Maxwell, 2006), CATHY(Camporese et al., 2010), and tRIBS (Ivanov et al., 2004), indicates thatthe coupled model stands closely alongside the array of model outputs
noted above. This underscores the general consistency of theFig. 4. Comparison of runoff outflow rates from a heterogeneous surface simulated using multiple models. The H1DK2-WS model (black-solid line) at the end of thedownhill plane (Plane 3) is compared with the hillslope-scale H1D-K2 code (brown-dashed lines) and HYDRUS-2D (red-dashed line). The HYDRUS-2D and H1D-K2(hillslope) models consider single hillslope domains with heterogeneous material properties. The H1DK2-WS domain includes three connected plane elements withvarying soil properties (see Fig. 2b and Table 1). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version ofthis article.)7M.B. Meles et al.

Journal of Hydrology 640 (2024) 131621

Fig. 5. Comparison of subsurface soil water content distributions over the simulation period between H1D-K2 (hillslope: broken line) and H1K2-WS (solid line). A
black solid line represents the initial soil moisture values at t = 0 min, while colored lines represent soil moisture profiles at various times throughout the simulation
period. Panels (a) and (c) illustrate soil moisture profiles under the uphill and downhill planes, respectively, with a Ks of 6.94 × 10-6 m/min, while panel (b)represents soil moisture profiles under the middle plane (Plane 2) with a Ks of 6.94 × 10-4 m/min.

calibrated saturated water content for the three events seems to be on
the higher end of the range for the soil type.
The calibrated model produced simulation results with highgoodness-of-fit, with a minimum Nash efficiency coefficient of 0.78 anda correlation coefficient of 0.91 for the smallest sub-watershed. The
computed root mean square error (RMSE) and mean error (ME) were
negligible, indicating the coupled model’s ability to simulate actual
watershed processes accurately. Some of the most significant discrepancies were related to the peak flow prediction in the larger flow events.
The calibrated model tends to overestimate (+ve peak runoff error) the
peak flow rate, especially in large events on 07/30/2015 and 07/26/
2010. The peak runoff error, computed as the difference between
simulated peak runoff and observation normalized by the observed peak
flow, showed less than 12 % on average, with the maximum discrepancy
on the larger flow events at the watershed outlet reaching 22 %. Additionally, comparing the time to peak error, a difference between simulated time minus observation expressed in minutes, demonstrated nearperfect prediction, with differences ranging between 2 min ahead (− vevalue) and 1 min behind the observed time to peak. This accuracy underscores the model’s potential applicability in flood early warning
prediction, highlighting its reliability in forecasting critical hydrologicalevents. It is also encouraging that the manual calibration implemented
in this study yielded predictions slightly superior to those of autocalibrated K2 models for the three events (e.g., Meles et al., 2021).
Previous research has underscored the importance of incorporating
spatial rainfall variability in both space and time to enhance the accuracy of runoff prediction in this study area (Goodrich et al., 1995; Faures
et al., 1995); to this end, the modeling exercise has included the two rain
gauges in Luck Hills representing the rainfall variability.
While acknowledging these positive outcomes, it is crucial toemphasize the potential for further enhancement in the predictive capabilities of the coupled model when applied to large watersheds. Future
efforts will focus on deepening our understanding of model parameter
sensitivity, identifiability, and importance across diverse scenarios andscales, requiring automatic optimization tools such as PEST (Doherty,2004) or GLUE (Beven and Binley, 1992). Optimization algorithms andadvanced sampling techniques for auto parameter identifiability andsensitivity analysis will be instrumental in improving predictive skills
and increasing the implementation of the model over a broader range of
scenarios. As we advance in these endeavors, the coupled model holds

Fig. 6. Comparison of hydrographs for the hypothetical tilted V-catchment
watershed simulated using various hydrologic models.

hydrographs produced by the coupled H1DK2-WS model with thosegenerated by multiple other physically based integrated or coupledsurface–subsurface flow models. The most significant deviations in thevisual comparison occurred during the rising phase of the hydrographrather than the recession limb phase, which is consistent with the observations in Maxwell et al. (2014).
4.4. Lucky Hills watershed
The calibrated coupled H1DK2-WS model accurately predicted theobserved hydrographs at the watershed outlet and the internal gaging
stations for three independent events in the Lucky Hills experimental
watershed (Fig. 7). We implemented a manual calibration utilizing theparameters outlined and parameter spaces specified in Table 2. In thisapproach, we independently calibrated model parameters for the subwatersheds leading to each internal gauging station, encompassing
both overland flow and corresponding H1D parameters. As shown in
Table 3, the model parameters are within acceptable ranges of theparameter spaces for the surface and subsurface properties of thewatershed with sandy loam soil type, as detailed in Section 2.4. The
8

M.B. Meles et al.

Journal of Hydrology 640 (2024) 131621

Table 3
Summary of statistical evaluation metrics evaluating the performance of thecalibrated H1DK2-WS model.
Error Type

Station 104

Station 102

Station 106

Correlation Coeff [-]RMSE [m3/s]
Mean Error [m3/s]
NSC [-]
Peak Runoff Error[%]Time to Peak Error [min]0.98 – 0.99
0.002 – 0.03
0.0007 – 0.0012
0.955 – 0.97
5.85 – 21.6
− 1.0 – 1.0

0.97 – 0.98
0.001 – 0.013− 0.0012
0.94 – 0.96
− 4.5–––20.1− 1.0 – 1.0

0.91 – 0.96
0.001––0.008− 0.0007
0.78 – 0.92
5.1–––15.4− 2.0–––0.0

substantial promise for even more accurate and reliable predictions infuture applications.
4.5. Water balance errors
The water balance errors in the coupled systems were minimal, witha maximum of 0.31 % for the tilted V-catchment (Fig. 8). The waterbalance error in this context refers to the difference between precipitation input, outflow (runoff), a storage change in the surface and subsurface, and outflow at the bottom of the soil column, ultimately
accounting for recharge.
4.6. Computational time
Utilizing the common exchange time of 1 min, we compared thecomputational time required to simulate both a single hillslope and aseries of heterogeneous hillslopes (3 connected uniform hillslopes)
against the computational time of the H1D-K2 and HYDRUS-2D models(Fig. 9). Similar to the enhancements demonstrated by the H1D-K2
model, the coupled watershed model (H1DK2-WS) exhibited a slight
improvement related to the implementation of the serialization technique (see Supplementary Material I). According to Chen et al. (2022), theshortest computation time for the hillslope model (H1D-K2) with threesoil profiles (refer to Section 4.2) was 21.36 s, followed by the HYDRUS2D model at 32.49 s. The analogous model configured with three connected hillslopes representing heterogeneous scenarios demonstrated
improvements in computation time even compared to the single hillslope H1D-K2 model. The H1DK2-WS model can easily integrate thechannel routing with three hillslopes with minimal additional computational demand. This improvement also extends to a single homogeneous hillslope domain, where the coupled watershed scale model(H1DK2-WS) required a mere 17.3 s compared to 21.36 s and 32.49 s for
the H1D-K2 and HYDRUS-2D models, respectively which is 14.3 % and
47 % improvement respectibely.
Despite substantial improvements in computational time achievedby the coupled model, particularly for simple scenarios, parallel
computing platforms present substantial opportunities for further
enhancement. This would enable the model not only to explicitly couple
various hydrologic components efficiently but also to create ensemble
modeling techniques and conduct long-term simulations for model
evaluation, assessment, and impact studies, such as those involving
evaluations of climate change and management effects.
Various efforts have demonstrated the utilization of distributed hydrologic models of varying capabilities on parallel computing platforms(e.g., Kollet et al., 2010; Kollet and Maxwell, 2008; Vivoni et al., 2011).
The utilization of an external exchange strategy and hierarchical simulations across plane and channel elements is well-suited for implementing various parallel computing strategies like functional
parallelism that conduct independent and simultaneous model simulations (e.g., Harada and Alba, 2020), as well as multifunctional parallelism (Vrugt et al., 2006; Liao et al., 2017).

Fig. 7. Comparison of three calibrated hydrographs at three gauging stations
simulated using the coupled H1DK2-WS model with observed data. The colors
show the gauging stations. The solid lines represent calibrated simulations, andthe dashed lines represent the observation data.
Table 2The initial H1K2-WS parameter values for the Luck Hills watershed estimatedfrom the existing databases and calibrated ranges.
ModelsParameter
description

Parameters

Unit

Databaseparameter range

K2

Roughness
Coefficient
Hydraulicconductivity
Residual watercontent
Saturated watercontentN

sec/m1/3

0.146 – 0.0125

Ks

8.14 – 4.54

θr

m/min
(*10-5)
m3/m3

θs

m3/m3

0.02 – 0.06

HYDRUS1D

0.51 – 0.55

5. Conclusion
In this study, we developed the H1DK2-WS model, which couples
9M.B. Meles et al.

Journal of Hydrology 640 (2024) 131621

Fig. 8. Total water balance errors of the coupled H1DK2-WS model for various simulation scenarios.

development and parallel computing implementation.
CRediT authorship contribution statement
Menberu B. Meles: Writing – original draft, Visualization, Validation, Software, Resources, Project administration, Methodology, Investigation, Formal analysis, Conceptualization. Lin Chen: Writing –review & editing, Methodology, Conceptualization. Carl Unkrich:
Writing – review & editing, Software, Conceptualization. Hoori Ajami:
Writing – review & editing, Methodology, Conceptualization. Scott A.
Bradford: Writing – review & editing, Supervision, Project administration, Methodology, Funding acquisition, Conceptualization. JiříŠimůnek: Writing – review & editing, Software, Methodology,
Conceptualization. Dave C. Goodrich: Writing – review & editing.

Fig. 9. Comparison of computational times for the reference simulations between the H1DK2-WS model (blue), the hillslope H1D-K2 code (red), andHYDRUS-2D (green). (For interpretation of the references to colour in thisfigure legend, the reader is referred to the web version of this article.)

Declaration of competing interest

H1D and K2 models at the watershed scale. This model is evaluatedagainst benchmark simulations and observed watershed data. The
integration encompasses a novel coupling approach with BC switchingto account for surface ponding and water exchange. The benchmark
simulations included a variety of test scenarios, ranging from hypothetical to real-world watershed conditions. Furthermore, we comparedthese benchmark simulations with models employing different structures and solution techniques, including HYDRUS-2D with an overland
flow routine, OpenGeoSys, and HydroGeoSphere.
The H1DK2-WS model demonstrated excellent consistency in termsof subsurface soil moisture redistribution and outflow rates at the hillslope scale for both homogeneous and heterogeneous subsurface scenarios consistent with the H1D-K2 (Chen et al., 2022), HYDRUS-2D(Šimůnek et al., 2006; Šimůnek, 2015). Despite slight differences due
to numerical methods, the model provided results consistent with multimodel outputs especially in the hypothetical benchmark simulations.
Applying combined water balances and BC switching in coupling
H1D and K2 has resulted in accurate model results. This approacheffectively accounts for surface ponding and water exchange between
the models, improving the accuracy in simulating hydrological fluxes
with minimal water balance errors.
The subdivision of hillslopes into a cascade of planes with varying
slopes, soil characteristics, and land cover conditions allows for incorporating complex spatial heterogeneity. This approach enhances themodel’s ability to capture the nuances of real-world landscapes,
including variability in topographic parameters (slope), roughness, land
use and land cover properties, and soil and subsurface properties,
thereby accurately representing hydrologic responses of domain complexities. Additionally, H1DK2-WS model demonstrates significant
computational efficiency, offering faster simulation times than
HYDRUS-2D. even without implementing any parallelization techniques. The findings affirm the model’s accuracy in simulating hydrological dynamics across scales and underscore its potential for further

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.
Acknowledgments
This research was supported by the United States Department ofAgriculture, Agricultural Research Service, National Program 211.
USDA is an equal opportunity provider, employer, and lender. Mentionof a proprietary product does not constitute endorsement by USDA or
the US government and does not imply its approval to the exclusion of
the other products that may also be suitable.
Appendix A
Supplementary Material I.
Appendix B
Supplementary Material II.
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. Menberu Meles: menberu.meles@usda.gov.
10M.B. Meles et al.

Journal of Hydrology 640 (2024) 131621

Appendix C. Supplementary data

Han, M., Zhao, C., Šimůnek, J., Feng, G., 2015. Evaluating the impact of groundwater on
cotton growth and root zone water balance using Hydrus-1D coupled with a crop
growth model. Agricultural Water Management 160, 64–75.
Harada, T., Alba, E., 2020. Parallel genetic algorithms: a useful survey. ACM Computing
Surveys (CSUR) 53 (4), 1–39.
Hartmann, A., 2016. Putting the cat in the box: why our models should consider
subsurface heterogeneity at all scales. Wiley Interdisciplinary Reviews: Water 3 (4),478–486.
Hernandez, Mariano, Scott N. Miller, David C. Goodrich, Bruce F. Goff, William G.
Kepner, Curtis M. Edmonds, and K. Bruce Jones. “Modeling runoff response to land
cover and rainfall spatial variability in semi-arid watersheds.” In Monitoring
Ecological Condition in the Western United States: Proceedings of the Fourth Symposiumon the Environmental Monitoring and Assessment Program (EMAP), San Franciso, CA,
April 6–8, 1999, pp. 285-298. Springer Netherlands, 2000.
Jajarmizad, M., Harun, S., Salarpour, M., 2012. A review on theoretical considerationand types of models in hydrology. Journal of Environmental Science and Technology
5 (5), 249–261.
Kan, G., He, X., Li, J., Ding, L., Hong, Y., Zhang, H., Liang, K., Zhang, M., 2019. Computer
aided numerical methods for hydrological model calibration: an overview and recent
development. Archives of Computational Methods in Engineering 26, 35–59.
Kennedy, J.R., Goodrich, D.C., Unkrich, C.L., 2013. Using the KINEROS2 modelingframework to evaluate the increase in storm runoff from residential development ina semiarid environment. Journal of Hydrologic Engineering 18 (6), 698–706.
Kollet, S.J., Maxwell, R.M., 2006. Integrated surface–groundwater flow modeling: A freesurface overland flow boundary condition in a parallel groundwater flow model.
Advances in Water Resources. 29 (7), 945–958. https://doi.org/10.1016/j.
advwatres.2005.08.006.
Korgaonkar, Y., Meles, M.B., Guertin, D.P., Goodrich, D.C., Unkrich, C., 2020. Global
sensitivity analysis of KINEROS2 hydrologic model parameters representing green
infrastructure using the STAR-VARS framework. Environmental Modelling &Software 132, 104814. https://doi.org/10.1016/j. envsoft.2020.104814.
Kornelsen, K.C., Coulibaly, P., 2014. Root-zone soil moisture estimation using datadriven methods. Water Resources Research 50 (4), 2946–2962.
Kurtz, W., Lapin, A., Schilling, O.S., Tang, Q., Schiller, E., Braun, T., Hunkeler, D.,
Vereecken, H., Sudicky, E., Kropf, P., Franssen, H.J.H., 2017. Integratinghydrological modelling, data assimilation and cloud computing for real-time
management of water resources. Environmental Modelling & Software 93, 418–435.
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 Resources Research. 50 (2),1531–1549. https://doi.org/10.1002/2013WR013725.
Meles, M.B., Goodrich, D.C., Gupta, H.V., Burns, I.S., Unkrich, C.L., Razavi, S.,
Guertin, D.P., 2021. Multi-criteria, time dependent sensitivity analysis of an eventoriented, physically-based, distributed sediment and runoff model. Journal of
Hydrology 598, 126268.
Miller, S.N., Semmens, D.J., Goodrich, D.C., Hernandez, M., Miller, R.C., Kepner, W.G.,
Guertin, D.P., 2007. The automated geospatial watershed assessment tool.
Environmental Modelling & Software 22 (3), 365–377.
Moradkhani, H., Hsu, K.L., Gupta, H., Sorooshian, S., 2005. Uncertainty assessment ofhydrologic model states and parameters: Sequential data assimilation using theparticle filter. Water Resources Research 41 (5).
Moran, M.S., Emmerich, W.E., Goodrich, D.C., Heilman, P., Holifield Collins, C.,
Keefer, T.O., Nearing, M.A., Nichols, M.H., Renard, K.G., Scott, R.L., Smith, J.R.,
Stone, J.J., Unkrich, C.L., Wong, J.K., 2008. Water Resources Research Vol. 44,
W05S01.
Panday, S., Huyakorn, P.S., 2004. A fully coupled physically-based spatially-distributedmodel for evaluating surface/subsurface flow. Advances in Water Resources 27 (4),361–382.
Parlange, J.Y., Lisle, I., Braddock, R.D., Smith, R.E., 1982. The three-parameter
infiltration equation. Soil Science 133 (6), 337–341.
Ponce, V.M., 1991. Kinematic wave controversy. Journal of Hydraulic Engineering 117(4), 511–525.
Richards, L.A., 1931. Capillary conduction of liquids through porous mediums. Physics 1,
318–333.
Schaap, M.G., Leij, F.J., Van Genuchten, M.T., 2001. Rosetta: A computer program forestimating soil hydraulic parameters with hierarchical pedotransfer functions.
Journal of Hydrology 251 (3–4), 163–176.
Semmens, D.J., Goodrich, D.C., Unkrich, C.L., Smith, R.E., Woolhiser, D.A. and Miller, S.
N., 2007. KINEROS2 and the AGWA modelling framework. Hydrological Modelling inArid and Semi-Arid Areas, p.49.
Shelia, V., Šimůnek, J., Boote, K., Hoogenbooom, G., 2018. Coupling DSSAT and
HYDRUS-1D for simulations of soil water dynamics in the soil-plant-atmospheresystem. Journal of Hydrology and Hydromechanics 66 (2), 232–245.
Shen, C., Phanikumar, M.S., 2010. A process-based, distributed hydrologic model basedon a large-scale method for surface–subsurface coupling. Advances in Water
Resources 33 (12), 1524–1541.
Shi, Y., Davis, K.J., Duffy, C.J., Yu, X., 2013. Development of a coupled land surfacehydrologic model and evaluation at a critical zone observatory. Journal of
Hydrometeorology 14 (5), 1401–1420.
Sidman, G., Guertin, D.P., Goodrich, D.C., Unkrich, C.L., Burns, I.S., 2015. Risk
assessment of post-wildfire hydrological response in semiarid basins: the effects ofvarying rainfall representations in the KINEROS2/AGWA model. International
Journal of Wildland Fire 25 (3), 268–278.
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.

Supplementary data to this article can be found online at https://doi.
org/10.1016/j.jhydrol.2024.131621.
References
Ajami, H., Khan, U., Tuteja, N.K., Sharma, A., 2016. Development of a computationallyefficient semi-distributed hydrologic modeling application for soil moisture, lateralflow and runoff simulation. Environmental Modelling & Software 85, 319–331.
Bicknell, B.R., Imhoff, J.C., Kittle Jr., J.L., Donigian, A.S., Johanson, R.C., 1993.
“hydrologic Simulation Program-FORTRAN (HSPF): Users Manual for Release, 10”.
Rep. No. EPA/600/R-93/174. U.S. EPA, Athens, Ga.
Bierkens, M.F., Bell, V.A., Burek, P., Chaney, N., Condon, L.E., David, C.H., de Roo, A.,Döll, P., Drost, N., Famiglietti, J.S., Flörke, M., 2015. Hyper-resolution globalhydrological modelling: what is next? “Everywhere and locally relevant”.
Hydrological Processes 29 (2), 310–320.
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.
Journal of Hydrology. 620, 129515 https://doi.org/10.1016/j.jhydrol.2023.129515.
Camporese, M., Paniconi, C., Putti, M., Orlandini, S., 2010. Surface-subsurface flow
modelling with path-based runoff routing, boundary condition-based coupling, andassimilation of multisource observation data. Water Resources Research 46,
W02512. https://doi.org/10.1029/2008WR007536.
Chen, Y., Li, J., Wang, H., Qin, J., Dong, L., 2017. Large-watershed flood forecasting withhigh-resolution distributed hydrological model. Hydrology and. Earth System
Science 21, 735–749. https://doi.org/10.5194/hess-21-735-2017.
Chen, L., Šimůnek, J., Bradford, S.A., Ajami, H., Meles, M.B., 2022. A computationallyefficient hydrologic modeling framework to simulate surface-subsurfacehydrological processes at the hillslope scale. Journal of Hydrology 614, 128539.
Chen, M., Willgoose, G.R., Saco, P.M., 2014. Spatial prediction of temporal soil moisturedynamics using HYDRUS-1D. Hydrological Processes 28 (2), 171–185.
de Saint-Venant, B., 1871. Theory of unsteady water flow, with application to riverfloods and to propagation of tides in river channels: French Academy of. Science v,
73.
Delfs, J.-O., Blumensaat, F., Wang, E.W., Krebs, P., Kolditz, O., 2012. Couplinghydrogeological with surface runoff model in a Poltva case study in Western
Ukraine. Environmental Earth Sciences 65, 1439–1457.
Delfs, J.-O., Wang, W., Kalbacher, W.T., Singh, A.K., Kolditz, O., 2013. A coupledsurface/subsurface flow model accounting for air entrapment and air pressurecounterflow. Environmental Earth Sciences 69 (2), 395–414. https://doi.org/
10.1007/s12665-013-2420-1.
Faures, J.-M., Goodrich, D.C., Woolhiser, D.A., Sorooshian, S., 1995. Impact of smallscale spatial rainfall variability on runoff modeling. Journal of Hydrology 173,
309–326.
Goodrich, D.C., Faures, J.-M., Woolhiser, D.A., Lane, L.J., Sorooshian, S., 1995.
Measurement and analysis of small-scale convective storm rainfall variability.
Journal of Hydrology 173, 283–308.
Goodrich, D.C., Lane, L.J., Shillito, R.M., Miller, S.N., Syed, K.H., Woolhiser, D.A., 1997.
Linearity of basin response as a function of scale in a semiarid watershed. Water
Resources Research 33 (12), 2951–2965.
Goodrich, D.C., Keefer, T.O., Unkrich, C.L., Nichols, M.H., Osborn, H.B., Stone, J.J.,
Smith, J.R., 2008. Long-term precipitation database, Walnut Gulch Experimental
Watershed, Arizona, United States. Water Resources Research 44, W05S04. https://doi.org/10.1029/2006WR005782.
Goodrich, D.C., Burns, I.S., Unkrich, C.L., Semmens, D.J., Guertin, D.P., Hernandez, M.,
Yatheendradas, S., Kennedy, J.R., Levick, L.R., 2012. KINEROS2/AGWA: Model use,
calibration, and validation. Trans. of the ASABE 55 (4), 1561–1574.
Goodrich, D.C., Bosch, D., Bryant, R., Cosh, M.H., Endale, D., Veith, T., Kleinman, P.,
Langendoen, E., McCarty, G., Pierson, F., Schonbetg, H., Smith, D., Starks, P.,
Strickland, T., Tsegaye, T., Awada, T., Swain, H., Derner, J., Bestelmeyer, B.,
Schmer, M., Baker, J., Carlson, B., Huggins, D., Aqrchera, D., Armendariz, G.A.,
2022. Long term agroecosystem research experimental watershed network.
Hydrological Processes 36(6), Article e14534.
Goodrich, D. C., Chehbouni, A., Goff, B., MacNish, B., Maddock, T., Moran, S.U.S.A.N.,
Shuttleworth, W.J., Williams, D.G., Watts, C., Hipps, L.H. and Cooper, D.I., 2000:
Preface paper to the Semi-Arid Land Surface-Atmosphere (SALSA) program special
issue, Agricultural and Forest Meteorology, 105(1–3), 3–20, doi:10.1016/ S0168-1923(00)00178-7.
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. Advances in Water Resources 16 (5), 277–284.
Guertin, D.P., Goodrich, D.C., Burns, I.S., Sheppard, B.S., Patel, J.B., Clifford, T.J.,
Unkrich, C., Kepner, W.G., Levick, L., 2019. Assessing the Impact of Wildland Fire on
Runoff and Erosion Using the Automated Geospatial Watershed Assessment Tool.
Research and Management for, Working Watersheds and Coastal Systems, p. 38.
Gupta, H., Razavi, S., 2017. Challenges and future outlook of sensitivity analysis.
Sensitivity Analysis in Earth Observation Modelling 397–415.

11

M.B. Meles et al.

Journal of Hydrology 640 (2024) 131621

Šimůnek, 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., G. Brunetti, D. Jacques, M. Th. van Genuchten, and M. Šejna, Recent
developments and applications of the HYDRUS computer software packages since2016, Vadose Zone Journal, 23(1), e20310, 29 p., doi: 10.1002/vzj2.20310, 2024.
Šimůnek, J., van Genuchten, M.T., Šejna, 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.
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, Version 1.0. PC
Progress, Prague, Czech Republic, p. 241.
Singh, V.P., Woolhiser, D.A., 2002. Mathematical modeling of watershed hydrology.
Journal of Hydrologic Engineering 7 (4), 270–292.
Skirvin, S., Kidwell, M., Biedenbender, S., Henley, J.P., King, D., Collins, C.H., Moran, S.,
Weltz, M., 2008. Vegetation data, walnut gulch experimental watershed, Arizona.
United States. Water Resources Research 44 (5).
Smith, R.E., Parlange, J.Y., 1978. A parameter-efficient hydrologic infiltration model.
Water Resources Research 14 (3), 533–538.
Smith, R.E., 1981. Supercritical flow flumes for measuring sediment-laden flow (No. 1655).
US Department of Agriculture, Agricultural Research Service, Tech. Bull. No. 1655,
70 p.
Stone, J.J., Paige, G.B., Hawkins, R.H., 2008. Rainfall intensity-dependent infiltration
rates on rangeland rainfall simulator plots. Transactions of the ASABE 51 (1), 45–53.
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. Advances in Water Resources. 33 (4), 456–467.
https://doi. org/10.1016/j.advwatres.2010.01.010.
Tetzlaff, D., Buttle, J., Carey, S.K., McGuire, K., Laudon, H., Soulsby, C., 2015. Tracerbased assessment of flow paths, storage and runoff generation in northerncatchments: A review. Hydrological Processes 29 (16), 3475–3490.

Therrien, R., Sudicky, E.A., Park, Y.J., McLaren, R.G., 2012. Hydrogeosphere: A threedimensional numerical modelling describing fully-integrated subsurface and surface
flow and transport. Ontario, Canada, Aquanty Inc, Waterloo.
Tonkul, S., Baba, A., Şimşek, C., Durukan, S., Demirkesen, A.C., Tayfur, G., 2019.
Groundwater recharge estimation using HYDRUS 1D model in Alaşehir sub-basin of
Gediz Basin in Turkey. Environmental Monitoring and Assessment 191, 1–19.
van Genuchten, M.T., 1980. A closed-form equation for predicting the hydraulicconductivity of unsaturated soils. Soil Science Society of America Journal 44 (5),892–898. https://doi. org/10.2136/sssaj1980.03615995004400050002x.
Vivoni, E.R., Mascaro, G., Mniszewski, S., Fasel, P., Springer, E.P., Ivanov, V.Y., Bras, R.
L., 2011. Real-world hydrologic assessment of a fully-distributed hydrological modelin a parallel computing environment. Journal of Hydrology 409 (1–2), 483–496.
Vrugt, J.A., Diks, C.G., Gupta, H.V., Bouten, W., Verstraten, J.M., 2005. Improved
treatment of uncertainty in hydrologic modeling: Combining the strengths of global
optimization and data assimilation. Water resources research 41 (1).
Wagner, S., Fersch, B., Yuan, F., Yu, Z., Kunstmann, H., 2016. Fully coupled atmospherichydrological modeling at regional and long-term scales: Development, application,
and analysis of WRF-HMS. Water Resources Research 52 (4), 3187–3211.
Wood, E.F., Roundy, J.K., Troy, T.J., van Beek, L.P.H., Bierkens, M.F.P., Blyth, E., deRoo, A., Doell, P., Ek, M., Famiglietti, J., Gochis, D., van de Giesen, N., Houser, P.,
Jaffe, P., Kollet, S., Lehner, B., Lettenmaier, D.P., Peters-Lidard, C., Sivapalan, M.,
Sheffield, J., Wade, A., Whitehead, P., 2011. Hyperresolution global land surface
modeling: meeting a grand challenge for monitoring Earth’s terrestrial water. Water
Resources Research 47, W05301. https://doi.org/10.1029/2010WR010090.
Woolhiser, D. A., Smith, R. E., and Goodrich, D. C. (1990). “KINEROS—A kinematic
runoff and erosion model: Documentation and user manual.” Rep. No. ARS-77,
USDA, Washington, D.C.
Yatheendradas, S., Wagener, T. and GOODRICH, D., 2005. A distributed real-timesemiarid flash-flood forecasting model utilizing radar data. IAHS-AISH publication,
pp.108-117.
Zhang, Y., Cheng, G., Li, X., Han, X., Wang, L., Li, H., Chang, X., Flerchinger, G.N., 2013.
Coupling of a simultaneous heat and water model with a distributed hydrologicalmodel and evaluation of the combined model in a cold region watershed.
Hydrological Processes 27 (25), 3762–3776.

12