Dam Break Simulation with HEC-RAS and OpenFOAM

A dam break is a natural disaster that can cause significant property damage and loss of life. It’s useful to identify potential flooding areas downstream in the event of a dam break. In this study both HEC-RAS [1] [2] and OpenFOAM [3] [4] are set up to simulate the inundation map downstream of the Dworshak dam in Idaho. Using the same topographical data from satellite observations, similar computational meshes are set up in both HEC-RAS and OpenFOAM. Where possible, identical or similar conditions are set up in HEC-RAS and OpenFOAM to model flooding patterns due to a dam break. The velocity of the water before reaching Ahsahka, the town located at the junction downstream from the dam, is 11.5% slower in HEC-RAS compared to OpenFOAM. The average velocity of water before reaching the end of the computational domain at Big Canyon Creek is about 20% slower in HEC-RAS compared to OpenFOAM. One notable discovery is that the water flow velocity in OpenFOAM appears to depend on the mesh resolution used in the simulation. A significant velocity difference is observed when water flows from one mesh refinement region to another mesh refinement region with a different resolution.


Introduction
Recent dam breaks in Brumadinho Brazil [8] [5], the iron range of Minnesota [6], and the state of Uttarakhand in India [7] remind us of the potential dangers of flooding due to dam breaks. HEC-RAS together with GIS [9] and satellite observations can be used to model such events. OpenFOAM provides an alternative approach to model such an event by performing a full 3 dimensional fluid dynamics simulation. In this study we apply both computational fluid 1 darrenjia@bernardsboe.com, Ridge High School, NJ USA dynamics (CFD) softwares to examine the potential flooding that could occur from a hypothetical dam break at Dworshak Dam.
Dworshak Dam ( Figure 1) is a concrete gravity dam in the United States, on the North Fork Clearwater River in Clearwater County, Idaho. With a height of 219 m, Dworshak is the third tallest dam in the United States and the tallest straight-axis concrete dam in the Western Hemisphere. On average the dam has a volumetric discharge rate about 169.9 m 3 /s. There are three downstream river stations monitoring water discharge (river station ids: 13340000, 13341550, 13341040) managed by USGS national water resources. The nominal water flow velocity downstream of the dam in Clearwater river is estimated to be less than 10 m/s on average using the river station data.
While the dam is currently intact, it's interesting to compare OpenFOAM with HEC-RAS for a hypothetical dam break simulation and evaluate the results from the solvers and the potential flooding that may occur downstream in an event where the dam breaks.
The digital elevation data is publicly available from USGS [10] at 1/3 arc second (a resolution around of 10 meters). This digital elevation data is used to generate a topography mesh in both HEC-RAS and OpenFOAM.

OpenFOAM Setup
OpenFOAM is a popular open source CFD package that supports a finite volume [3] approach to general computational fluid dynamics simulation. In addition, the finite volume of the mesh [13] is adaptive and can be refined near areas of interest to capture the details of the water flow near the rivers and populated areas.
To set up the dam break simulation, we first create a geometry of Dworshak Dam and the surrounding terrain using Blender (an open source CAD software, [14]) and the elevation data from USGS. Then a finite volume mesh with a refinement zone following the downstream terrain is created to be used in the OpenFOAM solver.

Dam Geometry
The digital elevation data of the dam and the terrain is converted into a geometry (Figure 2) through Blender. The reservoir portion of the geometry is lowered to simulate the depth of the lake. This allows for water to be placed inside the reservoir in OpenFOAM. Note that the flat geometry above the reservoir is artificially placed so that there will be a sufficient amount of water for the simulation. We only use a small portion of the digital elevation map (compare with Figure 1) near Dworshak reservoir in the simulation. The actual body of water above the dam is roughly 4,277,708,640 m 3 (by courtesy of USACE Dworshak Dam Reservoir data sheet) in volume. The artificial geometry allows more water to be stored and discharged in the Open-FOAM simulation. A trapezoidal breach ( Figure  3) with a width of 140m at the top, width of 60 meters at the bottom, and height of 30 meters is artificially created in the center of the dam.

