Application of Diffusion Hydrodynamic Model for Overland Flows

Theodore V. Hromadka II¹, Prasada Rao^2*

¹Department of Engineering-Mathematics, United States Military Academy, West Point, NY, USA.
²Department of Civil and Environmental Engineering, California State University, Fullerton, CA, USA.

DOI: 10.4236/ojfd.2019.94022                

Abstract

Modeling two-dimensional overland flow across complex real-world topography is a challenging problem. Predicting the overland flow variables for various whatif rainfall scenarios can facilitate designing water infrastructure components aimed at preventing inundation and urban flooding. Numerical models that are being used range from those that solve the simplified St. Venant equations to CFD models that solve the complete three dimensional Navier-Stokes equations. In this work, the performance of the USGS Diffusion Hydrodynamic Model (DHM) for a series of overland flow test problems was tested by comparing numerical solutions obtained for an event-driven simulation across various sensitive parameter combinations. The reliability of the model and its ability to incorporate various topographical characteristics in the domain are illustrated.

Keywords

Two-Dimensional Flow, Numerical Modeling, Natural Watercourse, Explicit Scheme, Hydrograph

1. Introduction

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 [1] [2] 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 [3] 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 [4]. 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 [5], TUFLOW [6], SWMM [7], HEC-HMS [8], MIKE URBAN [9], FLO-2D [10] and XPSWMM [11]. 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 [12]. DHM code is simple to use and as documented in [3], 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 [3]

$\frac{\partial q_x}{\partial x}+\frac{\partial q_y}{\partial y}+\frac{\partial H}{\partial t}=0$ (1)

$\frac{\partial q_x}{\partial t}+\frac{\partial }{\partial x}\left[\frac{q_x^2}{h}\right]+\frac{\partial }{\partial y}\left[\frac{q_x{q}_y}{h}\right]+gh\left[S_{fx}+\frac{\partial H}{\partial X}\right]=0$ (2)

$\frac{\partial q_y}{\partial t}+\frac{\partial }{\partial y}\left[\frac{q_y^2}{h}\right]+\frac{\partial }{\partial y}\left[\frac{q_x{q}_y}{h}\right]+gh\left[S_{fy}+\frac{\partial H}{\partial X}\right]=0$ (3)

in which $q_X$, $q_y$ are flow rates per unit width in the x, y—directions; $S_{fx}$, $S_{fy}$ 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

$m_Z+\left[S_{fz}+\frac{\partial H}{\partial Z}\right]=0,z=x,y$ (4)

where $m_Z$ 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:

$q_Z=\frac{1.486}{n}{h}^{5/3}{S}_{fz}^{1/2},z=x,y$ (5)

Equation (5) can be rewritten in the general case as

$q_Z=-K_Z\frac{\partial H}{\partial Z}-K_Z{m}_Z,z=x,y$ (6)

where

$K_Z=\frac{1.486}{n}\frac{h^{5/3}}{{\left|\frac{\partial H}{\partial S}+m_S\right|}^{1/2}},z=x,y$ (7)

The symbol S in Equation (7) indicates the flow direction which makes an angle of $\theta ={\tan }^{-1}\left(q_y/q_x\right)$ with the positive x-direction. By assuming the value of m to be negligible, the diffusion model can be expressed as:

$q_Z=-K_Z\frac{\partial H}{\partial Z},z=x,y$ (8)

Two-dimensional DHM is formulated by substituting Equation (8) into Equation (1)

$\frac{\partial }{\partial X}{K}_x\frac{\partial H}{\partial X}+\frac{\partial }{\partial y}{K}_y\frac{\partial H}{\partial y}=\frac{\partial H}{\partial t}$ (9)

If the momentum term groupings were retained, Equation (9) can be written as

$\frac{\partial }{\partial x}{K}_x\frac{\partial H}{\partial x}+\frac{\partial }{\partial y}{K}_y\frac{\partial H}{\partial y}+S=\frac{\partial H}{\partial t}$ (10)

where

$S=\frac{\partial }{\partial x}\left(K_x{m}_x\right)+\frac{\partial }{\partial y}\left(K_x{m}_y\right)$

and $K_x,K_y$ are also functions of $m_x,m_y$ 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 [3] 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 [13] 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 [1]. 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 [14] 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 [4]. 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

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.

5. Results

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

(a)(b)

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 [3] [13] [15].

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.

6. Conclusion

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.

(a)(b)

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

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.

References

[1]	Charles, H. (1997) Numerical Computation of Internal and External Flows. John Wiley and Sons, Chichester.
[2]	Ferziger, J.H., Peric, M. and Street, R.L. (2018) Computational Methods for Fluid Dynamics. Springer International, New York.
[3]	Hromadka, T.V. and Yen, C.-C. (1987) A Diffusion Hydrodynamic Model. U.S. Geological Survey, Denver, CO.
[4]	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. https://doi.org/10.1061/(ASCE)0733-9429(1993)119:2(244)
[5]	Flow Science. FLOW-3D User Manual. http://www.flow3d.com/
[6]	TUFLOW. https://www.tuflow.com/
[7]	Storm Water Management Model (SWMM). EPA. https://www.epa.gov/water-research/storm-water-management-model-swmm
[8]	“HMS.” HEC. https://www.hec.usace.army.mil/software/hec-hms/
[9]	“MIKE URBAN.” MIKE Powered by DHI. http://www.mikepoweredbydhi.com/products/mike-urban
[10]	“FLO-2D Software: Hydrologic and Hydraulic Modeling Software.” FLO2D Software. http://www.flo-2d.com/
[11]	“XPSWMM.” Innovyze. http://www.innovyze.com/en-us/products/xpswmm
[12]	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. https://doi.org/10.1002/hyp.6189
[13]	Hromadka, T.V. and Guymon, G.L. (1982) Nodal Domain Integration Model of One-Dimensional Advection—Diffusion. Advances in Water Resources, 5, 9-16. https://doi.org/10.1016/0309-1708(82)90022-7
[14]	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. http://openfoam.com/
[15]	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. https://doi.org/10.1080/02286203.2016.1261218