On well-balanced finite volume methods for non-conservative non-homogeneous hyperbolic systems

On well-balanced Finite Volume Methods for non-conservative non-homogeneous hyperbolic systems. ∗ Tom´as Chac´on Rebollo † , Enrique D. Fern´andez-Nieto ‡ Manuel J. Castro D´ıaz §& Carlos Par´es February 2, 2006 Abstract In this work we introduce a general family of finite volume methods for non-homogeneous hyperbolic systems with non-conservative terms. We prove that all of them are “asymptotically well-balanced”: They preserve all smooth stationary solutions in all the domain but a set whose measure tends to zero as ∆xtends to zero. This theory is applied to solve the bilayer Shallow-Water equations with arbitrary cross-section. Finally, some numerical tests are presented for simplified but meaningful geometries, comparing the computed solution with approximated asymptotic analytical solutions. Short title : Well-balanced Finite volume solvers. Keywords : Well-balanced, Finite Volume Method, upwinding, shallow water, source terms, two-layer flows. Subject Classifications : AMS (MOS) : 65N06, 76B15, 76M20, 76N99. 1 Introduction. This work deals with finite volume solvers for non-homogeneous hyperbolic systems with non-conservative terms. Our goal is to design numerical schemes that calculate with high accuracy the stationary solutions of the system. The accurate computation of stationary solutions of hyperbolic systems with source terms has been found in the past years as closely related to the accurate computation of transient solutions. Numerical schemes that do not solve steady solutions up to second order, at least, yield transient solutions that present large errors, unacceptable from the physical point of view. To avoid this difficulty, in [1] Berm´udez and V´azquez introduced the so-called “C-property” in the context of Shallow-Water equations. They ask the numerical schemes to exactly calculate at grid nodes the stationary solution corresponding to water at rest. They also present some numerical evidence of the large errors presented by schemes that do not verify this property. In that work, the scheme of Roe is generalized to non homogeneous hyperbolic systems, in such a way that it verifies the C-property. This property has given fruitful extensions to more general situations: numerical models, hyperbolic equations, and extensions of the condition itself. ∗This research was partially supported by Spanish Government Research Projects BFM2003-07530-C02-01 and BFM2003-07530-C02-02. †Departamento de Ecuaciones Diferenciales y An´alisis Num´erico, Universidad de Sevilla. C/ Tarfia, s/n. 41080 Sevilla, Spain ([email protected]). ‡Departamento de Matem´atica Aplicada I, Universidad de Sevilla. 41012 Sevilla, Spain ([email protected]). §Dpto. An´alisis Matem´atico, Universidad de M´alaga, Campus Teatinos s/n, M´alaga, Spain ([email protected]),([email protected]). 1 For instance, an extension to kinetic schemes was introduced by Perthame and Simeoni in [16]. This type of schemes is based upon the discretization of an equation for the density of particles. Also, Kurganov and Levy in [14] introduced the central upwind schemes for Shallow-Water equations, characterized by only including scalar numerical diffusion. No upwinding matrix is needed. Both schemes ensure the positiveness of the water height, and verify the C-property. The relevance of this property as a consistency condition was pointed out in the work of Perthame and Simeoni [17]. They prove the convergence, in the sense of distributions, of the numerical solution computed by a finite volume scheme for a non-homogeneous equation under some hypothesis of consistency of the numerical scheme. This hypothesis is fulfilled if the numerical scheme preserves certain equilibrium states of the system. An extension of the C-condition to a more general condition was introduced in Greenberg and Leroux [13]: the concept of “well-balanced”scheme, as a scheme that preserves all equilibria of the system at grid nodes. In this work, a numerical scheme for non-homogeneous scalar equations satisfying this condition is introduced. An extension to a class of non-conservative hyperbolic system of equations is introduced in [2] by Castro, Mac´ıas and Par´es by means of a generalized Roe’s scheme. The systems considered include a non-conservative product with the structure B(W)∂W ∂x where Bis a matrix and Wis the unknown. The numerical stability of this scheme is achieved by using upwinding matrices defined not only from the Jacobian of the flux, but also from the matrix Binvolved by the non-conservative terms. The scheme introduced is applied to the 1d two-layer shallow water system for channels with constant width. The scheme exactly computes water at rest. More recently, in [4] this Roe’s scheme has been extended to the more general case of channels with irregular geometry and not necessarily rectangular cross-section. This extended numerical scheme was applied to the simulation of the Strait of Gibraltar. The derivation of systematic techniques to build numerical schemes satisfying the C-property is also a relevant issue, as this is far from being a straightforward condition to obtain. In [5] Chac´on, Dom´ınguez and Fern´andez-Nieto introduced the idea of separately compensate the centered and decentered components of the numerical flux by centered and decentered discretizations of the source terms. In [6] the same authors introduced the idea of adding an additional source term to the parabolic equivalent equation. This term is specifically designed to compensate for the first-order error generated by the numerical diffusion term on steady solutions. Then, centered discretizations automatically provide second-order accurate approximations of all steady solutions. In this work, we step on this process by considering a general non-homogeneous hyperbolic system of equations with non-conservative terms, which also depends upon the space variable. We derive a general class of numerical schemes that we prove to be “asymptotically well-balanced” in the following sense: they calculate, with at least order two, all stationary solutions of the system in all the domain but on a set whose measure tends to zero as ∆xtends to zero. One of the free parameters of this class is the diffusion matrix, so we include both flux-splitting and flux-difference or Roe’s methods as particular cases. We also give general sufficient conditions under which the schemes exactly calculate a given stationary solution. These conditions may be read as meaning on one hand that the numerical flux of the scheme must compensate the centered part of the numerical source term, and on another hand that the numerical diffusion must compensate the decentered part of the numerical source term. This paper is organized as follows: In Section 2, a family of numerical schemes for a general non-homogeneous hyperbolic systems of equations with non-conservative terms is introduced. In Section 3, we introduce the concept of “asymptotically well-balanced” scheme and we prove that the introduced numerical schemes satisfy this property. Section 5 is devoted to the construction of numerical diffusion matrices that ensure the stability of our schemes for linear systems, while ensuring the asymptotically well-balanced computation of all steady solutions. In Section 5, we apply the theory developed to the bilayer Shallow-Water equations with arbitrary cross-sections to construct a family of asymptotically well-balanced numerical schemes. The Roe-type scheme introduced in [4] is a particular case. In Section 6, we present some numerical tests for simplified geometries, where exact 2 solutions can be determined. The purpose of these experiments is to compare the properties of the here introduced numerical schemes for the two-layer shallow water system. More precisely, we focus on the well-balance property and the quality of the numerical approximations for smooth transitions and shocks. We present our conclusions in Section 7. In particular we conclude that the use of pointwise adapted discretizations of the centered part of the numerical flux helps to construct accurate solvers for non-conservative equations. 2 Non-homogeneous non-conservative hyperbolic systems. We shall consider the following system:          ∂W (x, t) ∂t +∂ ∂x[F(x, W(x, t))] = G(x, W(x, t)) + B(x, W (x, t))∂W(x, t) ∂x , x∈(0, L), t ∈(0, T); W(x, 0) = W0(x)x∈(0, L); (1) where by F: [0, L]×RN→RNwe denote the physical flux function, by G: [0, L]×RN→RNthe source term, and B: [0, L]×RN→RN×RNis a regular matrix function of xand W. We will suppose that F,Gand Bare C2functions. By Lwe denote the length of the domain and by Tthe final time. We will complete (1) with appropriate boundary conditions, that we set in each problem. System (1) includes the non-conservative product B(x, W)∂W ∂x , which is not a distribution if Wis discontinuous. Its presence requires an adaptation of the usual concept of weak solutions: Across a discontinuity of W, the non conservative product can be defined along paths connecting the two states, following DalMasso-LeFloch & Murat [8]. System (1) also includes a dependency upon the space variable of the flux F. Both characteristics appear in two layer Shallow Water equations in channels of variable breadth; as we shall see later on. However, these characteristics also appear in other flow models such as boiling flows and two-phase flows among other (See Fowler [11]). We suppose that the system (1) is hyperbolic, in the sense that ∀x∈[0, L] and ∀W∈Rdthe matrix M(x, W) = A(x, W)−B(x, W) (2) can be diagonalized and all its eigenvalues are real and different, where A(x, W) = ∂F ∂W (x, W ) is the Jacobian matrix of the flux function F. We consider a partition {xi}M+1 i=0 of the interval [0, L], and a partition {tn}m+1 n=0 of [0, T]. By Wn i we define an approximation to the average of the solution in the control volume, (xi−1/2, xi+1/2) at t=tn. In this work, for the space discretization, we develop finite volume schemes that can be written in viscous form. For the homogeneous conservative system, these schemes are second order approximations in space, of a parabolic equivalent system: ∂W ∂t +∂ ∂xF(x, W)−ν∂ ∂x(D(x, W)∂ ∂x W) = 0, where νis equal to the half of the space step (ν= ∆x/2) and D(x, W ) is the numerical diffusion matrix: The term −ν∂ ∂x(D(x, W)∂W ∂x ),(3) 3 represents the numerical diffusion introduced by the scheme, stemming from the upwinding of the convection terms. The diffusion matrix D(x, W ) must take into account that the effective transport term in system (1) is M(x, W )∂W ∂x , where Mis the matrix defined by (2). If a centered discretization of the remaining terms in the non-homogeneous and non-conservative system (1) is performed, then the equivalent parabolic system is ∂W ∂t +∂ ∂xF(x, W)−ν∂ ∂x(D(x, W)∂ ∂x W) = G(x, W ) + B(x, W )∂W ∂x .(4) Notice that, if ¯ Wis a stationary solution of (1), then it verifies ∂ ∂x[F(x, ¯ W(x))] = G(x, ¯ W(x)) + B(x, ¯ W(x))∂¯ W(x) ∂x .(5) So, it is clear that ¯ Wis not a stationary solution of (4). There is a first order error due to the numerical diffusion term (3) that is not compensated. In order to obtain numerical schemes that calculate up to second order of accuracy all smooth stationary solutions, we propose to compensate the numerical diffusion term by an additional source term, that we denote by C(x, W ) (Correction term), that is, the numerical scheme must be a second order approximation of ∂W ∂t +∂ ∂xF(x, W)−ν∂ ∂xD(x, W)∂W ∂x =G(x, W ) + B(x, W)∂W ∂x +C(x, W ). In order to define the term C(x, W), we consider a smooth stationary solution ¯ Wof (1). We want the identity −ν∂ ∂xD(x, ¯ W)∂¯ W ∂x =C(x, ¯ W) (6) to hold. As ¯ Wverifies (5), then ∂F ∂x (x, ¯ W(x)) + A(x, ¯ W(x))∂¯ W ∂x (x) = ∂ ∂x[F(x, ¯ W(x))] = G(x, ¯ W(x)) + B(x, ¯ W(x))∂¯ W ∂x . Thus, M(x, ¯ W(x))∂¯ W ∂x (x) = G(x, ¯ W(x)) −Fx(x, ¯ W(x)),(7) where Fx(x, W) = ∂F ∂x (x, W ). Moreover, if M(x, ¯ W) is non-singular, we have ∂¯ W(x) ∂x =M−1(x, ¯ W(x))S(x, ¯ W(x)), where S(x, W) = G(x, W)−Fx(x, W). Therefore, −ν∂ ∂xD(x, ¯ W(x))∂¯ W ∂x (x)=−ν∂ ∂xD(x, ¯ W(x))M−1(x, ¯ W(x))S(x, ¯ W(x)).(8) Comparing (8) with (6), we propose to define C(x, W) = −ν∂ ∂xD(x, W )M−1(x, W)S(x, W);∀x∈[0, L], 4 ∀W∈RNsuch that M(x, W ) is non-singular. The equivalent system for the hyperbolic system with source term and non-conservative terms is the following, ∂W ∂t +∂ ∂xF(x, W)−ν∂ ∂xD(x, W)∂W ∂x = =G(x, W) + B(x, W )∂W ∂x −ν∂ ∂xD(x, W )M−1(x, W)S(x, W).(9) Every stationary solution of (1) is also a stationary solution of this system at points where M(x, W (x)) is non-singular. We conclude then that asymptotically well-balanced numerical schemes for the hyperbolic system (1) can be built by the use of approximations of second order in space of (9). Notice that system (9) has a conservative structure when the original system (1) is conservative (i.e., B= 0). Indeed, it can be re-written as ∂W ∂t +∂ ∂xF(x, W)−νD(x, W )∂W ∂x −M−1S(x, W )=G(x, W ) + B(x, W )∂W ∂x .(10) Based upon this structure, we propose general three-points schemes as follows Wn+1 i−Wn i ∆t+φS R(xi, xi+1, W n i, Wn i+1)−φS L(xi−1, xi, Wn i−1, Wn i) ∆x= =GC(xi−1, xi, xi+1, W n i−1, Wn i, Wn i+1),(11) where φS R(xi, xi+1, Wi, Wi+1) = φS(xi, xi+1, Wi, Wi+1)−1 2BR(xi, xi+1, Wi, Wi+1) (Wi+1 −Wi),(12) φS L(xi, xi+1, Wi, Wi+1) = φS(xi, xi+1, Wi, Wi+1) + 1 2BL(xi, xi+1, Wi, Wi+1)(Wi+1 −Wi) (13) and φS(xi, xi+1, Wi, Wi+1) = FC(xi, xi+1, Wi, Wi+1)− −νD(xi, xi+1, Wi, Wi+1)Wi+1 −Wi ∆x− ] M−1(xi, xi+1, Wi, Wi+1)SD(xi, xi+1, Wi, Wi+1).(14) Here, BLand BRare approximation of Bat x=xi+1/2, respectively. Moreover, GCis an approximation of G(x, W) in x=xi. By FC,D, ] M−1and SDwe denote the approximations of F,D,M−1and Sin x=xi+1/2, respectively. Let us suppose for instance that the singularities of Mare concentrated on a continuous hypersurface γ⊂Rn+1, where, say the j-th eigenvalue λjvanishes. In this case, we define ] M−1(x, y, U, V ) as ] M−1(x, y, U, V ) = X((x+y)/2,f W)g Λ−1(x, y, U, V )X−1((x+y)/2,f W),(15) where f W=f W(x, y, U, V ) is an intermediate state between Uand V, by Xwe denote the matrix defined by the eigenvectors of M((x+y)/2,f W), and g Λ−1(x, y, U, V ) = Diag(g λ−1 i(x, y, U, V ), i = 1,...,N) with g λ−1 i=λ−1 i,if i6=j, and g λ−1 j=λ−1 jif L[(x, U),(y, V )] ⊂RN+1 \γ, 0 otherwise; (16) 5 where the λjare the eigenvalues of M((x+y)/2,f W) and L[(x, U),(y, V )] is the closed segment in RN+1 that connects (x, U) and (y, V ). It is straightforward to extend this definition to a more general case. This definition of ] M−1implies that, when an eigenvalue vanishes, no upwinding of the source term is performed at the direction defined by the corresponding eigenvector. For non-homogeneous scalar hyperbolic equations we may prove that this leads to stable schemes. This definition of ] M−1allows the schemes to be asymptotically well-balanced even for a class of stationary solutions for which matrix Mis singular and, for convenient choices of the diffusion matrix D, for every stationary solutions. We shall prove this in the next section. Remark 1 In [7] it is proved that for homogeneous conservative systems like (10) both flux-difference and flux-splitting schemes fit into the general structure of fluxes given by (12) and (13) with convenient choices of the centered numerical flux FCand the numerical diffusion matrix D. We point out that for non-conservative systems, it makes no sense to consider flux-splitting schemes, as these are intrinsically based upon the conservative formulation. Indeed, these schemes use a splitting of the flux into (we assume B= 0) F(W) = A∗(W)W. for some matrix A∗. However, in general we may use flux-difference schemes, as these are based upon the use of the velocities and speeds of propagation of the total flux, given by the matrix M. In addition, other kind of schemes fit into the general framework set by the general structure (11). We shall introduce some of them in Section 5. 3 Analysis of the “asymptotically well-balanced” property The purpose of this section is to analyze the well-balanced property of the numerical schemes of the form (10)-(14). This analysis is performed by considering an arbitrary stationary solution and using Taylor expansions. As consequence, the results are valid for every stationary solution (and not only for particular cases as ’water at rest’ solutions for the shallow-water equations) but we shall only be able to analyze the well-balanced property in terms of orders of accuracy, and not exactly at grid nodes. Due to the lack of regularity of M−1(and also to the low regularity of D, usually only Lipschitz continuous) we cannot hope to balance steady solutions of system (1) up to second order of accuracy in the whole domain (0, L). However, we shall prove that this property holds on a subset of (0, L) except a subset whose measure tends to zero as ∆xtends to zero for a large class of stationary solutions. We will say in this case that the numerical scheme is “asymptotically well-balanced”. Based upon these observations, we introduce the following definitions: Definition 1(Second order well-balanced schemes) We say that the numerical scheme (11)-(14) is well-balanced for a regular stationary solution W(x)of the hyperbolic system (1) if the corresponding consistency error is second order in all nodes xiinside the domain (0, L). The following definition was introduced in [6]: Definition 2(Asymptotically well-balanced schemes) We say that the scheme (11)-(14) is asymptotically well-balanced for a regular stationary solution W(x)of the hyperbolic system (1) if there exists an increasing sequence of compacts {Kn}nof [0, L] such that 1) µ((0, L)\∪nKn) = 0, where by µwe denote the Lebesgue measure in R. 2) For all nthere exists a δn>0such that if 0<∆x < δn, then the scheme balances system (1) up to second order in Kn. 6 To perform our analysis, we will assume some hypothesis on the different elements defining the numerical scheme (11)-(14). Although these are somewhat cumbersome, this will allow to extend the analysis to a rather large class of schemes: Hypothesis 1 a)There exists a C12-tensor F1: [0, L]×RN→ L(Rn+1,Rn)such that ∀U, V ∈RN, ∀x, y ∈[0, L] FC(x, y, U, V ) = Fx+y 2,U+V 2+ +F1x+y 2,U+V 2 y−x V−U+o(|x−y|2+|U−V|2). b)There exists two C13-tensors B1,L and B1,R : [0, L]×RN→ L(Rn+1,L(Rn,Rn)) such that ∀U, V ∈RN,∀x, y ∈[0, L] BL/R(x, y, U, V ) = Bx+y 2,U+V 2+ +B1,L/R x+y 2,U+V 2 y−x V−U+o(|x−y|2+|U−V|2). c)∀x, y, z, ∈[0, L],∀U, V, W ∈RN, GC(x, y, z, U, V, W) = G(y, V ) + O(|y−x|2+|z−y|2+|V−U|2+|W−V|2). d)There exists a C12-tensor S1: [0, L]×RN→ L(Rn+1,Rn)such that ∀x, y ∈[0, L],∀U, V, ∈RN, SD(x, y, U, V ) = Sx+y 2,U+V 2+ +S1x+y 2,U+V 2 y−x V−U+O(|x−y|2+|V−U|2). e)There exists a continuous hypersurface γ⊂RN+1 such that Mis non-singular and C2in RN+1 \γ. Moreover, there exists a 3-tensor M1:RN+1 \γ→ L(Rn+1,L(Rn,Rn)) which is C1in RN+1 \γ, that verifies ] M−1(x, y, U, V ) = M−1x+y 2,U+V 2+ +M1x+y 2,U+V 2 y−x V−U+O(|x−y|2+|V−U|2). for close enough points (x, U),(y, V )such that the closed segment L[(x, U),(y, V )] ⊂RN+1 \γ. f)There exists a continuous hypersurface eγ⊂RN+1 and a 3-tensor D1:RN+1\eγ→ L(Rn+1,L(Rn,Rn)) which is C1in RN+1 \eγsuch that for all points (x, U),(y, V )with L[(x, U),(y, V )] ⊂RN+1 \eγ, D(x, y, U, V ) = Dx+y 2,U+V 2+ +D1x+y 2,U+V 2 y−x V−U+O(|y−x|2+|V−U|2). where, given two Banach spaces E,F,L(E, F)denotes the vector space of continuous linear maps. Hypothesis a), b) and c) state that FC, (BL+BR)/2∆xand GCrespectively are second order approximations of F,B∂W/∂x and G. Hypothesis e) formalizes the fact that M−1presents singularities at points where one of its eigenvalues vanishes. Hypothesis f) takes into account the fact that usually matrix Dwill only present Lipschitz regularity when one of its eigenvalues vanishes. A typical choice is, for example, D=|M|. To state our main result, we need another technical definition: 7 Definition 3Given a continuous function W: [0, L]→RNwe call “sets of regular and singular points” of Wwith respect to the scheme (11)-(14) to be the subsets of [0, L]defined respectively by ω(W) = {x∈(0, L)such that W(x)∈(RN\γ)∩(RN\eγ)} σ(W) = (0, L)\ω(W) Then we have Theorem 1We consider W: [0, L]→RNa stationary solution of class C2of the system (1). Then, a) If M(x, W)does not have any singular point, the scheme (11)-(14) balances the system (1) for W up to the second order. b) If the set of singular points of W,σ(W), has zero measure, then the scheme (11)-(14) asymptotically balances the system (1) for W. c) If Dhas the same eigenvectors as Mand the eigenvalues of Dvanish when the eigenvalues of M vanish, then the scheme (11)-(14) asymptotically balances the system (1) for W. The proof of this result is given in the Appendix, as it essentially uses Taylor expansions and finitedimensional compactness arguments. We also have the following result that yields sufficient conditions for the exact calculation of stationary solutions at grid points: Theorem 2Let Wbe a stationary solution of (1). Assume that FC,GC,SD,BLand BRevaluated on Wverify FC(xi, xi+1, Wi, Wi+1)−FC(xi−1, xi, Wi−1, Wi) ∆x=GC(xi−1, xi, xi+1, Wi−1, Wi, Wi+1)+ +1 2BL(xi−1, xi, Wi−1, Wi)Wi−Wi−1 ∆x+BR(xi, xi+1, Wi, Wi+1)Wi+1 −Wi ∆x,(17) Wi+1 −Wi ∆x= ] M−1(xi, xi+1, Wi, Wi+1)SD(xi, xi+1, Wi, Wi+1),(18) then, the scheme defined by (11)-(14) balances exactly system (1) for the stationary solution at grid points. And this occurs independently of the choice of the upwind matrix D(x, y, U, V ). To prove this Theorem it is enough to verify that under these hypothesis the numerical scheme is reduced to Wn+1 i=Wn i. These results mean that the general structure of our schemes allows to separately balance the numerical flux and the numerical diffusion terms, in a way largely independent of the discretization parameters appearing in each of them. In fact, the role of the matrix Dis mainly to give some stability to the scheme and not to contribute to the balancing of source terms. It remains, however, to construct matrices Dthat satisfy Hypothesis H1 f) of Theorem 1 in order to ensure the well-balance property for all stationary solutions. The construction of such matrices is discussed in Section 4 Remark 2 The definition of ] M−1given by (15) satisfies Hypothesis e)if the intermediate state f W satisfies for (x, U)and (y, V )close enough f W(U, V ) = U+V 2+ρU+V 2(V−U) + O(|V−U|2),(19) where ρ:RN→Ris a C1function. Indeed, if L[(x, U),(y, V )] ⊂RN+1 \γ, as (x+y 2,f W(U, V )) ∈ RN+1 \γand M−1is C2in RN+1 \γ, we deduce ] M−1(x, y, U, V ) = M−1x+y 2,f W(U, V )=M−1x+y 2,U+V 2+ 8 +∂W(M−1)x+y 2,U+V 2f W(U, V )−U+V 2+O f W(U, V )−U+V 2 2!= =M−1x+y 2,U+V 2+∂W(M−1)x+y 2,U+V 2ρU+V 2(V−U) + O(|V−U|2). We thus obtain Hypothesis e). Remark 3 Hypothesis f)is also recovered if the diffusion matrix is defined similarly to ] M−1: D(x, y, U, V ) = X((x+y)/2,f W)e ΛD(x, y, U, V )X−1((x+y)/2,f W), where Xdenotes the matrix composed by the eigenvectors of D, and e ΛD(x, y, U, V ) = Diag(e dj(x, y, U, V ), j = 1,...,N) with e dj(x, y, U, V ) = djif L[(x, U),(y, V )] ⊂RN\γ; 0if not; (20) where djdenote the eigenvalues of D((x+y)/2,f W). Notice that this is an usual way to define the diffusion matrix. Remark 4 Property (19) is satisfied by intermediate states of usual solvers such as Van-Leer’s or Roe’s for Shallow Water equations. Remark 5 The equalities (17) and (18) are discrete versions of the identities ∂ ∂xF(x, W(x)) = G(x, W(x)) + B(W(x))∂W ∂x , and ∂W ∂x (x) = M−1(x, W(x)) ·S(x, W(x)). 4 Some possible choices of matrix D. In this section we discuss some possible choices of the diffusion matrix Dwith a double purpose: To ensure the asymptotically well-balanced property for all steady smooth solutions, and the stability of the numerical scheme at least for the linear homogeneous case. To start with this last purpose, we consider the linear system, ∂W ∂t +A∂W ∂x =B∂W ∂x ,(21) where Aand Bare constant matrices. We suppose that the system is hyperbolic, that is, the matrix M=A−Bhas Nreal different eigenvalues mj∈R,j= 1,···, m. We study under which conditions the family of schemes previously introduced is L2stable. To do this, we observe that these schemes can be rewritten in viscous form, so the L2stability condition is given from the following result (Cf. Godlewsky and Raviart [12]) Theorem 3Let us denote by dj∈Rj= 1,...,N the eigenvalues of D. If Dand Mhave the same eigenvectors and their eigenvalues verify ∆t ∆xmj2 ≤∆t ∆xdj≤1, j = 1,...,N, (22) then the scheme is L2stable. 9 Theorem 5The numerical scheme defined by (29)-(44) exactly balances system (24) on the stationary solution (45) for any value of the parameter αand for any matrix D. To prove this theorem it is enough to verify conditions (17) and (18) on the stationary solution (45). Remark 7 As we previously mentioned, solving the one-layer Shallow-Water system with centered numerical fluxes FCas (29) with a good choice of the parameter α(α= 1/8) allows to correctly approximate sonic rarefaction waves, without the need of entropy-correction techniques. Also, in the case of stationary shocks the best choice is given by α= 0. This suggests that a better method can be obtained if the parameter αis replaced by a function α(x)taking the value 0at inter-cells where a stationary shock develops and 1/8elsewhere. Taking in mind that transitions depend on the sign of the internal eigenvalues, a suitable approximation αj+1/2to α(xj+1/2)can be given by αj+1/2=min(α− j+1/2, α+ j+1/2); with α± j+1/2=¯αif sgn(λi± j) = sgn(λi± j+1) ¯α(1 −0.5(1 + sgn(λi± j−λi± j+1))) otherwise,(46) where ¯αis a constant value (¯α= 1/8) and λi± jare the internal eigenvalues of the matrix M(xj, Wj). In most applications, this numerically-defined function α(x)will be constant in all the domain but a set whose measure tends to zero (a bounded number of cells), so that the balance properties are preserved. 6 Numerical tests In this section, some experiments illustrating the properties of the proposed family of numerical schemes are presented. The purposes of these experiments are mainly to validate the numerical schemes and to investigate the influence of the parameter αwhen the numerical viscosity vanishes. For this reason in all our tests we have chosen a unique numerical diffusion matrix D, given by the simplest choice in Section 4, i. e., D=|M|. The first experiment has been designed to verify whether the still water solution is exactly computed in practice, even for very irregular geometries. The second one provides a numerical study of the order of the numerical scheme. This study confirms the theoretical results. In section 6.3 the influence of the parameter αin smooth transitions is investigated: we observe that any positive value yields good results. Finally, the last test has been designed to analyze the influence of the parameter αfor a non smooth solution. In this case we see that the best results are obtained if the constant parameter αis replaced by the function αdefined in Remark 6. 6.1 Water at rest in a channel with irregular geometry. In this example, a numerical test is presented to verify how accurately the solution (45) is, in practice, computed for a channel with very irregular geometry. To do this, we have considered a 4m. long channel with rectangular cross-section and for which depth and breadth were randomly generated. The numerical scheme was applied for the values of the parameter α= 0, 1/8, 1/4 and 3/8 (only the values for α= 1/8 are shown, because all results are very similar), taking as initial condition a steady solution representing water at rest. The CFL parameter is set to 0.95 and ∆x= 0.05. Figures 2(a) and 2(b) show the bottom topography and the channel breadth, respectively. Figures 2(c) and 2(d) depict the free surface and the water interface at the initial time and after 15 seconds of simulation. Finally, figures 2(e) and 2(f) show the value of the discharge at each layer at the initial state and after 15 seconds of simulation, respectively. As predicted in Theorem 5, for any value of the parameter, the water elevation at the free surface and the interface as well as the discharges at each layer are accurately computed: only round-off errors are present and they do not grow with time. 16 Table 1: Test case 6.2 solved with α= 1/8. h1and h2. N. cells L1error h1L1order h1L1error h2L1order h2 25 5.20 ×10−3– 5.50 ×10−3– 50 1.66 ×10−31.65 1.87 ×10−31.56 100 4.20 ×10−41.98 4.83 ×10−41.95 200 1.04 ×10−42.02 1.15 ×10−42.07 400 2.38 ×10−52.12 2.72 ×10−52.08 800 5.67 ×10−62.07 6.48 ×10−62.07 6.2 Well-balancing test In this section, we test the order of the numerical scheme when it is applied to approach a stationary solution of the two-layer shallow water system with constant breadth function. The regular stationary solutions of this system satisfy:                      Q1= constant, v2 1 2−v2 2 2+g0h1= constant, Q2= constant, v2 1 2+g(h1+h2+zb) = constant, (47) where g0= (1 −r)gis the reduced gravity. In this numerical test, we consider σ(x, z) = 1.0, r= 0.98 and the depth function zb(x) = 0.5e−x2−2.0, x ∈[−3,3]. In order to obtain a particular stationary solution, we impose first its values at x=−3: h1(−3,0) = 0.5, h2(−3,0) = zb(−3) + 1.5, Q1(−3,0) = 0.15 and Q2(−3,0) = −0.15. From these values, we compute the constants in (47). Once the constants are known, the values of the solution at any point are obtained by solving (47). In the numerical experiment, this stationary solution is taken as initial condition. As the stationary solution chosen does not correspond to water at rest, the numerical scheme should be second order accurate. We have applied the numerical scheme for different values of the parameter α= 0, 1/8, 1/4 and 3/8 (only the values for α= 1/8 are shown, because all results are very similar). The CFL coefficient has been taken as 0.9. The results obtained for α= 1/8 are shown in Tables 1 and 2. As expected, the order obtained is 2. 6.3 Rectangular channel with a bump The validation of numerical schemes for solving the two-layer shallow water system with transcritical solutions is not a simple task, as exact solutions with smooth transitions cannot be easily obtained for this system. Nevertheless, Farmer and Armi in [9] introduced, in the context of the study of exchange flows through channels, an asymptotic technique to obtain steady solutions for channels with simplified geometries. In this particular case, they consider a channel of uniform width with a sill on the floor, separating two broad and deep reservoirs. Let us first briefly recall the main aspects of the theory developed by these authors. In their works, smooth steady state solutions are studied. The equations of the model are obtained by replacing the 17 Table 2: Test case 6.2 solved with α= 1/8. Q1and Q2. N. cells L1error q1L1order q1L1error q2L1order q2 25 8.42 ×10−4– 6.31 ×10−4– 50 2.30 ×10−41.87 1.84 ×10−41.78 100 5.96 ×10−51.95 4.89 ×10−51.91 200 1.38 ×10−52.11 1.23 ×10−51.99 400 3.11 ×10−62.14 2.97 ×10−62.05 800 6.91 ×10−72.17 6.79 ×10−72.13 fourth equation in (47) by the rigid lid assumption: h1+h2+zb= constant. (48) The equations are first written in non-dimensional form using some dimensionless variables Q0 i,h0 i, z0 bthat are such that h0 1+h0 2+z0 b= 1 at the section of minimal depth. Then, all the variables are expressed in terms of the internal Froude numbers, Fi=vi/p(g0hi), i= 1,2. Using the approximation r= 1 the following 2 ×2 nonilinear system is obtained from the second equation of (47) and (48): F−2 3 1+1 2F 4 3 1−1 2q−2 3 rF 4 3 2= ∆H0(Q0 1)−2 3,(49) F−2 3 1+q−2 3 rF−2 3 2="Q0 1 (1 −z0 b)3 2#−2 3 ,(50) where ∆H0is a constant corresponding to the dimensionless energy difference between the two layers and qr=Q1/Q2. Here, only the case |qr|= 1 will be considered. The solutions of the equations can be identified with the family of curves in the (F2 1,F2 2)-plane defined by (49). These curves haven been plotted in detail in [9]. If the reservoirs are sufficiently deep, except in the neighborhood of the sill, the speed of the thicker layer is negligible away from the influence of the sill. In this case, the Froude numbers F1and F2 have to be near of zero far from the sill at the reservoir corresponding to water of density ρ1and ρ2, respectively. A close inspection of the family of curves (49) shows that there is only a family of curves for which these two requirements are satisfied. Among the curves of this family, there is one for which the value of the exchanged flow Q0 1is maximal. This curve corresponds to ∆H0= 1.5 and is called the maximal two-way exchange: it represents the steady-state limit that occurs when the fluids are free to move in opposite directions across the sill. If the reservoir corresponding to water of density ρ1 is placed to the left, this solution is critical at the section of minimal depth, subcritical to the right of this section and supercritical to the left. If the condition F2= 0 is imposed at the exit section of the channel (that is, if the depth of the reservoir to the left is assumed to be infinity) the flow is critical again at this section. In the experiment, we consider a rectangular channel of 1 meter width and 6 meters long. Its bottom topography is defined by the function zb(x) = exp(−x2)−2, where x∈[−3,3]. The ratio between fluxes is set to 1 and the ratio between densities is set to r=ρ1/ρ2= 0.98. In order to reproduce the maximal two-way exchange with the numerical model, we solve the classical lock exchange problem. This problem consists in taking as initial state the two fluids separated by a vertical “artificial barrier” that, in this case, is located at x= 2 (see figure 4(a)). This barrier is dropped out at time t= 0 and the fluids are let to evolve until a stationary state is reached (see figure 4(b)). The CFL parameter is set to 0.9 and ∆x= 1/15. The boundary conditions in this test have been designed to simulate the presence of the infinite reservoirs at both ends of the channel and the fact that qr= 1: the two layers must be allowed 18 to go out of the channel as freely as possible, with the only constraint Q1=−Q2. To do this, we have modified the usual numerical treatment based on the introduction of two ghost cells and the duplication of states, which is used very often in practice to simulate open boundaries. The idea is as follows: let us suppose that, at the n-th time iteration, the approximation obtained at the first cell is Wn 1. If the corresponding value of Q2,Qn 2,1, is equal to zero, the lower layer has not yet reached the left boundary. In this case, the state is duplicated at the ghost cell: Wn 0= Wn 1, and the first cell is treated as an internal one during the calculation of the approximations at time tn+1. If, instead, Qn 2,1is not equal to zero, then the state at the ghost cell is modified so that the relation Q1=−Q2holds. More precisely, we define: Wn 0= Wn 1−∆Q |Qn 1,1|+|Qn 2,1|    0 |Qn 1,1| 0 |Qn 2,1|     being ∆Q=Qn 1,1+Qn 2,1. Again, the first cell is treated as an internal one. A similar treatment is performed at the right boundary. During the computations, complex eigenvalues appear. Following [10], when the appearance of complex eigenvalues is detected, a linear friction term is added to the momentum equations of both layers: −CQk, k = 1,2, where the local friction coefficient Cis great enough to make the numerical scheme linearly stable. These coefficients are function of the imaginary part of the eigenvalues. This extra amount of friction is intended to simulate the loss of momentum induced by the instabilities in real flows. Figures 4(b) and 4(c) show the interface at the steady state reached for the different values of the parameter α= 0, 1/8, 1/4, and 3/8. These steady states are compared with the interface corresponding to the maximal two-way exchange obtained in [9] (in the sequel, this solution will be called the F&A solution). To do this, once the value Q0 1has been calculated, the 2 ×2 nonlinear system (49)-(50) is solved numerically for every value of z0 bcorresponding to the intercells of the mesh (the graphical study of the solutions performed in [9] gives the necessary hints about the adequate choice of the initial guess for the numerical algorithms). Notice that all the numerical solutions differ slightly from the F&A solution at the right part of the channel (see 4(b)). These differences are due mainly to the rigid lid hypothesis. The main differences among the numerical solutions are located near the critical section, as expected. The parameter ¯α has practically no influence in regular zones, but it is important in transitions (regular or not). So, we focus only in a neighborhood of x= 0 ( x∈[−0.4,0.4]). Figure 4(c) shows the computed interfaces and the F&A interface in this zone. Notice that, if the parameter αis set to 0, the numerical scheme corresponds to a Roe-type scheme and it produces a non-entropy solution, as expected. Table 6.3 shows the L∞and L1errors corresponding to the comparison of the numerical interfaces with the analytical one. Observe that the values α= 1/8, α= 1/4, α= 3/8 produce similar results. The discharges computed by the model are Q1=−Q2= 9.38 10−2m3/s for α= 1/8, 1/4 and 3/8, and Q1=−Q2= 9.33 10−2m3/s for α= 0, while the F&A model gives Q1=−Q2= 9.213 10−2 m3/s. 6.4 Non smooth solution The purpose of this numerical test is to compare the behavior of the numerical schemes in the neighborhood of a discontinuity. We consider the two-layer shallow water system for channels with rectangular 19 α L∞-error L1-error 0 0.0331 0.0074 1/8 0.0085 0.0045 1/4 0.0085 0.0045 3/8 0.0085 0.0046 Table 3: L∞and L1errors (F&A interface against numerical computed ones) (x= [−0.4,0.4]) cross-section with constant breadth (σ= 1) and a flat bottom topography. In this case, the system reduces to: ∂W ∂t +A(W)∂W ∂x = 0,(51) where Ais the 4 ×4 matrix whose expression is given by: A(W) =            0 1 0 0 −Q2 1 A2 1 +gA12Q1 A1 gA10 0 0 0 1 rgA2 20−Q2 2 A2 2 +gA22Q2 A2            . The difficulty in defining a discontinuous solution of this system comes from the fact that the non-conservative products do not make sense within the distribution theory for such a solution. Nevertheless, there are several mathematical theories allowing to give a sense to these products in a weaker manner. In particular we consider here the definition of non-conservative product as Borel measures introduced by Volpert in [21]. This definition can be considered as a particular case of the concept of non-conservative product developed in [8]. In this latter work, the definition of the non-conservative product depends on the choice of a family of viscous paths. Volpert’s definition is equivalent to the choice of the family of segments. According to this definition, across a discontinuity with speed ξa weak solution of the system must satisfy the generalized Rankine-Hugoniot condition Z1 0AW−+s(W+−W−)·(W+−W−)ds =ξ(W+−W−),(52) where W−, W+are the left and right limits of the solution at the discontinuity. In [15], a Roe-type numerical scheme for the two-layer shallow water system was studied which is the particular case of the here introduced family that corresponds to the choice α= 0. In this reference it was shown that this numerical scheme is consistent with the notion of weak solutions associated to Volpert’s definition of nonconservative product in the following sense: let us suppose that the approximations Wn iand Wn i+1 obtained by the numerical scheme in two neighbor cells can be linked by a shock of speed ξsatisfying the jump condition (52). Then, ξis an eigenvalue of the intermediate matrix of the scheme, Ai+1/2, and Wn i+1 −Wn iis an associated eigenvector. As a consequence, the solution of the linear approximate Riemann solver used by the numerical scheme to advance in time is equal to that of the exact Riemann problem. In particular, this implies that the numerical scheme solves exactly stationary shocks. In order to numerically verify this property, in the cited reference the authors proposed a test consisting in solving the system with values g= 10 and r= 0.02, and the initial condition: W(x, t) = WLif x < 0 WRif x > 0, 20 with WL=    1 √0.1 1 √20    , WR=    0.3961560 √0.1 1.5820186 √20    . These states satisfy the generalized jump condition (52) with an error of about 10−5. We have applied to this test problem the numerical scheme with values of the parameter α= 0, 1/8, 1/4, and 3/8. If the numerical schemes were consistent with the jump condition (52), they should obtain a stationary discontinuity placed at x= 0 with values of the unknowns close to the initial conditions. Figure 5 depicts the steady state solutions (water depths and discharges) for the different values of the parameter ¯αobtained with ∆x= 0.01 and CFL= 0.90 for a rectangular channel of 1m long. The errors are shown in tables 2 and 3. Notice that in all the cases the numerical scheme behaves as expected: a stationary discontinuity placed at x= 0 and close to the initial condition is captured. h1 α L∞-error L1-error 0 0.0111 2.24·10−4 1/8 0.0963 0.0017 1/4 0.4693 0.0059 3/8 0.4817 0.0061 h2 α L∞-error L1-error 0 0.0114 1.14·10−4 1/8 0.0972 0.0015 1/4 0.4434 0.0055 3/8 0.4574 0.0056 Table 4: L∞and L1errors (water depths) Q1 α L∞-error L1-error 0 0.0114 4.33·10−4 1/8 0.0995 0.0014 1/4 0.1273 0.0018 3/8 0.1406 0.0020 Q2 α L∞-error L1-error 0 0.0127 1.28·10−4 1/8 0.1467 0.0021 1/4 0.2098 0.0031 3/8 0.2487 0.0039 Table 5: L∞and L1errors (discharges) From Tables 2 and 3 and Figure 5 it is clear that the smallest errors correspond to α= 0. For the other values an over-shooting near the stationary shock occurs (see figure 5), becoming larger for increasing values of ¯α. We have applied finally to the test problem the numerical scheme with the parameter αreplaced by the function α(x) defined in Remark 7 with ¯α= 1/8. Figure 6 depicts the steady state solution (water depths and discharges) obtained with again ∆x= 0.01 and CFL= 0.90. It coincides with that obtained with α= 0. In this case, once the stationary solution is reached the corresponding function α(x) takes the value 0 only at the inter-cell where the discontinuity stands and 1/8 at the others: figure 7 shows the values of α(x) close to the discontinuity. 7 Conclusion We have addressed in this paper the accurate numerical solution of non-conservative non-homogeneous hyperbolic systems in one space dimension by finite volume solvers. We have introduced a general class of schemes that we have proved to be “asymptotically well-balanced” in the sense that they 21 calculate, with at least order two, all stationary solutions of the system in all the domain but on a set whose measure tends to zero as ∆xtends to zero. One of the free parameters of this class is the diffusion matrix, so we include both flux-splitting (for conservative equations) and flux-difference solvers as particular cases. We also give general sufficient conditions under which the schemes exactly calculate a given stationary solution. These conditions may be read as meaning on one hand that the numerical flux of the scheme must compensate the centered part of the numerical source term, and on another hand that the numerical diffusion must compensate the decentered part of the numerical source term. In the case of two-layer Shallow Water flows, we finally have introduced a specific solver to deal with the numerical difficulties associated to non-conservative products combined with nonlinear convection and source terms effects. On one hand, this solver is able to exactly compute water at rest. On another hand, it yields accurate computation of subcritical to supercritical transitions without need of entropy corrections, and also of non smooth solutions due to the presence of the non-conservative products. Our main technical innovation is an adaptive discretization of the centered part of the numerical flux which, in our opinion, may be extended to construct accurate solvers for more general non-conservative non-homogeneous equations. Appendix Proof of Theorem 1: We consider first a regular stationary solution Wof (1) and a compact set, K⊂ω(W). Let Tbe the continuous operator appearing in system (1), T(W) = ∂ ∂x[F(x, W)] −G(x, W )−B(x, W)∂W ∂x and Ththe discrete operator appearing in the numerical scheme (11)-(14), Th(W) = φS R−φS L ∆x−GC (where we omit the arguments for brevity). We must prove that Th(W)(xi) = T(W)(xi) + O((∆x)2),∀xi∈K. The proof is based on identity (8): e T(W) = ∂ ∂xF(x, W)−νD(x, W )(∂W ∂x −M−1S(x, W ))−G(x, W)−B(x, W)∂W ∂x . It is enough to prove that every term appearing in Th: αi=FC(xi, xi+1, Wi, Wi+1)−FC(xi−1, xi, Wi−1, Wi) ∆x, βi=νD(xi, xi+1, Wi, Wi+1)(Wi+1 −Wi)−D(xi−1, xi, Wi−1, Wi)(Wi−Wi−1) (∆x)2(53) ηi=ν1 ∆xD(xi, xi+1, Wi, Wi+1) ] M−1(xi, xi+1, Wi, Wi+1)SD(xi, xi+1, Wi, Wi+1)− −D(xi−1, xi, Wi−1, Wi) ] M−1(xi−1, xi, Wi−1, Wi)SD(xi−1, xi, Wi−1, Wi)(54) δi=GC(xi−1, xi, xi+1, Wi−1, Wi, Wi+1), 22 λi=1 2∆xBR(xi, xi+1, Wi, Wi+1) (Wi+1 −Wi)−BL(xi−1, xi, Wi−1, Wi) (Wi−Wi−1) balances up to the second order the corresponding term in T(W). We prove the case of the more complex term ηi. The proof is similar for the other terms. We define R(x, U) = D(x, U)M−1(x, U)S(x, U), RC(x, y, U, V ) = D(x, y, U, V ) ] M−1(x, y, U, V )SD(x, y, U, V ), and we want to prove that ∆RC(xi) = RC(xi, xi+1, Wi, Wi+1)−RC(xi−1, xi, Wi−1, Wi) = = ∆x∂ ∂x[R(x, W (x))]|x=xi+O((∆x)2).(55) As Kis a compact set and RN+1 \γis an open set, there exists δK>0, such that if |U−W(x)|+ ∆x < δK, then (x+ ∆x, U)⊂RN+1 \γ, for all x∈K. On another hand, if Wi+θ= (1 −θ)Wi+θWi, then using the Mean Value Theorem we prove |Wi+θ−W(xi+θ)| ≤ ∆x 2 ∂W ∂x L∞(K) ;∀xi∈K. So, if ∆x < δK/(1 + k∂W ∂x k/2), we obtain that (xi+θ, Wi+θ)⊂RN+1 \γ, for θ∈[0,1]. Then if ∆xverifies the previous inequality, we obtain that L[(xi, Wi),(xi+1, Wi+1)] ⊂RN+1 \γfor all xi∈K. Then, we can apply the hypothesis H1 c), d) and f) and conclude that RC(xi, xi+1, Wi, Wi+1) = Rxi+1/2,Wi+Wi+1 2+ +R1xi+1/2,Wi+Wi+1 2 xi+1 −xi Wi+1 −Wi+O(|xi+1 −xi|2+|Wi+1 −Wi|2); for some C1matrix function R1. As Wis of class C2, we have that R1xi+1/2,Wi+Wi+1 2 xi+1 −xi Wi+1 −Wi−R1xi−1/2,Wi−1+Wi 2 xi−xi−1 Wi−Wi−1= = ∆xR1(xi+1/2, W(xi+1/2)) 1 ∂W ∂x (xi+1/2)!− −R1(xi−1/2, W(xi−1/2)) 1 ∂W ∂x (xi−1/2)!+O((∆x)2) = = (∆x)2∂ ∂xR1(x, W(x)) 1 ∂W ∂x !x=xi +O((∆x)2) = O((∆x)2).(56) 23 Moreover, Rxi+1/2,Wi+Wi+1 2−Rxi−1/2,Wi−1+Wi 2= ∆x∂ ∂xR(x, W(x))x=xi +O((∆x)2).(57) From (56) and (57) we conclude that (55) is verified. We can now prove the statements a) and b) of Theorem 1: a) If Whas no singular points, then ω(W) = (0, L) and it is enough to use the preceding argument for some compact set K⊂ω(W) containing all the internal nodes of the discretization. b) If meas(σ(W)) = 0, as ω(W) is an open set, we may find an increasing sequence of compact sets {Kn}n≥1such that ∪n≥1Kn=ω(W). Then, meas((0, L)\∪n≥1Kn) = 0. It is enough to use the preceding argument for each Kn. To prove statement c) let us remark that, by the previous argument, the scheme is asymptotically well-balanced in ω(W). Thus, it only remains to prove this property in σ(W) when this set has non-zero measure. Let us denote by int(σ(W)) the interior of σ(W). As meas(σ(W)−int(σ(W))) = 0, it is enough to prove the property in int(σ(W)). And, as we shall see this occurs because in int(σ(W)) Dand M−1coincide with smooth functions. For sake of simplicity, we shall assume N= 2 and that the first eigenvalue of Mvanishes in σ(W). Let us define the operator e T(W) = ∂ ∂xF(x, W)−H(x, W )∂W ∂x −C(x, W )S(x, W)− −G(x, W)−B(x, W )∂W ∂x ; where H(x, W) = X(x, W)0 0 0d2(x, W)X−1(x, W ); C(x, W) = X(x, W )0 0 0λ−1 2(x, W)X−1(x, W ). By d2and λ2we denote the second eigenvalue of Dand M, respectively. As Wis a steady solution of (1), M(x, W)∂W ∂x (x) = S(x, W). Then, C(x, W)S(x, W ) = C(x, W)M(x, W )∂W ∂x (x). But C(x, W)M(x, W )∂W ∂x (x) = X(x, W)0 0 0 1 X−1(x, W)∂W ∂x (x) = =X(x, W)1 0 0 1 −10 0 0X−1(x, W )∂W ∂x (x) = =∂W ∂x −X−1(x, W )∂W ∂x (x)1 X1(x, W), where by [ ·]1we denote the first component of a vector and by X1the first eigenvector of M. 24 Then, as X1is also the first eigenvector of H, H(x, W)∂W ∂x (x)−C(x, W )S(x, W)=λH,1(x, W)X−1(x, W )∂W ∂x (x)1 X1(x, W) = 0, as by construction, the first eigenvalue λH,1of His zero. We deduce that e T(W) = T(W). It is then enough to prove that for any compact set K⊂int(σ(W)) there exists δK>0 such that if 0<∆x < δK, then Th(W)(xi) = e T(W)(xi) + O((∆x)2),∀xi∈K. (58) Les us assume δK≤dist(K, ∂σ(W)). If xi∈K, then xi−1, xi+1 ∈int(σ(W)) and we may ensure that L[(xi, Wi),(xi+1, Wi+1)] cuts γ=eγ. Then, using our definition of ] M−1and D, we have ] M−1(xi, xi+1, Wi, Wi+1) = =X(xi+1/2,f Wi+1/2)0 0 0λ−1 2(xi+1/2,f Wi+1/2)X−1(xi+1/2,f Wi+1/2) = C(xi+1/2,f Wi+1/2), where we denote f Wi+1/2=f W(xi, xi+1, Wi, Wi+1). Also, D(xi, xi+1, Wi, Wi+1) = =X(xi+1/2,f Wi+1/2)0 0 0d2(xi+1/2,f Wi+1/2)X−1(xi+1/2,f Wi+1/2) = H(xi+1/2,f Wi+1/2). Now, Thcan be written as in the previous case by replacing Dby Hand ] M−1by Cin the definitions of βi,ηiand the proof is similar.  References [1] A. Berm´udez, M. E. V´azquez Cend´on, Upwind Methods for Hyperbolic Conservation Laws with Source Terms. Computers & Fluids 23(8), pp. 1049-1071 (1994). [2] M. J. Castro, J. Mac´ıas, C. Par´es, A Q-scheme for a class of systems of coupled conservation laws with source term. Application to a two-layer 1-d shallow water system. M2AN 35(1), pp. 107-127 (2001). [3] M. J. Castro, J. Mac´ıas, C. Par´es, J. A. Rubal, M. E. V´azquez, Two-layer numerical model for solving exchange flows through channels with irregular geometry. Proceedings of “ECCOMAS 2001”, Swansea (2001). [4] M. J. Castro, J. A. Garc´ıa, J. M. Gonz´alez, J. Mac´ıas, C. Par´es, M. E. V´azquez, Numerical simulation of two-layer Shallow Water flows in channels with irregular geometry. Journal of Comp. Phys. 195, pp. 202-235 (2004). [5] T. Chac´on Rebollo, A. Dom´ınguez Delgado, E. D. Fern´andez Nieto, A family of stable numerical solvers for Shallow Water equations with source terms. Comput. Methods Appl. Mech. Eng. 192, pp. 203-225 (2003). 25 −0.2 −0.15 −0.1 −0.05 0 0.05 0.1 0.15 0.2 0 0.02 0.04 0.06 0.08 0.1 0.12 Figure 7: α(x) near the discontinuity placed at x= 0. 32