To Split or Not to Split, That Is the Question in Some Shallow Water Equations
aa r X i v : . [ m a t h . NA ] N ov To Split or Not to Split, That Is the Question in SomeShallow Water Equations
Vicente Mart´ınez ∗ Departamento de Matem´aticas & Instituto de Matem´aticas y sus Aplicaciones de Castell´on (IMAC)Universitat Jaume I, Campus de Riu Sec, 12071 Castell´o, Spain
November 17, 2012
Abstract
In this paper we analyze the use of time splitting techniques for solving shallowwater equation. We discuss some properties that these schemes should satisfy so thatinteractions between the source term and the shock waves are controlled. This papershows that these schemes must be well balanced in the meaning expressed by Greenbergand Leroux [5]. More specifically, we analyze in what cases it is enough to verify an
Approximate C-property and in which cases it is required to verify an
Exact C-property (see [1], [2]). We also include some numerical tests in order to justify our reasoning.
Key words. splitting schemes, source terms, shallow water equations.
AMS subject classifications.
In this paper, our interest is to analyze time splitting schemes on conservation lawswith source terms, also called balance laws. A prototype, in one space dimension and undercertain regularity hypotheses, is given by the following system of partial differential equations ( W ( x, t ) t + F ( W ( x, t )) x = G ( x, W ( x, t )) , ( x, t ) ∈ R × R + ,W ( x,
0) = W ( x ) , x ∈ R , (1)where W : R × R + → R m is the vector of conserved variables, F : R m → R m is the vectorof fluxes and G : R m +1 → R m is the source term.Recently, there has been some controversy related to the application of time splittingtechniques on hyperbolic equations involving solutions with shock waves, as is the case of ∗ E-mail address: [email protected]. W ( x, t ) t + F ( W ( x, t )) x = 0 , (2)and the ordinary differential equation W ( x, t ) t = G ( x, W ( x, t )) . (3)LeVeque notices in [9] that such schemes can easily fail by the presence of shock wavesin solving (2). These shock waves involve large changes in the solution which can not becaptured in solving (3).On the other hand, some authors such as Ma, Sun and Yin (see [12]) use a time in-tegrating scheme with two-step predictor-corrector sequence quite successfully. Striba usesplitting techniques (see [13]) in meteorology models on the term that represent de Coriolisacceleration. In addition, Wicker and Skamarock (see [16]) use time-splitting methods forintegrating the elastic equations.In 1994, Berm´udez and V´azquez (see [1], [2]) introduce the concept of Exact C-property and
Approximate C-property in order to identify numerical schemes with an acceptable levelof accuracy in the resolution of shallow water equations (well-balanced scheme). Goingdeeper into well-balanced scheme idea, we can find the work of Greenberg and Leroux [5], inwhich they propose a numerical scheme that preserves a balance between the source termsand internal forces due to the presence of shock waves.Another point of view, more recently, it is provided by Lubich (see [11]), who givesan error analysis of Strang-type splitting integrators for nonlinear Schr¨odinger equations.Holdahl, Holden and Lie use an adaptive grid refinement and a shock tracking techniqueto construct a front-tracking method for hyperbolic conservation laws. They combine theoperator splitting to study the shallow water equations (see [6]). Holden, Karlsen, Risebroand Tao (see[8]) show that the Godunov and Strang splitting methods converge with theexpected rates if the initial data are sufficiently regular. Finally, in this way the reader canfind a deep study of splitting methods for partial differential equations in [7], where someanalysis of conservation and balance laws are included.We base our analysis on the ideas presented in [1], [2] and [5]. From these studies, itfollows that the numerical scheme used in solving (2) and the numerical scheme used insolving (3) cannot be whatever, even though they have a high degree of accuracy. Thesemust be balanced so that interactions between the source term in (3) and the shock waves in(2) are controlled. In this framework, we will analyze conditions to be verified by splittingschemes in order to avoid spurious oscillations, which are created in this type of equations.More specifically, in which cases it is enough to verify an
Approximate C-property and inwhich cases it is required to verify an
Exact C-property .The rest of this paper is organized as follows: in Section 2, we introduce the governingequations for the one dimensional shallow water model. We also analyze two kinds of timesplitting schemes to identify which conditions must satisfy a scheme for solving this typeof equations. In Section 3, we show the four test problems that we will use to test the2erformance of the schemes described above. In Section 4, we present the numerical resultswe have obtained. Finally, in Section 5 we reason the final conclusions we have reached byusing the results obtained with the presented schemes.
In this section we consider Eq. (1) with W = (cid:18) hq (cid:19) , F ( W ) = qq h + 12 gh , G ( x, W ) = (cid:18) − ghb ′ ( x ) (cid:19) , (4)where the unknowns of the problem are: the water height is h and the flow per unit lengthis q = hu. Here u is the average vertical speed in the direction of the axis x (see Fig.1). F isthe flux of conservative variables and g = 9 . ms − is the acceleration due to gravity. Thesource term G models the bottom variation given by the function b ( x ) . free surface reference level ✲ u ( x, t ) • A b ( x ) h ( x, t ) x Fig. 1: Shallow water variables.Let us consider numerical solvers based upon the decomposition F ( W ) x = A ( W ) W x , where A ( W ) = − q h + gh qh is the Jacobian matrix of F ( W ) . (5)3o that, system (1)-(4) is hyperbolic ( h > , then A = X Λ X − , whereΛ = (cid:18) λ λ (cid:19) , X = (cid:18) λ λ (cid:19) , X − = λ − λ (cid:18) λ − − λ (cid:19) , (6)where λ = qh + p gh and λ = qh − p gh. Our analysis needs the definition of some conservation properties given by Berm´udez andV´azquez in [2]. Since only the source term involves the bed slope, an inadequate choice ofthe numerical schemes can give spurious oscillations. So, these conservation properties try toidentify which kind of schemes have good behavior in equations with source term. Berm´udezand V´azquez characterize the good behavior of the numerical scheme in the manner in whichthe scheme approximates a steady solution representing the state of water at rest. Theyintroduce the stationary problem (
Problem SP ) given by q ( x, t ) = 0 and h ( x, t ) = H ( x ) anddefine the following conservation properties: Definition . We say that a scheme satisfies the Exact C-Propertyif it is exact when applied to the Problem SP.
Definition . We say that a scheme satisfies the Approxi-mate C-Property if it is accurate to the order
Θ(∆ x ) when applied to Problem SP. When a numerical scheme does not satisfy any of these conservation properties then thepropagation of spurious oscillations is also present in non stationary problems.
The differential formulation of the homogeneous equation given in (2) does not admit dis-continuous solutions. So, since these solution are physically relevant in this context, we needa suitable formulation of the problem to support discontinuous solutions. In this sense, it isusual to consider the following integral formulation Z ( W dx − F ( W ) dt ) = 0 . (7)The numerical schemes use usually (7) in order to approximate (2). To do that, it isintroduced a control volume in the space ( x, t ) of dimensions ∆ x × ∆ t . Next, it is evaluatedthe integral (7) in this volume control Z x j +1 / x j − / ( W ( x, t n +1 ) − W ( x, t n )) dx + Z t n +1 t n ( F ( W ( x j +1 / , t )) − F ( W ( x j − / , t ))) dt = 0 . x we obtain1∆ x Z x j +1 / x j − / W ( x, t n +1 ) dx = 1∆ x Z x j +1 / x j − / W ( x, t n ) dx − ∆ t ∆ x h t Z t n +1 t n F ( W ( x j +1 / , t )) dt − t Z t n +1 t n F ( W ( x j − / , t )) dt i . Thus, we deduce the conservation formula W n +1 j = W nj − ∆ t ∆ x [ F j +1 / − F j − / ] , (8)where W nj is an average W nj = 1∆ x Z x j +1 / x j − / W ( x, t n ) dx at time t = t n inside the interval I j = [ x j − / , x j +1 / ]whose length is ∆ x = x j +1 / − x j − / . The flux in (8) can be interpreted as the average in time of the physical flux, i.e., F j +1 / ≈ t Z t n +1 t n F ( W ( x j +1 / , t )) dt. (9)Conservative numerical methods for (2) are based in (8), and they are determined by theexpression of the numerical flux F j +1 / .The time splitting numerical schemes (see [14]), in each time step, act as follows: takinginto account the initial condition W n , we solve (2) and we obtain c W n +1 . Then, by using theinitial condition c W n +1 , we solve (3) and obtain W n +1 . More specifically, if we use the operators • A ( W n ) = W n +1 such as W n → ( W t + F ( W ) x = 0 , ( x, t ) ∈ R × [ t n , t n +1 ] ,W ( x, t n ) = W n , x ∈ R . ) → c W n +1 . • S ( W n ) = W n +1 such as W n → ( W t = G ( x, W ) , ( x, t ) ∈ R × [ t n , t n +1 ] ,W ( x, t n ) = c W n +1 , x ∈ R . ) → W n +1 . S (∆ t ) · A (∆ t ) ( W n ) = W n +1 . In the first step of the splitting procedure, we solve (2) at each time step. To do thatand for the sake of simplicity, we consider the Q -scheme of van Leer (see [15]) which uses amatrix Q satisfying some properties and the numerical fluxes F j − / = φ ( W nj − , W nj ) , F j +1 / = φ ( W nj , W nj +1 ) , (10)where the numerical flux function φ is given by φ ( U, V ) = F ( U ) + F ( V )2 − | Q ( U, V ) | ( V − U ) . (11)A possible choice of the matrix Q can be the Jacobian (5) of the system (4) evaluated atthe arithmetic mean, i.e. Q ( U, V ) = A (cid:18) U + V (cid:19) . (12)So that, we have (cid:12)(cid:12)(cid:12) Q ( W nj ± ) (cid:12)(cid:12)(cid:12) = X ( W nj ± , W nj ) (cid:12)(cid:12) Λ( W nj ± , W nj ) (cid:12)(cid:12) X − ( W nj ± , W nj ) , (13)where X ( W nj ± , W nj ) , (cid:12)(cid:12) Λ( W nj ± , W nj ) (cid:12)(cid:12) and X − ( W nj ± , W nj ) are evaluated at W nj ± + W nj . So, we obtain the numerical fluxes F j ± / = F ( W nj ± ) + F ( W nj )2 − (cid:12)(cid:12)(cid:12) Q ( W nj ± ) (cid:12)(cid:12)(cid:12) ( ∓ W nj ± W nj ± ) . (14)The second step of the procedure is to solve (3). To do that, we get the solutions ofthe homogeneous equation (2): ˆ h ( x j , t n +1 ) and ˆ q ( x j , t n +1 ) and we solve the following initialvalue ODE problem for each x j ddt h ( x j , t ) q ( x j , t ) ! = − gh ( x j , t ) b ′ ( x j ) ! ; t ∈ [ t n , t n +1 ] , h ( x j , t n ) q ( x j , t n ) ! = ˆ h ( x j , t n +1 )ˆ q ( x j , t n +1 ) ! . (15)The first equation has the solution: h ( x j , t n +1 ) = ˆ h ( x j , t n +1 ) . (16)6o solve the second one, we calculate dq ( x j , t ) dt = − gh ( x j , t ) b ′ ( x j ) ,q ( x j , t ) = − Z tt n gh ( x j , s ) b ′ ( x j ) ds ; t ∈ ( t n , t n +1 ] . Therefore, q ( x j , t n +1 ) = ˆ q ( x j , t n +1 ) − g · b ′ ( x j ) · Z t n +1 t n h ( x j , s ) ds. Finally, by using a trapezoidal rule, we obtain q ( x j , t n +1 ) = ˆ q ( x j , t n +1 ) − g · b ′ ( x j ) · ∆ t h ( x j , t n ) + ˆ h ( x j , t n +1 )) . (17) Remark 1
We denote the time splitting scheme given by (14)-(16)-(17) as
Q-tra1 . Proposition The numerical scheme
Q-tra1 satisfies an approximate C-property.
Proof
Problem SP satisfies λ = √ gh and λ = −√ gh. In addition (see Fig. 1), b ( x ) + h ( x, t ) = b ( x ) + H ( x ) = A. Therefore, b ( x j ) = A − H ( x j ) . (18)In order to solve (2) with (14), we compute (cid:12)(cid:12)(cid:12) Q ( W nj + ) (cid:12)(cid:12)(cid:12) = 1 λ − λ | λ | λ − | λ | λ −| λ | + | λ | λ λ ( | λ | − | λ | ) − λ | λ | + λ | λ | ! = r g h j +1 + h j r g h j +1 + h j . Then, F nj + = 12 ( h j +1 − h j ) r g h j +1 + h j
212 ( h j +1 + h j ) . Analogously, we obtain F nj − = 12 ( h j − h j − ) r g h j + h j −
212 ( h j + h j − ) . c W n +1 j = ˆ h j ˆ q j = h j + ∆ t x ( h j +1 − h j ) r g h j +1 + h j − ( h j − h j − ) r g h j + h j − ! − g ∆ t x (cid:0) h j +1 − h j − (cid:1) Next, to have the numerical approximation of Eq. (3), we use (16) and (17). By substi-tuting, we obtain h n +1 j = h j + ∆ t x ( h j +1 − h j ) r g h j +1 + h j − ( h j − h j − ) r g h j + h j − ! (19)and q n +1 j = − g ∆ t x (cid:0) h j +1 − h j − (cid:1) − g · b ′ ( x j ) · ∆ t h ( x j , t n ) + ˆ h ( x j , t n +1 )) . (20)From (19), we deduce that h n +1 j = h j is accurate to the order Θ(∆ x ) . And from (20)and (18) we have that q n +1 j = 0 is exact when we take b ′ ( x j ) = b ( x j +1 ) − b ( x j − )2 ∆ x . Thus, the proof ends. (cid:3)
Remark 2
If we use another convergent quadrature rule to approximate the integral Z t n +1 t n h ( x j , s ) ds, then Proposition
It is well-known that conservation laws with source terms can be solved with high accuracyand discontinuities are well captured by using upwind schemes. For instance, you can see theworks of LeVeque and Yee [10] or V´azquez-Cend´on [15]. Furthermore, this kind of schemeshave bigger stability regions than schemes using centered approximations of the source term8see [2]). In this section, we will use some approaches in order to build an upwind schemeusing some ideas contained in [15] and making changes to solve (3).In order to solve (3), we need to introduce two numerical source functions, G L on theleft and G R on the right, to upwind the source term. We also use two matrices D L and D R with the aim of making clear the upwind process. Below we explain the process in detail.We will use, in each time step, the average of solutions of the following equations: ( W t = G L ( x, W ) , ( x, t ) ∈ R × [ t n , t n +1 ] ,W ( x, t n ) = c W n +1 , x ∈ R . (21)and ( W t = G R ( x, W ) , ( x, t ) ∈ R × [ t n , t n +1 ] ,W ( x, t n ) = c W n +1 , x ∈ R , (22)where G L ( x, W ) = D L ( W L ) G ( x, W ) and G R ( x, W ) = D R ( W R ) G ( x, W );and functions W L and W R change in each computacional cell in the discretization problem W L = W j + W j − W R = W j +1 + W j . Thus, we obtain the solution of (21)-(22), W n +1 L and W n +1 R respectively. Then, in eachtime step, we take as solution of (3) the following average W n +1 = W n +1 L + W n +1 R . In other words, for the discretization of the problem, we can substitute the ODE in (15)by ddt h ( x j , t ) q ( x j , t ) ! = D L ( W L j ) G ( x j , W j ) = d L d L d L d L ! · − gh ( x j , t ) b ′ ( x j ) ! = − d L gh j b ′ j − d L gh j b ′ j ! . Therefore, instead of (15), we have two ODE ddt h j q j ! = − d L gh j b ′ j − d L gh j b ′ j ! ; t ∈ [ t n , t n +1 ] , h ( x j , t n ) q ( x j , t n ) ! = ˆ h ( x j , t n +1 )ˆ q ( x j , t n +1 ) ! . (23)9nd analogously from (22) we obtain ddt h j q j ! = − d R gh j b ′ j − d R gh j b ′ j ! ; t ∈ [ t n , t n +1 ] , h ( x j , t n ) q ( x j , t n ) ! = ˆ h ( x j , t n +1 )ˆ q ( x j , t n +1 ) ! . (24)where d R and d R are coefficients corresponding to matrix D R . Finally, we compute the solutions of (23) and (24), W n +1 Lj and W n +1 Rj . Remark 3
Matrices D L and D R must be chosen in coordination with the numerical methodused in solving (2) to get a well balanced scheme. For example, regarding (14), we can take D L = ( I + | Q | Q − ) and D R = ( I − | Q | Q − ) . (25) In this way we obtain consistency, since D L + D R = 2 I, we have G L + G R D L G + D R G D L + D R ) G G. Remark 4
We denote the time splitting scheme given by (14)-(23)-(24)-(25) as
Q-tra2 . Proposition The numerical scheme
Q-tra2 satisfies an exact C-property.
Proof
Acting in a similar way to
Proposition
1, we have c W n +1 j = ˆ h j ˆ q j = h j + ∆ t x ( h j +1 − h j ) r g h j +1 + h j − ( h j − h j − ) r g h j + h j − ! − g ∆ t x (cid:0) h j +1 − h j − (cid:1) . On the other hand, D ( W L j ) = q g hj + hj − q g h j + h j − and D ( W R j ) = − q g hj +1+ hj − q g h j +1 + h j . Now, taking in (23) b ′ ( x j ) = b ( x j ) − b ( x j − )∆ x and h ( x j , t ) = h j + h j − h L j = ˆ h j − ∆ t ∆ x ( b j − b j − ) r g h j + h j − q L j = ˆ q j − ∆ t x g ( b j − b j − )( h j + h j − ) (27)Acting in a similar way to (24), we have h R j = ˆ h j − ∆ t ∆ x ( b j +1 − b j ) r g h j +1 + h j q R j = ˆ q j − ∆ t x g ( b j +1 − b j )( h j +1 + h j ) (29)In addition, from (18) we have b j = A − h j , then h n +1 j = h L j + h R j h j and q n +1 j = q L j + q R j . This concludes the proof. (cid:3)
In this section, we will discuss four test problems for shallow water equations. The first testproblem is the dam-break problem with variable topography, which is one of the most basicproblems with source term. The second one is a stationary problem, also with source termthat represents a smooth bottom. The third one represents a tidal wave flow propagatingover an irregular topography. This test was discussed by Berm´udez and V´azquez [2] and itis one of the most popular test to check the performance of a numerical scheme which triesto be effective in solving shallow water equations. It represents a severe test regarding toirregular topography and the long time over which is applied. Finally, the last test representsa numerical simulation of a tidal wave on the shoreline with friction effects, which introducesa wet/dry front. 11 est 1: Dam-break with variable topography
We consider the one dimensional shallow water equations (4) on the domain ( x, t ) ∈ [0 , × R + . The initial data are b ( x ) =
18 cos (cid:18) π (cid:18) x − (cid:19)(cid:19) + 1 , < x < , , otherwise; W ( x,
0) = − b ( x )0 ! , x < , . − b ( x )0 ! , x > . Fig. 2 shows numerical results of Test. 1 for schemes Q-tra1 and Q-tra2. In this test, wehave chosen cf l = 0 . t = 0 . Initial condition t = 0 . ◦◦◦◦◦◦ Q-tra2 b ( x ) Fig. 2: Test 1. Initial condition and numerical result for t = 0 . . Test 2: Stationary problem with smooth bottom
We consider also Eq. (4) on the domain ( x, t ) ∈ [0 , × R + . The bottom function b ( x ) is thesame of test 1 and the initial data are W ( x,
0) = − b ( x )0 ! , ≤ x ≤ . This problem has a stationary solution W ( x, t ) = W ( x, , ≤ x ≤ t > .
12n this test, we have chosen cf l = 0 . t = 0 .
25 on 50 computational cells, it showsspurious oscillations for Q-tra1 (see Fig. 3). However, Q-tra2 scheme has a good behavior. ——– Analytical surface ◦◦◦◦◦◦
Q-tra1 ⋄⋄⋄⋄⋄⋄
Q-tra2 b ( x ) Fig. 3: Test 2. Numerical result for t = 0 . . Test 3: Tidal wave flow with irregular topography
We consider Eq. (4) on the domain ( x, t ) ∈ [0 , × [0 , , where length is in metersand time in seconds. The bottom function b ( x ) is shown in Fig. 4, and the initial data are W ( x,
0) = H ( x )0 ! , H ( x ) = H (0) − b ( x ) , ≤ x ≤ H (0) = 16 . Moreover, we have the following boundary conditions h (0 , t ) = H (0) + ϕ ( t ) ,ϕ ( t ) = 4 + 4 sin (cid:18) π (cid:18) t − (cid:19)(cid:19) ,q (1500 , t ) = 0 , ≤ t ≤ . Function ϕ ( t ) simulates a tidal wave of 4 m amplitud.13or this problem, we can obtain an asymptotic analytical solution (see [2]) given by h ( x, t ) = H ( x ) + ϕ ( t ) and it satisfies h ( x, , ≤ x ≤ . In order to check the behaviour of Q-tra2 scheme, we have chosen cf l = 0 . t = 10800 seconds and acomplex topography, however an excellent performance of this scheme is demonstrated (seeFig. 4). ◦◦◦◦◦◦ Q-tra2——– Analytical surface b ( x ) Fig. 4: Test 3. Tidal wave flow over an irregular topography using Q-tra2.
Test 4: Tidal wave on the shoreline with friction effects
In this test we take in account the bottom friction. Manning law is used to model frictionbetween fluid and bottom, so a new element appears in the source term. As a consequence,the source term is given by G ( x, W ) = − ghb ′ ( x ) − gqM (cid:12)(cid:12)(cid:12) qh (cid:12)(cid:12)(cid:12) h − , where M is the Manning coefficient.We consider the domain ( x, t ) ∈ [0 , × R + and the bottom function b ( x ) = ( . x + 0 . , ≤ x ≤ , . x −
3) + 0 . , ≤ x ≤ . Here, length is in meters, time in seconds and the reference level is located at 40 cm.The tidal wave is simulated by introducing a discharge q = 0 . m /s in the boundary x = 014uring the first 0 . x = 6 is also introduced a vertical wall condition to preserve the mass conservationin the domain.On the other hand, this test involves a wet/dry front, which complicates the stability ofthe scheme. Specific treatment is required to remove the spurious oscillations, the readercan find all detailed in [3] and [4]. For brevity, we summarize here two main actions of thistreatment:1.- Redefinition of the discretized bottom function to avoid the appearance of spuriouspressure forces. • If I j is dry, I j − wet and h j − + b j − < h j + b j , then b j = b j − + h j − . • If I j is wet, I j − dry and h j − + b j − > h j + b j , then b j = b j − − h j . • If I j is dry, then q n +1 j = 0 . • If I j is wet, estimate of q n +1 j < I j − is dry, then q n +1 j = 0 . • If I j is wet, estimate of q n +1 j > I j +1 is dry, then q n +1 j = 0 . In addition, friction causes spurious oscillations in q nj . However splitting schemes canavoid these effects by using a semi-implicit discretization. Below, we summarize the details.In order to compute W n +1 Lj , we consider Eq. (23) with friction ddt h j q j ! = − d gh j b ′ j − d gq j M (cid:12)(cid:12)(cid:12)(cid:12) q j h j (cid:12)(cid:12)(cid:12)(cid:12) ( h j ) − − d gh j b ′ j − d gq j M (cid:12)(cid:12)(cid:12)(cid:12) q j h j (cid:12)(cid:12)(cid:12)(cid:12) ( h j ) − ; t ∈ [ t n , t n +1 ] , h ( x j , t n ) q ( x j , t n ) ! = ˆ h ( x j , t n +1 )ˆ q ( x j , t n +1 ) ! . (30)Solving, we have h L j = ˆ h j − d ∆ tg h j + h j − b ′ j − d ∆ tg (cid:18) h j + h j − (cid:19) − q j + q j − M (cid:12)(cid:12)(cid:12)(cid:12) q j + q j − h j + h j − (cid:12)(cid:12)(cid:12)(cid:12) (31)15nd q L j = ˆ q j − d ∆ tg h j + h j − b ′ j − d ∆ tg (cid:18) h j + h j − (cid:19) − q L j M (cid:12)(cid:12)(cid:12)(cid:12) q j + q j − h j + h j − (cid:12)(cid:12)(cid:12)(cid:12) . (32)Now, from Eq. (32) we obtain that q L j = ˆ q j − d ∆ tg h j + h j − b ′ j d ∆ tg (cid:18) h j + h j − (cid:19) − M (cid:12)(cid:12)(cid:12)(cid:12) q j + q j − h j + h j − (cid:12)(cid:12)(cid:12)(cid:12) . (33)On the other hand, in order to compute W n +1 Rj , we act analogously and obtain h R j = ˆ h j − d ∆ tg h j +1 + h j b ′ j +1 − d ∆ tg (cid:18) h j +1 + h j (cid:19) − q j +1 + q j M (cid:12)(cid:12)(cid:12)(cid:12) q j +1 + q j h j +1 + h j (cid:12)(cid:12)(cid:12)(cid:12) (34)and q R j = ˆ q j − d ∆ tg h j +1 + h j b ′ j +1 d ∆ tg (cid:18) h j +1 + h j (cid:19) − M (cid:12)(cid:12)(cid:12)(cid:12) q j +1 + q j h j +1 + h j (cid:12)(cid:12)(cid:12)(cid:12) . (35) Remark 5
We denote the time splitting scheme given by (14)-(31)-(33)-(34) and (35) as
Q-tra3 . Test. 4 results using Q-tra3 scheme are shown in Fig. 5. This test has all the ingredientsof a real problem in shallow water: irregular topography, friction effects between bottomand water, and advancing front wet/dry. In Fig. 5, we can see the progress of the tidal wavefor different time values: for t = 1 the wave goes to head, in t = 2 the wave reaches theshoreline, in t = 3 the wave hits the wall and in t = 4 , it returns toward x = 0 . In t = 5 , wecan see how the wave returns back with negative velocity and it crashes with the wave from x = 0 with positive velocity. However, in front of this critical situation, the behavior of thescheme is stable.In this test, we have used a Manning coefficient M = 0 .
015 and cf l = 0 . t = 1 t = 2free surface wall b ( x ) free surface wall b ( x ) −1 0 1 2 3 4 5 6 7−0.200.20.40.60.81 −1 0 1 2 3 4 5 6 7−0.200.20.40.60.81 t = 3 t = 5free surface wall b ( x ) free surface wall b ( x ) Fig. 5: Test 4. Tidal wave evolution for different times using Q-tra3.
Some upwind time splitting schemes have been developed for solving 1-D shallow waterequations. In some test problems, our analysis shows that it is not enough to verify anapproximate C-property to obtain good approximations; we need that an exact C-propertyto be verified. Perhaps, for some undemanding tests, it is not necessary, but when thetest is demanding, then schemes that only satisfy an approximate C-property can bring tonumerical instabilities.We have shown that time splitting schemes for solving shallow water equations must bewell balanced in the following sense: the scheme used to solve the homogeneous equationEq. (2) and the scheme used to integrate the time-dependent ODE Eq. (3) cannot bewhatever. They must be linked so that an exact C-property is satisfied. Apparently, Q-tra2is computationally very expensive. But this is not true, since carried out calculations toobtain the numerical approximation of Eq. (2) are used to obtain D L and D R matrices.17oreover, time splitting schemes could give us the possibility to choose a semi-implicitor implicit solver in order to obtain stability when Eq. (3) is integrated. In Test. 4, we haveused a semi-implicit scheme that allowed us to obtain stability when shallow water modeltakes into account friction effects between fluid and bottom.In summary, we have presented three time splitting schemes with different properties,the second of which is the best suited to treat severe test on shallow water equations. Thesemi-implicit Q-tra3 scheme shows the versatility that can provide splitting techniques forsolving challenges of stiffness that present some partial differential equations. Acknowledgements