**Abstract**

The Diffusion Hydrodynamic Model or DHM is a coupled one- and two-dimensional (2-D) surface flow model based upon a diffusion formulation of the well-known Navier–Stokes equations, developed by research hydrologists of the USGS (United States Geological Survey) for use in modeling floodplains and dambreak situations. The Fortran 77 source code and various applications were published in 1987 by the USGS as a Technical Report authored by Hromadka and Yen. The DHM program led to the development of several subsequent computational programs such as the FLO-2D computational model and other similar programs. The original DHM program had a limit of applications to problems with no more than 250 nodes and modeling grids. That limitation was recently removed by a program version named EDHM (Extended DHM), which provides for 9999 nodes and grids. However, the computational code is preserved in order that the baseline code algorithmic procedures are untouched. In this paper, the DHM and EDHM are rigorously compared and examined to identify any variations between the two Fortran codes. It is concluded from this investigation that the two sets of algorithm codes are identical, and outcomes from either program are similar for appropriately sized applications.

**Keywords:** legacy Fortran codes, computational economy, large scale application, overland flow, flow through a constriction

### **1. Introduction**

Legacy Fortran 77 codes that have been developed in the 1980s continue to have a wide spread audience (for both research and commercial applications) across all the Computational Fluid Dynamics disciplines. Their popularity can be largely attributed to the extensive validation that these models have been subjected to with analytical and experimental data for single and multi-dimensional flows. The trust that the audience have in the end results from these legacy codes, and their ability to meet the user goals are other driving factors for making them popular among the modeling community. Few of the legacy Fortran 77 codes, developed by various groups include MFIX [1] (Open source software for simulating Multiphase Flow with Interphase eXchanges), VOF 2D [2] (Two- dimensional, transient, freesurface incompressible fluid dynamics program), FUN3D [3] (fully unstructuredgrid fluid dynamic simulations spanning incompressible flow to transonic flow), LAURA [4] (Langley Aerothermodynamic Upwind Relaxation Algorithm for

structured, multi-block, computational aerothermodynamic simulations) and INS3D [5] (incompressible Navier–Stokes equations in three dimensions for steady and transient flows).

Although theoretically, these legacy codes are on a firm footing, computationally, they are uneconomical. When applied over a large scale application, characterized by thousands of nodes, these codes are constrained by varying degrees. The limitations arise partly from (a) the lack of object-oriented tools and the absence of abstract modeling capabilities in Fortran 77 and (b) the required CPU time when these models are applied across large domains. Balancing the accuracy of simulation with acceptable CPU is a crucial element that the current modelers are looking for. Capturing the physics of some flow phenomena necessitates that the equations be applied across the small spatial grid and temporal scales, which can be an issue in applying legacy codes. Since the codes were written (at that time) for the then available computational resources, applying them, as they are, for large scale domains might not be feasible either due to the large CPU time that they need to complete the simulation or because of the array limitations or modifications that need to be made so that the codes can be compiled using the currently available Fortran 77 compilers. Addressing these limitations will result in a 'modernized' version of the legacy codes. Modernization does not mean a better numerical formulation or a more accurate code. It only refers to a computationally efficient code with perhaps a better user interface for input and for visualizing the results through colorful multi-dimensional graphs and tables. In fact, using any modernized code without extensive benchmark testing can lead to erroneous solutions.

The above first limitation was addressed by researchers either by rewriting the entire code from scratch using an object-oriented programming language or by using an incremental approach in which the computationally intensive modules in legacy codes are identified and replaced with their computationally efficient counterparts. For instance, in an implicit finite difference or finite element formulation, where the system of equations are assembled in the form Ax = B, much of the computational time is spent on solving for 'x' [6]. While for the application audience, a solution module merely represents a means to an end of solving the flow equations, for a solver developer, the application is a source of sparse equations to be solved. Using an appropriate solver and a preconditioner can significantly cut down the simulation time, thus partly addressing the limitation.

