Advances made in numerical algorithms used to solve the Navier-Stokes equations (which describe the fluid motion) and its simplified variants are pushing the frontiers of knowledge in all engineering disciplines. When these mathematical advances are integrated with the increased computational abilities, researchers can now attempt to solve large scale problems at small spatial grid size that hitherto have been difficult, primarily because of the lack of computational resources by which the modeling tasks can be accomplished in a reasonable amount of CPU time. While the popular algorithms used to solve the flow governing equations are from finite difference, finite element and finite volume methods, over the last few decades new algorithms from the family of spectral methods, Total variation diminishing, Non Oscillatory and Boundary element methods   although have not been rigorously applied are gaining popularity among the modeling audience.
The numerical codes that were written in the early years of computers for solving the flow equations were primarily in Fortran language. While some of these codes have been re-written using an object-oriented language, the rest have been integrated with advanced graphical user interface applications to make their execution and output visualization easier. Since hundreds of years of testing have gone into these codes by the audience from around the world, their results tend to act as a benchmark for codes that have evolved over the last two decades. As experimental investigations for complex fluid flow problems is not possible and as theoretical solutions for these problems do not exist, the application of numerical models will continue to become more popular in the coming decades.
In this work, we apply the DHM  for a complex natural overland surface flow domain. Overland flow is the surface runoff, which normally occurs as a sheet flow on the ground without any concentration in channels. DHM (written in Fortran) has been first developed for the United States Geological Survey, and the mathematical underpinnings in DHM have been used in some of the currently popular models . The DHM was originally based upon the groundwater program MODFLOW and then evolved into the control volume modeling approach. The chosen flow domain (Figure 1(a)) because of its topography and flow characteristics is complex by any benchmark. An intense rainfall event can bring a wave of flow and mud along the sloping domain. Their end combined impact on moving vehicles and properties has the potential to cause significant damage. Analyzing these flows is important from both the technical and the environmental point of view. Reliably predicting the transient flow variables (i.e., depth, velocity, and Froude number) for different what if flow scenarios and using them in engineering designs, can translate to reduced damage to the properties and can save human lives.
There are many codes that are being used for modeling overland flows. These include FLOW-3D , TUFLOW , SWMM , HEC-HMS , MIKE URBAN , FLO-2D  and XPSWMM . The respective model backgrounds and sample applications can be found in the given references. Identifying the most suitable code for a large scale application can be challenging because of the required optimal input variables and the ability of the codes to incorporate all the features in the physical domain. There is no universal model that can be reliably applied for all the domains. Depending on the site, the data requirements, and the flow variables of interest, the performance of one of the codes may tend to be better than others. The sophistication in the core engine across these codes depends on the assumptions made in the flow equations. Their complexity can vary from linear lumped black-box models to nonlinear
(a) (b) (c)
Figure 1. (a) Google map of the physical domain. The water rushes down the north direction primarily along the natural water course path before it encounters the asphalt roadway that crosses the study watercourse, after which most of the flow continues further north into the empty land; (b) The DHM domain with one row of grids, the 30 ft spaced cells are squares; (c) The DHM domain with two rows of 30 ft spaced square grids. Input hydrograph for (b) is at node 1 and the outflow hydrograph is at node 14. For (c), the inflow hydrograph can be either at node 1 or distributed across nodes 1 and 15. The outflow is at node 28.
physical-based distribution models . DHM code is simple to use and as documented in , can be used across a wide range of overland flow applications. The source code and related literature can be downloaded from http://www.diffusionhydrodynamicmodel.com and further enhanced by researchers.
The outline of this paper is as follows. The flow equations that describe the two-dimensional unsteady overland flow solved in DHM are described in the next section. The DHM numerical algorithm characteristics are briefly discussed after which the test problem, the computational grid, boundary flow conditions are detailed. The model results for different grid and roughness conditions are presented, and its performance analyzed.
2. Governing Equations
The two-dimensional flow continuity and momentum equations along the X and Y axis (assuming a constant fluid density without sources or sinks in the flow field and hydrostatic pressure distribution) can be written as 
in which , are flow rates per unit width in the x, y—directions; , represents friction slopes in x, y—directions; H, h, g stand for water surface elevation, flow depth, and gravitational acceleration, respectively; and x, y, t are spatial and temporal coordinates.
The local and convective acceleration terms can be grouped and Equations (2) and (3) are rewritten as
where represents the sum of the first three terms in Equations (1) and (2) divided by gh. Assuming the friction slope to be approximated by the Manning’s formulae, the flow equation in the U.S. customary units for flow in the x or y directions:
Equation (5) can be rewritten in the general case as
The symbol S in Equation (7) indicates the flow direction which makes an angle of with the positive x-direction. By assuming the value of m to be negligible, the diffusion model can be expressed as:
Two-dimensional DHM is formulated by substituting Equation (8) into Equation (1)
If the momentum term groupings were retained, Equation (9) can be written as
and are also functions of respectively.
3. Numerical Algorithm
To maintain continuity in the discussion, while salient aspects of the DHM numerical algorithm are presented here and readers are referred to  for a complete description of the model algorithm. Sample data input files and the executable model codes are also present at the model companion website, http://www.diffusionhydrodynamicmodel.com. For uniform grid elements, the integrated finite difference version of the nodal domain integration (NDI) method  is used to solve the flow equations. The model requires the computational domain to have uniform spaced square grids, and the grid connectors for each grid node need to be specified. The latter is accomplished at all the interior nodes by specifying the adjacent grid numbers along with the North, East, South, and West directions. Should an adjacent node be absent in any of these directions (i.e., along with the boundary nodes in the four directions), “0” is specified. The elevation data at the center of the grid is an input variable which can be obtained from the topographical map of the model area. The time step (Δt) in the explicit DHM formulation is governed by the Courant-Friedrich-Lewy stability condition . The model can self-adjust and refine the chosen time step, within the given tolerance, so that the stability condition is met.
DHM facilitates transient simulation. Predicting transient flow velocities and Froude numbers are important for any engineering designs. The model requires a boundary variable to be specified at the upstream and downstream end. In this application, at the upstream, the inflow hydrograph was specified (discussed in the next section). An unsteady flow model, like DHM, will enable to follow the peak hydrograph flow as it propagates into the flow domain.
Some of the features that were first presented in the DHM model but later carried over to current popular CFD models include: One of the long-recognized challenges in using explicit formulations is in using a constant time step as the solution progresses towards a steady-state solution or to the desired transient period. A constant time step can be computationally intensive and require a significant amount of CPU time. The self-adjustable or self-adaptive variable time step algorithms are powerful tools that are inbuilt into the current codes, which help in choosing an acceptable time step to ensure that the Courant-Friedrich-Lewy (CFL) stability criteria is met. Devoid of such a tool, the user has to enter a best guess, time step (dt) and manually iterate its value until the solution converges. The presence of source terms, which are an important component in the equations to address the physics of flow, coupled with varying bottom roughness coefficient in the domain, will prevent the user from theoretically calculating the optimal time step. Although most of the available commercial software  have this option, in early 1980s DHM is perhaps one of the few codes that had this paradigm inbuilt into it. DHM, based on user-defined time step, time increment and decrement values, and allowable tolerance, can compute the optimal variable time step value. The computational engine in the popular FLO-2D model evolved from the DHM . In FLO-2D, the DHM algorithm was expanded and revised to make it versatile across diverse flooding scenarios. The core features in DHM, i.e., discretizing the domain into uniform square grids, allocating a location number to each cell, assigning an elevation along with roughness value to each cell and routing the flood through the grids based on depth estimates form the background engine in FLO-2D solver.
4. Test Problem
The above described DHM model was applied to the flow domain shown in Figure 1(a). The overland flow direction in the domain is from Southwest to Northeast, and the slope falls at a steep rate. The width of the natural watercourse pathway at the southwest corner is about 30 ft, and it widens to around 60 ft before the water meets the road. The road is 60 ft wide and is relatively flat. The Manning roughness value for the earth was taken as 0.03 and for the road was 0.015. Figure 1(b) shows the computational domain with 30 ft square grids. The ground elevation at the center of grid 1 is 2690 ft, and the ground elevation at center of node 14 is 2668 ft. To take into account the wider natural watercourse pathway near the road, the DHM domain was reconfigured with two rows of 30 ft square grids shown in Figure 1(c). Results for these both computational domains are presented, later. Grid connections used in the model at a sample node for both these domains are shown in Figure 2. While node 11 in Figure 1(b) has grids in North and South directions, for the domain in Figure 1(c), it has an additional node in the East direction. In directions where there is no node, a zero is assigned, as shown in Figure 2. The cross-sectional view of the ground elevation along with the road location is shown in Figure 3.
Boundary conditions are a required component in all numerical models, and the boundary condition used should not only be mathematically consistent, but it should also represent the physics of flow accurately. In this application, for the domain shown in Figure 1(b), at the upstream node 1, the discharge hydrograph
Figure 2. (a) Grid connections for node 11 in Figure 1(b); (b) Grid connections for node 11 in Figure 1(c).
Figure 3. Cross-sectional view along the natural watercourse pathway. The transient water level at the outflow node is at time = 0.5 hours.
(Figure 4) was specified. The peak discharge of 756 cfs occurs at 1800 s and tapers to 22 cfs at 7200 s. The runoff hydrograph was calculated using a small-area hydrograph approach based upon Rational Method flow rate estimates. For the domain in Figure 1(c), the inflow hydrograph (Figure 4) can be applied either at node 1, or it is distributed across nodes 1 and 15 (Figure 5). Results from both of these combinations are presented next. At the downstream end, critical depth condition was specified for all the runs.
The DHM was applied to predict the maximum flow depth at the north end of the roadway. Models like DHM that can perform transient simulations give more flexibility to researchers for testing various flow combinations. Figure 6 is the plot of the transient depth and discharge profiles at the outlet for the domain in Figure 1(b). The total simulation was done until time = 2 hours. DHM requires the specification of flow resistance or “roughness” value, individually for
Figure 4. Inflow hydrograph at node 1 for the domain shown in Figure 1(b).
Figure 5. Inflow hydrograph at nodes 1 and 15 for the domain shown in Figure 1(c) [The inflow values from Figure 4 at each time step are halved, as the input hydrograph is distributed across two nodes].
each computational cell. A manning’s roughness value of 0.03 and 0.015 were used for earth and road. The output was recorded at 0.1 hour period. Overland flow models require smaller cell sizes to reliably predict the preferential flow paths (in the presence of any obstacles) and within the natural watercourse path. In this investigation, we chose a cell size of 30 ft, as the width of the natural watercourse pathway at the upstream is in the order of 30 ft. The optimal cell size is typically influenced by the detailed level that is being predicted and the associated CPU time. An ideal grid size should be small enough to capture the physics of flow and not large so that no additional errors are introduced. Although the currently available computational resources (both memory and processor speed) are enabling researchers to use grids in the order of millions even across a square feet area of the computational domain, the demands of researchers have
Figure 6. (a) Transient depth and (b) Outflow hydrographs at node 14 for domain shown in Figure 1(b).
also steadily risen. It is not uncommon for large scale application to use days of CPU time across multi-process parallel computers. DHM requires that the grids be uniform size square and it does not have tools like grid stretching or compression.
The results of numerical models tend to be sensitive across different parameters, including grid size, roughness coefficient, and time step. Several parameter sensitivity tests were carried out across the time step, grid spacing and roughness values (for earth and road). Figure 7 is the plot of depth profile for various roughness combinations. The flow depth prediction (Figure 7) is not strongly influenced by the roughness coefficient. This is an important observation as there are no universal set of constants that can be used in real-life modeling flows which are accompanied by complex topography and surface characteristics. Additional sensitivity tests on the time step and grid spacing did not change the trend of the results. Since the flow equations in DHM do not contain any turbulence terms, it is very likely that the value of the above parameters or its variants will have an impact on the variables of interest in equations that contain turbulence terms.
Figure 7. Effect of varying roughness values on the depth profile at node 14 in Figure 1(b). The roughness values shown are those of (earth, concrete).
A key component of CFD modeling algorithms is the way they incorporate turbulence terms. These terms are primarily, the fluctuations in pressure and velocity. Both Direct Numerical Simulation (DNS) and Large Eddy Simulation (LES) algorithms can well capture these as long as the mesh size is very fine, which will require more computational resources. To address this, Reynolds Averaged Navier-Stokes (RANS) approach, which is based on averaging the flow equations, is being widely used for engineering applications. However, flow turbulence effects can be reasonably modeled by calibration of the friction factor. Additionally, other flow effects can be similarly accommodated, such as flow debris among other factors. In the absence of experimental data for similar applications, the DHM results were compared across various time step, roughness and spatial steps. Results comparing the DHM output with analytical and experimental data for other publications can be found in literature   .
Figure 8 is a plot of the depth and flow hydrographs at node 14 for the domain shown in Figure 1(c). The cell size is 30 ft. As the natural watercourse path near the road widens to 60 ft, we used two rows of cells. The output data for the hydrographs at node 1 and nodes 1.15 are shown. The magnitude of the flow peak value and its timing are similar for both the input conditions, to underscore the reliability of the DHM model.
In this paper, the performance of the USGS Diffusion Hydrodynamic Model for overland flow application was analyzed. The model has been applied for overland flow over a sloping domain and the focus is on predicting the transient and peak flow variables at the end of the road. The results validate the theoretical basis of the model and the confidence in the values of the predicted variables. Although DHM solves reduced flow equations, it can be reliably used for flood inundation studies. The dominant transient flow processes that govern floodplain inundation are well captured by DHM.
Figure 8. (a) Transient depth and (b) Outflow hydrographs at node 28 for domain shown in Figure 1(c) [—Inflow hydrograph only at node 1;—Inflow hydrographs at nodes 1 and 15].
 Obrien, J.S., Julien, P.Y. and Fullerton, W.T. (1993) Two-Dimensional Water Flood and Mudflow Simulation. Journal of Hydraulic Engineering, 119, 244-261.
 Antti, T. and Bruen, M. (2006) Incremental Distributed Modelling Investigation in a Small Agricultural Catchment: 1. Overland Flow with Comparison with the Unit Hydrograph Model. Hydrological Processes, 21, 80-91.
 OpenCFD-OpenFOAM-Official Home of the Open Source Computational Fluid Dynamics (CFD) Toolbox. OpenFOAM-Official Home of the Open Source Computational Fluid Dynamics (CFD) Toolbox.
 Rao, P., Hromadka, T.V., Huxley, C., Souders, D., Jordan, N., Yen, C.C., Bristow, E., Biering, C., Horton, S. and Espinosa, B. (2017) Assessment of Computer Modeling Accuracy in Floodplain Hydraulics. International Journal of Modelling and Simulation, 37, 88-95.