Stability Analysis of an Upwind Difference Splitting Scheme for Two-Dimensional Saint–Venant Equations

: The paper is devoted to the construction and study of a numerical method for solving two-dimensional Saint–Venant equations. These equations have important applied signiﬁcance in modern hydraulic engineering and are suitable for describing waves in the atmosphere, rivers and oceans, and for modeling tides. The issues of formulation of the mixed problems for these equations are studied. The system of equations is reduced to a symmetrical form by transforming dependent variables. Then, the matrices of coefﬁcients are represented as the sums of two symmetric semideﬁnite matrices. This transformation allows constructing an upwind difference scheme in spatial directions to determine the numerical solution of the initial boundary value problem. The stability of the proposed difference scheme in energy norms is rigorously proved. The results of numerical experiments conducted for a model problem are provided to conﬁrm the stability of the proposed method.


Introduction
Numerical modeling of surface run-off, which is the main component of the water cycle, is known to be of critical importance in hydraulic and hydrological designs all over the world. For land surface and open channel currents, 1D, 2D, or coupled 1D and 2D depth-averaged Saint-Venant equations are widely used because of their simplicity and efficiency compared to the 3D Navier-Stokes equations [1][2][3].
Numerous scientific and practical studies in the field of hydrology are focused on the construction and study of difference schemes for the numerical solutions of mixed problems for the Saint-Venant system of equations.
Situations where the depth of the water is significantly smaller than its horizontal dimensions are quite common, so the two-dimensional system of Saint-Venant equations is widely used. It is applied by taking into account the Coriolis forces in the modeling of the atmosphere and the ocean as a simplification of the system of primitive equations describing the flows in the atmosphere. In this regard, one of the urgent problems is the solution of a multidimensional problem for the Saint-Venant system of equations. In this case, it is important to construct difference schemes for a mixed problem for the system of Saint-Venant equations and prove their stability. The one-dimensional theory of unsteady flows without discontinuities, i.e., flows not accompanied by the formation of discontinuous waves, can be considered the most developed to date (this includes the problem of a natural flood). This type of fluid motion is described by the classical Saint-Venant equations. The solutions of these systems of equations have been and remains one of the most important problems.
It must be noted that the application of the method of characteristics for complexshaped channels, especially rivers, was associated with significant difficulties and led to an unreasonably large number of computations. The disadvantage of explicit difference schemes, along with the method of characteristics, was the relatively large usage of computer time, since the computer technology of that time worked much slower than modern one. In [4], an implicit difference grid scheme was developed which makes it possible to perform calculations with a large time step. S. M. Shugrin [5] proposed a difference scheme that makes it possible to compute unsteady motion in complex channel systems.
One of the most difficult problems of unsteady motion in open channels is the problem of wave motion that occurs when a dam breaks. The use of computers makes it possible to solve the problem of breaking a dam in a fairly precise formulation. J. B. Whitham [6] and R. F. Dressler [7] proposed special methods for taking into account friction while considering the problem of the destruction of a dam with an outflow into a dry channel based on the shallow water equations. S. M. Shugrin [5] proposed a new mathematical statement of the problem of wave motion along a dry channel which is simpler than the previous ones. As in the works of C. Montuori [8], Shugrin believed that the Saint-Venant equations are valid for the entire flow region, and the boundary conditions at the wavefront are obtained by equating the front propagation velocity and the flow velocity near the wave head. The calculations performed in [9] in this formulation by the difference method using a moving grid showed satisfactory agreement with the experimental results.
The study of properties of unsteady flows in open flows described by the Saint-Venant system of equations has been studied by many scientists [10][11][12][13][14].
For the numerical solution of various problems for the system of Saint-Venant equations, scientists have proposed several methods [31][32][33][34][35][36]. The book [37] describes in detail approaches to discontinuity detection and methods of pass-through computation. The authors of this book focus on Godunov-type methods, which are based on exact and approximate solutions to the Riemann problem of the decay of an arbitrary discontinuity. Various studies are aimed at studying essential properties of difference schemes for these types of equations, namely, their conservatism [38], being well-balanced [39,40], the possibility of pass-through calculation and the choice of an explicit/implicit scheme [41].
An important property of difference schemes for hyperbolic equations is their stability, that is, the ability of the scheme to not accumulate computational perturbations. In the case of conditionally stable schemes, a restriction is imposed on the choice of the time integration step. A slight violation of the stability condition can lead to an increase in the error and non-physical solutions. For hyperbolic systems, the time step is determined from the Courant-Friedrichs-Lewy (CFL) condition [41].
Taking into account the practical significance of the problem, we devoted our study to the construction of an efficient numerical method for solving the system of two-dimensional Saint-Venant equations. The idea behind this method is based on the partition of matrices depending on the flow velocity into the sum of two symmetrical semidefinite matrices. This transformation makes it possible to construct a spatial splitting upwind difference scheme.
In this study, the main attention is devoted to a rigorous proof of the stability of the proposed difference scheme in energy norms. Specifically, sufficient conditions for obtaining a stable approximate solution are derived. The results of the theoretical analysis are verified on the example of a test problem.
The paper is structured as follows. In the next section, we formulate the problem statement and reduce it to a symmetrical form by transforming dependent variables. In Section 2.2, the matrices dependent on the flow velocity are presented as a sum of symmetrical matrices. Then, Section 2.3 provides an implicit upwind difference scheme for solving the obtained initial boundary value problem. In Section 3.1, we discuss the stability of the proposed difference scheme. Sections 3.2 and 3.3 present the results of the numerical experiments carried out for model problems, confirming the stability of the proposed method. Finally, we present a discussion of the results obtained in Section 4.

