Simulator for Hydrologic Unstructured Domains (SHUD v1.0): numerical modeling of watershed hydrology with the finite volume method

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 disconnect in the characteristic time and space scales of overland and groundwater flow. However, significant difficulties remain for the development of the efficient implementation of extensible modeling systems that operate robustly across complex regions. This paper introduces the Simulator for Hydro5 logic Unstructured Domains (SHUD), an integrated, multi-process, multi-scale, flexible-timestep model, in which hydrologic processes are fully coupled using the Finite Volume Method. SHUD integrates overland flow, snow accumulation/melt, evapotranspiration, subsurface and groundwater flow, and river routing, which realistically captures the physical processes in a watershed. SHUD incorporates one-dimension unsaturated flow, two-dimension groundwater flow, and river channel network fully connected with hillslopes via overland flow and baseflow. 10 The paper introduces the design of SHUD, from the conceptual and mathematical description of hydrologic processes in a watershed to computational structures. To demonstrate and validate the model performance, we employ three hydrologic experiments: the V-Catchment experiment, Vauclin’s experiment, and a model study of the Cache Creek Watershed in northern California, USA. Ongoing applications of the SHUD model include hydrologic analyses of the hillslope to regional scales (1m to 10km), water resource and stormwater management, and interdisciplinary research with related fields in limnology, 15 agriculture, geochemistry, geomorphology, water quality, ecology, climate and landuse change. The strength of SHUD is its’ flexibility as a scientific or resource evaluation tool where modeling and simulation are required. Copyright statement. TEXT

. The family tree of PIHM and SHUD. PIHM and SHUD share the same fundamental conceptual model, but use different realization.
The PIHMgis and SHUDtoolbox are GIS-tools for pre-and post-processing.
The 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 states-processesconstitutive relations depending on the objectives of the numerical experiment for research purpose.
As a distributed hydrologic model, the computational domain of the SHUD model is discretized using an unstructured 90 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). We developed the R package 95 SHUDtoolbox ( https://github.com/shud-system/SHUDtoolbox), helping users to generate the triangular domain. 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 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 100 the conceptual design and equations used in the model. In section 3, we employ three hydrologic experiments to demonstrate the simulation and capacity of the model. The three applications presented here are (1) the V-Catchment 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.  We begin our introduction to the SHUD model with a conceptual description of water movement in a watershed (Fig. 2) describing 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 plant-canopy interception 110 and through-fall. When precipitation exceeds the interception capacity, through-fall or excess precipitation falls to the land surface. Snowfall accumulation and melting is an important land-surface 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 115 layer) is constrained to 1-D vertical flow and the saturated groundwater layer admits 2-D flow. These layers overlie an impermeable or low-permeable 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 120 also accounts for upward capillary flow from a shallow water table depending on soil moisture and vegetation conditions. Lateral (2-D) 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 2-D 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 125 stage rises above bankfull storage during flooding events or in an arid region where the river recharges the local groundwater.
Evaporation generally represents the major 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.

130
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/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 ∆S = P − E − Q, where ∆S, P , E, Q are total water storage, precipitation, evapotranspiration and discharge, respectively. Modifications to these conditions can be applied to realize additional 135 flows such as pumping wells, irrigation, basin diversions, etc.. The model is able to specify boundary conditions based on the various research areas and purposes.
-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).

140
-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 any distributed models, as the various models still need discretized domains 145 instead of continuous space.
-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.

Mathematical structure
The notation used in this section is summarized in 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).
The algorithm of triangulation is not limited to Delaunay, so any GIS software and advanced programming language (R, Matlab 155 and Python) are capable of generating 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  relationship between river channels and hillslope cells, and (c) water fluxes between river segments and hillslope cells. elevation (z sf ) and groundwater table (z gw ) above datum, i.e. D us = z sf − z gw . When the groundwater table reaches the land 160 surface (z gw > z sf ), the thickness of unsaturated zone → 0. The hydrologic model uses the Method of Moments to reduce partial differential equations (PDE's) into ordinary differential equations (ODE's) and solves the ODE system using a global 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 Unsaturated zone Richards Equation where the discrete state vector is denoted by Y , Y 0 are the initial conditions and f (t, Y ) denotes the equations governing the hydrologic flow, which are described in this section.

