Modeling of Shallow-Water Equations by Using Implicit Higher-Order Compact Scheme with Application to Dam-Break Problem

To model the shallow water equations (SWE), various numerical methods such as explicit and implicit finite-difference methods with higher-order compact (HOC) scheme have been used by a number of investigators. The compact scheme is a class of high order finite difference methods, which provides an effective way of joining spectral method for accuracy and robust characteristics of finite difference schemes. The compact scheme uses smaller stencil to approximate leading truncation error terms of the governing equation and to provide better resolution. Such approximation provides a better representation at shorter length scales and is well applicable for the simulation of waves with high frequency. The application of compact scheme for various fluid flow problems can be obtained from the investigations [1-10]. Navon and Riphagen [1] developed compact fourth order algorithm by using alternating-direction implicit finite-difference scheme to solve non-linear shallow-water equations, expressed in conservationlaw form. Abarbanel and Kumar [2] obtained an unconditionally stable HOC scheme for the Euler Equations. Gupta [3] extended the high accuracy approximation method to solve Navier-Stokes equation for lead driven cavity problem. Lele [4] has investigated higher order compact scheme with spectral-like resolution with the approximation of the first and second derivative up to tenth order on a uniform grid. Spotz and Carey [5] developed higher-order difference scheme on compact stencil for boundary value problems and in particular to solve steady stream-function vorticity equations. Li and Tang [6] extended the previous work of Li et al. [11] for steady Navier-Stokes equation to unsteady case by using fourth order accurate compact scheme. Kalita et al. [7,9] proposed a transformation-free HOC for steady convectivediffusion equation and an efficient transient Navier-Stokes solver by using non-uniform grids. Pandit et al. [8] proposed an implicit HOC finite-difference scheme for two-dimensional unsteady Navier-Stokes equations in irregular geometries by using orthogonal grids. Shah et al. [10] developed a time accurate upwind higher order accurate compact finite difference scheme to Navier-Stokes equations by using artificial viscosity approach. The analysis of dam-break flow is important to capture spatial and temporal evolution of flood event and safety analysis. The dam-break is basically catastrophic failure of dam, leading to uncontrolled release of water causing flood in the downstream region. Several researchers have studied the flood wave caused due to the dam-break either experimentally or numerically. The experimental investigations can be obtained from the literature [12-15]. The successful representation of dam-break flow using numerical simulation have been reported by [16-26]. Townson and Al-Salihi [19] developed dam-break flow model for parallel, converging and diverging boundaries using the method of characteristics. Bellos et al. [20] developed a two-dimensional numerical method for dam-break problem by using the combination of finite-element and finite-difference methods. Mohapatra and Bhallamundi [21] simulated dam-break flow in channel transitions by using MacCormack [27] scheme. Zoppou and Roberts [22] obtained numerical solution for unsteady dam-break flow problem. Arico et al. [23] used shallow water approach for the representation of dambreak flow. The literature on open channel and dam-break flows can be obtained from Chaudhry [24]. Recently, three-dimensional numerical simulations for free surface flows induced by dam-break have been reported by Biscarini et al. [25]. Bellos and Hrissanthou [26] developed two numerical models for one-dimensional Shallow Water Equation (SWE) or Saint-Venants’s equation by using Lax-Wendroff [28] and MacCormack [27] schemes.


