the Creative Commons Attribution 4.0 License.
the Creative Commons Attribution 4.0 License.
Simulator for Hydrologic Unstructured Domains (SHUD v1.0): numerical modeling of watershed hydrology with the finite volume method
Paul A. Ullrich
Christopher J. Duffy
Hydrologic modeling is an essential strategy for understanding and predicting natural flows, particularly where observations are lacking in either space or time or where complex terrain leads to a disconnect in the characteristic time and space scales of overland and groundwater flow. However, significant difficulties remain for the development of efficient and extensible modeling systems that operate robustly across complex regions. This paper introduces the Simulator for Hydrologic Unstructured Domains (SHUD), an integrated, multiprocess, multiscale, flexibletimestep model, in which hydrologic processes are fully coupled using the finite volume method. SHUD integrates overland flow, snow accumulation/melt, evapotranspiration, subsurface flow, groundwater flow, and river routing, thus allowing physical processes in general watersheds to be realistically captured. SHUD incorporates onedimensional unsaturated flow, twodimensional groundwater flow, and a fully connected river channel network with hillslopes supporting overland flow and baseflow.
The paper introduces the design of SHUD, from the conceptual and mathematical description of hydrologic processes in a watershed to the model's computational structures. To demonstrate and validate the model performance, we employ three hydrologic experiments: the Vcatchment experiment, Vauclin's experiment, and a model study of the Cache Creek Watershed in northern California. Ongoing applications of the SHUD model include hydrologic analyses of hillslope to regional scales (1 m^{2} to 10^{6} km^{2}), water resource and stormwater management, and interdisciplinary research for questions in limnology, agriculture, geochemistry, geomorphology, water quality, ecology, climate and landuse change. The strength of SHUD is its flexibility as a scientific and resource evaluation tool where modeling and simulation are required.
The complexity of today's environmental issues, the multidisciplinary nature of scientific and resource management questions, and the diversity and incompleteness of available observational data have all led to the need for models as a means of synthesis. When models are computationally efficient and physically consistent, they become important tools for extrapolation across observations and systems. They help us better understand the physical history of a given system and make decisions about the future in light of socioeconomic, hydrologic, or climatological change. The datasets produced through modeling can assist with decisions for infrastructural planning, water resource management, flood protection, contamination mitigation, and other relevant concerns. Nonetheless, models of varying complexity are available, and the required model complexity (resolution, scale, and coupled states/fluxes) for a given problem depends on the particular research or management purpose, the questions to be answered, and the data available.
Nonetheless, environmental managers, policymakers, and stakeholders have a growing demand for highresolution and detailed information about hydrologic flows at fine temporal–spatial resolution across watersheds. This need reflects the everincreasing importance of detailed longterm predictions and projections for ecological systems, agricultural development, and food security under future climate change. Global climate modeling, typically performed with a general circulation model, also requires information on soil moisture and groundwater fluctuations, which are related to streamflow and reservoir management (Hrachowitz and Clark, 2017; Blöschl et al., 2019).
In hydrology, lumped models (Hawkins et al., 1985; Fleming, 2010; Bergström, 1992) have proven to be fast and stable tools for estimating the streamflow in rivers, requiring simplified meteorological data and limited observed flow data. Lumped models disregard the spatial heterogeneity of terrestrial characteristics while including the basic watershed features (e.g., contributing area, overall relief, average land use, soil conditions, etc.). Their purpose is inputoutput analysis without internal structure (Moradkhani and Sorooshian, 2008), and their parameters may lack precise physical meaning, which makes it challenging to interpret watershed characteristics or transfer parameters to other regions. On the other hand, distributed models (Beven, 2012; Lin et al., 2018; Gochis et al., 2015; Santhi et al., 2006; Liang et al., 1996; Vivoni et al., 2011; Refsgaard et al., 1998; Shen and Phanikumar, 2010) also have their limitations. One challenge for multiprocess distributed models is addressing uncertainties in spatial parameters (soils, hydrogeology, landsurface processes, etc.) and limited predictive skill for large, highresolution catchments, even if the estimated model parameters (e.g., soil properties, surface characteristics, and aquifer properties) and atmospheric inputs have incommensurate resolutions. The latter is particularly important to watershed calibration and has led to a major source of uncertainty (Beven, 2012; Blöschl et al., 2019). Nonetheless, at the continental and global scale, progress has been made in higherresolution elevation data, refined soil surveys (Soil Survey Staff, 2015), satellite land cover (Homer and Fry, 2012), and much higherresolution atmospheric inputs. This has created new opportunities for the development of new hydrologic models that leverage advances in data and mathematical/computational strategies.
Communities have become more sophisticated about their choice of models (Beven, 2012) and tailoring models to explicit purpose or applications, with models ranging from the simplest lumped models (HECHMS, Fleming, 2010, HBV, Bergström, 1992) to semidistributed models (Beven, 1989; Beven and Germann, 1982; Beven and Kirkby, 1979), complex distributed hydrologic models (WRFHydro, Lin et al., 2018; Gochis et al., 2015, PRMS, Leavesley et al., 1983, SWAT, Santhi et al., 2006, VIC, Liang et al., 1996, MIKESHE, Abbott and Refsgaard, 1996; Refsgaard et al., 1998, inHM, VanderKwaak, 1999, tRIBS, Vivoni et al., 2011, 2004, 2005 and PAWS, Shen and Phanikumar, 2010), and even cuttingedge hydrologic models based on machinelearning methods (Rasouli et al., 2012; Petty and Dhingra, 2018; Shen et al., 2018). In each case the choice of model requires the assessment of factors related to prediction variables, performance, flexibility, and availability of data.
The Simulator for Hydrologic Unstructured Domains (SHUD) is a multiprocess, multiscale model where major hydrologic processes are fully coupled using the finite volume method (FVM). SHUD encapsulates the strategy for the synthesis of multistate distributed hydrologic models using the integral representation of the underlying physical process equations and state variables. As an intellectual descendant of Penn State Integrated Hydrologic Model (PIHM), the SHUD model is a continuation of 16 years of PIHM model development in hydrology and related fields since the release of its first PIHM version (Qu, 2004).
The conceptual structure of the twostate integralbalance model for soil moisture and groundwater dynamics was originally devised by Duffy (1996), in which the partial volumes occupied by unsaturated and saturated moisture storage were integrated directly into a local conservation equation. This twostate integralbalance structure simplified the hydrologic dynamics while preserving the natural spatial and temporal scales contributing to runoff response. Brandes et al. (1998) used FEMWATER to capture the inflow/outflow behavior within a hillslope–stream scheme. In 2004, Qu (2004) embedded evapotranspiration and river networks and released the Penn State Integrated Hydrologic Model (PIHM) v1.0, which was the most important milestone of the twostate integralbalance model. Since PIHM v1.0 (Qu, 2004), the PIHM has been a generic hydrologic model applicable to general watersheds or basins. After that, PIHM v2.0 (Kumar et al., 2009; Kumar and Duffy, 2009) enhanced PIHM's land surface modeling component. A GIStool, PIHMgis (Bhatt et al., 2014) and the Essential Terrestrial Variables Data Server (HydroTerre, Leonard and Duffy, 2013) accelerated model deployment and motivated additional applications with PIHM. Because of the sophisticated hydrologic modeling capability and efficient spatial representation in PIHM, various model coupling projects were later initiated. For example, FluxPIHM coupled the NOAH Land Surface Model with PIHM to provide a more detailed calculation of energy balance and evapotranspiration (Shi et al., 2015, 2014). Zhang et al. (2016) coupled landscaping evolution with PIHM (LEPIHM). Bao (Bao, 2016; Bao et al., 2017) added a reactiontransport module to PIHM (RTPIHM, RTFluxPIHM). FluxPIHMBGC (Shi et al., 2018) provided biogeochemistry into FluxPIHM. The MultiModule PIHM (MMPIHM) project (https://github.com/PSUmodeling/MMPIHM, last access: 19 December 2019) was then devised to build a uniform repository for all coupled modules. Still, more PIHM coupling projects are ongoing, targeting sediments, lakes, crops, and other systems. Simultaneously, a Finite volumebased Integrated Hydrologic Modeling (FIHM) was also developed (Kumar et al., 2009), which used secondorder accurate finite volumes and solved for both 2D unsteady overland flow and 3D subsurface flow. Figure 1 shows the family tree of PIHM and SHUD. Within this tree, essentially every revision/branch received crosspollination from others, supporting the ever growing complexity of these systems. The PIHMgis (Bhatt et al., 2014) and rSHUD (Shu, 2020) are GIS, pre and postprocessing tools for PIHM and SHUD, respectively. Although PIHM and SHUD share the same fundamental conceptual twostate integral model, both the inputs and outputs are incompatible. Details of differences between them are summarized in Sect. 4.1 of this paper.
SHUD's design is based on a concise representation of a watershed and river basin's hydrodynamics, which allows for interactions among major physical processes operating simultaneously but with the flexibility to add or drop stateprocessconstitutive relations depending on the objectives of the numerical experiment.
As a distributed hydrologic model, the computational domain of the SHUD model is discretized using an unstructured triangular irregular network (e.g., Delaunay triangles) generated with constraints (geometric and parametric). A local prismatic control volume is formed by the vertical projection of the Delaunay triangles forming each layer of the model. Given a set of constraints (river network, watershed boundary, elevation, and hydraulic properties), an optimized mesh is generated. The optimized mesh allows the underlying coupled hydrologic processes for each finite volume to be calculated numerical efficiently and stability (Farthing and Ogden, 2017; Vanderstraeten and Keunings, 1995; Shewchuk, 1996). River volume cells are also prismatic, with trapezoidal or rectangular cross section, and maintain the topological relation with the Delaunay triangles. The local control volumes encapsulate all equations to be solved and are herein referred to as the model kernel. The R package rSHUD (Shu, 2020) has been developed to help users generate these triangular domains.
The objective of this paper is to introduce the design of SHUD, from the fundamental conceptual model of hydrology to governing hydrologic equations in a watershed to computational structures describing hydrologic processes. Section 2 describes the conceptual design and equations used in the model. In Sect. 3, we employ three hydrologic experiments to demonstrate the simulation and capacity of the model. The three applications presented here are (1) the Vcatchment experiment, (2) the Vauclin experiment (Vauclin et al., 1979), and (3) the Cache Creek Watershed (CCW), a headwater catchment in northern California. Section 4 summarizes the differences between SHUD and PIHM, then proposes possible applications of the SHUD model.
2.1 Conceptual description of hydrologic system
We begin our introduction to the SHUD model with a conceptual description of water movement in a watershed. Figure 2 depicts processes that are incorporated in the model. Catchment hydrology is driven by atmospheric inputs, including rainfall, snowfall, and irrigation. The land surface hydrology first interacts with vegetation through plantcanopy interception and throughfall. When precipitation exceeds the interception capacity, throughfall or excess precipitation falls to the land surface. Snowfall accumulation and melting is an important landsurface process and a major source of soil moisture and recharge in many temperate or mountainous climate regions. SHUD captures the partitioning of excess precipitation into surface runoff and infiltration into the soil.
Vertically, the aquifer is divided into two coupled layers based on its saturation status: the top unsaturated layer (or vadose layer) is constrained to 1D vertical flow, and the saturated groundwater layer admits 2D flow. These layers overlie an impermeable or lowpermeable layer such as bedrock or an effective depth of circulation where deeper flows are small and unlikely to contribute to baseflow. The vertical fluxes within the unsaturated zone include infiltration and exfiltration. Deep percolation or recharge to groundwater is fully coupled to soil moisture dynamics. SHUD accounts for conditions when the groundwater table reaches or exceeds the land surface, decreasing infiltration, allowing exfiltration as local ponding and runoff. The model also accounts for upward capillary flow from a shallow water table depending on soil moisture and vegetation conditions. Lateral (2D) groundwater flow represents the basic mechanism for baseflow to streams or rivers.
Surface runoff or overland flow is generated by excess rainfall and ponding and is represented as 2D shallow water flow in SHUD. Surface runoff complements baseflow runoff as the dominant source of streams or rivers. SHUD also allows reversing flows from the channel to the hillslope as surface inundation or channel losses to groundwater. This may occur when the river stage rises above bankfull storage during flooding events or in an arid region where the river recharges the local groundwater.
Evaporation generally represents the largest water loss from the catchment with four components: evapotranspiration (ET) from interception storage, surface ponding, soil moisture, and shallow groundwater. Transpiration occurs only when vegetation is present and could draw from the saturated groundwater when the groundwater level is high enough. Direct evaporation occurs from interception, ponding water, and soil moisture.
Following the above description, several assumptions and simplifications are made in the SHUD model:

