A computational investigation of the hydrodynamics of the Badush dam in northern Iraq

A computational analysis of the hydrodynamics of the Badush dam in Iraq is presented, which is planned to be reconstructed as a repulse dam, to prevent the Mosul city, in case of a failure of the Mosul dam. Computational Fluid Dynamics (CFD) is applied in combination with Geometric Information System (GIS) and Digital Elevation Model (DEM). In the first part of the study, a hydrologic study of a possible Mosul dam failure is performed, predicting the important parameters for a possible flooding of Mosul city. Here, a two-dimensional, depth-averaged shallow water equations are used to formulate the flow. Based on GIS and DEM, the required reservoir size and the water level of the Badush dam are predicted, for its acting as a repulse dam. Subsequently, a computational model of the reconstructed Badush dam is developed, combining the proposed construction with the local geographic topology to achieve a perfect fit. Finally, the water flow through the bottom outlets and stilling basin of the proposed dam is calculated by an unsteady, three-dimensional CFD analysis of the turbulent, free-surface flow. The CFD model is validated by comparing the predictions with measurements obtained on a physical model, where a quite satisfactory agreement is observed.


Introduction
In eighties, Mosul dam, which is a large multi-purpose dam located on the Tigris river at about 60 km Northwest of the Mosul city was constructed. Its total storage at maximum operation pool is 11.1 billion cubic meters and the maximum discharge of spillway at maximum pool is 13,000 m 3 /s. It is constructed on a foundation of soluble soils that are continuously dissolving, resulting in the formation of underground cavities and voids that place the dam under some continuing risk and thus requires a continuous grouting program to mitigate. The most significant feature regarding the dam foundation is the presence of anhydrite and gypsum bedding and the perpetual grouting program undertaken.
Since 1989, Mosul dam is considered as a floating dam because of the above-mentioned problems at its foundation. Thus, the dam is continuously under an intensive monitoring by the local and international expertise. It is obvious that a failure of the Mosul dam would have severe consequences.
For the case of a failure of the Mosul dam, it is of utmost importance to have a repulse dam between the Mosul dam and the Mosul city that can absorb the flooding water and protect the city from flooding. Such a dam should have the capability of storing a huge amount of water that would be released by the breaking Mosul dam in its upstream. As already identified by the previous studies, the already existing, partially constructed (less than 50%) Badush dam, which is situated 40 km downstream the Mosul dam (16 km northwest of Mosul) has the potential of acting as a repulse dam, provided that it is reconstructed for this purpose. The Badush dam currently has a design flood level of 250 m.a.s.l. corresponding to a storage capacity of 0.39x10 9 m 3 , which is too small. A reconstruction of the dam to a larger capacity to take the role of a repulse dam was already planned and various investigations for its reconstruction were performed [1,2]. However, a reconstruction of the dam could not be performed until now.
In the present work, the hydrodynamics of the bottom outlets and the stilling basin of the proposed construction [2] are computationally investigated. The present goal of the analysis is the validation of the computational models by comparisons of the predictions with the detailed measurements performed on the downscaled physical model [2]. A further goal, for the future studies, is to contribute to the improvements of the design, applying the validated computational models.