Meshing
To utilize the 1/3 arc second resolution from the digital elevation map, different levels of mesh refinement are applied to create a mesh that matches the terrain for the simulation ( Figure  4). A refined mesh in the shape of a box that surrounds the reservoir and the dam can be seen in Figure 4 near the right border. Rivers and areas of high population are further refined with visibly higher resolutions. Finally, the dam breach and portion of river immediately downstream from the dam are given the greatest refinement. As discussed later in Section 4.3, the resolution difference between different refinement regions causes a change in water velocity in Open-FOAM's simulation.
The coarse resolution mesh covering the entire computational domain has a cell resolution of 232.04 × 211.23 × 55.5 in meters. To better capture the flooding pattern, areas of interest, e.g. rivers and potential flooding regions have increased resolutions with higher mesh refinement levels. Each refinement level increases the coarse mesh resolution by a factor of 2. The area  surrounding the reservoir with water to be discharged has a resolution of 116.02×105.62×27.8 meters. Areas with high population and the portions of river far from the dam break have a resolution of 58.01 × 52.81 × 13.9 meters. The dam breach and the portion of river immediately downstream from the dam have a resolution of 29.01 × 26.40 × 6.9 meters. The entire Open-FOAM mesh has 198139 cells in total. The resulting mesh is shown in Figure 4 where the outlines of the downstream river and flooding regions are clearly visible from the cell resolution Figure 3: Geometry of the Dam: The top has a width of around 140m, the bottom has a width of around 60 meters and the height is around 30 meters. Due to the 10m resolution of the TIF file, the curved shape of the dam, along with the inaccuracies associated with meshing in OpenFOAM, it is difficult to get more precise measurements of the size of the breach unlike in HEC-RAS which can specify the breach geometry precisely. contrast.
To model water flow on an open terrain, Open-FOAM's Inter solver is used. The Inter solver works with 2 different phases of fluids: water and air. The shape of the reservoir is approximated using 3 boxes and a cylinder as shown in figure 5. Cells within the geometries are then set to water whereas all other cells are filled with air.
The fluid in the event of a dam break would likely be muddy as opposed to pure water. To simulate the muddy conditions of the water the parameters for activated sludge [16] [17] [18] [19] are used along with the Casson transport model [22] [23] instead of the default Newtonian model, allowing us to give the water a time dependent viscosity. The water is given a density of 1000 kg/m 3 which is identical to what's used in HEC-RAS.

HEC-RAS Setup
HEC-RAS is developed by the US Army Corps of Engineers (USACE) for hydrology modelling. One of the strengths of HEC-RAS is flooding simulations associated with dam breaks. Setting up the GIS projection of the digital elevation data, meshing, and configuring the solver are streamlined in HEC-RAS, making the simulation setup much easier.
A 2D storage area (SA) is created following the shape of the reservoir ( Figure 6). The 2D storage area is used to store water to be discharged and is connected to a downstream flow area through the dam. A 2D mesh is created downstream from the dam with a resolution of 30 × 30 meters which matches the resolution of the finest resolution cells in OpenFOAM. An SA connection breach ( Figure 7) is created between the storage area and the 2D mesh in order to allow for water to flow between the 2 areas. This SA connection breach is then breached with a trapezoidal dam break with a width of 140 meters and 60 meters on the top and bottom respectively and a height of 30 meters. The dam forms in 0 seconds in order to mimic the setup in OpenFOAM since there is presently no way to specify a breach formation time in OpenFOAM such that the breach can develop gradually. The 2D mesh (Figure 8) follows the shape of the river with outlets placed once the river exceeds the bounds meshed in OpenFOAM. The simulation uses default HEC-RAS settings with the exception of using the SWE-EM (stricter momentum) equation instead of Diffusion Wave and a Manning's number of 0.03, which is within the range of Manning's number for firm soil [21].