165
The system of ODEs describing the hydrologic processes are fully coupled and solved simultaneously at each time step (∆t = t n −t n−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 dY dt | tn−1 . The technical description of the CVODE solver can be found in the literature (Hindmarsh et al., 2019(Hindmarsh et al., , 2005Cohen and Hindmarsh, 1996). The governing equations in SHUD are provided in table 1. The change of fluxes of ET and interception within a short period (such as one hour) are relatively slow so that full coupling of the ET with soil water is not necessary for this model. Instead, the interception, ET and snow calculations solved explicitly at MTS, while the calculation of Y sf , Y us , Y gw and Y riv use the implicit time step (ITS).

175
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, dY s 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, 180 unsaturated layer, saturated layer, and river channel.

Vegetation and evapotranspiration
The interception is a 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 185 S * ic = C ic LAI, where LAI represents the coverage of vegetation canopy over the land area (area of leaves over the area of land, m 2 /m 2 ), and C ic is interception coefficient [m]. The default C ic is taken to be 0.2kg/m 2 as suggested in Dickinson (1984).
The interception is equal to the deficit of interception -the difference between interception capacity (S * ic ) and existing interception storage (S ic ). If precipitation is less than the deficit, interception is equal to the precipitation rate (see Fig. 6).
(2) Figure 6. The three conditions for the interception calculation within the imaginary canopy bucket. The through-fall or excess precipitation after interception is the water gaining on the land surface.
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,195 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, factored by varying albedo, LAI, and roughness length.

200
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 210 evaporation from soil moisture (E sm ), E s = E sp + E 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 sf /∆t > E 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 sm = E s , E sp = 0): Transpiration also has two potential sources: soil moisture and groundwater from the groundwater table and root depth for the land-use class. Once the groundwater table is higher than the root zone depth, vegetation uses groundwater, and soil moisture stress for transpiration is equal to zero (β s = 0).
Water balance associated with snow accumulation is quantified via Snow melt rate is determined with snow melt factor (m f ), air temperature (T ) and temperature threshold (T 0 ) at which snow melt occurs. This formulation is often referred to as the degree-day method, in which the values of the snow melt factor and temperature threshold are empirical (Maidment, 1993;Beven, 2012). The water from snow melt is considered as a direct water contribution to the land surface.

Water on the land surface
Water balance on the land surface is given by:

230
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 j s in direction j is calculated with Manning's equation (13). Here y sf is effective water height, determined by the gradient between two cells, Estimating infiltration utilizes 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 10cm, 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 245 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.

Unsaturated zone 250
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 the equation 16. The soil moisture content to field capacity controls the recharge rate.
Because of the simplification of two-layer description of vertical aquifer profile, we use 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.

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

265
The water balance of groundwater is controlled by the equation 18.
The calculation of horizontal groundwater flow uses the Boussinesq equation. When the bottom of the root zone is lower than The horizontal groundwater flux Q j g is determined by the Darcy equation and Dupuit-Forchheimer assumption. Above z b is the elevation of impervious bedrock, z j b is the bedrock elevation of its jth neighbor cell and d j is distance between the centroids of two adjcent cells, so the gradient between the two cells is (y The effective conductivity for the groundwater flow is the mean value of the effective horizontal conductivity over the two cells. The cross-sectional area along 275 the groundwater flux is equal to L j × y gw . In equation 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 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). table is below the level of the macropores, K eg is equal to 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.

Water in streams 290
The water balance in river channels is described by The mass balance in each river channel consists of four parts: Q j s , the overland flow from cells (1 to N c cells) that intersect with the river channel; Q j g , the lateral groundwater flux from intersection with the jth cell; Q j up , 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 295 upstream channels is nonnegative, but only one downstream channel is permitted. We assume 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. 4(b). 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.
where A cs is the cross-section area of the river reach, and P and s 0 are average wet perimeter and average slope of a river 305 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 pre-computation 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.

310
The overland flow between river segment and associated hillslope cell (Q sr ) is calculated as follows: H riv = y riv + z rb , H riv − z bank , H riv > z bank and H csf < z bank , H csf − z bank , H riv < z bank and H csf > z bank .

315
Here z bank is the elevation of riverbank or levee, 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

Applications
In this section, we present the results of applying SHUD to three applications: first, we use the V-catchment experiment to validate the calculation of overland flow and river routing in an idealized catchment; second, we use Vauclin's experiment 325 (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, a headwater catchment in Sacramento Watershed of Northern California.

V-Catchment
The V-Catchment (VC) experiment is a standard test case for numerical hydrologic models to validate their performance for 330 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 × 1000m 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 18mm/hr and stops after 90 min, producing 27 mm of accu-335 mulated 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 340 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/). 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 345 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 check the numerical method, we verify the bias of mass-balance in the model and assess the differences among input, output and storage change in the system (equation 35). The bias in the model result is ∼ 0.2%.

