Examination of Hydrologic Computer Programs DHM and EDHM

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, free-surface incompressible fluid dynamics program), FUN3D [3] (fully unstructured-grid 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 pagehttp://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 two-dimensional 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

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

in which qX,qyare the unit flow rates along the spatial directions; Sfx,Sfyrepresents 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

where mZrepresents 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

Eq. 5 can be rewritten in the general case as

where

The symbol S in Eq. 7 indicates the flow direction which makes an angle of θ= tan−1qy/qxwith 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 Eq. 8 into Eq. 1

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

where

and Kx,Kyare also functions of mx,myrespectively.

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

3. Computational limitations in DHM

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

1. The maximum number of cells (nodes) that can be accommodated in DHM is limited to 250

2. Inflow and outflow boundary nodes are limited to 10

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 domain, warranting modifications to DHM, as discussed next.

4. Features in extended DHM (EDHM)

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 output files.

5. Application of EDHM

EDHM is qualitatively and quantitatively evaluated by comparing its results with the output from DHM, for two test simulations. The focus was to check if the EDHM solution resembles DHM output for these two cases for varying inflow and other model parameters.

Case 1: Flow in a transition.

Open channel flow through a linear contraction under the framework of two-dimensional flow is a common phenomenon and has drawn the attention of many experimental and numerical studies. The flow characteristics in the contraction depend on the Froude number at the upstream end. Flow in a contraction has acted as a benchmark simulation in investigations [12] that compared the performance of various CFD and hydraulic models. Figure 1 is the definition sketch of the test problem. The rectangular channel is 380 ft. long and 260 ft. wide. The constriction portion of the channel is 60 ft. x 60 ft. The channel length before and after constriction is 120 ft. and 200 ft., respectively. The cell size in the domain is 20 ft., and the total number of cells are 239. Figure 2 illustrates the cell numbers in the domain. The elevation of cells along the north and south boundaries was assigned a high value to physically denote that they are walls. The flow is confined within these boundaries. At the upstream end, cells 3-11 (nine cells) were specified with a constant inflow of 33.33 cfs. At the downstream end, a free outfall boundary was specified. The transient simulation was carried out until time = 1 hour.

Figure 3 plots the water depth profile for the two models along the channel centerline at time = 0.9 hours for a channel bottom Manning’s roughness value of 0.015. The close agreement of results gives confidence that the changes made to arrive at EDHM did not lead to different output data. The analysis over varying inflow discharge and bottom roughness values did not change the trend of the results. Figures 4 and 5 compare the depth profile for the roughness coefficient of 0.024 and 0.04. A similar agreement in output data (including velocity components) was observed across other longitudinal sections.

Case 2: Overland flow.

Overland flow over a hill slope generated by a rainfall event is characterized by varying hydraulic properties, roughness values, topography, and physical features in the domain. For predicting the hydraulic and hydrologic properties of flow, various models that solve a range of equations from a one-dimensional hydrodynamic equation for homogeneous place surfaces [13, 14] to 2D full non-linear shallow water Equations [15] have been applied. Figure 6 shows the flow domain. The number of cells in the domain is 56 and are 30 ft. in size. The roughness value ranged between 0.015 to 0.03Inflow hydrograph (Figure 7) was applied at cells 1 and 15. The transient model was run until time = 2 hours. The elevation drop along the north and south end of the domain is 26 ft., and the drop along the west and east boundaries is 5 ft. Cells 42 to 56 were specified as critical depth outflow nodes. Figures 8 and 9 compare the transient depth profiles from DHM and EDHM (time = 0 to 2 hours) at cells 36 and 26, respectively. The trend of the depth and velocity values at other cells also indicated the close agreement of flow data from the two models.

6. Conclusions

USGS Diffusion Hydrodynamic Model is a legacy Fortran 77 that has been widely applied for multiple one and two-dimensional flow scenarios. The computational limitations in the model which prevent its application over large domains have been addressed in this paper. The enhanced model can accommodate 9999 cells and 99 inflow and outflow nodes. Based on the analysis that was carried out and the close agreement of the results between the two models gives confidence in the reliability of the extended model. Current findings encourage future development of parallel EDHM in order to reduce the computational time.