The advent of high-performance computing tools and their availability to all the audience (based on commodity processors) saw the evolution of legacy Fortran codes to their full or semi-parallel versions, thus addressing the second limitation in legacy codes. There are many advantages to using parallel codes. For an application modeler using a well-written parallel code, the primary advantage is the reduction in the computational time. This reduction is on the order of the number of processors used in the communicator. Parallelization allows a program to be executed by as many processors available within the sub-complex simultaneously, thus facilitating a steep reduction in the computational time required for the solution to converge to the desired simulation time. This reduction enables the code to be applied over domains with millions of nodes or cells, to better capture the physics of the flow. Using standard parallel libraries that use MPI protocols like PETSc [7], ScaLAPACK [8], and BLAS [9], existing serial codes can be converted to parallel codes with reduced effort. Alternatively, highly optimized parallel versions of legacy codes can be written from scratch, which involves higher costs for developing and testing the code.

Diffusion Hydrodynamic Model (DHM) is a legacy model developed in the 1980s for USGS by the first author and his colleague [10]. The report and the Fortran 77 source code are also available at the DHM companion web page http://

**41**

*Examination of Hydrologic Computer Programs DHM and EDHM*

are discussed in Section 5. Conclusions are presented in Section 6.

field and hydrostatic pressure distribution) can be written as [10].

2

2

h h

h h

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

> <sup>x</sup> <sup>y</sup> <sup>0</sup> <sup>∂</sup> <sup>∂</sup> <sup>∂</sup> ++= ∂∂ ∂ *H q q txy*

x x <sup>q</sup> q qx y S 0

y y q qqx y S 0

<sup>∂</sup> ∂∂ ∂ + + + += ∂∂ ∂ <sup>∂</sup> *fy*

in which , *qX qy* are the unit flow rates along the spatial directions; , *Sfx Sfy* represents friction slopes; and h, H, h, g denote flow depth, water surface elevation,

The local and convective acceleration terms can be grouped and Eqs. 2 and 3 are

m S 0,z x,y <sup>Z</sup> ∂ + + == <sup>∂</sup> *fz H*

where mZ represents the sum of the first three terms in Eqs. 1,2 divided by gh. Using Manning's formula to calculate the frictional slope, the flow equation can be

*<sup>q</sup> <sup>H</sup> gh tx y <sup>X</sup>* (2)

*<sup>q</sup> <sup>H</sup> gh ty y <sup>X</sup>* (3)

*<sup>Z</sup>* (4)

<sup>∂</sup> ∂∂ ∂ + + + += ∂∂ ∂ <sup>∂</sup> *fx*

(1)

