Byunghyun Kim1 , Taehyung Kim2 , Jinho Kim3 , Kunyeon Han4
1UC Center for Hydrologic Modeling, Department of Civil and Environmental Engineering University of California, Irvine, CA 92697, United States
2, 4School of Architecture and Civil Engineering, Kyungpook National University, Daegu, 702-701, Korea
3Department of Mechanical Engineering, Yeungnam University, Gyeongsan, 712-749, Korea
Journal of Vibroengineering, Vol. 16, Issue 3, 2014, p. 1574-1589.
Received 24 February 2014; received in revised form 5 April 2014; accepted 11 April 2014; published 15 May 2014
Copyright © 2014 JVE International Ltd. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
A two-dimensional finite volume model is developed for the unsteady, and shallow water equations on arbitrary topography. The equations are discretized on quadrilateral control volumes in an unstructured arrangement. The HLLC Riemann approximate solver is used to compute the interface fluxes and the MUSCL-Hancock scheme with the surface gradient method is employed for second-order accuracy. This study presents a new method for translation of discretization technique from a structured grid description based on the traditional (, ) duplet to an unstructured grid arrangement based on a single index, and efficiency of proposed technique for unsplit finite volume method. In addition, a simple but robust well-balanced technique between fluxes and source terms is suggested. The model is validated by comparing the predictions with analytical solutions, experimental data and field data including the following cases: steady transcritical flow over a bump, dam-break flow in an adverse slope channel and the Malpasset dam-break in France.
Keywords: shallow water equations, unsplit scheme, finite volume method, well-balanced, unstructured grid, Riemann solver, dam-break.
A sudden increase in discharge from dam or levee breaks can generate shock waves of discontinuous flow properties. Discharge discontinuity ultimately leads to water surface discontinuity. It is difficult to find exact solutions to these discontinuous flows through the linearization of non-linear equations. With the introduction of the total variation diminishing (TVD) method and an approximate Riemann solver, there has been an innovative advance in second or higher order accurate schemes that can be applied to analyzing discontinuous shock waves such as dam-break waves. However, the high resolution schemes face challenges such as an imbalance between fluxes and source terms in steady flows over an irregular bed. The imbalance is caused by the inaccurate discretization of source terms, such as the bed slope in shallow water equations. Accordingly, various numerical techniques for treating source terms in an effective way, so-called well-balanced schemes, which would make it possible to maintain the balance between fluxes and source terms under steady-state conditions, have been proposed.
Bermudez and Vázquez  introduced the C-property, which shows validity for conservation laws in still water by dividing the vector and scalar terms in the shallow water equations and discretizing the source terms by upwind schemes. LeVeque  proposed a quasi-steady wave algorithm that recomposes the fluxes by introducing an artificial discontinuity within each computational cell. However, this algorithm cannot provide accurate solutions for a steady transcritical flow with a shock. Zhou et al.  proposed a surface gradient method (SGM) for accurately treating source terms such as the bed slope. The method was extended to the SGMS , which is capable of treating bed topographies including a vertical step, and was validated by laboratory data with various bed slopes. Bradford and Sanders  introduced a new technique for handling numerical truncation errors induced by the pressure and bed slope terms in an arbitrary topography with moving lateral boundaries. Burfau et al.  solved the 1D and 2D unsteady shallow water equations using the cell-centered Roe finite volume scheme and an upwind discretization of the bed slope. Rogers et al.  proposed an algebraic technique using Roe’s approximate Riemann solver for keeping the balance between flux gradients and source terms. Valiani and Begnudelli  introduced the divergence form for bed slope source term (DFB) as a possible basis of any numerical integration method for the 1D or 2D SWE. Caleffi and Valiani  presented a new method called a finite volume well-balanced weighted essentially nonoscillatory (WENO) scheme to manage bed discontinuities. Pu et al.  proposed a surface gradient upwind method (SGUM) to combine the merit of SGM and depth gradient method (DGM) by including the source terms into the inviscid flux discretization.
This study develops a well-balanced finite volume model applicable to natural rivers using the unsplit method. The unsplit method considers the conserved variables in all directions of space simultaneously. Hence, this approach is very effective and flexible in dealing with unstructured meshes which are robust in generating computational domains with highly complex boundary shapes. To calculate the flux at the interfaces between the cells through the use of the unsplit method, information for both computational and neighboring cells should be known. To achieve this, this study introduces an efficient algorithm to automatically find neighboring cells of computational cells and to determine the order of the neighboring cells even in complex calculation domains composed of unstructured meshes. Also, this work proposes a simple technique to preserve a well-balanced solution in a dry bed and a new boundary-handling technique that can enable the effective application of an unsplit method.
The proposed model is validated through the comparison of numerical solutions and analytical solutions to steady transcritical flow over a bump. It is also applied to the dam-break flow in an adverse slope channel using laboratory data. Finally, the accuracy and reliability of the numerical model based on a Godunov-type scheme are verified with a natural topography through the simulation of the Malpasset dam-break by comparing simulated results, field surveyed data and experimental measurements.
The numerical model solves the two-dimensional shallow water equations as governing equations. The two-dimensional shallow water equations can be written in matrix form as:
where is the conserved variable vector, and are the fluxes in the and directions, respectively, represents a source term, and are the velocity in the and directions, respectively, is the water depth and and are the bed slope and friction slope terms, respectively.
The governing equations are discretized based on the cell centered finite volume method. By integrating Eq. (1) over an arbitrary computational cell, the basic equation of finite volume method can be obtained as:
where is the boundary of the cell, is the outward unit vector that is normal to the boundary , and is the normal flux vector and can be determined by using the rotational invariance property :
where is the angle between the outward unit vector and the -axis, represent the transformed conserved variables, is the transformed normal flux, and and are the transformation matrix and its reverse, respectively:
therefore, Eq. (3) can be discretized as :
where is the area of an integration cell, is the number of cell interfaces, is the index of the cell interface, and is the width of each interface.
The HLL approximate Riemann solver  cannot take into account intermediate waves, such as shear waves or contact discontinuities . To correct the defects of the HLL approach, Toro et al.  proposed a modified method called the HLLC scheme, in which C stands for Contact. The HLLC scheme provides results that satisfy the entropy condition . If wave speeds are properly calculated and applied, the HLLC scheme can easily handle dry beds as well . It also provides successful approximate solutions to practical applications .
Fig. 1. Wave structure of the HLLC approximate Riemann problem
Fig. 1 shows the structure of waves assumed in the HLLC approximate Riemann solver. The conserved variable is divided into four parts according to wave speeds, and . If the conserved variables in each region have been determined, the HLLC numerical fluxes are evaluated from :
where , , and and are the left and right Riemann states at each interface, respectively. and are the numerical fluxes in the left and right parts of the star region of the Riemann solution separated by a contact wave, respectively, and , and are the speeds of the left, contact and right waves, respectively.
To properly apply the HLLC approach to shallow water equations, it is important to calculate correct wave speeds , and . In this study, a computational cell and a neighboring cell are set as the left and right hand sides of the Riemann interface for correctly calculating wave speeds using the unsplit method (Fig. 2(b)). At the same time, the wave speeds , and are estimated by taking into account the wet and dry bed conditions. The combinations of existing relations among wave speeds and are shown in Table 1.
Table 1. Wave speeds for wet and dry beds
Computational cell ()
Neighboring cell ()
In Table 1, , , and are the left and right initial values for a local Riemann problem. The contact wave speed , water depth and () are calculated from :
The MUSCL-Hancock scheme is a second-order extension of the Godunov upwind method. This scheme reconstructs conserved variables at the interface between cells prior to the calculation of fluxes for second-order accuracy in space and predicts fluxes at for second-order accuracy in time. To achieve second-order accuracy in space and time, the two-step MUSCL-Hancock method consisting of prediction and correction procedures is employed.
In the prediction step, the reconstruction of conserved variables is implemented at the interface between cells by using the data of both computational cells and neighboring cells. The simultaneous reconstruction of conserved variables at each interface is required to predict fluxes using the unsplit method. Thus, information about a computational cell and its neighboring cells, including cell numbers and primitive variables, must be known, and the order of the neighboring cells needs to be determined for two-dimensional calculations. To obtain the information necessary for calculation, the following algorithm is introduced in this study:
1) The angles produced by the centroid coordinate of a computational cell and the coordinates of nodes consisting of a computational cell are calculated using Eq. (11) for all cells in the computational domain:
2) As shown in Fig. 2(a), the index of nodes is determined in ascending order based on the angles calculated in the previous step.
3) Two nodes from the lowest number are read counter-clock wise for a computational cell. At the same time, the node numbers for the other cells, except the computational cell, are read in the same manner. The process is repeated to find the nodes matched with the two previously read nodes.
4) The cells with two shared nodes become the computational cell and one of its neighboring cells. The order of neighboring cells is determined in the order of angles produced between the centroid and the nodes of the computational cell.
The proposed algorithm can be easily applied without any limits to the configuration of cells, even in complex flow domains such as natural rivers.
Fig. 2. Schematic diagram of topological relationships between cells a) The order of neighboring cells and b) Linear reconstruction at interface for unsplit scheme
The SGM  is coupled with the MUSCL-Hancock scheme so as to keep the balance between flux gradients and source terms. The SGM uses the water level instead of the water depth of depth gradient method (DGM) for reconstruction of conserved variables at interfaces in the prediction procedure. Accordingly, , which is the variation of , instead of the variation of is estimated as:
where is a computational cell and 1, 2, 3 and 4 represent the index of neighboring cells composing cell (Fig. 2(a)). To prevent numerical oscillations that may occur near discontinuities with the use of the second-order-accurate MUSCL-Hancock method, the TVD scheme, which limits the slopes of conserved variables to be reconstructed at the interface, is applied in this study. The Superbee limiter  is extended to the and directions for unsplit scheme application as Eq. (13):
For linear reconstruction of flow variables at interface, the method which uses the distance between the center of a computational cell and midpoint of edges is employed (Fig. 2(b)). This technique is efficient to treat uniform or nor-uniform mesh since it can capture the correct behavior of linear function within a cell . Therefore, for each interface between cells is extrapolated as:
where (1, 2, 3 and 4) is the distance from the center of cell to each midpoint of interface (Fig. 2(b)). To estimate the conserved variable needed to calculate the fluxes at the interfaces, and , both of which are conserved variables of momentum equation calculated in Eq. (14), and , which is the conserved variable of the continuity equation calculated in Eq. (15), are used :
where , and refer to the water depth, water level and average bed elevation at interfaces between a computational cell () and neighboring cells (). Audusse et al.  suggested instead of Eq. (15) to preserve nonnegative water depth. However, this suggestion can cause imbalance between fluxes and sources terms in dry bed . In this study, a simple technique propose to preserve nonnegative water depth and the well-balanced solutions in a dry bed application as Eq. (16):
where is the counterpart number of (e.g., 1, 3).
The prediction fluxes and are obtained by substituting the conserved variables at interfaces into Eq. (2). For second-order accuracy in time, a conserved variable in is computed as:
where is the area of the cell , is the computational time step, is the number of the neighboring cells of cell , and are the and components of the normal outward unit vector of the interface ( and is the width of the interface.
In the correction step, a conserved value at the time level is achieved by solving a series of Riemann problems at the interface based on data from the prediction step [3, 19]:
where the correction fluxes and in the and directions are computed using the HLLC approximate Riemann solver.
In developing a numerical model to compute shallow water flows with an irregular bed topography, the numerical treatment of source terms largely affects the accuracy of the calculated results [17, 22-23]. The SGM is accurate and effective in treating source terms and can be easily coupled with the MUSCL-Hancock scheme . Accordingly, this study improved the accuracy and stability of the developed model by combining the two schemes. The source terms consist of two parts, including the bed slope and the friction slope.
For the discretization of the bed slope terms, the application of the divergence theorem yields the rate of change of bottom elevation in the and directions of cell as follows :
where subscripts 1, 2, 3 and 4 represent the index of nodes composing the computational cell, as shown in Fig. 2(b). The friction slopes in the and directions are estimated using Manning’s formula with the wide channel assumption :
where is Manning coefficient, and the primitive variables , and are estimated using the average values over the computational cell .
Fig. 3. Assignment of ghost cell numbers for boundary condition treatment
The finite volume method requires ghost cells for the treatment of the boundary condition. Fig. 3 depicts the computational domain and ghost cells. In this study, the numbers to be assigned to each ghost cell are determined using boundary cell numbers as shown in Fig. 3. This technique delivers an efficient treatment of complex topographies at the boundary. For general shallow water problems, two types of boundary conditions including the transmissive boundary condition and the reflective boundary condition are taken into account without considering the flow regime . The Eq. (21) and (22) are used for a transmissive, and a reflective boundary condition, respectively:
where the subscripts and refer to the left hand side and right hand side of the Riemann states at the interface between boundary cells and ghost cells, respectively.
The developed model is applied to various test cases to assess the efficiency and accuracy of the algorithm proposed in this study. A steady flow over a bump and a dam-break flow in an adverse slope channel with flow conditions including steady transcritical flow, hydraulic jump, wet and dry fronts and reverse flow are considered, and the calculated results are compared with analytical solutions and experimental data. Additionally, the numerical model is applied to the Malpasset dam-break in France, and the application results are compared with field surveyed data and experimental data.
If boundary conditions remain constant in a computational domain in which stationary flow is taken as the initial state, the flow will become steady-state despite the presence of a bed slope. However, nonphysical oscillations occur around the discontinuity of the bed slope if the source terms are not treated appropriately. Simulations investigating the convergence to a steady state are generally used to verify the balance between the flux gradients and source terms of a numerical model. In order to demonstrate convergence to a steady state, this study takes into account a channel proposed by Goutal and Manual . The channel is frictionless and rectangular, 25 m in length and 1 m in width. The bed elevation is defined as:
According to the initial and boundary conditions, the steady flow over a bump can be subcritical or transcritical flow with or without a shock . Of the three cases, subcritical flow and steady transcritical flow with a shock are selected for model verification, because the two cases are considered most challenging for a numerical model. We refer to Goutal and Maurel  for analytical solutions.
For the case of a subcritical flow, a constant discharge per unit width of 4.42 m2/s is specified at the upstream boundary and a fixed water level of 2 m is imposed at the downstream boundary. 2 m is also used as an initial value for the entire domain. For the computational domain, 100 cells with 0.25 m are composed. The Courant number is 0.9, and the simulation time is 360 seconds. Fig. 4 indicates that the model predictions are in good agreement with the analytical solutions for water level and discharge per unit width. Deviations between the simulation results and the analytical solutions are similar to other researchers’ results [3, 26-27].
For the transcritical flow with a shock, a discharge per unit width of 0.18 m2/s is imposed at the upstream boundary and a water level of 0.33 m is fixed at the downstream boundary. The initial water level is 0.33 m throughout the entire computational domain, in which case the hydraulic jump occurs at the end of the bump.
Fig. 4. Steady subcritical flow over a bump a) Water level and b) Discharge per unit width
To investigate accuracy dependence on cell size, simulations for 100 and 200 cells were carried out, respectively. Fig. 5 shows a comparison between the numerical predictions and analytical solutions for the water level and discharge per unit width estimated after 360 seconds. In both cases, numerical errors occur at the end of the bump. When the number of cells is 100, discharge near the point where the hydraulic jump occurs cannot meet steady-state conditions. On the other hand, when the number of cells is increased to 200, the numerical prediction is shown to have good agreement with analytical solutions. The maximum deviation between the discharge predictions and the analytical solutions is 15 % and 2 % for 100 and 200 cells, respectively. It is clearly shown that the error between the numerical solution and analytical solutions decreases as the number of cells increases.
Fig. 5. Steady transcritical flow over a bump with a shock a) Water level and b) Discharge per unit width
The simulated results of the steady flow over a bump show close agreement with analytical solutions without spurious oscillations. This suggests that the unsplit method applied to the proposed algorithm in this study and the source-term treatment technique combined with the SGM maintains a numerical balance between fluxes and source terms.
In order to assess the efficiency and accuracy of the source term treatment with the SGM in the developed model, a case with the DGM is employed for comparison. The SGM is the same as the DGM in procedure, except that it uses the water level instead of the water depth for the reconstruction of conserved variables at the interfaces between cells. The SGM and DGM are applied to the previous two test cases of steady flow over a bump, respectively. In those applications, 200 cells are used. Fig. 6 and 7 show a comparison of the unit discharge predictions of the SGM and DGM to subcritical flow and transcritical flow with a shock, respectively. In the comparison results, the numerical treatment of source terms using the DGM resulted in far larger numerical errors than that of source terms using SGM. This occurs because the DGM does not reflect the real variation of water depth . As a result, the DGM cannot accurately reconstruct water depth at the interface. Thus, the incorrect reconstruction of water depth affected the entire calculation by causing an inaccurate estimation of fluxes at the interface. This demonstrates that the application of the SGM for the treatment of source terms can improve the accuracy and efficiency of the developed model.
Fig. 6. Comparison of the DGM and SGM using the proposed model for steady subcritical flow over a bump a) DGM and b) SGM
Fig. 7. Comparison of the DGM and SGM using the proposed mode for steady transcritical flow over a bump with a shock a) DGM and b) SGM
In this section, the accuracy of the proposed algorithm is verified through simulating dam-break experiments designed by Aureli et al. . Laboratory experiments set up with different upstream depths and adverse slopes were intended to induce shock origination and propagation, a reverse flow and wetting and drying conditions. As shown in Fig. 8, the experimental domain was composed of a rectangular reservoir and a channel with an adverse slope connected to it. The length and width of this channel were set at 7 m and 1 m, respectively. To apply the numerical model, a computational domain consisting of 140 cells with 0.05 m and 1.0 m was composed. The Manning roughness coefficient is set to 0.010 m-1/3s, and the reservoir has an initial water depth of 0.25 m. The initial water depth downstream of the dam is zero for the dry bed test, and the total simulation time is set to 40 seconds. The side walls and upstream region of the channel have closed boundary conditions, and the downstream end has transmissive boundary conditions. Under these same conditions, Aureli et al.  and Tseng  validated a numerical model using the TVD-Maccormack scheme and the flux-limited TVD scheme, respectively. Fig. 8 shows the experimental configurations and simulation conditions. The bed topography is defined as:
Fig. 8. Schematic view of the dam-break experiment in an adverse slope channel
Fig. 9 shows a comparison of the calculated water depth and measured water depth at four different points. The simulated results for predicting the propagation time of the forward and reverse waves, wet/dry fronts and water depth are in close agreement with the observed values. From these simulated results, it can be demonstrated that a new method for a linear reconstruction of conserved variables in effectively applying an unsplit scheme and a simple technique to preserve balance between fluxed and source terms in a dry bed can handle reverse flow and wetting and drying beds very well.
Fig. 9. Dam-break flow in an adverse slope channel: comparison of water depths at a) G1, b) G2, c) G3 and d) G4
To verify the applicability of the developed model to a natural river with an irregular bed elevation and dry terrain, it is applied to the Malpasset dam-break simulation. The dam was located in the valley of the Reyran River about 12 km upstream of the Frejus city in France and was constructed to provide irrigation and drinking water for that city in 1954. It was an arch dam with a maximum height of 66.5 m, a crest length of 223 m and a storage capacity of 55×106 m3. It collapsed due to a sudden rise of water level in the reservoir resulting from exceptional rainfall on December 2, 1959, and the resulting flood waves were propagated along the Reyran valley to Frejus, leaving 421 casualties . This dam-break case was simulated by other researchers [30-33] to verify their numerical model.
The high water marks of flood waves resulting from the dam-break were surveyed by the police along the left and right banks of the Reyran River. Fig. 10 shows the surveyed points, including P1-P17. In 1964, a laboratory experiment on the Malpasset dam-break was conducted using a physical model scaled at 1:400, made by LNH-EDF. In the experiment, the arrival time of the flood wave and the maximum free surface elevations were measured in S6-S14 gauges as shown in Fig. 10 . These data are used to validate the numerical model in this study. The boundary condition for the upstream region is an imposed reflective boundary condition without any inflow into the reservoir because it was very small in comparison to the discharge resulting from the dam-break . The downstream boundary conditions are set as a free flow which flood waves run into the sea. For the initial conditions, the water level of the reservoir is imposed at 100 m, and the downstream region of the dam is considered as a dry bed. The 13,550 quadrilateral meshes composing the computational domain are used and Manning coefficient is specified at 0.033 m-1/3s for the entire area.
In Fig. 11, the predicted arrival time of the flood front and the maximum water surface elevation are compared with the measurements of the physical model. Fig. 12 shows a comparison of the calculated maximum water surface elevation with field data surveyed in the left and right banks. The predictions in the proposed model demonstrate relatively good agreement with the field and experimental data. However, there are small discrepancies between the simulation results and high water marks at surveyed points, P4, P5, P7, P8, P13 and P16. It should be noted that, in the simulation, bottom elevation data digitized from maps dating back to 1931 are used . However, the field data were collected after the dam-break, which was accompanied by abrupt changes in topography. This may be the reason for differences between the numerical results and the high water marks.
Fig. 10. Topography and locations of field surveyed points (P1-P17) and gauges at physical model (S6-S14)
Fig. 11. Comparison between predicted results and experimental data a) Propagation time at gauges, b) Maximum water level at gauges
Fig. 12. Comparison between predicted results and field data a) Maximum water level at right banks, b) Maximum water level at left banks
Fig. 13 illustrates the three-dimensional flood inundation extent and water depth distributions at 0, 3, 12, 21, 30 and 35 minutes after the dam-break. As shown in Fig. 13, the main channel in the center of the river valley is included in the topography and flood waves due to dam-failure are rapidly propagated along the steep and narrow valley, starting at the dam.
Fig. 13. 3D inundation water depth and extent at a) 0, b) 3, c) 12, d) 21, e) 30 and f) 35 min
However, as the flood waves travel downstream, where the floodplain becomes wider, they are propagated more slowly and the water depth decrease. Additionally, the wave front is also propagated over a wider area. The water level, which is initially set to 100 m in the reservoir, dramatically drops for 35 minutes after the dam-break. Fig. 13 shows that the main characteristics of the flood resulting from the Malpasset dam-break can be well reproduced by the proposed model. It also reveals that the techniques presented in this study can be simply and efficiently applied to irregular topographies such as natural channels.
This study developed a well-balanced finite volume model that can be applied to irregular boundaries and complex bed topographies. The proposed model is based on the conservative form of the two-dimensional shallow water equations, and the HLLC approximate Riemann solver with MUSCL-Hancock scheme solves the equations by applying the unsplit method to achieve geometrical flexibility. An algorithm for finding neighboring cells to a computational cell for the efficient application of the unsplit method and a new scheme for handling ghost cells for the treatment of complex boundaries are proposed. In addition, a simple and robust technique to preserve nonnegative water depth and well-balanced solutions in a dry bed is suggested. The proposed model is validated through the comparison of predictions with analytical solutions, experimental data and field data. In this study, three test cases are considered. One is a steady flow over a bump, another is dam-break flow with an adverse slope and the last one is the Malpasset dam-break in France. Simulation results for the first two cases show that a balance between the flux gradient and source terms can be reached with the developed model. In addition, the proposed model gives accurate solutions to complex flows, including the interactions of multiple flows, reverse flow, dry/wet shock fronts and hydraulic jumps. In the Malpasset dam-break simulation, the developed model is capable of providing accurate predictions for the arrival time and the maximum water surface elevation of floods occurring in natural rivers. From the reasonable results of several test cases chosen for a numerical verification, it is demonstrated that the proposed model can provide high accuracy, flexibility and efficiency in handling various types of shallow water flows with arbitrary topography.