Formulation of the Problem
The system of Saint-Venant equations describing flows without turbulent diffusion components has the following form [35]: where H = h + ζ and x, y are directions of orthogonal coordinate axes, t is time, u, v are velocity components, H is the water depth, ρ is the water density, g is the gravitational acceleration, τ x , τ y are components of wind stress on the water's surface, h is the distance from the formation mark to the reference horizontal surface, ζ is a deviation of the free surface from the reference, C s is Chézy's coefficient and l is the Coriolis parameter. System (1) can be represented in a vector form: where By transforming dependent variables, we obtain the following symmetrical form of the system (2): where c is a real and non-negative number. Systems (1)-(3) are systems of quasi-linear hyperbolic equations. In this paper, the issues of mixed problem statement for a two-dimensional system of Saint-Venant equations are investigated. Depending on the decomposition of the matrix coefficients of the equation system, problem statements are formulated in 25 possible cases with dissipative boundary conditions. For the numerical solution of mixed problems for a system of equations, a new upwind difference scheme of splitting is constructed. Note that the approaches of constructing difference splitting and upwind schemes (separately) are known. The paper proposes a combination of these two approaches. Namely, splitting is carried out in x and y directions. An explicit upwind scheme is used in y-direction and an implicit upwind scheme is used in x-direction. As a result, we get a new upwind difference scheme of splitting. To the best of our knowledge, this form of the difference scheme is considered and investigated for the first time here. The stability of the proposed difference scheme in the energy norm is proved.

A Special Form of the Saint-Venant Equations
We represent the matrix A depending on u as the sum of two matrices: a positive semidefinite matrix A + and a negative semi-definite matrix A − . Consider the following cases.
(1) If u < −c, then Here and below, λ(A) denotes the eigenvalue of the matrix A.
(2) If −c ≤ u < 0, then (5) If c ≤ u, then Similarly, the matrix B depending on v can be represented as the sum of two matrices: a positive semi-definite matrix B + and a negative semi-definite matrix B − .
The matrix A in system (3) can be represented as Similarly, the matrix B in system (3) can be represented as In this case, system (3) takes the following notation: Partitioning (4) of matrices A and B into the sum of two semidefinite matrices allows us to construct an upwind difference scheme. Note that the following inequalities hold for the eigenvalues of symmetric matrices A + i , B + i and i = 1, 5: It is well known [35,36] that it is necessary to specify the initial condition and the following boundary conditions for the system of Saint-Venant equations: At the border of the inflow: (1) For a sub-critical flow, −c ≤ w n < 0, it is necessary to set two conditions, where w n is the projection of the velocity vector onto the outer normal to the boundary.
(2) Three conditions are necessary for a supercritical flow, w n < −c.
At the border of the outflow: (1) Sub-critical flow requires only one condition.
(2) No boundary conditions are required for supercritical flow, c ≤ w n . On a solid boundary, w n = 0, only one condition is required. The authors of [36] established some additional equations, based on the characteristic form of the system of Equation (3), which gives a closed system for computing the solution on these boundaries in combination with the given boundary conditions.