Introduction
To model the shallow water equations (SWE), various numerical methods such as explicit and implicit finite-difference methods with higher-order compact (HOC) scheme have been used by a number of investigators. The compact scheme is a class of high order finite difference methods, which provides an effective way of joining spectral method for accuracy and robust characteristics of finite difference schemes. The compact scheme uses smaller stencil to approximate leading truncation error terms of the governing equation and to provide better resolution. Such approximation provides a better representation at shorter length scales and is well applicable for the simulation of waves with high frequency. The application of compact scheme for various fluid flow problems can be obtained from the investigations [1][2][3][4][5][6][7][8][9][10]. Navon and Riphagen [1] developed compact fourth order algorithm by using alternating-direction implicit finite-difference scheme to solve non-linear shallow-water equations, expressed in conservationlaw form. Abarbanel and Kumar [2] obtained an unconditionally stable HOC scheme for the Euler Equations. Gupta [3] extended the high accuracy approximation method to solve Navier-Stokes equation for lead driven cavity problem. Lele [4] has investigated higher order compact scheme with spectral-like resolution with the approximation of the first and second derivative up to tenth order on a uniform grid. Spotz and Carey [5] developed higher-order difference scheme on compact stencil for boundary value problems and in particular to solve steady stream-function vorticity equations. Li and Tang [6] extended the previous work of Li et al. [11] for steady Navier-Stokes equation to unsteady case by using fourth order accurate compact scheme. Kalita et al. [7,9] proposed a transformation-free HOC for steady convectivediffusion equation and an efficient transient Navier-Stokes solver by using non-uniform grids. Pandit et al. [8] proposed an implicit HOC finite-difference scheme for two-dimensional unsteady Navier-Stokes equations in irregular geometries by using orthogonal grids. Shah et al. [10] developed a time accurate upwind higher order accurate compact finite difference scheme to Navier-Stokes equations by using artificial viscosity approach.
The analysis of dam-break flow is important to capture spatial and temporal evolution of flood event and safety analysis. The dam-break is basically catastrophic failure of dam, leading to uncontrolled release of water causing flood in the downstream region. Several researchers have studied the flood wave caused due to the dam-break either experimentally or numerically. The experimental investigations can be obtained from the literature [12][13][14][15]. The successful representation of dam-break flow using numerical simulation have been reported by [16][17][18][19][20][21][22][23][24][25][26]. Townson and Al-Salihi [19] developed dam-break flow model for parallel, converging and diverging boundaries using the method of characteristics. Bellos et al. [20] developed a two-dimensional numerical method for dam-break problem by using the combination of finite-element and finite-difference methods. Mohapatra and Bhallamundi [21] simulated dam-break flow in channel transitions by using MacCormack [27] scheme. Zoppou and Roberts [22] obtained numerical solution for unsteady dam-break flow problem. Arico et al. [23] used shallow water approach for the representation of dambreak flow. The literature on open channel and dam-break flows can be obtained from Chaudhry [24]. Recently, three-dimensional numerical simulations for free surface flows induced by dam-break have been reported by Biscarini et al. [25]. Bellos and Hrissanthou [26] developed two numerical models for one-dimensional Shallow Water Equation (SWE) or Saint-Venants's equation by using Lax-Wendroff [28] and MacCormack [27] schemes.
In the present study, we propose HOC scheme to solve transient two-dimensional SWE. The distinct advantage of this scheme is high order accuracy associated with compact stencils which provides accurate numerical solution even for relatively coarser grids. To demonstrate the application of higher order compact scheme, classical dam-break problem is considered to simulate the flood wave in channel transition. We consider parallel-parallel channel with complete breach [21] as shown in Figure 1a, 1b. The flow is considered to be inviscid and incompressible and accordingly the non-linear SWE are expressed as These above equations are described in primitive variable form are obtained from Navier-Stokes (NS) equations by integrating over the depth and by assuming hydrostatic pressure distribution. Further, vertical velocity component and corresponding shear stress is also considered insignificant. The set of partial differential equations (PDEs) are nonlinear first-order and hyperbolic in nature. Here, h is water depth; u is depth averaged velocity in the x-direction; v is depth averaged velocity in the y-direction; g is acceleration due to gravity; S is bottom friction in the y-direction. Generally, the bottom friction can be estimated by using the Manning's formula: Where n is Manning roughness coefficient; and c 0 is a dimensional constant. The paper is organized as follows: Section-2 provides the mathematical formulation and discretization procedure. Section-3 discusses computation procedure and Section-4 discusses boundary conditions. In Section-5, linearized stability analysis is derived. The numerical results of test cases are given in Section-6. Finally, Section-7 provides results and discussion.