Boundary Condition
The bottom boundary that follows the shape of the terrain is set to a rough wall. In HEC-RAS, A Manning's number of 0.03 is used to simulate wall roughness for the terrain. In OpenFOAM the wall roughness boundary condition is set to a nutURoughWallFunction, giving friction to the Figure 8: HEC-RAS is setup such that the mesh covers the river and mirrors the terrain in Open-FOAM terrain. We have chosen a roughness height of K s = 0.05, a roughness constant of C s = 1, and a roughness factor 100. The rough wall can also be set to a nutkRoughWallFunction which is not exercised in this study. Examining the documentation and source code of OpenFOAM, it's found that the roughness factor is applied to the roughness height in the boundary layer friction calculation (Equation 1 by courtesy of Open-FOAM documentation and source code). Having a roughness factor of 100 is equivalent to multiplying the roughness height by 100.
dKs + dy + = roughness factor * roughness height No transformation exists currently from Manning's number n used in HEC-RAS to roughness parameters used in OpenFOAM. The final roughness wall boundary condition for the terrain is chosen such that OpenFOAM most closely matches the distance traveled by the water modelled in HEC-RAS in the same amount of time.

Initial Condition
In order to match the simulation from HEC-RAS and OpenFOAM, both simulations start when the breach is fully formed. In HEC-RAS, this means the breach formation time is reduced to 0. Then the extension of water flow and the flooding inundation map are compared between HEC-RAS and OpenFOAM.

Solvers
OpenFOAM's solver is a full 3D solver for the transient incompressible Navier-Stokes Equations. The thermodynamics aspects of the fluid flow are not modelled by the OpenFOAM solver chosen. HEC-RAS's solver is based on the 2D shallow water equations which does not simulate the thermodynamics aspects either. Both solvers allow for a dynamic time stepping constraint by prescribed Courant numbers. This feature allows the solvers to choose the largest time step possible that makes working with the variable resolution mesh possible and efficient throughout the simulation. Both simulations are calculation demanding, which means fixed time stepping can be extremely inefficient if the fixed time steps are set to values too small.

OpenFOAM Solver
The incompressible fluid solver in OpenFOAM uses a finite volume [3] approach to solve the incompressible viscous fluid flow [11] [12]. The governing equations of the incompressible vis-cous fluids are the typical incompressible steady Navier-Stokes equations: During each time step, the pressure-velocity coupled equation is solved by decoupling the pressure and momentum fields through predictor-corrector steps. During the momentum predictor step, H matrix is solved from the momentum equation on the finite volume mesh.
Now we can start the iterative process solving for pressure P and velocity V . Start with the momentum equation, Substituting the V equation into continuity equation leads to a Poisson equation of pressure that can be solved by under relaxation method.
In the predictor step, the momentum equation is solved using initial pressure P and velocity boundary condition to find velocity V The velocity field V is then used in the Poisson equation to solve pressure P The pressure field is then used in equation ( 9) to correct V at the boundary. This is the corrector step. This process is iterated until a solution of V and P converges for the computational domain.
In this particular simulation, the k− [15] turbulence mode is enabled. The Reynolds number Re = vL ν for the water rushing down the dam far exceeds the laminar regime of air flow.

HEC-RAS Solver
There are two options for solvers in HEC-RAS, the simplified diffusion equation versus the full momentum shallow water equation. The faster simplified diffusion equation is used initially but it produced unreasonably high water depth leading the water flow downstream. The full momentum equation 14 is computationally more expensive but does produce more realistic water depth results compared with the simplified diffusion equations.
The bottom friction coefficient c f can be calculated from the Manning formula c f = n 2 g √ u 2 +v 2 R and is dependent on the Manning's coefficient n [20].

Results
Both models are able to model the downstream water flow and inundation map successfully. In this section, we discuss the main results simulated from both models.

HEC-RAS Results
The development of the flooding pattern can be seen in Figures 9 after 4 minutes, 10 after 11 minutes, and 11 after one hour. Table 1 shows the water flow extension per minute as a function of time. For example the table indicates 1079.2 meters of flooding occurs from minute 0 to minute 1 and from minute 1 to 2, the length of flooding increases by 663.8 meters. Since the flooding diverges into east and west directions, the table measures only the river's expansion in the west direction. The extension from minute 0 to 1 is measured from the base of the dam to the furthest point the water reaches along the river. The flow extension between minute 1 to minute 2 is measured from the furthest point the water reaches at minute 1 to the furthest point the water reaches at minute 2. This is repeated for the next 20 minutes as shown in the table. As expected, the water flows rapidly upon exiting the breach, gaining kinetic energy as it falls. The extension then slows as the water interacts with the terrain and flows further away from the dam, attaining a relatively constant speed of around 60 m/s as shown in Figure  19. Overall, the water flows slower compared to OpenFOAM and has less variation in velocity.