Vauclin's experiment
Vauclin's laboratory experiment (Vauclin et al., 1979)  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 hours at several locations along the length of the box. (Vauclin et al., 1979) also uses a 2-D (vertical and horizontal) numeric model to simulate the soil moisture and groundwater table. The maximum bias between measurement and simulation was 0.52m, according to the digitalized value of Vauclin et al., 1979, Fig. 10.

365
Besides the parameters specified in (Vauclin et al., 1979), additional information is needed in the SHUD model, including the α and β in the van Genutchen 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.001m 3 /m 3 , α = 0.3 and β = 5.2. The simulated results in our 370 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) the validity of the vertical and horizontal flow assumptions in the SHUD model. SHUD simulated the groundwater table at all four measurement points (see Fig. 10(b)). The maximum bias between simu-375 lation and Vauclin's observations is 5.5cm, with R 2 = 0.99, which is comparable to the bias 5.2cm of numerical simulation in (Vauclin et al., 1979). By adding more layer structure, the bias in simulation decreases to 3cm. Indeed, the simplifications employed by SHUD for the unsaturated and saturated zone benefits the computation efficiency while limiting the applicability where fine-scale soil profile information is required.
The simulations, compared against Vauclin's experiment, generally validate the algorithm for infiltration, recharge, and 380 lateral groundwater flow. A more reliable vertical flow within the unsaturated layer requires multiple layers, which is planned in the next version of SHUD.

Cache Creek Watershed
The Cache Creek Watershed (CCW) is a headwater catchment with area 196.4km 2 in the Sacramento Watershed in Northern California (Figures 11 (a), (b) and (c)). The elevation ranges from 450m to 1800m, with a 0.38m/m average slope, where the 385 steep implies a particular challenge for numeric models.
Based on NLDAS-2 from 2000 to 2017, the mean temperature and precipitation were 12.8 • C and ∼ 817mm, respectively, in this catchment. Precipitation varies seasonally; Winter and spring are the predominant wet seasons. (Fig. 12).

