Extension of WAF type methods to non-homogeneous Shallow Water Equations with pollutant

Extension of WAF type methods to non-homogeneous Shallow Water Equations with pollutant E. D. Fern´andez-Nieto, G. Narbona-Reina ∗ July 27, 2015 Abstract This paper deals with the extension of the WAF method to discretize Shallow Water Equations with pollutants. We consider two different versions of the WAF method, by approximating the intermediate waves using the flux of HLL or the direct approach of HLLC solver. It is seen that both versions can be written under the same form with different definitions for the approximation of the velocity waves. We also propose an extension of the method to non-homogeneous systems. In the case of homogeneous systems it is seen that we can rewrite the third component of the numerical flux in terms of an intermediate wave speed approximation. We conclude that – in order to have the same relation for non-homogeneous systems – the approximation of the intermediate wave speed must be modified. The proposed extension of the WAF method preserves all stationary solutions, up to second order accuracy, and water at rest in an exact way, even with arbitrary pollutant concentration. Finally, we perform several numerical tests, by comparing it with HLLC solver, reference solutions and analytical solutions. Short title : the WAF method for non-homogeneous SWE with pollutant. Keywords : Finite Volume Method, well-balanced, upwinding, shallow water, source terms, WAF, HLLC, pollutant. Subject Classifications : AMS (MOS) : 65N06, 76B15, 76M20, 76N99. 1 Introduction In this paper we study the extension of the WAF (Weighted Average Flux) method to Shallow Water Equations (SWE in what follows) with topography and transport of pollutants. ∗Departamento de Matem´atica Aplicada I, E.T.S. Arquitectura. Universidad de Sevilla. Avda. Reina Mercedes 2. 41012 Sevilla, Spain ([email protected], [email protected]). 1 This method was introduced by Toro in [29]. WAF is a one-step method of Godunovtype to solve hyperbolic conservation laws that achieves second order accuracy by averaging the solution of the conventional Riemann problem with piece-wise constant initial data. There are many works related to the applications of the WAF method, for example to the Shallow Water and Euler equation (see [30, 31]). The extension to the multidimensional case was performed by Billet and Toro in [3]. In [25] the WAF method is combined with WENO reconstructions. As it is well known, due to Godunov’s theorem, linear schemes with high order accuracy generate spurious oscillations near large gradients of the solution. To avoid this problem, the WAF method is used with a flux limiter function, getting a non-linear TVD (Total Variation Diminishing) scheme of second order accuracy. This property, and the fact that it does not need characteristic decomposition, make the WAF method a quick solver with a good order of accuracy for solving Shallow Water equations with pollutant. With respect to high order numerical schemes, we can find the ADER (Advection- Diffusion-Reaction) approach, developed by Toro et al. (See for example [33]). The ADER is a finite-volume Godunov-type approach that allows to obtain schemes with arbitrary order of accuracy, involving the use of data reconstruction and a first order flux solver. We can also find the extension of ADER to the multidimensional case [24, 37, 27, 36] and to the non-homogeneous case [36, 39, 34]. Regarding the WAF method, we can find in [35] the introduction of the ADER-TVD approach. In this version of ADER methods a TVD flux is used to construct the ADER scheme instead of the conventional use of a first-order monotone flux. Specifically, the WAF method is taken as the only direct method of second order accuracy. In the case of no reconstruction, the ADER-WAF scheme provides a second order TVD scheme, obtaining a superior scheme to the conventional ADER. There are different ways of extending a numerical scheme to non-homogeneous hyperbolic systems. The pioneering work of Roe [23] relates the choice of a quadrature formula for the source term with the property of preserving stationary solutions. Different extensions have been done: see for instance Bermudez-V´azquez [2], Greenberg-Leroux [16], LeVeque [19], Zhou et al. [41], Chac´on et al. [8], Perthame-Simeoni [22]. These techniques focus on the problem of a hyperbolic system with source term (or balance law). Another different form to tackle the problem is to rewrite the balance law as a hyperbolic system in non-conservative form by adding a new equation. For example, for SWE with topography, if zbdenotes the function corresponding to the bed of the channel, the equation ∂tzb= 0 is added. Then a numerical scheme for the non-conservative system is applied and finally the scheme is rewritten in function of the original unknowns. See for instance LeFloch [18], Par´es-Castro [21], Gosse [15]. In this paper we adopt the strategy of studying the problem in form of a balance law. We follow the technique introduced in [12] where the concept of ‘asymptotically wellbalanced method’ is introduced. The purpose is for the method to preserve all stationary solutions up to second order over a set whose measure of the complementary tends to zero, when the space step also tends to zero (see Definition 3.1). To perform the extension of the method we must localize the corresponding numerical diffusion that must be neglected 2 for stationary solutions. The extension of the WAF method put forward in this paper is actually asymptotically well-balanced. Moreover it exactly preserves water at rest. The case of water at rest is the only stationary solution that allows to have an arbitrary pollutant concentration, it also introduces an additional difficulty in order for this solution to be preserved. In [17] Kurganov and Petrova present a central-upwing scheme using a second order state reconstruction verifying that the gradient of the height of the water column at intercells coincides with the gradient of the free water surface. Using this reconstruction it is proved that the proposed method is positivity preserving. Other positivity preserving techniques applied to VFROE with hydrostatic reconstruction [1] are introduced in [20] and [5]. We cannot prove the method positivity preserving for general situations. We have not found negative values in the numerical tests. Specifically in Tests 4, 5 and 6 dry/wet fronts and vacuum generation appear. Using a dry/wet technique it is possible to prove positivity of the height in this front, in some cases by using a reduced CFL condition. For the tests to be developed in Section 4 we use the dry/wet treatment proposed in [6]. It consists of a local redefinition of the topography in the intercells corresponding to dry/wet transition. Another dry/wet technique is presented in [7] considering the Godunov flux only in the intercell corresponding to the transition. The nonlinear Rieman problem whose solution depends on the tail of the jump in the bottom is solved. In [13] this technique is generalized to be used with high-order finite volume methods through state reconstructions. Another important issue is that although the scheme is well-balanced, a bad approximation of the pollutant concentration may be obtained, as it is shown in Test 1. We observe that the problem arises from the definition of the approximation of the intermediate wave speed. In fact in this work we show that the third component of the numerical flux for the WAF method can be rewritten in terms of an approximation of the intermediate wave speed. In order to have the same property for the extended version of the WAF method for non-homogeneous systems we conclude that the approximation of the intermediate wave speed must be modified. In Test 1 we can actually see that we obtain a correct profile for pollutant concentration with this modified approximation of the intermediate wave speed. After this introduction, the paper continues with Section 2, where we summarize how the WAF method is obtained for the homogeneous case. We consider two different versions corresponding to two different approaches for the intermediate flux and we rewrite them under an unique formulation. Later, in Section 3 we perform the extension of the method for SW equations with pollutant and topography. In Subsection 3.1 we summarize the technique introduced in [10] which is based in the study of the numerical diffusion introduced by the numerical scheme. In 3.2 we propose the extension of the WAF method for balance laws, in 3.3 the well-balanced properties of the proposed scheme are studied. Later, in 3.4 we present the intermediate wave speed depending on the definition of the source term. Finally in Section 4 several numerical tests are presented: a numerical study of the well-balanced properties of the method, three tests that involve the treatment of dry/wet fronts, another one to study the order of accuracy and a last one where we make a comparison with a reference solution when a periodic time-varying pollutant concentration 3 is imposed as a boundary condition. 2 The WAF method In this section we summarize the WAF method for the homogeneous SWE with pollutant. We begin by considering the homogeneous SWE given by the system      ∂W ∂t +∂F(W) ∂x = 0, x ∈[0, L], t ∈[0, T], W(x, 0) = W0x∈[0, L]; (1) where W=  h q r  , F(W) =       q q2 h+1 2gh2 qr h       , the unknowns are h, the height of the water column, q, the discharge; if we denote by ψ the pollutant concentration, or another passive scalar, then r=hψ.W0is the initial data, Lis the length of the domain, Tthe final time and gis the constant gravity. We consider a partition of the domain {xi}i={i∆x}iwhere, by simplicity, we take ∆x a constant space step, and we denote by tn=tn−1+ ∆tthe time values, with ∆tthe time step. If we use a finite volume method in conservative form to approximate the solution of this problem, we have Wn+1 i=Wn i−∆t ∆x(φn i+1/2−φn i−1/2),(2) where we denote by Wn ian approximation of the mean value of the solution on the control volume (xi−1/2, xi+1/2) at time t=tn: Wn i≈1 ∆xZxi+1/2 xi−1/2 W(x, tn)dx, and by φn i+1/2=φ(Wn i, Wn i+1) the numerical flux function that characterize each method. As we have mentioned before, WAF type methods have second order accuracy. This is obtained by averaging the solution of the associated Riemann problem in time and space. We summarize the form in which it is obtained. We consider the Riemann problem associated to (1) with initial data Wn iand Wn i+1:          ∂W ∂t +∂F(W) ∂x = 0; W(x, 0) = Wix < 0 Wi+1 x > 0 (3) 4 0 t x S S ∆t ∆t / 2 S 1 23 4 3 2 1 ω ωω ω −∆ ∆x /2 x /2 Figure 1: Computational grid and intermediate waves to compute the WAF method. We remove nsuperindex for simplycity of notation. We denote by Sifor 1 ≤i≤3 the approximation of the characteristic velocities. We consider a computational grid [x1, x2]×[t1, t2] in the x−tplane, and take t1= 0, t2= ∆t, x1=−∆x/2 and x2= ∆x/2 (see Figure 1). Then, the numerical flux of the WAF method is obtained by integrating in this volume. By using the midpoint rule for the time integral, that is, by evaluating the flux at t= ∆t/2, we get: φWAF i+1/2=1 ∆xZ∆x/2 −∆x/2 F(˜ W(x, ∆t 2))dx. Being ˜ Wthe solution of the Riemann problem. We denote by ωkthe measure of the interval kbetween the waves Sk−1and Skdivided by ∆x(see Figure 1), given by: ωk=1 2(ck−ck−1), c0=−1, cN+1 = 1,with ck=∆t ∆xSk, where Nis the number of waves. Concretely, N= 3 for SWE with pollutant. So, we can write the numerical flux as follows: φWAF i+1/2= N+1 X k=1 wkF(k) i+1/2,(4) with F(k) i+1/2the value of the flux function in the interval k. If we expand the expression (4), having into account the values of ωk, we obtain the next expression for the WAF flux: φWAF i+1/2=1 2(Fi+Fi+1)−1 2 N X k=1 ck∆F(k) i+1/2,(5) 5 where we have denoted by ∆F(k) i+1/2=F(k+1) i+1/2−F(k) i+1/2, and Fiis the flux function at point xi. The WAF scheme is a second order accuracy method in time and space, so according to Godunov’s theorem, it produces spurious oscillations near great gradients of the solution. To overcome this fact, a TVD stabilization is taken. If we denote by B(v) a flux limiter function, then the limiter function is defined by A(v, c)=1−(1 − |c|)B(v). So the TVD-WAF flux function becomes as follows: φWAF i+1/2=1 2(Fi+Fi+1)−1 2 N X k=1 Mk∆F(k) i+1/2,(6) where Mk=sign(Sk)Ak. Some suitable choices for Akcan be found in [32]. We have used in particular the Van Albada’s limiter: Ak=A(v(k), ck) =    1 if v(k)≤0 1−(1 − |ck|)v(k)(1 + v(k)) 1 + v(k)2if v(k)≥0,(7) where v(k)=            p(k) i−p(k) i−1 p(k) i+1 −p(k) i if ck>0 p(k) i+2 −p(k) i+1 p(k) i+1 −p(k) i if ck<0 , being pa scalar value. In our case we set p=h+zbfor k= 1,3 (waves S1and S3) and p=ψfor k= 2 (wave S2). the WAF method depends on the choice of an approximation for the intermediate waves. The usual choice is the HLLC flux (see [32]). This method is constructed from the HLL method. The improvement of HLLC vs. HLL is that in flux approximation it has into account the intermediate eigenvalue (corresponding to the wave S2) which is associated to the third unknown, so it gives a better approximation for this variable. The numerical flux function for HLLC scheme can be written as follows: [φHLLC i+1/2](j)= [φHLL i+1/2](j)j= 1,2; [φHLLC i+1/2](3) = [φHLL i+1/2](1)ψ∗where ψ∗=ψiif S2≥0 ψi+1 if S2<0. (8) Where the HLL flux is defined as: φHLL i+1/2=       Fiif S1≥0; FHLL =S3Fi−S1Fi+1 +S3S1(Wi+1 −Wi) S3−S1 if S1≤0≤S3; Fi+1 if 0 ≥S3. (9) 6 We can find different possibilities for S1and S3(see [32]), for instance we can set: S1=ui−pghiS3=ui+1 +pghi+1. There are also several choices for S2, we take the following value: S2=S1qi+1 −S3qi−S1S3(hi+1 −hi) qi+1 −qi−hi+1S3+hiS1 .(10) This approximation of the intermediate speed S2is proposed by Toro in [32] for HLLC solver. Remark 2.1 There are several possibilities to estimate S2. The definition of S2must be an approximation of the velocity in the usually called ‘star region’, that is, the region defined between S1and S3waves in Figure 1. The approximation of the velocity at the star region can be done using a Riemann solver. We can find different possiblilities based on several properties. For example, the exact depth positivity preservation or the presence of two-shock or two-rarefaction waves, see [32] for details. The choice of S2in equation (10) is presented in [32] for the homogeneous HLLC solver. This one is obtained by imposing that h(2) =h(3), being h(i),i= 2,3the approximations of the height in regions w2and w3respectively (see Figure 1). Following [32] we can see that this value of S2has the property of being exact when one of the states is a dry bed state, because in this case S2is the speed of the front. To compute the flux, since F(2) =F(3) for the first and second components, we can use FHLL, defined by equation (9). For the third component of the flux, following the HLLC approach, we can use FHLLψwith ψ=ψiat the left of S2and ψ=ψi+1 at the right of S2. Then, we can write the WAF flux as follows: [φWAF ](j)=Fi+Fi+1 2−1 2M1(FHLL −Fi) + M3(Fi+1 −FHLL)(j) for j= 1,2.(11) [φWAF ](3) =Fi+Fi+1 2(3) −1 2M1([FHLL](1)ψi−[Fi](3)) + M2[FHLL](1)(ψi+1 −ψi) + +M3([Fi+1](3) −[FHLL](1)ψi+1).(12) Another variant that we consider in this paper is to set F(2) and F(3) equal to the definition of φHLL instead of FHLL, obtaining a method with smaller dependence on the choice of the flux limiters. Although it is not possible to prove second order accuracy for this version 7 of the WAF method, numerically we have obtained similar results for both versions. We obtain the following numerical flux in this case [φWAF−2](j)=                  hFi+Fi+1 2−1 2 ∆t ∆xM3(Fi+1 −Fi)i(j)0≤S1; [φWAF ](j)S1≤0≤S3; hFi+Fi+1 2−1 2 ∆t ∆xM1(Fi+1 −Fi)i(j)0≥S3. for j= 1,2.(13) and for the third component: [φWAF−2](3) =                  hFi+Fi+1 2i(3) −1 2 ∆t ∆x[M3ψi+1(Fi+1 −Fi) + M1Fi(ψi+1 −ψi)](1) 0≤S1; [φWAF ](3) S1≤0≤S3; hFi+Fi+1 2i(3) −1 2 ∆t ∆x[M1ψi(Fi+1 −Fi) + M3Fi+1(ψi+1 −ψi)](1) 0≥S3. (14) 2.1 Rewriting the WAF method We have presented two different versions of the numerical flux for the WAF method, corresponding to definitions (11)-(12) and (13)-(14). In this section we rewrite both versions of the WAF method under a compact formulation by introducing a generalization of the definitions of Si,i= 1,2,3. From now on, with the purpose of clarity, we denote SL=S1, S∗=S2, SR=S3,and analogously ML=M1, M∗=M2MR=M3. (15) We introduce the values SL,S∗,SR. For the first version (11)-(12) these values coincide with their original values SL=SL, S∗=S∗, SR=SR. We can also see that the flux defined by (13)-(14) only varies in this definition. In this case we define the extended version of these values as follows: SL= min{SL,0};S∗=   SL0≤SL; S∗SL≤0≤SR; SRSR≤0. SR= max{SR,0}.(16) These definitions allow us to write both versions of the WAF methods under the same compact expression φWAF (j)=Fi+Fi+1 2−1 2 1 SR−SL(MRSR−MLSL)(Fi+1 −Fi)− −SRSL(MR−ML)(Wi+1 −Wi)(j); for j= 1,2.(17) 8 φWAF (3) =Fi+Fi+1 2(3) −1 2 1 SR−SL(SRMRψi+1 −SLMLψi)(Fi+1 −Fi) + +SLSR(MLψi−MRψi+1)(Wi+1 −Wi) + +M∗(ψi+1 −ψi)SRFi−SLFi+1 +SLSR(Wi+1 −Wi)(1). (18) In the following section we shall use this formulation to perform the extension to nonhomogeneous case. 3 Extension to non-homogeneous system We consider now the non-homogeneous shallow-water equations with pollutant:      ∂W ∂t +∂F(W) ∂x =G(x, W), x ∈[0, L], t ∈[0, T], W(x, 0) = W0x∈[0, L]. (19) If we take the topography source term and we denote by zb(x) the height of the topography at point x, we have: G(x, W) =   0 −ghz0 b(x) 0  . For the treatment of the non-homogeneous system, the usual technique is to approximate the numerical flux function by a numerical flux for the homogeneous system and to consider an upwind discretization of the source term (see [2]). Nevertheless, in order to obtain the structure of the numerical scheme in conservative form, we must integrate equation (19) on the control volume (xi−1/2, xi+1/2). Thus, we obtain Wn i+1 −Wn i ∆t+Fi+1/2−Fi−1/2 ∆x=1 ∆xZxi+1/2 xi−1/2 G(x, W(x))dx. Where, if we consider a Godunov solver, Fi+1/2denotes an approximation of F(˜ W(xi+1/2)) with ˜ Wthe solution of the Riemann problem associated to the non-homogeneous system. The key point is that the Riemann problem depends on the original system, so it depends on the source term. Then, Fi+1/2depends on the definition of the source term. So we use the following structure for the scheme to approximate (19): Wn i+1 −Wn i ∆t+φn G,i+1/2−φn G,i−1/2 ∆x=Gn C i.(20) 9 (a) 01234 1 0.5 0 0.5 1 1.5 2 Velocity S S *G * (b) 0 0.5 1 1.5 2 2.5 3 3.5 4 0 2 4 6 8 10 12 14 16 Conc. WAF− S*G Conc. WAF− S* Bottom Figure 2: (Test 1). (a) Velocity approximations, S∗and S∗G. (b) Pollutant concentration for the WAF method taking S∗or S∗G. a great error near the bump. So, from now on we will only consider the approximation of the intermediate wave speed given by S∗G. In Figure 3 we present the bottom function, water surface and pollutant concentration provided by the WAF solver and HLLC (see [11]). We observe a less numerical diffusion of the WAF method. In Table 1, we make a comparison of the runtime (seconds, in an old Pentium II machine) of two methods, WAF and HLLC, in relation to the error for the pollutant concentration in L1norm. For this aim we have solved the problem with meshes ∆x= 0.08,0.04,0.02,0.01 and we compare it with a reference solution corresponding to the WAF method with 1000 points. We observe that the runtime is very similar for both methods. However to get the same order of accuracy, for instance an error of 0.02, the the WAF method spends 1.3 seconds whereas HLLC spends 74.29 seconds. So, in this sense, the WAF method is again better than HLLC. ∆x L1err WAF Runtime WAF L1err HLLC Runtime HLLC 0.08 0.0225 1.3 0.0834 1.21 0.04 0.011 4.98 0.0584 4.87 0.02 0.0047 18.51 0.0406 18.72 0.01 0.0014 74.5 0.0256 74.29 Table 1: (Test 1). Comparison of runtime and error. Test 2: Well-balanced tests We consider two different well-balanced tests. In the first one (Test 2.a) we study the 16 (a) 0 0.5 1 1.5 2 2.5 3 3.5 4 0 2 4 6 8 10 12 14 16 18 Water surface Concentration (HLLC) Concentration (WAF) Bottom (b) 0 0.5 1 1.5 2 2.5 3 3.5 4 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Concentration (HLLC) Concentration (WAF) Figure 3: (Test 1). (a) Bottom, water surface and pollutant concentration for HLLC and the WAF method using S∗G. (b) Pollutant concentrations. exact well-balanced property of the method for water at rest. In the second one (Test 2.b) we study the second order well-balanced accuracy of the method for a subcritical stationary solution. Test 2.a: Exact well-balanced property We consider a test of water at rest using a random bottom (See Figure 4 (a)). The case of water at rest is the only one of a stationary solution for which the pollutant concentration is also stationary for arbitrary profiles (otherwise, if q6= 0, ψmust be constant to obtain a stationary solution). Then, we also consider a random profile of ψ(See Figure 4 (b)). We set a domain of L= 1 m. discretized with 100 points and a CFL condition equal to 0.8. In Table 2 we present the error in L1norm. The solution is preserved up to the machine accuracy. h q r,ϕ=r/h L1error 3.3827 ×10−16 1.1725 ×10−15 1.7493 ×10−15 1.1328 ×10−15 Table 2: (Test 2.a). Errors. Exact well-balanced property. Test 2.b: Approximated well-balanced property In this test we verify numerically the approximated well-balanced property of the WAF method (Theorem 3.1 i)). We consider a subcritical stationary solution. The length of the domain is 20 m. The discharge is constant, q= 4.42, and the bottom function zb(x) = 0.2e−0.16(x−10)2, x ∈[0,20]. 17 (a) 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 −2 −1.8 −1.6 −1.4 −1.2 −1 −0.8 −0.6 −0.4 −0.2 0 Water surface Bottom (b) 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Concentration Figure 4: (Test 2.a). (a) Random bottom, (b) Random concentration of pollutant The analytical solution for h(x) can be obtained as follows. We have that the equation for the velocity u=q/h is ∂t(u) + ∂x(u2 2+g(h+zb)) = 0. Then, for a regular stationary solution, u2/2 + g(h+zb) = k, with ka constant value. If we rewrite this equation in terms of hand q, we obtain that his a root of P(h) = h3+ (zb−k g)h2+q2 2g. To define k, we impose at x= 0 that h= 2 −zb(0). The subcritical solution correspond with the bigger root of P(h). In Table 3 we observe that effectively we achived the second order accuracy. Nodes L1err h L1order h L1err q L1order q 20 1.58 ×10−3- 5.02 ×10−3- 40 5.07 ×10−41.646 1.21 ×10−32.0513 80 1.3×10−41.967 3.04 ×10−41.9965 160 3.2×10−51.995 7.6×10−51.9985 320 9 ×10−61.904 1.9×10−51.9987 Table 3: (Test 2.b). Errors and order. Subcritical stationary solution. Test 3: Order of accuracy We consider a test proposed in [40] to study the order of accuracy of the method. Although we also introduce a non constant profile of the concentration of pollutant. We compute up to t= 0.1 where any shock still appears in any of the unknowns. The domain is [0,1], we impose periodic boundary conditions and initial conditions h(x, 0) = 5 + ecos(2πx), q(x, 0) = sin(cos(2πx)), ψ(x, 0) = 1 2(cos(4πx) + 1). 18 The bottom functions is defined by zb(x) = sin2(πx). As in this case we do not know an analytical solution for the problem, we compute a reference solution with the WAF method using 12800 cells to compute the errors. In all cases we fix the CFL condition to 0.8. In Table 4 we present the results for the errors in norm L1and the corresponding order. The second order of accuracy is achieved. Nodes L1err h L1order h L1err q L1order q L1err r L1order r 25 2.802 ×10−2- 3.143 ×10−1- 2.112 ×10−2- 50 1.021 ×10−21.456 9.702 ×10−21.695 7.464 ×10−31.501 100 3.228 ×10−31.661 2.677 ×10−21.857 2.533 ×10−31.559 200 9.15 ×10−41.819 6.594 ×10−32.021 7.51 ×10−41.753 400 2.53 ×10−41.851 1.553 ×10−32.085 2.07 ×10−41.858 Table 4: (Test 3). Errors and order. Test 4: Oscillating lake In this test we study an oscillating lake over a non flat bottom. This test has been proposed by Audusse et al in [1]. A shore line variable in time appears in the problem. The initial conditions are h(x, 0) = max(0,0.4−zb(x)+0.04 sin((x−0.5)/0.25) + 0.04 max(0., −0.4 + zb(x))), q(x, 0) = 0, zb(x) = 1 2(1 −1 2(cos(π(x−0.5) 0.5) + 1)). The length of the domain is L= 1 m., we consider 200 points and a CFL condition equal to 0.8. In Figure 5 the water surface is plotted at the initial condition and at t= 19.87 s., corresponding to a period where the wave is located at the opposite side of the initial condition. Test 5: Small perturbation of a steady-state solution This test has been proposed by Kurganov and Petrova in [17]. It consists of studying a small perturbation in the water surface for a steady-state solution, where the height of the water column is zero in some isolated points. The domain is [0,1], we compare the approximations obtained with 200 points and 3200 points. The CFL condition is set to 0.8. The initial conditions are h(x, 0) = 1 + ε−zb(x) 0.1<x<0.2 1−zb(x) otherwise,q(x, 0) = 0, 19 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 t=0 t=19.87 bottom Figure 5: (Test 4). Oscillating lake. Water surface at t= 0 and t= 19.87. where ε= 0.001 is the small perturbation on the water surface located between 0.1 and 0.2. The bottom function is defined by zb(x) =        10(x−0.3),0.3≤x≤0.4, 1−0.0025 sin2(25π(x−0.4)),0.4≤x≤0.6, −10(x−0.7),0.6≤x≤0.7, 0 otherwise. In Figure 6 we present the initial conditions. In Figure 7 we present a zoom on [0,1] × [0.995,1.005] of the solution for t= 1 s. computed with 200 points. In Figure 8 we make a comparison with the solution for 3200 points. In Figure 8 (a) we present a zoom on the first point of jump of the bottom found by the wave. We have a vertical wall at x= 0.4, but in discret we have a plane with slope 1/∆x. That explains the difference in this area between both values of ∆x(∆x= 1/200 and ∆x= 1/3200). Otherwise, in Figure 8 (b) we present a zoom on [0.25,0.75] ×[0.9994,1.0005] of the water surface. We observe that the results obtained with 200 and 3200 points are not so different in the intermediate areas between the two vertical walls. And as we have previously explained, it is greater close to the vertical walls at x= 0.4 and x= 0.6. Test 6: Dry zone by a double rarefaction wave This test was proposed by Gallou¨et in [14]. It consists of a double rarefaction wave (see [32]) including an obstacle at the bottom. The numerical treatment of dry/wet areas has been done following the technique described in [6]. The length of the domain is 20 meters and the space step ∆x= 0.1. We set the CFL condition equal to 0.8. The initial conditions are h0= 10 −zb;q0=350 x > 50/3 −350 otherwise ;r0=h025/3<x<12.5 0otherwise ; 20 (a) 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 0.2 0.4 0.6 0.8 1 Water surface (t=0) bottom (b) 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0.995 0.996 0.997 0.998 0.999 1 1.001 1.002 1.003 1.004 1.005 Water surface (t=0) bottom Figure 6: (Test 5). Initial condition. Water surface and bottom. (b) zoom on [0,1] × [0.995,1005] 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0.995 0.996 0.997 0.998 0.999 1 1.001 1.002 1.003 1.004 1.005 Water surface (t=1) bottom Figure 7: (Test 5). Solution for t= 1 s. Zoom on [0,1] ×[0.995,1005] 21 (a) 0.39 0.4 0.41 0.42 0.43 0.44 0.45 0.9975 0.998 0.9985 0.999 0.9995 1 1.0005 200 points 3200 points bottom (b) 0.25 0.3 0.35 0.4 0.45 0.5 0.55 0.6 0.65 0.7 0.75 0.9994 0.9995 0.9996 0.9997 0.9998 0.9999 1 1.0001 1.0002 1.0003 1.0004 1.0005 200 points 3200 points Figure 8: (Test 5). Line with circles: 200 points. Continuous line: 3200 points. (a) Zoom on [0.39,0.45] ×[0.9975,1.0005], Water surface and bottom. (b) Zoom on [0,25,0,75] × [0.9994,1.0005], Water surface with zb(x) = 1 25/3<x<12.5 0otherwise . In Figure 9 we compare the water surface obtained by WAF and HLLC with a reference solution computed with the WAF method with a mesh of 3000 points. We observe how the WAF method better captures the shock at x= 25/3, and the bump on the water surface around x= 6.5. In Figure 9 (b) we show the solution for the pollutant concentration at time T= 0.18, again we see that the WAF method introduces a less numerical diffusion. Test 7: Stationary transcritical flux with a shock and a periodic time-varying pollutant concentration We consider a classical test for a stationary solution for hand qbut including in this case a periodic time pollutant source as a boundary condition. The domain length is 20 meters and the bottom function is defined by a bump at the middle of the domain: zb(x) = 0.2−0.05(x−10)28<x<12 0otherwise . We consider a space step ∆x= 0.1 and a CFL condition equal to 0.9. The initial conditions are: h0= 0.33 −zb;q0= 0; r0= 0, and the boundary conditions q= 0.18 and r=rp(t) at x=0; h= 0.33 at x= 20. The function rpis a periodical function in time that represents a contribution of pollutant of three seconds every ten seconds; and it is given by: rp(t) = h t ∈[10i, 10i+ 3] 0otherwise ,for i∈N. 22 (a) 0 5 10 15 20 25 0 2 4 6 8 10 HLLC WAF Reference solution Bottom (b) 0 5 10 15 20 25 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 HLLC WAF Reference solution Figure 9: (Test 6). (a) Water surface at T= 0.25. (b) Pollutant concentration at T= 0.18. For this test the analytical solution for hand qcan be obtained. qis constant and equals to the boundary condition, so q= 0.18. To obtain h, it is necessary to solve the third order polynomial P(h) = h3+ (zb−k g)h2+q2 2g, where kis the constant energy value and gis the gravity. This solution has a shock that connects two different constant values of the energy. The position of the shock is obtained by the Rankine-Hugoniot condition. With this data we can calculate the velocity of the fluid u=q/h. So we can predict the position of the pollutant fronts by approximating the trajectory of the points, by solving the characteristic curves problem ˙x(t) = u(x(t), t); x(0) = 0; (38) To approximate x(t) we have used a Runge-Kutta method of third order. Once we have x(t), we can build the profile of the pollutant concentration, that is a passive scalar transported with the velocity of the fluid. In Figure 10, we compare the exact solution for hand qwith the numerical results obtained with the proposed the WAF method. In Figure 10 (a) we compare the water surface. This is a transcritical problem with a critical point at x= 10. We observe that the method does not produce the phenomenon known as ‘dog leg’ as presented in Roe method. In Figure 10 (b) the velocity of the fluid is compared. Usually the discharge value is presented for this test but we are interested in the velocity, which is the responsible for the transport of the pollutant concentration. We see that the approximated velocity suits the analytical solution, even close to the shock. 23 (a) 0 2 4 6 8 10 12 14 16 18 20 0 0.2 0.4 Exact sol. WAF Bottom (b) 0 2 4 6 8 10 12 14 16 18 20 0 0.5 1 1.5 2 Exact sol. WAF Bottom Figure 10: (Test 7). (a) Water surface and topography. (b) Velocity. To compare the pollutant concentration obtained by the WAF and HLLC solver with the exact solution, we begin at t= 189, a time that ensures the solution for hand qto be stationary. In Figure 11 we compare the pollutant concentration for the WAF and HLLC solver for different times values. We observe the great numerical diffusion introduced by the HLLC solver in comparison with the WAF method. In Figure 12 we compare the WAF method with the approximated solution computing solving (38). The considered times show the transition of a column of pollutant through the bump where the velocity of the fluid grows. One of the characteristics of the solution for the pollutant concentration is that the width of the rectangles of pollutant that is introduced by the right boundary is not the same in all the domain. That is because the velocity is not constant in all the domain, then the velocity of the particle at the front of the rectangle is different to those at the back. If we look at Figure 10 (b) where the velocity is drawn, we observe two regions of constant velocity and a fast increase of the velocity over the bump until the shock. Then, according to this profile of velocity we can see how the rectangle of pollutant expands in this zone in Figures 11 and 12. In Figure 12 we observe that the biggest difference between the numerical and the exact solution just corresponds to time t= 195, over the area of bigger velocity previous to the shock (check the velocity profile of Figure 10 (b)). This difference at the back of the rectangle after the bump is motivated by the numerical diffusion of the scheme because some points have not still pass through the shock. We observe that after the shock, in the zone of constant velocity, this difference between exact and numerical solution is reduced. 24 0 2 4 6 8 10 12 14 16 18 20 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 t = 189s HLLC WAF Bottom 0 2 4 6 8 10 12 14 16 18 20 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 t = 191s 0 2 4 6 8 10 12 14 16 18 20 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 t = 193s 0 2 4 6 8 10 12 14 16 18 20 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 t = 194s 0 2 4 6 8 10 12 14 16 18 20 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 t = 195s 0 2 4 6 8 10 12 14 16 18 20 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 t = 197s Figure 11: (Test 7). Concentration for WAF and HLLC methods. 25