Modelling
The present work is a part of a study with a broader perspective, which also addresses the Mosul dam brake and the resulting flood scenarios with an accurate consideration of the terrain topology.
Remote sensing and Geographic Information Systems (GIS) are being effectively used as important tools to determine the quantitative description of morphometry of a basin [3], which provide an integrated platform for using Digital Elevation Model (DEM) data. In the present work, the package ArcGIS 10.5 [4] is used as the GIS software tool in order to process 10 m resolution of DEM for extracting related hydrological feature classes and information of the regions under study.
Although the simulation of flood in case of Mosul dam break is not the main scope of the present paper, the modelling framework of the flood simulation shall shortly be outlined here: For the flood simulation, a twodimensional approach is used, based on a depth-averaged shallow water formulation. This is done within the framework of the hydrologic prediction software HEC-RAS 5.0.3 [5], while adapting commonly employed [6] Manning's roughness coefficients [7].
The unsteady, three-dimensional flow simulations for the Badush dam are performed by the general-purpose finite-difference/finite-volume based CFD code FLOW-3D [8], which uses a special technique, based on open area and volume fractions, for modelling complex geometries on a rectangular grid. The prevailing flow is a two-phase flow of water and air. Unlike entrained twophase flows, where, e.g. one phase is dispersed in the other [9], the air and water phases are essentially separated by an interface (free-surface). In such cases, the rather state-of-the art approach would be to define a single velocity field, which is governed by the Navier-Stokes equations in the water and air regions, both, based on the local material properties in the corresponding domains, respectively, along with an adequate interface treatment.
Besides this type of modelling, which belongs to the category of volume-of-fluid (vof) [10] model, a simpler model, the so-called "one-fluid" model, is additionally provided in FLOW-3D [8] as a reasonable and efficient one, adequate for the present class of flows, which leads to a much sharper representation of the free-surface [8].
Here, the governing equations are solved only in the liquid domain, but not in the gas one. Although the governing equations are not solved for the gas (in our case air), the model has the capability of taking air entrainment (e.g. by turbulent mixing) empirically into account [11].
In turbulence modelling, advanced, high-resolution approaches such as Large Eddy Simulation (LES) [12] or Detached Eddy Simulation (DES) are being increasingly used in different classes of flow problems [13][14][15]. However, in the present work, which presents an initial effort in the field, a RANS (Reynolds Averaged Numerical Simulation) [16] approach is adopted as the turbulence modelling framework, which is being applied in the majority of engineering problems in a large spectrum of applications [17][18][19][20][21][22]. As turbulence model, the RNG (ReNormalization Group) k-ε model [23] is used, which is amended by the so-called wall-functions [8,16].
The velocity-pressure coupling is treated by a pressure-correction method [8], along with a staggeredgrid arrangement for the velocity and pressure nodes on a structured grid. The convective terms are discretized by a second-order accurate monotonicity-preserving upwind scheme [24]. In time, a second-order implicit scheme is used [8].