In the default configuration, the watershed boundary is generally handled as a closed domain, in which precipitation and evapotranspiration are the major vertical fluxes, and river flow is the major lateral flow into and out of the domain. It is a common water balance assumption of $\mathrm{\Delta}S=PEQ$, where ΔS, P, E, and Q are total water storage, precipitation, evapotranspiration, and discharge, respectively. Modifications to these conditions can be applied to realize additional flows such as pumping wells, irrigation, basin diversions, and so forth. The model is able to further support boundary conditions as needed for various research questions.

The hydraulic gradient is vertical within the soil column and is controlled by gravity and capillary potential. This assumption is invalid for microscale soil water movement but useful when the model grid spacing ranges from meters to kilometers (Beven, 2012).

The evaporative fluxes that occur due to ET from rivers are assumed to be small and can be approximated by the evapotranspiration from the riparian or hyporheic area of the model where the shallow water table and high soil moisture conditions exist.

The hydrologic characteristics, including all physical parameters in soil, landuse, and terrain, are homogeneous within each cell. This is a typical assumption in distributed models, as these models still need discretized domains.

At present SHUD requires all geographic, topographic, and hydraulic parameters do not change in time.

Finally, SHUD uses a simplified representation of the geometry of the river networks due to the limitation of such data. This assumption is due to the inherent challenges in measuring the geometry of the river cross section everywhere along with the stream network.
2.2 Mathematical structure
The notation used in this section is summarized in the list of symbols in the Appendix.
Figure 3 depicts the geometric structure of the discrete cells in SHUD. The watershed domain is discretized using an irregular unstructured triangular network (Delaunay triangles) generated with imposed spatial constraints by triangle (Shewchuk, 1996). NonDelaunay triangulations are also permitted, so any GIS software and advanced programming language (R, MATLAB, and Python) could potentially generate the triangular domain. A prismatic control volume is formed by the vertical extension of the Delaunay triangles to produce three layers: land surface layer, unsaturated zone, and groundwater layer. The modeler is responsible for defining the aquifer depth (from the land surface to the impervious bedrock) based on measurements or terrestrial characteristics. The thickness of the unsaturated zone (D_{us}) is determined by the difference between the land surface elevation (z_{sf}) and groundwater table (z_{gw}) above datum, i.e., ${D}_{\mathrm{us}}={z}_{\mathrm{sf}}{z}_{\mathrm{gw}}$. When the groundwater table reaches the land surface (z_{gw}>z_{sf}), the thickness of unsaturated zone →0.
Figure 4 depicts the exchange of water between the rivers and hillslope cells. Within each river channel, there are two longitudinal fluxes and two lateral fluxes: upstream (Q_{up}), downstream (Q_{dn}), overland (Q_{sf}), and groundwater (Q_{sub}).
The hydrologic model uses the method of moments to reduce the partial differential equations (PDEs) into ordinary differential equations (ODEs) and solve the ODE system using a globally implicit numerical solver. The state variables include water height on the land surface (Y_{sf}), soil moisture storage (Y_{us}), groundwater depth (Y_{gw}), and river stage (Y_{riv}). The initial value problem for these ODEs is formulated as
where the discrete state vector is denoted by Y,
Y_{0} contains the initial conditions, and f(t,Y) denotes the equations governing the hydrologic flow, which are described in this section.
The system of ODEs describing the hydrologic processes are fully coupled and solved simultaneously at each time step ($\mathrm{\Delta}t={t}_{n}{t}_{n\mathrm{1}}$) using CVODE, a stiff solver based on Newton–Krylov iteration (Hindmarsh et al., 2019). In brief, the CVODE solver calculates Y(t_{n}), given Y(t_{n−1}) and $\frac{\mathrm{d}\mathit{Y}}{\mathrm{d}t}{\mathrm{}}_{{t}_{n\mathrm{1}}}$. The technical description of the CVODE solver can be found in the literature (Hindmarsh et al., 2019, 2005; Cohen and Hindmarsh, 1996). The governing equations in SHUD are provided in Table 1.
Figure 5 depicts the workflow within the SHUD model. The explicit model time step (MTS) $\mathrm{\Delta}t={t}_{n}{t}_{n\mathrm{1}}$ is user specified, typically varying from 1 min to 1 h. Within the MTS, the relevant gradient law determines the fluxes and hence all state values.
The fluxes of ET and interception change relatively slowly within short periods (such as 1 h) so that full coupling of ET with soil water is not necessary for this model. Instead, the interception, ET, and snow calculations are solved explicitly at the MTS, while the calculation of Y_{sf}, Y_{us}, Y_{gw}, and Y_{riv} uses the implicit time step (ITS).
The CVODE solver determines the ITS automatically based on both the specified tolerances and the error function of Y and dY. The initial ITS equates to the explicit MTS. Within the ITS, dYs is calculated based on Y from the last MTS. If the CVODE solver converges with the current value of the ITS, it returns the updated Y. Otherwise, a convergence failure occurs that forces an ITS reduction.
The mathematical model underlying SHUD consists of five components: vegetation and evapotranspiration, land surface, unsaturated layer, saturated layer, and river channel. These are described in the following subsections.
2.2.1 Vegetation and evapotranspiration
Interception refers to the direct water loss of precipitation when vegetation cover exists and is treated as a simple storage bucket – namely, precipitation cannot reach the land surface until the interception storage capacity is satisfied. The capacity of this storage is the maximum interception volume, which is a function of the vegetation leaf area index (LAI) and satisfies equation ${S}_{\mathrm{ic}}^{*}={C}_{\mathrm{ic}}\mathrm{LAI}$, where LAI represents the coverage of vegetation canopy over the land area (area of leaves over the area of land, in m^{2} m^{−2}), and C_{ic} is interception coefficient (in m). The default C_{ic} is taken to be 0.2 kg m^{−2} as suggested in Dickinson (1984).
The interception is equal to the deficit of interception – the difference between interception capacity (${S}_{\mathrm{ic}}^{*}$) and existing interception storage (S_{ic}). If precipitation is less than the deficit, interception is equal to the precipitation rate (see Fig. 6).
Potential evapotranspiration (PET) is the quantity of water that would evaporate and transpire from an ideal surface if extensive free water was available to meet the demand (Maidment, 1993; Kirkham, 2014). As such, PET is a practical and rapid estimation of water flux from land to atmosphere. The PET (E_{0}) is governed by Penman–Monteith equation (Penman, 1948):
Here we do not elaborate on this equation, as it is common among different hydrologic models (Allen, 1998; Maidment, 1993). At each ET step, the model calculates PET in terms of the prescribed forcing data. PET values are conditioned on the parameters from various land cover types, including factors for albedo, LAI, and roughness length.
The total actual evapotranspiration (AET) consists of three parts: evaporation from interception (E_{c}), transpiration from vegetation canopy (E_{t}), and direct evaporation of soil (E_{s}). The calculation of AET for these three components follows from the equations below:
Here, E_{c} is controlled by PET from interception storage. Both E_{s} and E_{t} are affected by soil water stress (β_{s}) and impervious area fraction (α_{imp}). Impervious area is also considered a barrier of evapotranspiration in the model. E_{s} is referred to as the demand water evaporation from soil and emerges from two sources, namely the evaporation from ponding water (E_{sp}) and evaporation from soil moisture (E_{sm}), ${E}_{\mathrm{s}}={E}_{\mathrm{sp}}+{E}_{\mathrm{sm}}$. Ponded water has higher priority to evaporate and direct evaporation only uses the water in the surface when ponding water is able to meet the E_{s} demand, i.e., ${y}_{\mathrm{sf}}/\mathrm{\Delta}t>{E}_{\mathrm{s}}$. When ponding water is insufficient to meet E_{s}, soil water balances the difference between demand and available water in the surface; when ponding water does not exist, direct evaporation extracts water from the soil profile (${E}_{\mathrm{sm}}={E}_{\mathrm{s}},{E}_{\mathrm{sp}}=\mathrm{0}$):
Transpiration also has two potential sources: soil moisture and groundwater from the groundwater table and root depth for the landuse class. Once the groundwater table is higher than the root zone depth, vegetation uses groundwater, and soil moisture stress for transpiration is equal to one (β_{s}=1). Water balance associated with snow accumulation is quantified via
Snow melt rate is determined by the snow melting factor (m_{f}), air temperature (T), and temperature threshold (T_{0}) at which snow melt occurs. This formulation is often referred to as the degreeday method, in which the values of the snow melting factor and temperature threshold are empirical (Maidment, 1993; Beven, 2012). The water from snowmelt is considered as a direct water contribution to the land surface.
2.2.2 Water on the land surface
Water balance on the land surface is given by
The water balance of net precipitation (P_{n}), infiltration (q_{i}), evaporation from the ponding layer (E_{sp}), and horizontal overland flow (Q_{j}) determine the storage of water on the land surface. Net precipitation (P_{n}) is the total residual water after adjusting for rainfall/snow interception and snowmelt. The overland flow ${Q}_{\mathrm{s}}^{j}$ in direction j is calculated with Manning's equation (13). Here $\stackrel{\mathrm{\u203e}}{{y}_{\mathrm{sf}}}$ is effective water height, determined by the gradient between two cells,
Infiltration is estimated using Richards equation:
The infiltration rate is a function of soil saturation ratio (Θ), soil properties (k_{x}, k_{m} α, β, and α_{h}), and ponding water height (existing ponding water plus precipitation/irrigation). Infiltration occurs in the top soil layer (D_{inf}), and the infiltration rate results from the ponding water height and soil moisture. The default value of D_{inf} is 10 cm, which is adjustable in calibration files. The application rate y_{s}∕Δt combines ponding water, irrigation, and precipitation together, and that determines the hydraulic gradient applied on the top soil layer. Finally, K_{max} is the infiltration capacity determined by both soil matrix and macropore characteristics. When application rate is less than the maximum infiltration capacity, the infiltration is controlled by soil matrix flow; when application rate is larger than K_{max}, effective conductivity is a function of the soil matrix and macropores (Chen and Wagenet, 1992). The infiltration equation takes the macropore effect into account, so the algorithm allows faster infiltration under heavy rainfall events and enables the soil to hold water for vegetation under dry conditions.
2.2.3 Unsaturated zone
As discussed above, the horizontal flow in the vadose zone is neglected compared to the dominant vertical flow. There are three processes controlling the water in vadose zone: infiltration (q_{i}), ET in soil moisture (E_{sm}), and recharge to groundwater (q_{r}). The calculation of infiltration and ET is explained in the previous subsection. Recharge to groundwater is calculated with Eq. (16). The soil moisture content to field capacity controls the recharge rate.
Because of the simplifications required for the twolayer description of the vertical aquifer profile, SHUD uses the relationship between soil moisture and field capacity as the gradient to drive the recharge instead of the hydraulic gradient. K_{er} is the effective conductivity for recharge and is equal to the arithmetic mean of the conductivity of the unsaturated zone and saturated zone.
When the bottom of the vegetation root zone is below the groundwater table, then E_{tg}>0 and vegetation extracts water from the saturated zone. Otherwise E_{tg}=0, meaning that transpiration uses soil moisture. When ponding water exists on the land surface, direct evaporation extracts water from ponding water first; when ponding water is depleted via evaporation, then the remainder of evaporation (E_{sm}) uses water from soil moisture based on the water stress.
2.2.4 Groundwater
The water balance of groundwater is controlled by Eq. (18).
The calculation of horizontal groundwater flow uses the Boussinesq equation. When the bottom of the root zone is lower than the groundwater table, then E_{t}>0, otherwise, E_{t}=0.
The horizontal groundwater flux ${Q}_{\mathrm{g}}^{j}$ is determined by the Darcy equation and Dupuit–Forchheimer assumption. Above z_{b} is the elevation of impervious bedrock; ${z}_{\mathrm{b}}^{j}$ is the bedrock elevation of its jth neighbor cell, and d_{j} is distance between the centroids of two adjacent cells, so the gradient between the two cells is $\left[\left({y}_{\mathrm{gw}}+{z}_{\mathrm{b}}\right)\left({y}_{\mathrm{gw}}^{j}+{z}_{\mathrm{b}}^{j}\right)\right]{d}_{j}^{\mathrm{1}}$. The effective conductivity for the groundwater flow is the mean value of the effective horizontal conductivity over the two cells. The crosssectional area along the groundwater flux is equal to ${L}_{j}\times \stackrel{\mathrm{\u203e}}{{y}_{\mathrm{gw}}}$.
In Eq. (20), the effective horizontal conductivity (K_{eg}) is a function of the groundwater table and characteristics of the macropores. The calculation of effective horizontal conductivity of each cell is given by
where k_{g} and k_{m} are the saturated hydraulic conductivity of soil matrix and macropores, z_{m}, z_{gw}, and z_{cb} are elevations of macropore, groundwater table, and bedrock, and α_{v} is the vertical areal macropore fraction (in m^{2} m^{−2}).
The effective horizontal conductivity captures the effect of spatially varying conductivity on the saturated flow when the groundwater level rises (Jiang et al., 2009; Bobo et al., 2012; Chen et al., 2018; Cheema, 2015; Taylor, 1960; Lin et al., 2007). Figure 7 reveals the effective horizontal conductivity changes along with different groundwater levels. When the groundwater table is below the level of the macropores, K_{eg} is equal to the saturated conductivity. When the groundwater level is above the macropore level, the effective conductivity increases with the groundwater level, taking into consideration the conductivity and area fraction of macropores in the soil profile. The effective conductivity reaches its maximum value once the groundwater table level reaches the land surface.
2.2.5 Water in streams
The water balance in river channels is described by
The mass balance in each river channel consists of four parts: ${Q}_{\mathrm{s}}^{j}$, the overland flow from cells (1 to N_{c} cells) that intersect with the river channel; ${Q}_{\mathrm{g}}^{j}$, the lateral groundwater flux from intersection with the jth cell; ${Q}_{\mathrm{up}}^{j}$, the longitudinal flow from upstream channels; and Q_{dn}, the flux to the downstream channel. N_{u} is the number of upstream channels; in the model, the number of upstream channels is nonnegative, but only one downstream channel is permitted. SHUD assumes river channels can converge to a single downstream channel. The convergence rule does not affect the topological relationship between river channels and cells.
The topological relationship between cells and river channels is shown in Fig. 4b. As depicted, the river consists of a series of river reaches which intersect prismatic elements. A single reach is a sorted collection of multiple river segments, while each segment lies within hillslope cells. Surface and groundwater exchanges occur between each segment and the underlying cell. The sum of overland flow from multiple cells contributes to the net storage in a river reach.
The downstream channel flux Q_{dn} is based on the 1D diffusive wave equation that is simplified as Manning's equation for open channel:
where A_{cs} is the crosssection area of the river reach, and $\stackrel{\mathrm{\u203e}}{P}$ and $\stackrel{\mathrm{\u203e}}{{s}_{\mathrm{0}}}$ are average wet perimeter and average slope of a river reach and its downstream reach.
The upstream flux Q_{up} is equal to the sum of Q_{dn} from the multiple upstream reaches. The water balance equation in the river channel neglects evaporation and precipitation because the area of open water in the watershed is relatively small, and the area of open water is already included in precomputation for the cells. Therefore, the channel routing represents the water exchange between the river and hillslope and takes the overland flow and baseflow into account.
The overland flow between river segment and associated hillslope cell (Q_{sr}) is calculated as follows:
Here z_{bank} is the elevation of riverbank or levee (Fig. 4c), implying that the relative height of the land surface or river stage controls the direction of water exchange between the land surface and river segment.
The groundwater exchange between river segment and hillslope cell is described by Q_{gr}, which is calculated as
In this section, we present the results of using SHUD for three applications: first, we use the Vcatchment experiment (Shu, 2019e) to validate the calculation of overland flow and river routing in an idealized catchment; second, we use Vauclin's experiment (Vauclin et al., 1979) to assess the calculation of infiltration, unsaturated flow in the vadose zone and horizontal saturated flow; finally, we apply the model to a hydrologic simulation in the Cache Creek Watershed (Shu, 2019a), a headwater catchment in the Sacramento Watershed of northern California.
3.1 Vcatchment
The Vcatchment (VC) experiment is a standard test case for numerical hydrologic models to validate their performance for overland flow along a hillslope and in the presence of a river channel (Shen and Phanikumar, 2010). The VC domain consists of two inclined planes draining into a sloping channel (Fig. 8). Both hillslopes are 800 m×1000 m with Manning's roughness n=0.015. The river channel between the hillslopes is 20 m wide and 1000 m in length with n=0.15. The slope from the ridge to the river channel is 0.05 (in the x direction), and the longitudinal slope (in the y direction) is 0.02.
Rainfall in the VC begins at time zero at a constant rate of 18 mm h^{−1} and stops after 90 min, producing 27 mm of accumulated precipitation. Since this experiment excludes the evaporation and infiltration, the total outflow from lateral boundaries and the river outlet must be the same as the total precipitation (following conservation of mass).
Figure 9 illustrates the discharge from the side plane to the river channel and at the river outlet. The specific discharge (the volume discharge divided by the total area of the catchment) increases with precipitation until it reaches the maximum discharge rate, which is equal to the precipitation rate. Discharges along lateral boundaries and from the river outlet reach the maximum discharge rate but at different times; namely, the discharge rate from the side plane reaches the maximum value earlier than in the river outlet. The dots are discharge digitalized from Shen and Phanikumar (2010) with WebPlotDigitizer (https://automeris.io/WebPlotDigitizer/, last access: 19 December 2019). The results suggest SHUD can correctly capture the processes in overland flow and channel routing, although flow from the river outlet occurs earlier than the prediction in Shen and Phanikumar (2010). Both the fluxes from side plane and outlet meet the maximum flow rate, which is the same magnitude of precipitation after a short period of rainfall. The flux rates start decreasing after precipitation stops. The accumulated volume of flux confirms the correct mass balance of both fluxes.
To verify correct performance of the numerical method, we calculate the bias of mass balance in the model and assess the differences among input, output, and storage change in the system (Eq. 35). For this simulation the bias ends up being ∼0.2 %.
3.2 Vauclin's experiment
Vauclin's laboratory experiment (Vauclin et al., 1979) is to assess groundwater table change and soil moisture in the unsaturated layer under precipitation or irrigation. The experiment was conducted in a sandbox with dimensions of 3 m long ×2 m deep ×0.05 m wide (see Fig. 10). The box was full of uniform sand particles with measured hydraulic parameters: the saturated hydraulic conductivity was 35 cm h^{−1}, and porosity was 0.33 m^{3} m^{−3}. The left and bottom of the sandbox were impervious layers, and the top and the right side were open. The hydraulic head was 0.65 m constantly on the right boundary. Constant irrigation (1.48 cm h^{−1}) was applied over the first 50 cm of the top left of the sandbox, while the rest of the top was covered to avoid water loss via evaporation.
The experiment's initial condition is an equilibrium water table with a uniform hydraulic head. Irrigation was initiated at t=0. The groundwater table was then measured at 2, 4, 6, and 8 h at several locations along the length of the box. (Vauclin et al., 1979) also uses a 2D (vertical and horizontal) numeric model to simulate the soil moisture and groundwater table. The maximum bias between measurement and simulation was 0.52 m, according to the digitalized value of Vauclin et al., 1979, Fig. 10.
Besides the parameters specified in Vauclin et al. (1979), additional information is needed in the SHUD model, including the α and β in the van Genuchten equation and residual water content (θ_{r}). Therefore, we use a calibration tool to estimate the representative values of these parameters. The use of calibration in this simulation is reasonable because the model, inevitably, simplifies the real hydraulic processes. The calibration thus nudges the parameters to representative values that approach or fit the true natural processes. The calibrated values are θ_{r}=0.001 m^{3} m^{−3}, α=0.3 ,and β=5.2. The simulated results in our modeling and literature (Vauclin et al., 1979; Shen and Phanikumar, 2010) show a modest error between the simulations and measurements.
This error is likely due to (1) the need for a detailed aquifer layer description of soil layers or (2) possible invalidity of the vertical and horizontal flow assumptions in the SHUD model.
SHUD simulated the groundwater table at all four measurement points (see Fig. 10b). The maximum bias between simulation and Vauclin's observations is 5.5 cm, with R^{2}=0.99, which is comparable to the bias 5.2 cm of numerical simulation in Vauclin et al. (1979). By adding more parameters into the calibration, the bias in the simulation decreased to 3 cm. Indeed, as discussed earlier, the simplifications employed by SHUD for the unsaturated and saturated zone mainly benefits computational efficiency while limiting applicability where finescale soil profile information is required.
These simulations, compared against Vauclin's experiment, generally validate the algorithm for infiltration, recharge, and lateral groundwater flow. A more reliable vertical flow within the unsaturated layer requires multiple layers, which is planned in the next version of SHUD.
3.3 Cache Creek Watershed
The Cache Creek Watershed (CCW) (Shu, 2019a) is a headwater catchment with area 196.4 km^{2} in the Sacramento Watershed in northern California (Fig. 11a, b and c). The elevation ranges from 450 to 1800 m, with a 0.38 m m^{−1} average slope, where the steep implies a particular challenge for numeric models.
Based on NLDAS2 (Xia et al., 2012) from 2000 to 2017, the mean temperature and precipitation were 12.8 ^{∘}C and ∼817 mm, respectively, in this catchment. Precipitation varies seasonally; winter and spring are the primary wet seasons (Fig. 12).
(McKay et al., 2012)(U.S. Geological Survey, 2016)(Soil Survey Staff, 2015)(Homer and Fry, 2012)(Xia et al., 2012)Table 2 lists the geospatial and forcing data supporting hydrologic modeling in the CCW. The DEM is 30 m resolution raster data from National Elevation Dataset (NED) (U.S. Geological Survey, 2016). Forcing data, including precipitation, temperature, relative humidity, wind speed, and net radiation, are from NLDAS2 (Xia et al., 2012, https://ldas.gsfc.nasa.gov/nldas/v2/forcing, last access: 19 December 2019). Our simulation in the CCW (Shu, 2019a) covers the period from 2000 to 2007. Because of the Mediterranean climate in this region, the simulation starts in summer to ensure adequate time before the October start to the water year. In our experiment, the first year (1 July 2000 to 30 June 2001) is the spinup period, the following 2 years (1 July 2001 to 30 June 2003) are the calibration period, and the period from 1 July 2003 to 1 July 2007 is for validation.
The unstructured domain of the CCW (Fig. 11d) is built with rSHUD (Shu, 2020). The domain consists of 1147 triangular cells, with a mean area of 0.17 km^{2}. The total length of the river network is 126.5 km and consists of 103 river reaches in which the highest order of stream is 4. With a calibrated parameter set, the SHUD model (Shu, 2019a) took 5 h to simulate 18 years (2000–2017) in the CCW, with a nonparallel configuration (OpenMP is disabled on Mac Pro 2013 Xeon 2.7 GHz, 32 GB RAM).
Figure 13 reveals the comparison of simulated discharge against the observed discharge at the gage station of USGS 11451100 (https://waterdata.usgs.gov/ca/nwis/uv/?site_no=11451100, last access: 19 December 2019). The calibration procedure exploits the covariance matrix adaptation–evolution strategy (CMAES) to calibrate automatically (Hansen, 2016). The calibration program assigns 72 children in each generation and keeps the best child as the seed for the next generation, with limited perturbations. The perturbation for the next generation is generated from the covariance matrix of the previous generation. After 23 generations, the calibration tool identifies a locally optimal parameter set.
Within the calibration period, Nash–Sutcliffe efficiency (NSE; Nash and Sutcliffe, 1970), KlingGupta efficiency (KGE; Gupta et al., 2009), and R^{2} are 0.72, 0.83, and 0.72, respectively (Fig. 13). The goodnessoffit in the validation period is less than calibration period (as expected), with NSE =0.66, KGE =0.67, and R^{2}=0.65. Although the SHUD model captures the flood peaks after rainfall events, the magnitude of high flow in the hydrograph is less than the gage data. There are two potential causes of this bias: (1) underestimated precipitation intensity from NLDAS2 data or (2) overfitting in the calibration, as the NSE tends to capture the mean value of the observational data rather than the extremes.
Figure 14 represents the monthly water balance in CCW, in which the PET is 3 times the annual precipitation, but the actual evapotranspiration (AET) is only 27 % of the precipitation. This result emerges because the summer is the peak of PET, while winter is the peak of precipitation and water availability. The AET is subjected to PET and water availability, so the maximum of AET occurs in early summer. The runoff ratio is about 73 %.
We use the groundwater distribution (Fig. 15) to demonstrate the spatial distribution of hydrologic metrics calculated from the SHUD model. Figure 15 illustrates the annual mean groundwater table in the validation period. Because the model fixes a 30 m aquifer, the results represent the groundwater within this 30 m aquifer only. The groundwater table and elevation along the green line on the map are extracted and plotted in Figure 15b. The gray ribbon is the 30 m aquifer, and the blue line is the location where groundwater storage is larger than zero. The light blue polygons (with the right axis for scale) are the groundwater storage along the cross section. The groundwater levels tend to follow the terrain, with groundwater accumulated in the valley or along relatively flat flood plains. In the CCW, groundwater does not stay on steep slopes suggesting the high conductivity of upland slopes.
The features of the SHUD model are as follows:

SHUD is a physically based process, spatially distributed catchment model. The model applies national geospatial data resources to simulate surface and subsurface flow in gaged or ungaged catchments. SHUD represents the spatial heterogeneity that influences the hydrology of the region based on national soil data and superficial geology. Several other groups have used PIHM, a SHUD ancestor to couple processes from biochemistry, reaction transport, landscape, geomorphology, limnology, and other related research areas.

SHUD is a fully coupled hydrologic model, where the conservative hydrologic fluxes are calculated within the same time step. The state variables are the height of ponding water on the land surface, soil moisture, groundwater level, and river stage, while fluxes are infiltration, overland flow, groundwater recharge, lateral groundwater flow, river discharge, and exchange between river and hillslope cells.

The global ODE system in SHUD is solved with a stateoftheart parallel ODE solver, known as CVODE (Hindmarsh et al., 2005) developed at Lawrence Livermore National Laboratory.

SHUD permits adaptable temporal and spatial resolution. The spatial resolution of the model varies from centimeters to kilometers based on modeling requirements computing resources. The internal time step of the iteration is adjustable and adaptive; it can export the status of a catchment at time intervals from minutes to days. The flexible spatial and temporal resolution of the model makes it valuable for coupling with other systems.

SHUD can estimate either a longterm hydrologic yield or a singleevent flood.

SHUD is an opensource model, available on GitHub (Shu, 2019b).
4.1 Differences from PIHM
As a descendant of PIHM, SHUD inherits many fundamental ideas and the conceptual structure from PIHM, including the solution of hydrologic variables using CVODE. The code has been completely rewritten in a new programming language, with a new discretization and corresponding improvements to the underlying algorithms, adapting new mathematical schemes and flexible input/output data formats. Although SHUD is forked from PIHM, SHUD still inherits the use of CVODE for solving the ODE system but modernizes and extends PIHM's technical and scientific capabilities. The major differences are as follows:

SHUD is written in C++, an objectoriented programming language with functionality to avoid memory leaks from C. Most of the functions in SHUD do not exist in PIHM, which are newly defined because of the brandnew design of the SHUD model. Only a few functions related to physical equations, such as Manning's equation and van Genuchten's equation, are shared in SHUD and PIHM. So although SHUD and PIHM follow a similar fundamental perceptual model, they follow different mathematical and computational strategies.

SHUD implements a redesign of the calculation of water exchange between hillslope and river. The PIHM defines the river channel as adjacent to bank cells – namely, the river channel shares the edges with bank cells. This design leads to sink problems in cells that share one node with a starting river channel that can reduce simulation performance.

Although the mathematical equations underlying SHUD are generally the same as PIHM, the numerical formulation of processes, the coupling strategy, and input/output structures have been greatly enhanced. Common computations in different functions within various processes are extracted and defined as shared inline functions, maintaining calculation consistency and facilitating code updates. The elimination of redundant variables and functions also advances consistency and efficiency.

SHUD adds massbalance control within the calculation of each layer of cells and river channels, important for accurate and efficient longterm and microscale hydrologic modeling.

The internal data structure and external input/output formats have been redesigned for efficiency and userfriendly formats supporting ASCII and binary. The binary format is particularly important for efficient writing and postprocessing over numerous model domains.
We now briefly summarize the technical model improvements and technical capabilities of the model compared to PIHM. This elaboration of the relevant technical features aims to assist future developers and advanced users with model coupling. Compared with PIHM, SHUD (Shu, 2019b) does the following:

supports compatibility with the implicit Sundial/CVODE solver version 5.0 (the most recent version at the time of writing) and above;

supports OpenMP parallel computation;

utilizes objectoriented programming (C++);

supports humanreadable input/output files and filenames;

exposes unified functions to handle the timeseries data (TSD) (in a standardized spreadsheet format), including forcing, leaf area index, roughness length, boundary conditions, and melting factor;

exports model initial condition at specific intervals that facilitate warm starts of continuous simulation;

automatically checks the range of physical parameters and forcing data;

adds a debug mode that monitors potential errors in parameters and memory operations; and

includes a series of R codes for pre and postprocessing data, visualization, and data analysis (that will be discussed in future work).
The Simulator for Hydrologic Unstructured Domains (SHUD) is a multiprocess, multiscale, and multitemporal model that integrates major hydrologic processes and solves the physical equations with the finite volume method. The governing equations are solved within an unstructured mesh domain consisting of triangular cells. The variables used for the surface, vadose layer, groundwater, and river routing are fully coupled together with a fine time step. The SHUD uses the 1D unsaturated flow and 2D groundwater flow. River channels connect with hillslope via overland flow and baseflow. The model, while using distributed terrestrial characteristics (from climate, land use, soil, and geology) and preserving their heterogeneity, supports efficient performance through parallel computation.
SHUD is a robust integrated modeling system that has the potential for providing scientists with new insights into their domains of interest and will benefit the development of coupling approaches and architectures that can incorporate scientific principles. The SHUD modeling system can be used for applications in (1) hydrologic studies from hillslope scale to regional scale; area of the model domain ranges from 1 m^{2} to 10^{6} km^{2}; (2) water resource and stormwater management; (3) coupling research with related fields, such as limnology, agriculture, geochemistry, geomorphology, water quality, and ecology; (4) climate change; and (5) landuse change. In summary, SHUD is a valuable scientific tool for any modeling task associating with hydrologic responses.
Evapotranspiration Calculation  
Δ  Slope vapor pressure curve (kPa ^{∘}C^{−1}) 
γ  Psychrometric constant (kPa ^{∘}C^{−1}) 
λ  Latent heat of vaporization (MJ kg^{−1}) 
ρ_{a}  Density of Air (kg m^{−3}) 
c_{p}  Specific heat at constant pressure (MJ kg^{−1} ^{∘}C^{−1}) 
e_{a}  Actual vapor pressure (kPa) 
e_{s}  Saturation vapor pressure (kPa) 
G  Soil heat flux density (MJ m^{−2} s^{−1}) 
R_{n}  Net radiation at the crop surface (MJ m^{−2} s^{−1}) 
r_{s}  Surface resistance of vegetation (s m^{−1}) 
r_{a}  Aerodynamic resistance (s m^{−1}) 
Hydrologic metrics  
α  van Genuchten soil parameter (m^{−1}) 
α_{h}  Horizontal macropore areal fraction (m^{2} m^{−2}) 
α_{imp}  Impervious area fraction (m^{2} m^{−2}) 
α_{v}  Vertical macropore areal fraction (m^{2} m^{−2}) 
α_{veg}  Vegetation fraction (m^{2} m^{−2}) 
$\overline{K}$  Average conductivity (m s^{−1}) 
$\overline{y}$  Effective height of overland flow between two adjacent cells (m) 
β  van Genuchten soil parameter (–) 
β_{s}  Soil moisture stress to evapotranspiration (–) 
Δt  Time interval between consequential time steps (m) 
$\stackrel{\mathrm{\u203e}}{{y}_{\mathrm{gw}}}$  Effective water height for groundwater flow calculation (m) 
$\stackrel{\mathrm{\u203e}}{{y}_{\mathrm{sf}}}$  Effective water height for overland flow calculation (m) 
ψ  Soil matrix potential head (m) 
θ  Soil moisture content (m^{3} m^{−3}) 
Θ  Relative saturation ratio (–) 
θ_{fc}  The soil moisture content of field capacity (m^{3} m^{−3}) 
θ_{r}  Residual soil moisture content (m^{3} m^{−3}) 
θ_{s}  Porosity of soil (m^{3} m^{−3}) 
A_{r}  Area of river open water (m^{2}) 
A_{c}  Area of a cell (m^{2}) 
b_{g}  Effective height of groundwater flow between the river segment and hillslope cell (m) 
b_{s}  Effective height of overland flow between the river segment and hillslope cell (m) 
C_{w}  Coefficient of discharge (m) 
C_{ic}  Coefficient of interception (m) 
D_{us}  The deficit of soil column; thickness of vadose layer (m) 
d_{j}  Distance between centroids of the current cell and neighbor j (m) 
d_{rb}  Thickness of river bed; for calculation of baseflow to rivers (m) 
E_{0}  Potential evapotranspiration (m s^{−1}) 
E_{c}  Evaporation from interception (m s^{−1}) 
E_{s}  Evaporation from soil (m s^{−1}) 
E_{sm}  Evaporation from the soil matrix (m s^{−1}) 
E_{sp}  Evaporation from ponding water on land surface (m s^{−1}) 
E_{t}  Transpiration (m s^{−1}) 
E_{tg}  Transpiration from saturated layer (m s^{−1}) 
E_{tg}  Transpiration from saturated layer (m s^{−1}) 
H_{cgw}  Hydraulic head of water in cell groundwater (m) 
H_{csf}  Hydraulic head of water on land surface (m) 
H_{riv}  Hydraulic head in a river channel (m) 
K_{ei}(Θ)  Effective infiltration conductivity (m s^{−1}) 
K_{er}  Effective recharge conductivity (m s^{−1}) 
K_{eg}  Effective horizontal conductivity (m s^{−1}) 
k_{g}  Saturated horizontal conductivity (m s^{−1}) 
k_{v}  Saturated vertical conductivity of saturated layer (m s^{−1}) 
k_{m}  Saturated conductivity of soil macropore (m s^{−1}) 
K_{r}(Θ)  Relative conductivity, which is a function of saturation ratio (–) 
k_{rb}  Saturated conductivity of the river bed (–) 
k_{x}  Saturated conductivity of the top soil (m s^{−1}) 
L_{j}  Length of edge j of a cell (m) 
L_{s}  Length of river segment that overlay with a cell (m) 
LAI  Leaf area index (m^{2} m^{−2}) 
m_{f}  Snow melting factor (m s^{−1} ^{∘}C^{−1}) 
n  Manning's roughness (s m${}^{\mathrm{1}/\mathrm{3}}$) 
N_{c}  Number of cells overlaying a river reach (–) 
N_{u}  Number of upstream reaches flowing to a river reach (–) 
P  Atmospheric precipitation or irrigation (m s^{−1}) 
P_{n}  Net precipitation (m s^{−1}) 
P_{sn}  Snowfall (m s^{−1}) 
Q_{dn}  Volume flux to the downstream river channel (m^{3} s^{−1}) 
Q_{g}  Groundwater flow between two cells (m^{3} s^{−1}) 
Q_{gr}  Volume flux between river and cells via groundwater flow (m^{3} s^{−1}) 
q_{i}  Infiltration rate, positive is downward (m s^{−1}) 
q_{r}  Recharge rate, positive is downward (m s^{−1}) 
Q_{s}  The overland flow between two cells (m^{3} s^{−1}) 
q_{sn}  Snow melting rate (m s^{−1}) 
Q_{sr}  Volume flux between river and hillslope cells via overland flow (m^{3} s^{−1}) 
Q_{up}  Volume flux from upstream river reaches (m^{3} s^{−1}) 
s_{0}  The slope of land surface (m m^{−1}) 
S_{ic}  Water storage of interception layer (canopy) (m) 
${S}_{\mathrm{ic}}^{*}$  Maximum interception capacity (m) 
S_{sn}  Snow storage (m) 
s_{y}  Specific yield (m m^{−1}) 
T  Air temperature (^{∘}C) 
y_{gw}  Groundwater head (above impervious bedrock) of a cell (m) 
y_{riv}  River stage in a river channel (m) 
y_{sf}  Surface water storage in a cell (m) 
y_{us}  Unsaturated storage equivalence of a cell (m) 
z_{b}  Elevation of impervious bedrock from the datum (m) 
z_{bank}  Elevation of the riverbank from the datum (m) 
z_{gw}  Elevation of groundwater table from the datum (m) 
z_{m}  Elevation of macropore from the datum (m) 
z_{sf}  Elevation of land surface from the datum (m) 
T_{0}  Temperature threshold for snow melt to occur (^{∘}C) 
Variables used in CVODE  
Y_{0}  The initial conditions to start the simulation (m) 
Y_{gw}  Vector of cell groundwater head (above impervious bedrock) (m) 
Y_{riv}  Vector of river stage in all river channels (m) 
Y_{sf}  Vector of surface water storage of all cells (m) 
Y_{us}  Vector of unsaturated storage equivalence of all cells (m) 
Y  Vector of conserved state variables in CVODE (m) 
t  Time (s) 
t_{n−1}  Previous time (s) 
t_{n}  Current time (s) 
The source code of the SHUD model (Shu, 2019b) is kept and updated at https://github.com/SHUDSystem/SHUD. The code and data used for this page are archived at Zenodo as follows: SHUD model (https://doi.org/10.5281/zenodo.3561293, Shu, 2019b), User manual (https://doi.org/10.5281/zenodo.3561295, Shu, 2019c), rSHUD (https://doi.org/10.5281/zenodo.3758097, Shu, 2020), Vcatchment (https://doi.org/10.5281/zenodo.3566022, Shu, 2019e), Vauclin (1979) experiment (https://doi.org/10.5281/zenodo.3566020, Shu, 2019d), and Cache Creek Watershed (https://doi.org/10.5281/zenodo.3566034, Shu, 2019a).
LS contributed to conceptualization, investigation, methodology, software, validation, visualization, writing of the original draft, and editing. PU contributed to supervision, investigation, writing of the original draft, and editing. CD contributed to supervision, investigation, writing of the original draft, and editing.
Paul Ullrich is a member of the editorial board of the journal.
Author Lele Shu was supported by the California Energy Commission grant “Advanced StatisticalDynamical Downscaling Methods and Products for California Electrical System” project (award no. EPC16063). Coauthor Paul Ullrich was supported by the U.S. Department of Energy Regional and Global Climate Modeling Program (RGCM) “An Integrated Evaluation of the Simulated Hydroclimate System of the Continental US” project (award no. DESC0016605) and the National Institute of Food and Agriculture, U.S. Department of Agriculture, California Agricultural Experiment Station hatch project accession no. 1016611.
This research has been supported by the U.S. Department of Energy, Office of Science (grant no. DESC0016605), and the California Energy Commission (grant no. EPC16063).
This paper was edited by Andrew Wickert and reviewed by two anonymous referees.
Abbott, M. B. and Refsgaard, J. C. (Eds.): Distributed Hydrological Modelling, vol. 22 of Water Science and Technology Library, Springer Netherlands, Dordrecht, https://doi.org/10.1007/9789400902572, 1996. a
Allen, R. G.: Crop evapotranspiration – Guidelines for computing crop water requirements – FAO Irrigation and drainage paper 56, FAO, 1998. a
Bao, C.: Understanding Hydrological And Geochemical Controls On Solute Concentrations At Large Scale, PhD thesis, Pennsylvania State University, 2016. a
Bao, C., Li, L., Shi, Y., and Duffy, C.: Understanding watershed hydrogeochemistry: 1. Development of RTFluxPIHM, Water Resour. Res., 53, 2328–2345, https://doi.org/10.1002/2016WR018934, 2017. a
Bergström, S.: The HBV model – its structure and applications, Tech. rep., 1992. a, b
Beven, K.: Changing ideas in hydrology – The case of physicallybased models, J. Hydrol., 105, 157–172, https://doi.org/10.1016/00221694(90)90161P, 1989. a
Beven, K.: RainfallRunoff Modelling, John Wiley & Sons, Ltd, Chichester, UK, https://doi.org/10.1002/9781119951001, 2012. a, b, c, d, e
Beven, K. and Germann, P. F.: Macropores and water flows in soils, Water Resour. Res., 18, 1311–1325, https://doi.org/10.1029/WR018i005p01311, 1982. a
Beven, K. J. and Kirkby, M. J.: A physically based, variable contributing area model of basin hydrology, Hydrolog. Sci. B., 24, 43–69, https://doi.org/10.1080/02626667909491834, 1979. a
Bhatt, G., Kumar, M., and Duffy, C. J.: A tightly coupled GIS and distributed hydrologic modeling framework, Environ. Modell. Softw., 62, 70–84, https://doi.org/10.1016/j.envsoft.2014.08.003, 2014. a, b
Blöschl, G., Bierkens, M. F., Chambel, A., et al.: Twentythree unsolved problems in hydrology (UPH) – a community perspective, Hydrolog. Sci. J., 64, 1141–1158, https://doi.org/10.1080/02626667.2019.1620507, 2019. a, b
Bobo, A. M., Khoury, N., Li, H., and Boufadel, M. C.: Groundwater flow in a tidally influenced gravel beach in prince william sound, alaska, J. Hydrol. Eng., 17, 478–494, https://doi.org/10.1061/(ASCE)HE.19435584.0000454, 2012. a
Brandes, D., Duffy, C. J., and Cusumano, J. P.: Stability and damping in a dynamical model of hillslope hydrology, Water Resour. Res., 34, 3303–3313, https://doi.org/10.1029/98WR02532, 1998. a
Cheema, T.: Depth dependent hydraulic conductivity in fractured sedimentary rocksa geomechanical approach, Arab. J. Geosci., 8, 6267–6278, https://doi.org/10.1007/s1251701416038, 2015. a
Chen, C. and Wagenet, R.: Simulation of water and chemicals in macropore soils Part 1. Representation of the equivalent macropore influence and its effect on soilwater flow, J. Hydrol., 130, 105–126, https://doi.org/10.1016/00221694(92)901066, 1992. a
Chen, Y. F., Ling, X. M., Liu, M. M., Hu, R., and Yang, Z.: Statistical distribution of hydraulic conductivity of rocks in deepincised valleys, Southwest China, J. Hydro., 566, 216–226, https://doi.org/10.1016/j.jhydrol.2018.09.016, 2018. a
Cohen, S. D. and Hindmarsh, A. C.: CVODE, a stiff/nonstiff ODE solver in C, Comput. Phys., 10, 138–143, https://doi.org/10.1063/1.4822377, 1996. a
Dickinson, R. E.: Modeling Evapotranspiration for ThreeDimensional Global Climate Models, Climate Processes and Climate Sensitivity, 29, 58–72, https://doi.org/10.1029/GM029p0058, 1984. a
Duffy, C. J.: A twostate integralbalance model for soil moisture and groundwater dynamics in complex terrain, Water Resour. Res., 32, 2421–2434, https://doi.org/10.1029/96WR01049, 1996. a
Farthing, M. W. and Ogden, F. L.: Numerical Solution of Richards' Equation: A Review of Advances and Challenges, Soil Sci. Soc. Am. J., 81, 1257, https://doi.org/10.2136/sssaj2017.02.0058, 2017. a
Fleming, M. J.: Hydrologic Modeling System HECHMS Quick Start Guide, U.S Army Corps of Engineers, 2010. a, b
Gochis, D., Yu, W., and Yates, D.: The NCAR WRFHydro Technical Description and User's Guide, version 3.0, NCAR Technical Document, Tech. Rep. May, NCAR, 2015. a, b
Gupta, H. V., Kling, H., Yilmaz, K. K., and Martinez, G. F.: Decomposition of the mean squared error and NSE performance criteria: Implications for improving hydrological modelling, J. Hydrol., 377, 80–91, https://doi.org/10.1016/j.jhydrol.2009.08.003, 2009. a
Hansen, N.: The CMA Evolution Strategy: A Comparing Review, in: Towards a New Evolutionary Computation, SpringerVerlag, Berlin/Heidelberg, 75–102, https://doi.org/10.1007/11007937_4, 2016. a
Hawkins, R. H., Hjelmfelt, A. T., and Zevenbergen, A. W.: Runoff probability, storm depth, and curve numbers, J. Irrig. Drain. E., 111, 330–340, https://doi.org/10.1061/(ASCE)07339437(1985)111:4(330), 1985. a
Hindmarsh, A. C., Brown, P. N., Grant, K. E., Lee, S. L., Serban, R., Shumaker, D. E., and Woodward, C. S.: SUNDIALS: Suite of nonlinear and differential/algebraic equation solvers, ACM T. Math. Software, 31, 363–396, 2005. a, b
Hindmarsh, A. C., Serban, R., and Reynolds, D. R.: User documentation for CVODE v5.0.0, Tech. rep., Center for Applied Scientific Computing Lawrence Livermore National Laboratory, 2019. a, b
Homer, C. and Fry, J.: The National Land Cover Database, US Geological Survey Fact Sheet, 2012. a, b
Hrachowitz, M. and Clark, M. P.: HESS Opinions: The complementary merits of competing modelling philosophies in hydrology, Hydrol. Earth Syst. Sci., 21, 3953–3973, https://doi.org/10.5194/hess2139532017, 2017. a
Jiang, X. W., Wan, L., Wang, X. S., Ge, S., and Liu, J.: Effect of exponential decay in hydraulic conductivity with depth on regional groundwater flow, Geophys. Res. Lett., 36, 3–6, https://doi.org/10.1029/2009GL041251, 2009. a
Kirkham, M.: Potential Evapotranspiration, in: Principles of Soil and Plant Water Relations, Elsevier, chap. 28, 501–514, https://doi.org/10.1016/B9780124200227.000288, 2014. a
Kumar, M. and Duffy, C. J.: Detecting hydroclimatic change using spatiotemporal analysis of time series in Colorado River Basin, J. Hydrol., 374, 1–15, https://doi.org/10.1016/j.jhydrol.2009.03.039, 2009. a
Kumar, M., Duffy, C. J., and Salvage, K. M.: A SecondOrder Accurate, Finite VolumeBased, Integrated Hydrologic Modeling (FIHM) Framework for Simulation of Surface and Subsurface Flow, Vadose Zone J., 8, 873–890, https://doi.org/10.2136/vzj2009.0014, 2009. a, b
Leavesley, G. H., Lichty, R. W., Troutman, B. M., and Saindon, L. G.: Precipitationrunoff modeling system; user's manual, Tech. rep., USGS, Denver, Colorado, https://doi.org/10.3133/wri834238, 1983. a
Leonard, L. and Duffy, C. J.: Essential Terrestrial Variable data workflows for distributed water resources modeling, Environ. Modell. Softw., 50, 85–96, https://doi.org/10.1016/j.envsoft.2013.09.003, 2013. a
Liang, X., Lettenmaier, D. P., and Wood, E. F.: Onedimensional statistical dynamic representation of subgrid spatial variability of precipitation in the twolayer variable infiltration capacity model, J. Geophys. Res., 101, 21403, https://doi.org/10.1029/96JD01448, 1996. a, b
Lin, L., Jia, H., and Xu, Y.: Fracture network characteristics of a deep borehole in the Table Mountain Group (TMG), South Africa, Hydrogeol. J., 15, 1419–1432, https://doi.org/10.1007/s100400070184y, 2007. a
Lin, P., Yang, Z. L., Gochis, D. J., Yu, W., Maidment, D. R., SomosValenzuela, M. A., and David, C. H.: Implementation of a vectorbased river network routing scheme in the community WRFHydro modeling framework for flood discharge simulation, Environ. Modell. Softw., 107, 1–11, https://doi.org/10.1016/j.envsoft.2018.05.018, 2018. a, b
Maidment, D. R.: Handbook of hydrology, vol. 9780070, McGrawHill New York, 1993. a, b, c
McKay, L., Bondelid, T., Dewald, T., Johnston, J., Moore, R., and Rea, A.: NHDPlus version 2: user guide, Tech. rep., US Environmental Protection Agency, 2012. a
Moradkhani, H. and Sorooshian, S.: General Review of RainfallRunoff Modeling: Model Calibration, Data Assimilation, and Uncertainty Analysis, in: Hydrological Modelling and the Water Cycle, Springer Berlin Heidelberg, Berlin, Heidelberg, 1–24, https://doi.org/10.1007/9783540778431_1, 2008. a
Nash, J. E. and Sutcliffe, J. V.: River Flow Forecasting Through Conceptual Models Part Ia Discussion of Principles, J. Hydrol., 10, 282–290, https://doi.org/10.1016/00221694(70)902556, 1970. a
Penman, H. L.: Natural evaporation from open water, hare soil and grass, P. Roy. Soc. Lond. A Mat., 193, 120–145, https://doi.org/10.1098/rspa.1948.0037, 1948. a
Petty, T. R. and Dhingra, P.: Streamflow Hydrology Estimate Using Machine Learning (SHEM), J. Am. Water Resour. As., 54, 55–68, https://doi.org/10.1111/17521688.12555, 2018. a
Qu, Y.: An integrated hydrologic model for multiprocess simulation using semidiscrete finite volume approach, PhD thesis, Pennsylvanis State University, 2004. a, b, c
Rasouli, K., Hsieh, W. W., and Cannon, A. J.: Daily streamflow forecasting by machine learning methods with weather and climate inputs, J. Hydrol., 414415, 284–293, https://doi.org/10.1016/J.JHYDROL.2011.10.039, 2012. a
Refsgaard, J. C., Sørensen, H. R., Mucha, I., Rodak, D., Hlavaty, Z., Bansky, L., Klucovska, J., Topolska, J., Takac, J., Kosc, V., Enggrob, H. G., Engesgaard, P., Jensen, J. K., Fiselier, J., Griffioen, J., and Hansen, S.: An integrated model for the Danubian Lowland – methodology and applications, Water Resour. Manag., 12, 433–465, https://doi.org/10.1023/A:1008088901770, 1998. a, b
Santhi, C., Srinivasan, R., Arnold, J. G., and Williams, J. R.: A modeling approach to evaluate the impacts of water quality management plans implemented in a watershed in Texas, Environ. Modell. Softw., 21, 1141–1157, https://doi.org/10.1016/j.envsoft.2005.05.013, 2006. a, b
Shen, C. and Phanikumar, M. S.: A processbased, distributed hydrologic model based on a largescale method for surfacesubsurface coupling, Adv. Water Resour., 33, 1524–1541, https://doi.org/10.1016/j.advwatres.2010.09.002, 2010. a, b, c, d, e, f, g
Shen, C., Laloy, E., Elshorbagy, A., Albert, A., Bales, J., Chang, F.J., Ganguly, S., Hsu, K.L., Kifer, D., Fang, Z., Fang, K., Li, D., Li, X., and Tsai, W.P.: HESS Opinions: Incubating deeplearningpowered hydrologic science advances as a community, Hydrol. Earth Syst. Sci., 22, 5639–5656, https://doi.org/10.5194/hess2256392018, 2018. a
Shewchuk, J. R.: Triangle: Engineering a 2D quality mesh generator and Delaunay triangulator, in: Applied Computational Geometry Towards Geometric Engineering, edited by: Lin, M. C. and Manocha, D., Springer Berlin Heidelberg, Berlin, Heidelberg, 203–222, 1996. a, b
Shi, Y., Davis, K. J., Zhang, F., Duffy, C. J., and Yu, X.: Parameter estimation of a physically based land surface hydrologic model using the ensemble Kalman filter: A synthetic experiment, Water Resour. Res., 50, 706–724, https://doi.org/10.1002/2013WR014070, 2014. a
Shi, Y., Baldwin, D. C., Davis, K. J., Yu, X., Duffy, C. J., and Lin, H.: Simulating highresolution soil moisture patterns in the Shale Hills watershed using a land surface hydrologic model, Hydrol. Process., 29, 4624–4637, https://doi.org/10.1002/hyp.10593, 2015. a
Shi, Y., Eissenstat, D. M., He, Y., and Davis, K. J.: Using a spatiallydistributed hydrologic biogeochemistry model with a nitrogen transport module to study the spatial variation of carbon processes in a Critical Zone Observatory, Ecol. Model., 380, 8–21, https://doi.org/10.1016/j.ecolmodel.2018.04.007, 2018. a
Shu, L.: ModelIntercomparisonDatasets/CacheCreek v0.0.1, Zenodo, https://doi.org/10.5281/zenodo.3566034, 2019a. a, b, c, d, e
Shu, L.: SHUDSystem/SHUD: v1.0, Zenodo, https://doi.org/10.5281/zenodo.3561293, 2019b. a, b, c, d
Shu, L.: SHUDSystem/SHUD_User_Guide: v1.0, Zenodo, https://doi.org/10.5281/zenodo.3561295, 2019c. a
Shu, L.: ModelIntercomparisonDatasets/Vauclin1979: v0.0.1, Zenodo, https://doi.org/10.5281/zenodo.3566020, 2019d. a
Shu, L.: ModelIntercomparisonDatasets/VCatchment: v0.0.1, Zenodo, https://doi.org/10.5281/zenodo.3566022, 2019e. a, b
Shu, L.: SHUDSystem/rSHUD: Ver 1.0, Zenodo, https://doi.org/10.5281/zenodo.3758097, 2020. a, b, c, d, e
Soil Survey Staff: Gridded Soil Survey Geographic (gSSURGO) Database for the Conterminous United States, Tech. rep., United States Department of Agriculture, 2015. a, b
Taylor, G. S.: Drainable porosity evaluation from outflow measurements and its use in drawdown equations, Soil Sci., 90, 338–343, https://doi.org/10.1097/0001069419601200000004, 1960. a
U.S. Geological Survey: USGS National Elevation Dataset (NED) 1 arcsecond Downloadable Data Collection from The National Map 3D Elevation Program (3DEP) – National Geospatial Data Asset (NGDA), Tech. rep., U.S. Geological Survey, 2016. a, b
VanderKwaak, J. E.: Numerical simulation of flow and chemical transport in integrated surfacesubsurface hydrologic systems, PhD thesis, University of Waterloo, 1999. a
Vanderstraeten, D. and Keunings, R.: Optimized partitioning of unstructured finite element meshes, International Journal for Numerical Methods in Engineering, 38, 433–450, https://doi.org/10.1002/nme.1620380306, 1995. a
Vauclin, M., Khanji, D., and Vachaud, G.: Experimental and numerical study of a transient, twodimensional unsaturatedsaturated water table recharge problem, Water Resour. Res., 15, 1089–1101, https://doi.org/10.1029/WR015i005p01089, 1979. a, b, c, d, e, f, g, h
Vivoni, E. R., Ivanov, V. Y., Bras, R. L., and Entekhabi, D.: Generation of Triangulated Irregular Networks Based on Hydrological Similarity, J. Hydrol. Eng., 9, 288–302, https://doi.org/10.1061/(asce)10840699(2004)9:4(288), 2004. a
Vivoni, E. R., Ivanov, V. Y., Bras, R. L., and Entekhabi, D.: On the effects of triangulated terrain resolution on distributed hydrologic model response, Hydrol. Process., 19, 2101–2122, https://doi.org/10.1002/hyp.5671, 2005. a
Vivoni, E. R., Mascaro, G., Mniszewski, S., Fasel, P., Springer, E. P., Ivanov, V. Y., and Bras, R. L.: Realworld hydrologic assessment of a fullydistributed hydrological model in a parallel computing environment, J. Hydrol., 409, 483–496, https://doi.org/10.1016/j.jhydrol.2011.08.053, 2011. a, b
Xia, Y., Mitchell, K., Ek, M., Sheffield, J., Cosgrove, B., Wood, E., Luo, L., Alonge, C., Wei, H., Meng, J., Livneh, B., Lettenmaier, D., Koren, V., Duan, Q., Mo, K., Fan, Y., and Mocko, D.: Continentalscale water and energy flux analysis and validation for the North American Land Data Assimilation System project phase 2 (NLDAS2): 1. Intercomparison and application of model products, J. Geophys. Res.Atmos., 117, D03109, https://doi.org/10.1029/2011JD016048, 2012. a, b, c
Zhang, Y., Slingerland, R., and Duffy, C.: Fullycoupled hydrologic processes for modeling landscape evolution, Environ. Modell. Softw., 82, 89–107, https://doi.org/10.1016/j.envsoft.2016.04.014, 2016. a