Upwind Implicit Difference Scheme
In the domain Q = {(t, x, y) : 0 ≤ t < +∞, 0 ≤ x ≤ X, 0 ≤ y ≤ Y}, a difference grid with steps ∆t, ∆x and ∆y in directions t, x and y, respectively, is constructed. Here J∆x = X and L∆y = Y.
Denote t κ = κ∆t, κ = 0, K, x j = j∆x, j = 0, J and y l = l∆y, l = 0, L, where K, J and L are some positive integers. The nodal points of the difference grid, i.e., the intersections of straight lines t = t κ , x = x j and y = y l , are denoted by t κ , x j , y l . The set of nodal points of the difference grid is denoted by Q h , where and the values of the numerical solution at the nodal points are denoted by To find a numerical solution to the mixed problem (1)-(3) over the difference grid Q h , we propose the following implicit upwind difference scheme: In the difference scheme (I) and (II), the value of W jl is an auxiliary intermediate value for determining the value of V k+1 jl through V k jl . In the algorithm of finding a numerical solution, the values of W jl are first determined from (I), and then substituted into (II), and thereby the values of V k+1 jl are determined. Consider the following system of linear hyperbolic equations: andū,v,c,f 1 ,f 2 ,f 3 are the given functions. Note that the system of Equation (5) is obtained by linearization from the system of Equation (3) with respect to the stationary solution.
The well-posedness of the mixed problem for a two-dimensional symmetric hyperbolic system with dissipative boundary conditions was proved in [41]. It is not difficult to verify that system (5) is a two-dimensional symmetric hyperbolic system. Further, to ensure the well-posedness of the mixed problem on the domain under consideration, we found the boundary conditions that satisfy dissipativity conditions (non-negative definiteness of some quadratic forms).
For the sake of simplicity, we investigate the stability of solutions in the stationary case (i.e., whenF = 0) with constant coefficients.
Let us dwell in more detail on the implementation of the proposed difference scheme. First, consider the first splitting step in y-direction.
(1) Ifv < −c, then consider the difference scheme as the first splitting step: Decompose the matrixB − 1 as follows: Substitute the decomposition of the matrixB − 1 into (6) and multiply the resulting equality on the left by the matrix Z * P. As a result, scheme (6) reduces to the following form: Next, we introduce the following unknown dependent vector function T = Z * PV. whereT Similarly, one can easily prove the equality The component-by-component notation of the scheme (7) is the following: When setting the boundary conditions, we require that they be dissipative (see [41]). For example, no boundary conditions are required on the left boundary l = 0, and we set the following boundary conditions on the right boundary l = L: The initial data were set at κ = 0 as below: where v 10 (x, y), v 20 (x, y) and v 30 (x, y) are given smooth functions. Note that From (9) at l = 0, it follows that the following inequality holds for the left boundary quadratic form: Similarly, by taking into account the boundary conditions (7b), it follows that the following inequality holds for the boundary quadratic form on the right boundary l = L: (2) If −c ≤v < 0, then we consider the difference scheme as the first splitting step: wherē Decompose the matricesB ± 2 as follows: Substitute the matrix expansionsB ± 2 into (10) and multiply the resulting equality on the left by the matrix Z * P. As a result, scheme (10) will look like the following: Let us rewrite this scheme in terms of a vector of functions T κ jl andT κ jl : The component-wise notation of the difference scheme (11) looks like the following: On the left boundary, at l = 0, it is required to set one boundary condition: and on the right boundary, at l = L, it is required to set two boundary conditions: The initial conditions are set in the form (8).
(3) Ifv = 0, then we consider the following difference scheme as the first splitting step: Decompose the matricesB ± 3 as follows: Substitute the matrixB ± 3 expansions into (12) and multiply the resulting equality by the matrix Z * P. As a result, scheme (12) will have the following form: Let us rewrite this scheme in terms of a vector of functions T κ jl andT κ jl : The component-by-component notation of the scheme (13) looks as follows: On the left boundary, at l = 0, it is required to set one boundary condition: and on the right boundary, at l = L, it is required to set two boundary conditions: The initial conditions are set in the form (8).
(4) If 0 <v <c, then we consider the following difference scheme as the first splitting step:V κ jl =V κ jl − r yB wherē Decompose the matricesB ± 4 as follows: Substitute the matrixB ± 4 expansions into (14) and multiply the resulting equality by the matrix Z * P on the left. As a result, scheme (14) will look as follows: Let us rewrite this scheme in terms of a vector of functions T κ jl andT κ jl : The component-by-component notation of the scheme (15) looks as follows: On the left boundary, at l = 0, it is required to set two boundary conditions: and on the right boundary, at l = L, it is required to specify one boundary condition: The initial conditions are set in the form (8).
(5) Ifc ≤v, then we consider the following difference scheme as the first splitting step: Decompose the matricesB ± 5 as follows: Substitute the matrixB ± 5 expansions into (16) and multiply the resulting equality on the left by the matrix Z * P. As a result, scheme (16) will look as follows: Let us rewrite this scheme in terms of a vector of functions T κ jl andT κ jl : The component-wise notation of this scheme looks as follows: On the left boundary, at l = 0, it is required to set three boundary conditions: and on the right boundary, at l = L, it is not necessary to set the boundary condition. The initial conditions are set in the form (8).
The second splitting step in x-direction is constructed in a similar way.