OpenFOAM Results
In order to obtain a reasonable roughness factor for the simulation, the first 60 seconds of flow are simulated for terrains using nutURough-  WallFunction in OpenFOAM (Section 2.3) with a roughness height of K s = 0.05, a roughness constant of C s = 1, and roughness factor of 1, 10, and 100, as well as a frictionless terrain, shown in Figures 12, 13, 14, and 15 respectively. The nutURoughWallFunctions of roughness factor 1 and 10 do not modify the flow significantly when compared to a frictionless terrain. The nutURoughWallFunction with a roughness factor of 100 reduces the distance traveled by 9.8% when compared to the frictionless terrain and matches HEC-RAS closely. Within the first minute, the water travels 1106 meters in Open-FOAM and 1079 meters in HEC-RAS.
The flow is simulated for 1 hour using a nutURoughWallFunction with roughness factor 100 for the terrain. Figures 16, 17, and 18  one hour respectively corresponding to the HEC-RAS Figures 9, 10, and 11. Table 2 displays the extension in the flow for each minute using the same method used for Table 1. The flow extension for the first minute is measured from the base of the dam to the furthest point the water reaches after 1 minute. The data values for the next 19 minutes measure the distance from the previous furthest point of flooding to the new furthest point of flooding.
After 13 minutes, the water slows down dramatically before speeding up again. The slowing of the water coincides with the water reach-  Figure 19.
Overall, OpenFOAM has a faster flow as shown in Figure 20.

Resolution Dependence
The flow of the water in HEC-RAS has negligible dependence on resolution. Changing the mesh resolution from 30 by 30 meters to 60 by 60 meters results in the water flowing 9.3 meters further during the first 60 seconds. This is  only 0.9% more than the length the water travelled at 30 by 30 meter resolution. In comparison, OpenFOAM's flow is very sensitive to resolution changes. In particular, water slows down or stops at the boundary of the cells. Simulating the first 60 seconds of the dam break shows that water also tends to flow faster in a lower resolution mesh compared to a higher resolution mesh. Table 3 shows the distance the water travels in the first 60 seconds for a mesh with a low reso-  lution of 67.68 × 39.76 × 13.9 meters along the river compared to a mesh with a high resolution of 29.01 × 26.40 × 6.9 for different roughness factors. For roughness factors of 1 and 10, the lower resolution mesh travels roughly 30% faster and for a roughness factor of 100, the low resolution mesh travels 14% faster.   1  1620  1230  10  1600  1230  100 1250 1100 Table 3: The distance the water travels (in meters) in the first 60 seconds for a low resolution mesh and a higher resolution mesh with varying roughness factors 2.3.

Conclusion and future work
Based on a 3D digital elevation map, the potential flooding area downstream of Dworshak dam is modelled with state of the art CFD softwares, HEC-RAS and OpenFOAM. The two simulations are set up as similar to each other as possible in initial condition, boundary condition, geometry and meshing. Despite differences in numerical schemes and the treatment of roughness boundary conditions, the overall results including flow velocity and flooding area from OpenFOAM are in good agreement with HEC-RAS. There are two no-table differences between HEC-RAS and Open-FOAM's results.
First, OpenFOAM's simulation shows significant dependences on mesh resolution: 1) The overall water flow extension and velocity increases as mesh resolution decreases (Table 3); 2) Modelled water velocity can increase or decrease when water flows through mesh regions of different resolution. Our results presented in Section 4 are based on a 30 by 30 meter resolution in OpenFOAM which shows good agreement with HEC-RAS. We suspect the resolution dependent behavior could be caused by an implementation issue in OpenFOAM.
Second, the modelled water depth in Open-FOAM significantly exceeds that modelled by HEC-RAS, evidenced by  Further work should look into the origin of OpenFOAM's dependence on mesh resolution and excessive water depth modelled. The dependence on mesh resolution in particular is con-