Novel Slope Source Term Treatment for Preservation of Quiescent Steady States in Shallow Water Flows

This paper proposes a robust method for modeling shallow-water flows and near shore tsunami propagation, applicable for both simple and complex geometries with uneven beds. The novel aspect of the model includes the introduction of a new method for slope source terms treatment to preserve quiescent equilibrium over uneven topographies, applicable to both structured and unstructured mesh systems with equal accuracy. Our model is based on the Godunov-type finite volume numerical approximation. Second-order spatial and temporal accuracy is achieved through high resolution gradient reconstruction and the predictor-corrector method, respectively. The approximate Riemann solver of Harten, Lax, and van Leer with contact wave restoration (HLLC) is used to compute fluxes. Comparisons of the model’s results with analytical, experimental, and published numerical solutions show that the proposed method is capable of accurately predicting experimental and real-time tsunami propagation/inundation, and dam-break flows over varying topographies.


Introduction
Shallow-water equations (SWEs) have a wide range of applications in hydraulic, coastal, and environmental engineering; these equations are extensively used to simulate various geophysical flows [1][2][3][4][5]. With appropriate numerical techniques, SWEs can model complex flows that involve sharp discontinuities and shocks. The most challenging issues in modeling shallow flows involve maintaining second-order accuracy, satisfying well-balanced condition, and avoiding numerical instabilities at wet-dry fronts [6][7][8].
A higher accuracy for any flow solver guarantees better results and reduces computational time by enabling the use of coarser grids. Some of the popular methods of establishing second-order accuracy are Monotonic Upwind Scheme for Conservation Laws (MUSCL) scheme [1,[7][8][9], the Weighted Essentially Non-Oscillatory (WENO) scheme [10][11][12], and a high-resolution reconstruction procedure developed by [13] which has been widely employed for shallow water flows [5,14,15]. Reconstruction by MUSCL scheme may result in negative water depths at the interface which consequently give unphysical high velocities because of discharge divided by small water depths [16]. Application of the WENO scheme give accurate results but it is expensive and complicated. Therefore, in this study, we establish second-order accuracy by using the multidimensional gradient reconstruction method proposed by [13]; because of its adaptability to both structured and unstructured grids and accuracy even on highly distorted grids due to dependence on a wide computational stencil.
Following the pioneering work of Bermudez and Vazquez [17] on the preservation of the conservation property (C-property) to achieve a well-balanced solution, various methods have been developed for discretization of slope source terms. The main concept of a well-balanced solution is to preserve a quiescent equilibrium over an uneven bed elevation (see details in Section 2.2). Failing to achieve balance between momentum flux and slope source terms results in spurious oscillations near the steady states. It is worth mentioning some of the classical well-balanced numerical approaches. Popular approaches include the upwind discretization of slope terms [18,19], the surface water gradient (SGM) methods in which the water surface level is used as a conservative variable instead of water depth [7,20,21]; the hydrostatic reconstruction method proposed by [22] is another popular method among researchers [16,[23][24][25]. The upwind method performs well for first-order schemes, but needs additional source term treatment for higher-order schemes, rendering this method both complicated and expensive. The surface gradient method is typically proposed for structured grids and does not offer the flexibility of being applicable to all mesh types [26]. The hydrostatic reconstruction method is applicable for both structured and unstructured meshes but for second-order schemes additional source term treatments are required [2,22]. Hou et al. [26] proposed an improvement to the existing slope source term treatment methods by offering a scheme that could be extended to high-order accuracies without complication. Their method was based on the splitting of a cell into sub-cells and approximating the flux by adding fluxes through the sub-cells. Hou et al. [16] further improved the scheme proposed in [26] through local modification of the bed elevation in slope source terms and through a better approximation of water depths at the interface by hydrostatic reconstruction.
Analysis of the existing techniques for obtaining a well-balanced solution shows that correct approximation of averaged water depth at a control volume interface is a difficult task. A compromise is often made on the accuracy, simplicity, and computational efficiency of the solution. Therefore, in this study, we propose a novel slope source term treatment method which possesses all the features of a robust numerical scheme without any compromise on accuracy, simplicity and computational efficiency. We convert the slope terms into fluxes at the interface as previously done in [16,26], but with a different mathematical formulation to average water depth in convective flux and slope source terms. Unlike [16] where the scheme needs to switch to first-order accuracy at wet-dry interfaces to avoid numerical instabilities, the proposed method, paired with the multidimensional reconstruction method, is able to maintain second-order accuracy throughout the domain including at the wet-dry interfaces. Instead of averaging water depths at the cell center and left side of the interface in the slope term discretization, this work proposes an averaging of the reconstructed water depth between the cell center and both sides of the interface. This method results in accurate flux approximation across the interface even in the case of a very abrupt bed level change between adjacent cells. This is because the discontinuity of the Riemann states at the left and right of an interface that can result during reconstruction is avoided by averaging the water depth between the cell center and both sides of interface. The proposed method of slope source term treatment results in a numerical scheme which is: (1) able to estimate accurate average water depth at cell interfaces on unstructured grids with complex geometries; (2) can maintain second-order accuracy throughout the domain; (3) suitable for both structured and unstructured mesh systems; (4) computationally inexpensive since it does not require any additional correction terms of slope source terms treatment; (5) simple and robust; and (6) preserves the quiescent equilibrium revealing that the proposed slope source terms treatment satisfies the well-balanced condition.
Wetting and drying may result in a non-conservation of mass, resulting in numerical instabilities which lead to negative or very small water depths at wet-dry interfaces and unphysically high velocities. Unphysical velocities can be satisfactorily avoided by employing a positivity preserving the hydrostatic reconstruction method to modify flow variables at the interface between two adjacent cells and setting a threshold limit on water depth [23,[26][27][28]. We set a threshold of 10 −5 m in this study to prevent occurrence of negative water depths; a cell is considered dry if water depth in a cell is less than threshold limit which sets the velocity at interfaces of a cell to zero. After the identification of dry cells, the variables at the wet-dry interface are modified by the hydrostatic reconstruction method proposed by Audusse et al. in [22].
Various approaches-each with its own merits and shortcomings-for satisfying second-order accuracy, well-balanced condition, and handling instabilities at the wet-dry interfaces are briefly discussed in this study. We integrate robust numerical techniques to address all of the above challenges and propose a new method for treating slope source terms. A number of test cases are modelled to test various aspects of the proposed scheme, including two-dimensional flow modeling abilities, shock capturing, the handling of numerical instabilities at wet-dry interfaces, wave propagation over uneven topographies, bore reflection, the flooding and draining of initially dry cells, and flow transition from subcritical to supercritical flows. Our proposed method models the benchmark test cases accurately while preserving well-balanced conditions.

Governing Equations and Numerical Model
Flow in rivers, open channels and floodplain systems can be satisfactorily represented through a set of two-dimensional SWEs. Non-linear SWEs over an arbitrary control volume can be expressed in the integral form as: where U is the vector of conserved variables, F is the flux vector, and S is a vector of source terms; Ω and Γ are the volume and boundary integrals of the control volume, respectively. In Equation (1): E and G are flux vectors in the xand y-directions, respectively, η is the water surface elevation, h is the water depth, g is the gravitational acceleration, and q x , q y and u, v are the discharge and velocity components in the xand y-axis directions, respectively. The source terms vector comprises the bed slope and friction terms: where S ox , S oy , and S f x , S f y are the bed slope and bottom frictional source terms in the xand y-directions, respectively. The source terms are estimated as follows: where z is the bed elevation, and C f = n 2 , n is the Manning's roughness coefficient. Flux in Equation (1) can also be expressed as: where n denotes the unit outward normal vector to the considered edge, defined as n x , n y T , and F(U).n is the flux vector normal to the interface. Equation (5) is summed over the interfaces of an arbitrary triangle: where M and l represent the index and length of the cell edge, respectively. To achieve second-order temporal accuracy, a two-step predictor-corrector method is used [29]. In the predictor step, the boundary-extrapolated values at each cell interface evolve by half-time steps in terms of the corresponding conserved variable [20]: where A is the cell area and L M is the cell side vector defined as the side length multiplied by the outward-pointing unit normal vector, U M is the vector of flow variables at midpoint of each interface of an element, and F(U M ) is the flux vector through the interface.
In the corrector step, the solution is updated over a full-time step based on data from the predictor step: where the flux vector F U L M , U R M is calculated by solving local Riemann problem at each cell interface. U L M and U R M are vectors of conserved variables on the left and right sides of the cell interface, respectively. The flux vector is calculated using the HLLC approximate Riemann solver which is suitable for wet-dry interface applications. Wave speeds and water depth needed to compute numerical fluxes by HLLC are approximated as [29]: if right side is dry max(u R + gh R , u * + gh * ) if right side is wet (9) and h * = 1 g where S L , S R are the respective speeds of the left-and right-moving waves, h L , h R are the respective water depths on the left and right sides of the cell interface, u L , u R are the respective velocities on the left and right sides of the cell interface, and h * , u * are the respective water depth and velocity values along the shear wave.

Variable Reconstruction and Limiting
Higher order numerical schemes guarantee increased accuracy and high resolution in regions of smooth and abrupt flows. One way to achieve a high-order scheme is through a gradient-reconstruction of the cell-center variables at the interface of the control volume. We employed a second-order piecewise linear reconstruction technique proposed by [13]; this method uses the information beyond the adjacent cells to compute left and right states for flux computation at interface thus utilizing a wide computational stencil, making it suitable for application on meshes with poor grid connectivity. The value of conserved variables stored at cell center is reconstructed at the interface by Taylor series expansion: where U i is the averaged value of conserved variables at a cell center, r i is the vector extending from the cell center to the center of the interface, and ∇U 1 i is the limited cell gradient. We used the Green-Gauss theorem to compute the cell gradient: where A is the cell area. A pseudo-Laplacian formula proposed Holmes and Connell in [30] is used to compute the values of conserved variables at vertices for gradient computation. The vertex values are given by the following equation: where v and e represent the vertex and number of cells surrounding the vertex, respectively, and ω is the weightage factor depending on the dimensions of a cell (readers are referred to reference [14] for details). The unlimited gradient at the interface between two cells is computed by taking the area-weighted average of gradients of two triangles formed by joining the ends (vertices) of the interface with the barycenter of the cell as shown in Figure 1. For interface M = 1, the gradient can be computed as: A similar procedure is adopted to find gradients at other interfaces of the cell i. The unlimited gradient of the triangular cell i is then computed by taking the area-weighted average of the gradients computed by Equation (15) at the three interfaces. The unlimited cell gradient is given as: In order to suppress the oscillations associated with high-order methods, we employ the multidimensional slope limiter proposed by Jawahar and Kamath [13]. For a two-dimensional triangular unstructured grid, the limited gradient within the cell i shown in Figure 1 is given as: where ∇U 1 , ∇U 2 and ∇U 3 are the unlimited gradients;α N , α R , and α P are the weights given as: where g 1 , g 2 , and g 3 are the functions of vertex gradients estimated by the square of L 2 norm: The explicit nature of the scheme requires the time step ∆t to be limited by the Courant-Friedrichs-Lewy (CFL) condition: where i and j represent the cell and its interface, respectively, R i is the minimum distance from the cell center to the cell vertex, and c is the wave celerity. Friction source terms are treated explicitly in the numerical scheme.
proposed source term treatment preserves the quiescent equilibrium perfectly. Equation (36) may seem similar to Equation (30) in [26], but is based on a different mathematical formulation and slope source flux computation through the interface. The proposed method of considering the reconstructed water depth at both sides of the interface in averaging water depths for source terms and convective flux serves to eliminate any chance for instabilities to arise due to a discontinuity of conserved variables at an interface caused by abrupt change in bed elevation. Our method is also able to obtain appropriate water depth averages at interfaces for both structured and unstructured meshes, which is verified through test case simulations. We demonstrated through several theoretical cases that the well-balanced conditions were perfectly preserved by the proposed slope source treatment.

Wet-Dry Front Treatment
To prevent unphysical high velocity and negative water depth at the wet-dry interface, we adopt a hydrostatic reconstruction procedure proposed in [22]. The bed elevation at the interface between two cells is chosen as: where , L R z z are the bed elevations at left and right of the interface, respectively. After selecting the higher bed elevation at the interface, the bed elevation is modified by the following step: Water depths and velocities at left and right side of the interface are then evaluated as:

Slope Term Treatment and Preservation of C-Property
Since the introduction of a well-balanced condition for the shallow-water equations in [17], a number of methods have been proposed with regard to the slope source term treatment. A scheme is regarded to be well-balanced if it satisfies the C-property or maintains quiescent equilibrium over uneven beds. The velocity and water levels for the quiescent state become: Corresponding water levels at the cell center and interface are then given as: where h o , h M , z o , z M are respective water depths and bed elevations at cell center and interface. Following [16] we take the slope source term integral over a triangular control volume shown in Figure 1 followed by expressing the integrals in form of fluxes. The slope source flux through the cell can be expressed as a summation of the flux through each interface: where F sM denotes the slope source term flux through the interface M and l M is the length of the interface M. The integral of the slope term over the cell abc can be expressed as: For demonstration, we estimate the integral over sub-cell abM; the same procedure can be applied to compute flux through other interfaces: where h 1 is the average water depth at interface M = 1, Equation (25) can then be expressed as: The slope source term flux through all the faces of cell abc can then be given as: For a two-dimensional arbitrary cell, slope source term flux can be expressed as: where h ML , h MR are the reconstructed water depths and z M is the bed elevation at interface M, estimated as: The next step is to verify whether the proposed slope source term treatment in Equation (29) balances the convective flux terms represented by E and G in Equation (2). For the motionless steady state given by Equation (21), the convective flux terms become: where F c represents the convective fluxes. We use Equation (26) to approximate h M in Equation (31): Making use of the quiescent equilibrium condition expressed in Equation (22), Equation (29) becomes: Subtracting the convective flux in Equation (32) from the slope source flux in Equation (33) yields: Since we have converted the slope source terms into flux through the interface, Equation (8) can also be written as Friction source terms are discretized separately using an explicit scheme; therefore, they are disregarded in Equation (35). Differences in the slope source and convective flux given by Equation (34) are introduced in Equation (35): (37) which is the primary requirement of a motionless steady flow. Equation (37) shows that the proposed source term treatment preserves the quiescent equilibrium perfectly. Equation (36) may seem similar to Equation (30) in [26], but is based on a different mathematical formulation and slope source flux computation through the interface. The proposed method of considering the reconstructed water depth at both sides of the interface in averaging water depths for source terms and convective flux serves to eliminate any chance for instabilities to arise due to a discontinuity of conserved variables at an interface caused by abrupt change in bed elevation. Our method is also able to obtain appropriate water depth averages at interfaces for both structured and unstructured meshes, which is verified through test case simulations. We demonstrated through several theoretical cases that the well-balanced conditions were perfectly preserved by the proposed slope source treatment.

Wet-Dry Front Treatment
To prevent unphysical high velocity and negative water depth at the wet-dry interface, we adopt a hydrostatic reconstruction procedure proposed in [22]. The bed elevation at the interface between two cells is chosen as: where z L , z R are the bed elevations at left and right of the interface, respectively. After selecting the higher bed elevation at the interface, the bed elevation is modified by the following step: Water depths and velocities at left and right side of the interface are then evaluated as:

Boundary Conditions
Boundary conditions are implemented in the ghost cells outside the computational domain at inlet and outlet for all test cases. Following the theory of characteristics, two boundary conditions among the three flow variables (η, q x , q y ) are specified for subcritical flows, and in case of supercritical flows all three flow quantities are specified in ghost cells. Depending upon the outflow conditions, transmissive and solid reflective boundary conditions are considered at the outlet.

Solitary Wave Interaction with Conical Island
Non-linear SWEs have been validated by several researchers in the coastal engineering community [4,31,32] for nearshore tsunami propagation, including the benchmark case of a solitary wave run-up over a conical island proposed by Briggs et al. [33]. The unexpectedly large tsunami heights observed at circular-shaped Babi Island and pear-shaped Okushiri Island (1993)

Boundary Conditions
Boundary conditions are implemented in the ghost cells outside the computational domain at inlet and outlet for all test cases. Following the theory of characteristics, two boundary conditions among the three flow variables ( , , x y q q  ) are specified for subcritical flows, and in case of supercritical flows all three flow quantities are specified in ghost cells. Depending upon the outflow conditions, transmissive and solid reflective boundary conditions are considered at the outlet.

Solitary Wave Interaction with Conical Island
Non-linear SWEs have been validated by several researchers in the coastal engineering community [4,31,32] for nearshore tsunami propagation, including the benchmark case of a solitary wave run-up over a conical island proposed by Briggs et al. [33]. The unexpectedly large tsunami heights observed at circular-shaped Babi Island and pear-shaped Okushiri Island (1993)    In this study, we modeled a solitary wave height of ζ = 0.032 m propagating over an initial water depth of h = 0.32 m (ζ/h = 0.1 m). The discretized domain consisted of 2577 elements; the CFL number was set to 0.6, and the Manning's roughness coefficient was taken as 0.016 for the concrete surface of the basin and island. To emulate the experimental conditions, radiation boundary conditions were applied at lateral boundaries to reduce wave reflection that could interact with the propagating wave. At x = 0, the solitary wave was generated by specifying the free surface elevation and velocity [4,34]: where c is the wave celerity equal to g(h + A). There were 27 gauges in the experimental set up to record wave positions with time. We compared the computed results with the measured results at four gauge locations in Figure 3. The model was run for 20 s because afterwards the wave motion becomes dispersive [4,33].
In this study, we modeled a solitary wave height of  = 0.032 m propagating over an initial water depth of h = 0.32 m (  /h = 0.1 m). The discretized domain consisted of 2577 elements; the CFL number was set to 0.6, and the Manning's roughness coefficient was taken as 0.016 for the concrete surface of the basin and island. To emulate the experimental conditions, radiation boundary conditions were applied at lateral boundaries to reduce wave reflection that could interact with the propagating wave. At x = 0, the solitary wave was generated by specifying the free surface elevation and velocity [4,34]: where c is the wave celerity equal to ( ) g h A  . There were 27 gauges in the experimental set up to record wave positions with time. We compared the computed results with the measured results at four gauge locations in Figure 3. The model was run for 20 s because afterwards the wave motion becomes dispersive [4,33].  A time series of water surface elevations obtained from our proposed model was compared with the measured values at four gauges; the comparisons can be seen in Figure 3. The wave arrival time was well predicted at all gauges. The predicted overall wave shape is in general agreement with the experimental observations; however, due to the non-dispersive nature of the non-linear SWEs a steeper wave is observed at gauges 6, 9 and 16. This wave front shape was also observed in prior studies that modeled tsunami propagation by SWEs [4,31]. The depression that followed the lead wave and the secondary depressions are underestimated at gauges 6, 9, and 22; this is consistent with results reported by [35]. Observations at gauges 9 and 22 show that the measured free surface elevation remains below the still water depth for a considerably longer time than the model simulation, predicting that the free surface elevation approaches the still water depth faster and maintains that profile. The leading wave at gauge 16 has a slightly wider crest that causes a delayed depression, but the overall waveform and amplitude values are in good agreement.
A comparison of the maximum predicted and measured run-ups around the island is shown in Figure 4. β is the angle measured counter-clockwise to the island, and R is the maximum run-up. We observe a trend with regard to the run-up being over-predicted for β between 0 • and 65 • (the front of the island); however, from 90 • to 120 • , the predicted results are in good agreement with the experiment until the wave reaches the back of the island where we observe a gradual increase in the under-estimation of run-up heights. The overall agreement between the computed and measured values shows that our proposed scheme could adequately model tsunami propagation and the various processes associated with it, including: wave breaking, shoaling, reflection, run-ups and run-downs. A time series of water surface elevations obtained from our proposed model was compared with the measured values at four gauges; the comparisons can be seen in Figure 3. The wave arrival time was well predicted at all gauges. The predicted overall wave shape is in general agreement with the experimental observations; however, due to the non-dispersive nature of the non-linear SWEs a steeper wave is observed at gauges 6, 9 and 16. This wave front shape was also observed in prior studies that modeled tsunami propagation by SWEs [4,31]. The depression that followed the lead wave and the secondary depressions are underestimated at gauges 6, 9, and 22; this is consistent with results reported by [35]. Observations at gauges 9 and 22 show that the measured free surface elevation remains below the still water depth for a considerably longer time than the model simulation, predicting that the free surface elevation approaches the still water depth faster and maintains that profile. The leading wave at gauge 16 has a slightly wider crest that causes a delayed depression, but the overall waveform and amplitude values are in good agreement.
A comparison of the maximum predicted and measured run-ups around the island is shown in Figure 4. β is the angle measured counter-clockwise to the island, and R is the maximum run-up. We observe a trend with regard to the run-up being over-predicted for β between 0° and 65° (the front of the island); however, from 90° to 120°, the predicted results are in good agreement with the experiment until the wave reaches the back of the island where we observe a gradual increase in the under-estimation of run-up heights. The overall agreement between the computed and measured values shows that our proposed scheme could adequately model tsunami propagation and the various processes associated with it, including: wave breaking, shoaling, reflection, run-ups and run-downs.

Dam-Break Flow in an L-Shaped Channel
We considered the flow behavior in an L-shaped channel with a 90° bend. Flow within the bend shows strong two-dimensional features, making it suitable for assessing the two-dimensional modeling capabilities of our numerical scheme. The channel was discretized into 3948 unstructured elements and 2132 nodes, shown in Figure 5. The reservoir bed was 0.33 m below the channel bed and a friction coefficient of 0.019 was assumed over the entire domain. Based on the experimental studies of [36], an initial water depth of 0.58 m was contained in the reservoir while the channel bed is dry. A dam at x = 2.39 m separated the reservoir from the channel.
At t = 0.0 s, the dam was suddenly removed and supercritical flow was induced; the simulations were run for 15 s. Propagation of the wave front within the channel, its reflection from the bend, and resulting bore propagation can be seen in Figure 6 at t = 6 and 15 s.
The wave reached the bend approximately 6 s after the dam-break, as shown in Figure 6a, followed by an upstream travelling reflection and a depression in the water surface elevation immediately downstream from the bend. The reflected bore smoothly propagated upstream without exhibiting any oscillatory behavior at t = 15 s after the dam-break, as shown in Figure 6b. Model prediction for the water depth in the channel was further compared with the experimental observations in [36] and the published numerical results of [37] in Figure 7, at t = 3, 5, 7, and 14 s. It is evident that our scheme achieves a very satisfactory agreement with the experimental observations,

Dam-Break Flow in an L-Shaped Channel
We considered the flow behavior in an L-shaped channel with a 90 • bend. Flow within the bend shows strong two-dimensional features, making it suitable for assessing the two-dimensional modeling capabilities of our numerical scheme. The channel was discretized into 3948 unstructured elements and 2132 nodes, shown in Figure 5. The reservoir bed was 0.33 m below the channel bed and a friction coefficient of 0.019 was assumed over the entire domain. Based on the experimental studies of [36], an initial water depth of 0.58 m was contained in the reservoir while the channel bed is dry. A dam at x = 2.39 m separated the reservoir from the channel.
At t = 0.0 s, the dam was suddenly removed and supercritical flow was induced; the simulations were run for 15 s. Propagation of the wave front within the channel, its reflection from the bend, and resulting bore propagation can be seen in Figure 6 at t = 6 and 15 s.
The wave reached the bend approximately 6 s after the dam-break, as shown in Figure 6a, followed by an upstream travelling reflection and a depression in the water surface elevation immediately downstream from the bend. The reflected bore smoothly propagated upstream without exhibiting any oscillatory behavior at t = 15 s after the dam-break, as shown in Figure 6b. Model prediction for the water depth in the channel was further compared with the experimental observations in [36] and the published numerical results of [37] in Figure 7, at t = 3, 5, 7, and 14 s. It is evident that our scheme achieves a very satisfactory agreement with the experimental observations, demonstrating that it could satisfactorily model two-dimensional supercritical flows involving bends and sudden directional changes. demonstrating that it could satisfactorily model two-dimensional supercritical flows involving bends and sudden directional changes.   demonstrating that it could satisfactorily model two-dimensional supercritical flows involving bends and sudden directional changes.

Flow over a Hump
This test problem is widely used to verify well-balanced property of a numerical scheme towards moving steady state solution [7,8,25]. The moving steady-state water equilibrium is given by:

Flow over a Hump
This test problem is widely used to verify well-balanced property of a numerical scheme towards moving steady state solution [7,8,25]. The moving steady-state water equilibrium is given by: where m and E are equilibrium variables. The bed profile of a 25 m-long and 0.1-m wide channel with a bump is defined as Boundary conditions are adjusted to generate subcritical flow with a shock, or a transcritical flow with or without a shock. The domain is treated as two-dimensional by discretizing with unstructured triangular cells and the CFL number of 0.6 was maintained for all cases.

Transcritical Flow with a Shock
We simulate a transcritical flow with a shock by setting a water depth (0.33 m) and discharge (q x = 0.18 ms −2 ) at the left and right boundaries, respectively. The initial conditions are fixed as: At upstream boundary a discharge of 1.534 ms −2 is imposed and a water depth of 0.66 m is imposed at downstream boundary.

Subcritical Flow
The initial conditions are fixed as: At upstream boundary a discharge of 4.42 ms −2 is imposed and a water depth of 2 m is imposed at the downstream boundary.
For steady states, the initial condition given in Equations (43), (45) and (46) should be exactly preserved. We run the simulations for 20.0 s; comparisons of the numerical results for the water surface elevation and discharge with the analytical solution were in excellent agreement, as shown in Figures 8-10. The comparison of numerical and analytical solutions shows excellent agreement; however, discharge for flow with shock in Figure 8 shows discrepancy with the analytical solution at x = 12 m, which is where shock is generated. This disagreement was also observed previously in [20,28,38]. We measured convergence of solution for transcritical flow with a shock by using global relative error: where h is water depth; i is cell number, N is the total number of cells; n is time level; and A is the area of cell i. To demonstrate that the steady state is exactly preserved we measure L 1 error for water height and discharge for all three cases at t = 20 s. The values in Table 1 show smaller deviation between numerical and exact solutions of water depth than for discharge which is due to discrepancy between exact and numerical discharge immediately downstream of the hump, however, the overall agreement between computed and exact solution is satisfactory. Comparison of convergence rates for transcritical flow with a shock in Figure 11 shows improvement by the proposed method where convergence is reached after 15,000 iterations. For SGM [20] and RSGM [7] methods, convergence is achieved after 16,000 and 15,600 iterations, respectively. between numerical and exact solutions of water depth than for discharge which is due to discrepancy between exact and numerical discharge immediately downstream of the hump, however, the overall agreement between computed and exact solution is satisfactory. Comparison of convergence rates for transcritical flow with a shock in Figure 11 shows improvement by the proposed method where convergence is reached after 15,000 iterations. For SGM [20] and RSGM [7] methods, convergence is achieved after 16,000 and 15,600 iterations, respectively.

Thacker's Planar Solution
Flow inside a parabolic basin is an important test case in hydraulics and coastal engineering, modelled by a number of researchers [4,24,26] for testing the ability of the model to handle wet-dry fronts over uneven beds. Simulating flow inside a parabola is challenging because of periodic wetting and drying involving run-up and run-down. We verify the wet-dry front handling ability of our model by considering a parabolic computational domain of 4 m × 4 m with a frictionless bottom.

Thacker's Planar Solution
Flow inside a parabolic basin is an important test case in hydraulics and coastal engineering, modelled by a number of researchers [4,24,26] for testing the ability of the model to handle wet-dry fronts over uneven beds. Simulating flow inside a parabola is challenging because of periodic wetting and drying involving run-up and run-down. We verify the wet-dry front handling ability of

Thacker's Planar Solution
Flow inside a parabolic basin is an important test case in hydraulics and coastal engineering, modelled by a number of researchers [4,24,26] for testing the ability of the model to handle wet-dry fronts over uneven beds. Simulating flow inside a parabola is challenging because of periodic wetting and drying involving run-up and run-down. We verify the wet-dry front handling ability of our model by considering a parabolic computational domain of 4 m × 4 m with a frictionless bottom. Geometry of the domain is given by the equation: where h o is the water depth at center of domain and a is the distance from center to zero elevation of shore. The analytical solution to the problem is given as [39]: where σ is the amplitude parameter fixed equal to 0.5; ω represents rotational frequency given by Water surface elevations and velocity results obtained from numerical scheme are compared with the analytical solution at time 4T, shown in Figure 12. The comparisons show very good agreement between numerical results and analytical solution for water levels at y = 2 m. The comparison of vector field provided for velocity show that very small perturbations exist near the wet-dry front, but importantly we observed that those perturbations do not affect the accuracy of the solution by growing in time; moreover, it has been observed in previous studies too, that small perturbations near the wet-dry front are very difficult to remove. The shoreline propagates smoothly between wet-dry surfaces with no spurious oscillations which indicate that water depth and velocity at wet-dry interfaces are estimated accurately with a hydrostatic reconstruction method and proposed treatment of slope source terms.
To assess the asymptotic order convergence for water depth h and velocity u, we measure the errors in terms of L 1 and L 2 norms for successively refined unstructured grids. The error norms are computed by the following equation: where p is the error norm being measured; T is the flow variable, in this case water depth h and velocity u; A is the area of the cell, and N is the total number of elements in the mesh. Table 2 shows the properties of successively refined grids and values of error norms for water depth and velocity. We observe behavior in values of error norms for both flow variables since only distorted unstructured grids are used to check the convergence. Mass conservation of the solution is verified by the following equation: where V o and V t are the volume of water initially and after time t, respectively, and δV t is the difference between initial and final volume of water. For flow inside the parabola δV t = 0, because there is no flow through the boundaries. If the water depth in cells, classified as dry by the threshold limit, is set to zero it will lead to subtraction of water volume from the total system which results in non-conservation of mass. Therefore, to preserve mass exactly, after identifying the cells with water depth 0 ≤ h ≤ 10 −5 m, we find the volume of water in those cells before setting the water depth to zero. This volume of water is then added to the volume in the neighboring wet cells at every time step. This step preserves the mass conservation exactly. The maximum relative mass error computed with Equation (53) is 4.9 × 10 −14 , which is accurate up to round-off error and confirms the mass conservation. Water surface elevations and velocity results obtained from numerical scheme are compared with the analytical solution at time 4T, shown in Figure 12. The comparisons show very good agreement between numerical results and analytical solution for water levels at y = 2 m. The comparison of vector field provided for velocity show that very small perturbations exist near the wet-dry front, but importantly we observed that those perturbations do not affect the accuracy of the solution by growing in time; moreover, it has been observed in previous studies too, that small perturbations near the wet-dry front are very difficult to remove. The shoreline propagates smoothly between wet-dry surfaces with no spurious oscillations which indicate that water depth and velocity at wet-dry interfaces are estimated accurately with a hydrostatic reconstruction method and proposed treatment of slope source terms.
To assess the asymptotic order convergence for water depth h and velocity u, we measure the errors in terms of 1 L and 2 L norms for successively refined unstructured grids. The error norms are computed by the following equation: where p is the error norm being measured; T is the flow variable, in this case water depth h and velocity u; A is the area of the cell, and N is the total number of elements in the mesh. Table 2 shows the properties of successively refined grids and values of error norms for water depth and velocity. We observe behavior in values of error norms for both flow variables since only distorted unstructured grids are used to check the convergence.

Dam-Break over a Triangular Obstacle
This test case models the complicated interactions of shock and rarefaction waves, to test the shock handling and propagation of the proposed scheme. This case was previously modeled in [1,7,24] to validate their numerical schemes. A thin wall acting as a hypothetical dam was positioned 15.5 m from the left boundary. To generate a shock wave, a triangular obstacle 0.4 m high and 6 m wide was installed 10 m downstream from the dam.
The computational domain discretized into 2590 elements is shown in Figure 13. The water depth was set to 0.75 m upstream, and 0 m downstream from the dam. Free outfall and solid boundary conditions were imposed on the left and right boundaries, respectively. The laboratory experiment was performed by EU CADAM in a 38 m-long channel with a Manning's roughness coefficient equal to 0.0125. The simulation results after 90 s were compared with the experimental observations in Figure 14 at seven gauges, located at 2 m, 4 m, 8 m, 10 m, 11 m, 13 m, and 20 m downstream from the dam.
Due to the relative closeness of gauges 1 and 2, almost similar water depth profiles were observed. We observed a gradual increase in water depth on an initially dry bed until approximately 12 s, after which there was an abrupt increase in depth caused by a propagating shock wave reflected from the triangular obstacle. However, the right moving rarefaction wave caused an almost linear draw down in water depth for up to 32 s for gauge 1 and 35 seconds for gauge 2. At that instant, the second shock wave (reflected from the triangular obstacle) again generated an increase in water depth at gauges 1 and 2. However, this time the maximum water depth was lower than before due to the dissipation of energy by friction. The interaction of a rarefaction wave again caused a drop in water level, until the third shock wave reached gauges 1 and 2 at around 45 s.
This wave interaction process continued until all the energy was dissipated and the mass losses due to transmissive flow at the downstream boundary made the flow steady. The same wave propagation processes and interactions were also observed at gauges 3-7. An observed trend in the numerical computations was the early arrival of the second-to-last shock for gauges 1 and 2, and the last shock for the other gauges. This could be due to inaccurate estimation of wave speeds in the HLLC flux computation, causing the waves to travel too fast. However, this was not certain, since the arrival times of the other shocks and rarefaction waves were accurately predicted. Another source for the discrepancy could be error in the experimental observations; there was no consistency in the observations plotted between 65 and 90 s by Liang and Marche [1], observations were both on the higher and lower side. In such a case, the choice of plotting higher or lower observations may affect the resultant comparisons.
Regardless of this small discrepancy, the overall agreement between the numerical results and experimental observations was satisfactory, showing that our scheme is capable of modeling complex wave interactions and wetting/drying. Mass conservation of the solution is verified by the following equation: where o V and t V are the volume of water initially and after time t, respectively, and t V  is the difference between initial and final volume of water. For flow inside the parabola 0 there is no flow through the boundaries. If the water depth in cells, classified as dry by the threshold limit, is set to zero it will lead to subtraction of water volume from the total system which results in non-conservation of mass. Therefore, to preserve mass exactly, after identifying the cells with water depth 5 0 10 h    m, we find the volume of water in those cells before setting the water depth to zero. This volume of water is then added to the volume in the neighboring wet cells at every time step. This step preserves the mass conservation exactly. The maximum relative mass error computed with Equation (53) is 14 4.9 10   , which is accurate up to round-off error and confirms the mass conservation.

Dam-Break over a Triangular Obstacle
This test case models the complicated interactions of shock and rarefaction waves, to test the shock handling and propagation of the proposed scheme. This case was previously modeled in [1,7,24] to validate their numerical schemes. A thin wall acting as a hypothetical dam was positioned 15.5 m from the left boundary. To generate a shock wave, a triangular obstacle 0.4 m high and 6 m wide was installed 10 m downstream from the dam.
The computational domain discretized into 2590 elements is shown in Figure 13. The water depth was set to 0.75 m upstream, and 0 m downstream from the dam. Free outfall and solid boundary conditions were imposed on the left and right boundaries, respectively. The laboratory experiment was performed by EU CADAM in a 38 m-long channel with a Manning's roughness coefficient equal to 0.0125. The simulation results after 90 s were compared with the experimental observations in Figure 14     Due to the relative closeness of gauges 1 and 2, almost similar water depth profiles were observed. We observed a gradual increase in water depth on an initially dry bed until approximately 12 s, after which there was an abrupt increase in depth caused by a propagating shock wave reflected from the triangular obstacle. However, the right moving rarefaction wave caused an almost linear draw down in water depth for up to 32 s for gauge 1 and 35 seconds for gauge 2. At that instant, the second shock wave (reflected from the triangular obstacle) again generated an increase in water depth at gauges 1 and 2. However, this time the maximum water depth was lower than before due to the dissipation of energy by friction. The interaction of a rarefaction wave again caused a drop in water level, until the third shock wave reached gauges 1 and 2 at around 45 s. This wave interaction process continued until all the energy was dissipated and the mass losses due to transmissive flow at the downstream boundary made the flow steady. The same wave propagation processes and interactions were also observed at gauges 3-7. An observed trend in the numerical computations was the early arrival of the second-to-last shock for gauges 1 and 2, and the last shock for the other gauges. This could be due to inaccurate estimation of wave speeds in the HLLC flux computation, causing the waves to travel too fast. However, this was not certain, since the arrival times of the other shocks and rarefaction waves were accurately predicted. Another source for the discrepancy could be error in the experimental observations; there was no consistency in the observations plotted between 65 and 90 s by Liang and Marche [1], observations were both on the higher and lower side. In such a case, the choice of plotting higher or lower observations may affect the resultant comparisons.
Regardless of this small discrepancy, the overall agreement between the numerical results and experimental observations was satisfactory, showing that our scheme is capable of modeling complex wave interactions and wetting/drying.

Conclusions
In this paper, we introduced a new, accurate, well-balanced numerical scheme for shallow-water flows using a high-resolution Godunov-type finite volume method. The proposed slope source term discretization is based on averaging water depths at the cell center and reconstructed water depths at

Conclusions
In this paper, we introduced a new, accurate, well-balanced numerical scheme for shallow-water flows using a high-resolution Godunov-type finite volume method. The proposed slope source term discretization is based on averaging water depths at the cell center and reconstructed water depths at both sides of the interface and multiplying this with the difference in bed elevation between the cell center and interface. Applying the model to numerical examples shows that the proposed interpretation of slope source terms can accurately maintain well-balanced conditions for still water. Results from wave run-up/run-down on the conical island and dam-break over triangular obstacle indicate that the proposed scheme can accurately handle wet-dry fronts. Flow over the hump demonstrates that the slope source term exactly balances convective fluxes maintaining a well-balanced condition. The proposed model is tested for various benchmark coastal and hydraulic engineering cases; comparisons between the experimental, analytical, and published numerical results show that the scheme is robust and has a wide range of applications with very satisfactory modeling capabilities for a wide variety of mesh types. Based on the results of several test cases, we can confidently predict that application of the proposed method with regard to slope source term discretization for complicated flow problems (real life and laboratory) will result in well-balanced and accurate solutions.