Stability of the Difference Scheme
To prove Theorem 1, let us first introduce the following notation: We define values q pw κ j , q p w κ j , q = 1, 5, p = 1, 3 for each q = 1, 5. Let us start with the case q = 1. For this purpose, we consider the difference initial boundary value problem (7), (7a), (7b), (8), and the following discrete functions as 1 pw κ j , 1 p w κ j and p = 1, 3: Lemma 1. The following inequality holds under the conditions of Theorem 1: Proof. Denoting the Courant number r yv by ρ 1 , rewrite the first difference equation of system (7a) as follows: Taking into account this equation, we obtain the following expression for 1 1w κ According to the CFL condition, the inequality 1 + ρ 1 > 0 is satisfied; therefore, the inequality −ρ 1 (1 + ρ 1 ) > 0 is valid. Using an algebraic inequality 2ab ≤ a 2 + b 2 , we have: Then, it follows from (18) that Lemma 1 is proved.

Corollary 1.
The following inequalities hold under the CFL condition:

Corollary 2. The following inequalities hold under the CFL condition:
1Wκ Corollary 3. If the condition of dissipativity of the boundary conditions (7b) is satisfied, then the following inequality holds: We now turn to the case q = 2. To this end, we consider the difference initial boundary value problem (11), (11a), (11b), (11c), (8), and we take the following discrete functions as 2 pw κ j , 2 p w κ j and p = 1, 3:

Lemma 3.
Let the CFL conditions of Theorem 1 be satisfied. Then, the solution of the difference initial boundary value problem (11), (11a), (11b), (11c), (8) satisfies the inequality Proof. When denoting the Courant number r y (v +c) for the second difference equation of system (11a) by ρ 2 , we have Taking into account this equation, we obtain the following expression for 2 2w κ Under the CFL condition, the inequality 1 − ρ 2 > 0 is valid. Then, Then, it follows from (19) that Lemma 3 is proved.

Corollary 4.
It follows from Lemmas 2 and 3 that the following inequality holds under the CFL condition: Proof. It follows from the first difference equation of the system (13a) that Square both parts of this equality, multiply by ∆y and sum over l from 1 to L − 1 to have which implies the assertion of the lemma.

Lemma 6.
Let the CFL conditions of Theorem 1 be satisfied. Then, the following inequalities hold for the solution of the difference initial boundary value problem (13), (13a), (13b), (13c), (8): The proof of Lemma 6 follows directly from Lemma 2.

Corollary 7.
It follows from Lemma 7 that the following inequality holds under the CFL condition: Corollary 8. If the dissipativity condition for the boundary conditions (15b) and (15c) is satisfied, then the following inequality holds according to Corollary 7: Finally, consider the case q = 5, the difference initial boundary value problem (16), (16a), (16b), (8) and the following discrete functions as 5 pw κ j , 5 p w κ j and p = 1, 3: The proof of Lemma 8 follows directly from Lemma 7.