Mosul dam break
Although this is not the main topic of the present paper, this issue is shortly addressed, as it provides the main motivation for the studies of the Badush dam, and also dictates the boundary conditions for the modelling of the Badush dam.
DEM data with 10 m resolution is used. Landsat satellite image with 15 m resolution is used for land use/land cover (LULC) classification. For the Manning friction coefficients, the following values are used: 0.03 for Tigris river channel, 0.04 for vegetation and river banks, 0.01 for Mosul city and for 0.05 soil texture.
It is assumed that the process of the flooding will start at elevation 330 m.a.s.l. as a worst scenario. The distance between the Mosul dam to the Mosul city is considered to be 69,000 m. The approximate total area of Mosul city is 209.5 km 2 , the elevation of Mosul city is about 228 m.a.s.l. and elevation of the Tigris river bed is 210 m.a.s.l..
The two-dimensional domain is meshed by an equidistant mesh. As a result of grid-independence study, a grid with 500x500 cells is used.
The predicted flooding situation after the reach of the flood wave to Mosul city is displayed in Figure 1, where the position of the Badush dam intended as the repulse dam is also indicated.  According to the prediction, the flood wave would reach Mosul city approximately in two hours after the Mosul dam failure, and the maximum flood wave height of approximately 25 meters would be reached in the city in approximately 8 hours.  Figure 2, based on a DEM with a resolution of 10m. The important parameters of the intended Badush dam re-construction are summarized in Table 2.  [1,2]. A novelty of the present model is in the advanced methodology applied to fit the dam structure to the topology utilizing DEM, whereas surveyance based topologic data was used in the previous studies as the basis [1,2]. To this purpose, both concrete part and earthfill part of the dam structure which are saved in STL format are investigated to ensure that the dam structures are locating on the exact location that should be. The exact location of the large scaled dam was already determined as ground control points (GCPs) by locally measured data which have been collected via the process of surveying using both GPS and total station instruments. During preprocessing of the collected data, the results showed that the dam structures are coming in fit with its location on the DEM, exactly as the case in the reality; this insured that there is a perfect compatibility between designed dam structure and the DEM.
The central parts of the created virtual dam model is depicted in Figure 3. In the figure, the bottom outlets and the stilling basin are indicated which are of main concern in the present study, as well as the spillways, which play a role especially for high values of the water level (e.g. for a situation that would occur in case of Mosul dam failure). Left to the block containing the bottom outlets, stilling basin and the spillways, the hydropower unit can be seen. The number of bottom outlets, which is planned to be eight in total in the considered design (Table 2), and their opening scenarios are important for preventing possible harmful effects on the downstream side of the dam.
It is necessary to define the number of bottom outlets that can be opened as a function of the headwater elevation. For headwater equal to or lower than the normal water level (NWL) of 250 m.a.s.l (meters above sea level), all the eight can be opened (Case 1), while for the headwater attaining a level higher than the NWL (exceptional, but critical operation in the case of Mosul dam failure) not all of the bottom outlets can be totally opened at the same time (Case 2). In that case, bottom outlets need to be operated possibly in different combinations with the spillways, depending on the prevailing water level in the reservoir.
The main role of the stilling basin is energy dissipation including the spillway and low level outlets operation, to minimize any undesirable effects in the downstream including scouring. The dissipation of the kinetic energy in the stilling basin can be affected by the basic layout of the stilling basin, with respect to the prescribed lengths and elevations of the apron part and sill part, which may be designed to trigger a hydraulic jump, as well as the introduced artificial structures, roughness elements such as rows and blocks. The definition of the solution domain is depicted in Figure 4. On the free lateral boundaries of the solution domain, symmetry conditions are applied. The boundary surfaces that enclose the domain on the top are also prescribed as symmetry boundaries. At the downstream, a pressure outlet boundary is applied, which is further adjusted for local flow being sub-or supercritical [8]. At the inlet a pressure distribution corresponding to the hydrostatic pressure field is prescribed. The remaining boundaries are solid walls, along which no-slip conditions are apply that are augmented by wall-functions in conjunction with the applied two-equation turbulence model [8].
For achieving an adequate grid resolution, a grid independence study is performed. A grid with approx. 2 million cells is found to deliver sufficient grid independence with respect to the monitored value, the flow rate, which is, then, used in the subsequent calculations. A view of the used computational grid is presented in Figure 5. As can be seen in Figure 5, a structured, rectangular grid, consisting of two main blocks is used. The curved/irregular boundaries that do not follow the straight grid lines, i.e. the solid walls of the dam structure and the ground, are approximated by fractional areas and fractional volumes of the cells on interface surfaces, by the feature FAVOR of FLOW-3D [8].
The shape of the free-surface at an instant of time is presented in Figure 6, where the color represents the velocity magnitude (m/s) in Fig. 6a, and the gage static pressure (Pa) in Fig. 6b.
The predicted fields of velocity magnitude (m/s), the velocity component in the horizontal direction, the gage static pressure (Pa), the fluid depth (m) as well as the Froude number based on the depth-average velocity are displayed in Figure 7, for an instant of time, in a plane extending in the main flow direction through the middle of a bottom outlet.
As one can see in the figure, the flow abruptly accelerates at the bottom outlet, filling up the initial part, but forming a free surface in the downstream. The thinning of the water layer in further downstream with a corresponding increase of the water velocity indicate a supercritical flow state, which is transformed into a subcritical one by a hydraulic jump in the stilling basin (Fig. 7).
Inspecting the horizontal velocity component (which is defined to be negative in the main flow direction, i.e. form right to left) recirculation patterns can be detected. Just downstream the hydraulic jump, a rather strong recirculation via a reverse flow on the surface is observed (Fig. 7b). Another major recirculation zone is observed downstream the sill.
Measurements were performed on a 1:60 physical model of the dam, which are documented in [2]. For the water level of 250 m.a.s.l., the measured (and upscaled) value of the discharge (Q), the depth-averaged velocity magnitudes at the apron (VAPRON) and sill (VSILL), as well as the water depths at the apron (hAPRON) and sill (hSILL) of the stilling basin are compared with the predicted value in Table 3.
One can see that the predictions show a reasonable agreement with the available experimental data obtained on a 1:60 scale physical model [2].

Conclusions
In the present study, the hydrodynamics characteristics of the Badush dam in northern Iraq is computationally investigated, the reconstruction of which as a repulse dam is being planned to protect the Mosul city from flooding, in case of a possible failure of the Mosul dam. Different computational tools such as CFD, GIS and DEM are used in combination. The study is performed in two stages. In the first stage, a hydrologic study of a probable Mosul dam break is performed, computationally, based on a two-dimensional, depthaveraged shallow water formulation. As result, besides the important parameters for the flooding of the Mosul city, the required reservoir size and water level of Badush dam are predicted for its reconstruction as a repulse dam.
Subsequently, a virtual model of the reconstructed Badush dam is developed, combining a previously proposed construction with the actual local geographic topology to achieve a perfect fit. The water flow through the bottom outlets and stilling basin of the proposed Badush dam are calculated by an unsteady, threedimensional CFD analysis of the turbulent, free-surface flow of water, for a given water level. A specialized version of the volume of fluid model is used to treat the two-phase flow that solves the transport equations in the water domain only, and the turbulence is modeled by the renormalization group theory version of the k-ε turbulence model. The applied CFD model is validated by comparing the predictions with the measurements obtained on a 1:60 physical model. A quite satisfactory agreement is observed between the computational predictions and the experimental analysis.
Application of the developed and validated model for design optimizations of the Badush dam is planned as future work.