diffusionhydrodynamicmodel.com/*.* The model was extensively tested in the 1980s for various free surface flow scenarios, and its back engine has laid foundation blocks for other popular models like FLO-2D [11]. The DHM was originally developed as the first (or one of the first) 3-D CFD computational programs but was subsequently revised into the form published due to computer limitations of the day. In the last 30 years, DHM has proven to be a practical and reliable tool for predicting two-dimensional surface flow characteristics associated with gradually varied flows and is popular among the hydraulic community. DHM solves a simplified twodimensional diffusion wave equation, the solution of which is sufficient for many free surface flows that are commonly encountered. For gradually varied flows, the model predicts the values of flow depth and velocity. The model does include any turbulence terms, and it cannot be used for rapidly varying flows. In this work, we address the computational limitations in DHM so that it can be applied over larger domains. The modified DHM, herein, is referred to as Enhanced DHM (EDHM). The layout of this document is as follows. In Section 2, the flow equations and other salient characteristics in DHM are briefly described. Section 3 lists the computational limitations in DHM. In Section 4, the modifications done in DHM to arrive at EDHM are detailed. Performance tests to validate the reliability of EDHM

*DOI: http://dx.doi.org/10.5772/intechopen.94283*

**2. Overview of DHM**

and gravity, respectively.

rewritten as

simplified to

#### *Examination of Hydrologic Computer Programs DHM and EDHM DOI: http://dx.doi.org/10.5772/intechopen.94283*

diffusionhydrodynamicmodel.com/*.* The model was extensively tested in the 1980s for various free surface flow scenarios, and its back engine has laid foundation blocks for other popular models like FLO-2D [11]. The DHM was originally developed as the first (or one of the first) 3-D CFD computational programs but was subsequently revised into the form published due to computer limitations of the day. In the last 30 years, DHM has proven to be a practical and reliable tool for predicting two-dimensional surface flow characteristics associated with gradually varied flows and is popular among the hydraulic community. DHM solves a simplified twodimensional diffusion wave equation, the solution of which is sufficient for many free surface flows that are commonly encountered. For gradually varied flows, the model predicts the values of flow depth and velocity. The model does include any turbulence terms, and it cannot be used for rapidly varying flows. In this work, we address the computational limitations in DHM so that it can be applied over larger domains. The modified DHM, herein, is referred to as Enhanced DHM (EDHM).

The layout of this document is as follows. In Section 2, the flow equations and other salient characteristics in DHM are briefly described. Section 3 lists the computational limitations in DHM. In Section 4, the modifications done in DHM to arrive at EDHM are detailed. Performance tests to validate the reliability of EDHM are discussed in Section 5. Conclusions are presented in Section 6.

### **2. Overview of DHM**

*Hydrology*

and transient flows).

structured, multi-block, computational aerothermodynamic simulations) and INS3D [5] (incompressible Navier–Stokes equations in three dimensions for steady

extensive benchmark testing can lead to erroneous solutions.

down the simulation time, thus partly addressing the limitation.

Although theoretically, these legacy codes are on a firm footing, computationally, they are uneconomical. When applied over a large scale application, characterized by thousands of nodes, these codes are constrained by varying degrees. The limitations arise partly from (a) the lack of object-oriented tools and the absence of abstract modeling capabilities in Fortran 77 and (b) the required CPU time when these models are applied across large domains. Balancing the accuracy of simulation with acceptable CPU is a crucial element that the current modelers are looking for. Capturing the physics of some flow phenomena necessitates that the equations be applied across the small spatial grid and temporal scales, which can be an issue in applying legacy codes. Since the codes were written (at that time) for the then available computational resources, applying them, as they are, for large scale domains might not be feasible either due to the large CPU time that they need to complete the simulation or because of the array limitations or modifications that need to be made so that the codes can be compiled using the currently available Fortran 77 compilers. Addressing these limitations will result in a 'modernized' version of the legacy codes. Modernization does not mean a better numerical formulation or a more accurate code. It only refers to a computationally efficient code with perhaps a better user interface for input and for visualizing the results through colorful multi-dimensional graphs and tables. In fact, using any modernized code without

The above first limitation was addressed by researchers either by rewriting the entire code from scratch using an object-oriented programming language or by using an incremental approach in which the computationally intensive modules in legacy codes are identified and replaced with their computationally efficient counterparts. For instance, in an implicit finite difference or finite element formulation, where the system of equations are assembled in the form Ax = B, much of the computational time is spent on solving for 'x' [6]. While for the application audience, a solution module merely represents a means to an end of solving the flow equations, for a solver developer, the application is a source of sparse equations to be solved. Using an appropriate solver and a preconditioner can significantly cut

The advent of high-performance computing tools and their availability to all the audience (based on commodity processors) saw the evolution of legacy Fortran codes to their full or semi-parallel versions, thus addressing the second limitation in legacy codes. There are many advantages to using parallel codes. For an application modeler using a well-written parallel code, the primary advantage is the reduction in the computational time. This reduction is on the order of the number of processors used in the communicator. Parallelization allows a program to be executed by as many processors available within the sub-complex simultaneously, thus facilitating a steep reduction in the computational time required for the solution to converge to the desired simulation time. This reduction enables the code to be applied over domains with millions of nodes or cells, to better capture the physics of the flow. Using standard parallel libraries that use MPI protocols like PETSc [7], ScaLAPACK [8], and BLAS [9], existing serial codes can be converted to parallel codes with reduced effort. Alternatively, highly optimized parallel versions of legacy codes can be written from scratch, which involves higher costs for devel-

Diffusion Hydrodynamic Model (DHM) is a legacy model developed in the 1980s for USGS by the first author and his colleague [10]. The report and the Fortran 77 source code are also available at the DHM companion web page http://

**40**

oping and testing the code.

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

$$\frac{\partial H}{\partial t} + \frac{\partial q\_x}{\partial x} + \frac{\partial q\_y}{\partial y} = 0 \tag{1}$$

$$\frac{\partial \mathbf{q}\_{\mathbf{x}}}{\partial t} + \frac{\partial}{\partial \mathbf{x}} \left[ \frac{\mathbf{q}\_{\mathbf{x}}}{\mathbf{h}}^2 \right] + \frac{\partial}{\partial \mathbf{y}} \left[ \frac{\mathbf{q}\_{\mathbf{x}} \mathbf{q}\_{\mathbf{y}}}{\mathbf{h}} \right] + gh \left[ \mathbf{S}\_{\delta \mathbf{x}} + \frac{\partial H}{\partial \mathbf{X}} \right] = \mathbf{0} \tag{2}$$

$$
\frac{\partial \mathbf{q}\_y}{\partial t} + \frac{\partial}{\partial \mathbf{y}} \left[ \frac{\mathbf{q}\_y^{-2}}{\mathbf{h}} \right] + \frac{\partial}{\partial \mathbf{y}} \left[ \frac{\mathbf{q}\_x \mathbf{q}\_y}{\mathbf{h}} \right] + gh \left[ \mathbf{S}\_{fi} + \frac{\partial H}{\partial \mathbf{X}} \right] = \mathbf{0} \tag{3}
$$

in which , *qX qy* are the unit flow rates along the spatial directions; , *Sfx Sfy* represents friction slopes; and h, H, h, g denote flow depth, water surface elevation, and gravity, respectively.

The local and convective acceleration terms can be grouped and Eqs. 2 and 3 are rewritten as

$$\mathbf{m}\_{\mathbf{z}} + \left[ \mathbf{S}\_{f\mathbf{z}} + \frac{\partial H}{\partial \mathbf{Z}} \right] = \mathbf{0}, \mathbf{z} = \mathbf{x}, \mathbf{y} \tag{4}$$

where mZ represents the sum of the first three terms in Eqs. 1,2 divided by gh. Using Manning's formula to calculate the frictional slope, the flow equation can be simplified to

*Hydrology*

$$\mathbf{q}\_{\mathbf{z}} = \frac{\mathbf{1}.486}{\mathbf{n}} \mathbf{h}^{5/3} \mathbf{S}\_{\text{fz}}^{1/2}, \ \mathbf{z} = \mathbf{x}, \mathbf{y} \tag{5}$$

Eq. 5 can be rewritten in the general case as

$$\mathbf{q}\_{\rm Z} = -\mathbf{K}\_{\rm Z} \frac{\partial H}{\partial \mathbf{Z}} - \mathbf{K}\_{\rm Z} \mathbf{m}\_{\rm Z}, \quad \mathbf{z} = \mathbf{x}, \mathbf{y} \tag{6}$$

where

$$\mathbf{K}\_{\mathbf{z}} = \frac{\mathbf{1}.486}{\mathbf{n}} \stackrel{\text{h}^{5/3}}{\underset{\beta\mathbf{S}}{\left|\beta\mathbf{S}\right|}} \stackrel{\text{h}^{5/3}}{\underset{\beta\mathbf{S}}{\left|\alpha\mathbf{S}\right|}} + \mathbf{m}\_{\mathbf{s}} \stackrel{\text{h}^{1/2}}{\underset{\beta\mathbf{S}}{\left|\alpha\mathbf{S}\right|}} \tag{7}$$

The symbol S in Eq. 7 indicates the flow direction which makes an angle of θ = ( ) <sup>1</sup> y x tan q / q <sup>−</sup> with the positive x-direction. By assuming the value of m to be negligible, the diffusion model can be expressed as,

$$\mathbf{q}\_{\rm Z} = -\mathbf{K}\_{\rm Z} \frac{\partial H}{\partial \mathbf{Z}}, \quad \mathbf{z} = \mathbf{x}, \mathbf{y} \tag{8}$$

Two-dimensional DHM is formulated by substituting Eq. 8 into Eq. 1

$$\frac{\partial \mathbf{H}}{\partial t} = \frac{\partial}{\partial \mathbf{X}} \mathbf{K}\_{\mathbf{x}} \frac{\partial \mathbf{H}}{\partial \mathbf{X}} + \frac{\partial}{\partial \mathbf{y}} \mathbf{K}\_{\mathbf{y}} \frac{\partial \mathbf{H}}{\partial \mathbf{y}}.\tag{9}$$

If the momentum term groupings were retained, Eq. 9 can be written as

$$\frac{\partial \mathbf{H}}{\partial t} = \frac{\partial}{\partial \mathbf{x}} \mathbf{K}\_{\mathbf{x}} \frac{\partial H}{\partial \mathbf{x}} + \frac{\partial}{\partial y} \mathbf{K}\_{\mathbf{y}} \frac{\partial H}{\partial y} + \mathbf{S} \tag{10}$$

**43**

*Examination of Hydrologic Computer Programs DHM and EDHM*

b.Inflow and outflow boundary nodes are limited to 10

domain, warranting modifications to DHM, as discussed next.

cells and the inflow hydrograph at each of them need to be inputted. The number and the outflow boundary cell numbers should be identified. Since the formulation is explicit, the choice of time step (Δt) is limited by the Courant-Friedrich-Lewy stability condition. Starting from time = 0, the explicit solution is marched in the direction of time, until the required transient time level is reached. DHM gives the option of printing the output variables at any time level. The output at the center of each cell includes the flow depth, elevation, and flow velocities along the four directions.

For modeling flows over large sizes domain, the two primary shortcomings in

a.The maximum number of cells (nodes) that can be accommodated in DHM is

Both these limitations were largely due to the computational resources that were available to the developers in the 1980's. Application of DHM over large flow domains would require using a higher number of nodes in the computational

The changes made in DHM can be grouped into three categories, as detailed below. No changes were made to the format or structure of variables input and

a.Increased the array size of variables (FP,FC) that stores the location of cells and the initial depth, average elevation, roughness coefficient values, velocity (VEL), maximum water depth, and the corresponding time at the cell (DMAX,

b.Increased the array size of variables that store the stage curve data for the

d.Increased the array size of variables relating to the inflow and outflow hydrograph nodes (KIN, KOUT), depth of the specified stage-discharge curve (HOUT), inflow boundary condition nodal points (KINP) and the inflow

e.Increased of the array which stores the nodal points where outflow hydro-

The changes in the above array sizes were done in the main code and the

c.Increased effective rainfall intensity data pairs from 10 to 99.

graphs are being printed (NODFX, NODCC) from 50 to 99

associated subroutines FLOODC, QFP, QFC, and CHANPL.

*DOI: http://dx.doi.org/10.5772/intechopen.94283*

**3. Computational limitations in DHM**

**4. Features in extended DHM (EDHM)**

DHM are

output files.

**4.1 Major enhancements**

TIMEX) from 250 to 9999.

channel (NOSTA, STA) from 10 to 99.

hydrograph details (HP) from 10 to 99

limited to 250

where

$$\mathbf{S} = \frac{\widehat{\boldsymbol{\mathcal{O}}}}{\widehat{\boldsymbol{\mathcal{O}}} \mathbf{x}} (\mathbf{K}\_{\mathbf{x}} \mathbf{m}\_{\mathbf{x}}) + \frac{\widehat{\boldsymbol{\mathcal{O}}}}{\widehat{\boldsymbol{\mathcal{O}}} \mathbf{y}} (\mathbf{K}\_{\mathbf{y}} \mathbf{m}\_{\mathbf{y}}),$$

and Kx , Ky are also functions of mx , my respectively.

To maintain continuity in the discussion, while salient aspects of the DHM numerical algorithm are presented here, readers are referred to [10] for a detailed description of the numerical formulation and the input file format. The domain is divided into uniform square grids or cells. For each interior grid, its connectivity with adjacent grids along the North, East, South, and West directions is specified. For grids that are on the boundaries, '0' is specified along the directions that do not have adjacent nodes. The flow equation is solved using the nodal domain integration method. Apart from the grid connectors, at the center of each cell, the required input variables are the roughness value, ground elevation, initial flow depth. The number of inflow

*Examination of Hydrologic Computer Programs DHM and EDHM DOI: http://dx.doi.org/10.5772/intechopen.94283*

cells and the inflow hydrograph at each of them need to be inputted. The number and the outflow boundary cell numbers should be identified. Since the formulation is explicit, the choice of time step (Δt) is limited by the Courant-Friedrich-Lewy stability condition. Starting from time = 0, the explicit solution is marched in the direction of time, until the required transient time level is reached. DHM gives the option of printing the output variables at any time level. The output at the center of each cell includes the flow depth, elevation, and flow velocities along the four directions.