Corollary 9.
It follows from Lemma 8 that the following inequality holds under the CFL condition:

Corollary 10.
If the dissipativity condition for the boundary condition (16b) is satisfied, then the following inequality holds according to Corollary 9: Corollary 11. If the condition of dissipativity of the boundary condition (16b) is satisfied, then the following inequality holds: In a similar way, one can prove that W κ+1 ≤W κ . Consequently, we arrive at the chain of inequalities which imply the stability of the proposed difference scheme in the norm √ W k . Thus, Theorem 1 is proved. Remark 1. It follows from Theorem 1 and Lax equivalence theorem by virtue of the linearity of the mixed problems for the system of Saint-Venant equations that the proposed difference scheme convergences with the first order in both the spatial and temporal variables.

Numerical Example 1
The numerical experiment was carried out on a computer with the following technical characteristics: Intel(R) Core(TM) i7-8700, 8Gb, Intel(R) UHD Graphics 630, Windows 7 Pro and PTC Mathcad Prime 7. The elapsed time for calculating numerical example 1 was 43 ms, and for numerical example 2 was 34 ms.
This section provides the results of a numerical experiment to validate the theoretical analysis for the difference scheme proposed in Section 2.3. The numerical experiment was carried out with the assumption thatū = u * ,v = v * ,c = c * under the fulfillment of the inequalitiesc ≤v andc ≤ū.
The numerical experiment was conducted for the following parameters of the problem: the length and width of the domain were X = 100 m and Y = 6 m, respectively; final time was T = 1 sec. The steps in the spatial variables and time were chosen to be ∆x = 0.5, ∆y = 0.2 and ∆t = 0.02; Coriolis parameter l = −10.3469869 × 10 −5 ; and Chézy's coefficient was C S = 35.5383100397. In addition, it was assumed that u = 0.3 sin(x + y) + 2.9,v = 0.2 sin(x + y) + 2.9, c = 0.1 sin(x + y) + 2.5.
One can verify that the selected input data satisfy the conditions of Theorem 1. As can be seen in Figures 1 and 2, the proposed difference scheme makes it possible to obtain a stable approximate solution of the problem under consideration.  According to the plots of Figures 1 and 2, it can be seen that with the well-posed statement of the initial boundary value problem, the obtained numerical solution comes to a stationary state in one time step.

Numerical Example 2
The numerical experiment was carried out in the case of −c ≤v ≤ 0 and −c ≤ū ≤ 0. The numerical test was conducted for the following parameters of the problem: the length and width of the domain were X = 50 m and Y = 6 m, respectively; final time was T = 1 s. The steps in the spatial variables and time were chosen to be ∆x = 0.5, ∆y = 0.2 and ∆t = 0.02; Coriolis parameter was l = −10.3469869 × 10 −5 ; and Chézy's coefficient was C S = 35.5383100397. In addition, it was assumed that u = −(0.9 sin(x + y) + 1),v = −(0.2 sin(x + y) + 0.9), c = 0.1 sin(x + y) + 2.
The results of the numerical experiments are depicted in Figures 3 and 4. (a) Velocity in x-direction, u (b) Velocity in y-direction, v

Discussion
Based on the theoretical analysis carried out and the results of the computational experiments, the following conclusions can be outlined.
(1) The issues of the well-posed statement of a mixed problem for a system of twodimensional Saint-Venant equations have been studied.
(2) To solve a mixed problem for the two-dimensional Saint-Venant equation, a new upwind splitting scheme was proposed. The stability of the proposed difference scheme according to energy norms was strictly proved.
(3) The well-posedness of the discrete analogue of the mixed problem for the twodimensional Saint-Venant equation was established.
(4) Computational experiments carried out on the examples of the test problems have shown that the spatial-splitting upwind-difference scheme proposed in this work allows obtaining a stable approximate solution to a system of two-dimensional Saint-Venant equations.
Since the work was mainly of a theoretical nature, the main attention was paid to obtaining sufficient conditions for the stability of the proposed upwind scheme. In subsequent studies, the authors intend to solve a more realistic example based on the constructed numerical scheme.