Mathematical Formulation and Discretization Procedure
The governing equations (1)-(3) can be written in conservative law form or divergence form as where φ, F(φ), G(φ) and S(φ) are the column matrices Where φ denotes a vector containing the primitive variables h, u, v; F(φ) and G(φ) are vectors written in flux form and are functions of φ; S(φ) is the source Term. The scheme developed here is in rectangular Cartesian coordinate system using nine point stencil for grid point (i, j) with two-time levels at n and n+1 as shown in Figure 2. We discretize equation (1) by adopting two procedures: (i) non-conservative form and (ii) conservative law form. Both the approaches are applied to solve shallow water equations (SWE). However, conservative law formis used to solve equation (5).

Discretization in non-conservative form
It is necessary for some numerical schemes to represent the governing equations in non-conservative form. Therefore, after rearranging equation (5) in non-conservation form, one can write x y (7) where A, B are the Jacobian of F(φ), G(φ) with respect to φ: First we develop HOC formulation for the steady-state form of equation (7) which is obtained when φ, A, B and S (φ) are independent of t. Accordingly, equation (7) becomes The rectangular domain is divided in to a uniform mesh size ∆x and ∆y along x and y directions respectively; φ ij refers to the approximation   Figure 2: Higher order compact stencil using (9,9) grid points.
of φ(x i , y j ) where (x i , y j ) is the coordinate of a typical node. In the computational domain on a grid point (i, j), the stencil is similar to the one at the n or (n+1) th time level as shown in Figure 2. A second-order central difference approximation is considered by including the leading truncation error terms at the (x i , y j ) th node for the derivatives appears in equation (9) ( ) ( ) where δ x and δ y are the first order central difference operators along the x and y directions respectively. To obtain a fourth order compact scheme, we approximate the derivatives by differentiating equation (9) and include them in the difference equation (10). The successive differentiation of (9) with respect to x yields ( ) Using central differencing and writing in operator form, equations (11) and (12) become ( ) Substituting equations (13), (14) into equation (10) we obtain the following approximation to (9) ( After rearranging, we obtain where the coefficients A ij , B ij , C ij , D ij , E ij , F ij , G ij and ̅ S ij are as follows: After rearranging, equation (23) can be finally written as ( ) (25) and I is the identity matrix, n represents the time level and the coefficients A ij , B ij , C ij , D ij , E ij , F ij , G ij and ̅ S ij to be evaluated at the time level n. We introduce a weighted average parameter μ for the approximation of time derivative with differencing ( ) This provides a class of integrators, i.e., for μ=0, p=1 gives forward Euler, μ=1, p=1 gives backward Euler and μ μ=1/2, p=2 gives Crank-Nicolson. Equation (26) can be written as where To derive a finite difference scheme for the unsteady case, we replace the source function S(φ) in equation (9) p q x y (28) and ( ) Equation (27) represents HOC finite difference approximation for the two-dimensional unsteady SWE. For uniform grids, the spatial accuracy of the proposed scheme becomes fourth order accurate.

Discretization in conservation form
From equation (5) the steady-state form becomes We consider a second-order central difference approximation by including the leading truncation error terms at the (x i , y j ) the node for the derivatives appears in (32) The successive differentiation of (32) with respect to x and y yields  (27), (28), (29) and (30) remain unchanged whereas equations (31) change along with the matrices C ij , D ij , E ij , F ij , G ij and these are

Computational Procedures and Algorithm
To solve the algebraic systems of equation (27), the following computational procedures are adopted. Using the proposed finite difference approximation, the system of equation (27) can be written in matrix form as Where M is a coefficient matrix. For the coefficient matrix M with grid size of (m×n), the structure is 9-band block matrix whose elements are also 3×3 matrices with (3×m×n) dimension. σ n+1 and f(σ n ) are (3×m×n) component vectors can be expressed as Where φ ij are three-element vectors of flow variables h, hu, hv. It can be mentioned here that after each iteration, the coefficient matrix M changes. We use bi-conjugate gradient stabilized method (BiCGStab) with preconditioning. It is also possible to partition the coefficient matrix M to a block tridiagonal matrix with its elements (3×n)×(3×n) matrix where m×n is the grid size. Equation (42) shows the block tridiagonal matrix for grid size 5×5 as an example with its elements 15×15 with Neumann boundary conditions at the side walls and Dirichlet boundary conditions at the inflow and outflow boundaries.

Solution algorithm
Step 1: Initialize the height and velocity fields Step 2: Evaluate the variable coefficient matrix M and variable coefficient f (σ n ) in (27) using equation (40) and boundary conditions (52) and (53). The source term in n and n+1 time level is evaluated with one step delay.
Step 4: Filter out short-wave noise from the u, v, and h using Wallington filter (70) Step 5: Store σ and time step in a file.
Step 6: Increase time step by ∆t, if t ≤ t max go to Step 2 .

Compact Boundary Conditions
To solve discretized equation (27), Dirichlet boundary conditions for inflow and outflow and Neumann boundary conditions for both the walls are considered. Prior to dam breach the upstream section acts as a reservoir. Therefore, Dirichlet boundary condition is applied at the downstream end after the dam breach. Neumann boundary conditions are applied in reflective or solid boundaries. The advantage of Dirichlet conditions is that they are simple and straightforward and relatively easy to implement in a HOC setting. We implement Adam [29] boundary conditions and Taylor series analysis of backward and forward difference to develop HOC wall boundary conditions for ∂φ/∂y. However, for Neumann type conditions, it is possible to use one order lower in accuracy at the boundary approximations although such approximations to the interior points remain unaffected [30]. We adopt two boundary conditions Adam [31] and Taylor series analysis for HOC approximations for wall boundaries.

Adam boundary condition
Using the Adam boundary conditions [29], HOC approximation at the walls and interior points in y-direction become Here j=1 and j=n indicate the left and right side walls respectively.

Taylor series analysis
Using Taylor's series analysis and using the forward difference for j=1 , we get Using similar approach as described in equation (27) with μ=0, we get ϕ ϕ δ δ ϕ δ δ ϕ The following linear extrapolation equations are used at both the boundaries And the amplification matrix  G is given by After simplification, we obtain For stability, all the eigenvalues of amplitude matrix  G must be on the unit circle. If ti, i=1,2,3 are the eigenvalues of T and  i g , i=1,2,3 are the eigenvalues of  G then  i g can be written as Henrici [31]  ( ) We consider two cases corresponding to non-conservative and conservative forms.

Analysis for non-conservative form
In this case, the eigen values of matrix T can be defined as real and imaginary parts as Accordingly, equation (60) can be written as This indicates that the method is stable provided the above condition is satisfied.

Conservative form
For conservative form the matrix T can be obtained by substituting equations (37)-(39) in equation (59) ( ) Since in equation (67) real part is zero (Appendix-1), t i =it" i . Accordingly, equation (66) becomes It means that the method is unconditionally stable for μ ≥ 1/2.

Numerical Experiment and Model Validation
To demonstrate the validity and effectiveness of the proposed scheme, unsteady 2D dam-break flow problem is considered. The schematic diagram of Dam-break flow in straight walled transition is shown in Figure 1. Numerical computations are made for the following input values which correspond to the experiments conducted by Townson and Al-Salihi [19]. The computational domain comprises of 4 m long and 0.1 m wide channel. The length of the channel upstream of the dam is 1.8 m and the channel downstream is 2.2 m. The initial flow depth on the upstream side of the dam is 0.1 m for both wet bed and dry bed whereas downstream water depth was assumed to be 0.00001 m and 0.006 m respectively. The velocities, u and v at inflow boundary are specified as zero for all times t and the flow depth at this boundary computed by using the following linear extrapolation.
The boundary conditions at the downstream end are the same as the initial condition in terms of depth and velocities. The symmetry boundaries are treated at the sidewalls. A higher order compact finite difference implicit scheme is used for numerical simulation on a 21×5 grid with ∆t=0.01. The corresponding grid sizes become ∆x=0.2 m and ∆y=0.025 m along x and y directions respectively. To reduce the higher order dissipation due to aliasing error introduced by the fourthorder compact scheme, Wallington filter [32] is applied. This consists of periodic successive application of the following two-point operators where one filtering is performed after eighteen iterations.
The numerical wave profiles are compared with the experimental values, 2.5 s after the instantaneous dam-break for the wet bed conditions whereas for dry bed condition 1.5 s is considered. Figure   3a shows very good agreement of water surface profiles between numerical and experimental results for wet bed condition. For dry bed condition as shown in figure 3b, the result shows very good agreement for upstream section and also closer to the downstream boundary. However, in both the cases compact formulation provides much better results as compared to the results of Mahapatra and Bhallamudi [21]. With such coarser discretization, desired accuracy has been achieved in comparison to the traditional schemes.

Results and Discussions
In this problem, dam breach takes place instantaneously with full breach condition. Due to discontinuous initial conditions, it enforces several restrictions and most of the numerical schemes fail to capture this phenomenon. The algorithm developed here is second order accurate in time and fourth order accurate in space O((∆t) p, h 4 , k 4 ) , on the nine-point stencil using third order non-centered difference at the wall boundaries. For solving the algebraic systems, we adopt biconjugate gradient stabilized method (BiCGStab) with preconditioning for 21×5 grid points along x and y directions. Interpolation has been performed along y direction while plotting. It is also possible to partition the coefficient matrix to a block tri-diagonal matrix then find Solution by simple block tri-diagonal solvers. We present our results for wet bed and dry bed conditions for downstream of the dam. The actual water depth at t=0 is specified as the initial condition in the downstream. The tail water/reservoir ratio for wet bed and dry bed is considered as 0.176 and 0.025 respectively. In the case of dry bed, downstream water level is fixed at 0.0025 m and below to this level would cause blow up due to not satisfying the continuity equation. Figure 3c shows the comparison of water depth profiles for wet bed and dry bed conditions at the downstream after 1.5 s of dam-break. During this time the wave due to dam-break has not reached the upstream boundary. It can be noticed that the water depth profiles for wet bed and dry bed condition remains same till mid location of the upstream section and subsequently water depth profile of wet bed increases gradually towards downstream direction. Figure 3d shows the comparison of water depth profiles for wet bed and dry bed conditions for time t=2.5 s after the dam-break. The wave front clearly indicates the effect of dam break closer to the Water surface profiles at t=0.5s, 1.0s after the dam-break for dry bed condition. c: Water surface profiles at t=1.5s after the dam-break for dry-bed condition.
downstream boundary for wet bed condition. Figures 4a-4c show the three-dimensional view of water surface profiles for dry bed taken at time t=0.5 s, 1.0 s and 1.5 s respectively. It can be noted that water surface profiles manifest smooth transition even with less number of grid points. At time t=0.5 s, (Figure 4a), the flood wave is yet to reach the downstream region and therefore water level remains at 0.0025 m. Figure 4b shows that after t=1.0 s, the flood wave almost reach downstream region barring couple of grid points close to the boundary. At time t=1.5 s, the flood wave reaches the downstream completely and drop in water level in the upstream boundary is noticed (Figure 4c). We consider four cases to illustrate the effect of flood wave for wet bed conditions for time t=0.5 s, 1.0 s, 1.5 s and 2.5 s and the corresponding water surface profiles are shown in Figres 5a-5d respectively. Similar to the dry bed conditions, the flood wave does not reach the downstream completely at time t=0.5 s and 1.0 s as shown in Figures 5a and 5b respectively. Nevertheless, corresponding rise in the downstream water levels is due to the imposition of higher water level during initial condition (wet bed). Figure 5c and 5d show that the wave due to dambreak has reached upstream end with significant decrease in flow depth and the movement of wave front through numerical simulation.

Conclusions
In this paper we have developed two types of fourth order compact schemes for steady as well as unsteady cases using non-conservation and conservative forms for two-dimensional shallow water equations. This fourth order compact scheme is applied for dam-break problem and the numerical results are compared well with the experimental results. Due to the implicit nature and less stringent stability criteria, the present model proved to be more versatile even with coarser grids and can provide a useful framework for spatial discretization. It is also possible to partition the coefficient matrix to block tridiagonal matrix and then obtain solution by using the simple tridiagonal solver. The present model could be a good alternative to the existing explicit Mac-Cormack [27] or implicit Beam and Warming [33] methods for solving hyperbolic systems in conservative law form, particularly for solving the shallow water equation (SWE).