395
The unstructured domain of the CCW (Fig. 11 (d)) is built with SHUDtoolbox. The number of triangular cells is 1147, with a mean area of 0.17km 2 . The total length of the river network is 126.5km and consists of 103 river reaches and in which the highest order of stream is 4. With a calibrated parameter set, the SHUD model took 5 hours to simulate 18 years (2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015)(2016)(2017) in the CCW, with a non-parallel configuration (OpenMP is disabled on Mac Pro 2013 Xeon 2.7GHz, 32GB RAM). Figure 13 reveals the comparison of simulated discharge against the observed discharge at the gage station of USGS 400 11451100 (https://waterdata.usgs.gov/ca/nwis/uv/?site_no=11451100). The calibration procedure exploits the Covariance Matrix Adaptation -Evolution Strategy (CMA-ES) to calibrate automatically (Hansen, 2016). The calibration program assigns 72 children in each generation and keeps the best child as the seed for 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.

405
In the calibration period, Nash-Sutcliffe efficiency (NSE Nash and Sutcliffe (1970)), Kling-Gupta Efficiency (KGE, Gupta et al. (2009)) and R 2 is 0.72, 0.83 and 0.72 respectively (Fig. 13. The goodness-of-fit 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 NLDAS-2 data, or (2) over-fitting in the calibration, as the 410 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 three 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%. 415 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 30m aquifer, the results represent the groundwater within this 30-meter aquifer only. The groundwater table and elevation along the green line on the map are extracted and plotted in the bottom figure. The gray ribbon is the 30-meter aquifer, and the blue line is the location where groundwater storage is larger than zero. The green polygons with the right axis are the 420 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.

Summary
The features of the SHUD model is summarized following. Bottom: fluxes of precipitation, actual evapotranspiration, potential evapotranspiration and discharge at the outlet.
-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 the national soil data and surficial geology. Several other groups have used PIHM, a SHUD ancestor to couple processes for biochemistry, reaction transport, landscape, geomorphology, limnology, and other related research areas.

430
-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 solved in SHUD solves with a state-of-the-art parallel ODE solver, known as CVODE (Hind- -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 is valuable for community model coupling.

440
-SHUD can estimate either a long-term hydrologic yield or a single-event flood.
-SHUD is an open-source model -available on GitHub.

Differences from PIHM
As a descendant of PIHM, SHUD inherits the fundamental idea of conceptual structure and solving hydrologic variables in CVODE. The code has been completely rewritten in a new programming language, with a new discretization and corresponding 445 improvements to the underlying algorithms, adapting new mathematical schemes and flexible input/output data formats. Although SHUD is forked from PIHM's track, 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 the following: 1. SHUD is written in C++, an object-oriented programming language with functionality to avoid risky memory leaks from C. Most of the functions in SHUD do not exist in PIHM, which are newly defined because of the brand-new design Equation, are shared in SHUD and PIHM. So SHUD and PIHM follow a similar fundamental perceptual model, but different mathematical and computational strategies.
2. SHUD implements a re-design 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 455 sink problems in cells that share one node with a starting river channel, and fatherly slow down the performance of the simulation.
3. Although the mathematical equations underlying SHUD are generally the same as PIHM, the numerical formulation of processes, the coupling strategy and input/output have been greatly enhanced. The common computations in different functions within various processes are extracted and defined as shared inline functions, which maintain the consistency adds a debug mode that monitors potential errors in parameters and memory operations.

480
includes a series of R codes for pre-and post-processing data, visualization, and data analysis (that will be discussed in future work).

Conclusion and future work
The Simulator for Hydrologic Unstructured Domains (SHUD) is a multi-process, multi-scale and multi-temporal model that integrates major hydrologic processes and solves the physical equations with the Finite Volume Method. The governing equa-485 tions are solved within an unstructured mesh domain -triangular cells. The variables in the surface, vadose layer, groundwater and river routing are fully coupled together with a rather fine time-step. The SHUD uses the one-dimensional unsaturated flow and two-dimensional 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.

490
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 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) 495 climate change, and (5) land-use change. In summary, SHUD is a valuable scientific tool for any modeling task associating with hydrologic responses.
Several tests and experiments have been left for the future due to the length of this paper. Future work concerns applications of the model to improve our understanding of hydrologic responses.
-This paper does not mention data preparation, even though a tutorial is available online. The data preparation covers 500 1) source and format of data, 2) processing the spatial data and time-series data, 3) how to build the optimal triangular mesh, and 4) result handling and visualization.
-A comprehensive comparison between the PIHM and SHUD models is necessary. Benchmarks of several data-rich watersheds could be used to demonstrate their comparative performance. Since they are not identical-input files cannot be converted between them directly; such a comparison would require building test cases on the same raw data (watershed 505 boundary, stream network, DEM, soil, geology, land use, and forcing data).
-It would be valuable to perform a comprehensive parameter sensitivity test within the model. The sensitivity experiments could focus on the impact of these parameters on soil moisture, groundwater level, ET ratio, as well as streamflow.
-Experiments examining model spin-up could be interesting. The spin-up issue exists in every fully coupled surfacesubsurface hydrologic model because of the relatively slow groundwater flow that influences the estimation of ground-510 water table, baseflow and soil moisture. Deeper groundwater in the simulation requires a more extended spin-up period.
The available data, however, limits the period of spin-up. So a model must perceive the optimal spin-up of the modeling domain with any initial conditions, or warm-start conditions, e.g., the conditions approaching the equilibrium.
-Modelers may consider examining various calibration strategies with the SHUD model. Besides existing calibration methods, modelers may calibrate on groundwater, soil moisture and baseflow, as well as streamflow -spatially and 515 temporally.
Code and data availability. The source code of SHUD model is kept updating at https://github.com/SHUD-System/SHUD. The code and data used for this page is archived at ZENODO: