Coupling superposed 1D and 2D shallow-water models: Source terms in finite volume schemes

Coupling Superposed 1D and 2D Shallow-Water Models: Source Terms in Finite Volume Schemes E.D. Fern´andez-Nieto aJ. Marin bJ. Monnier c,∗ aDepartamento de Matem´atica Aplicada I, Universidad de Sevilla. E.T.S. Arquitectura. Avda, Reina Mercedes, s/n. 41012 Sevilla, Spain bUniv. Grenoble & INRIA, Lab. LJK, B.P. 53, F-38041 Grenoble Cedex 9, France cUniv. Toulouse - INSA, Institut de Math´ematiques, F-31077 Toulouse cedex 4, France Abstract We study the superposition of 1D and 2D shallow-water equations with non-flat topographies, in the context of river - flood modeling. Since we superpose both models in the bi-dimensional areas, we focus on the definition of the coupling term required in the 1D equations. Using explicit finite volume schemes, we propose a definition of the discrete coupling term leading to schemes globally well-balanced (the global scheme preserves water at rest whatever if overflowing or not). For both equations (1D and 2D), we can consider independent finite volume schemes based on well-balanced Roe, HLL, Rusanov or other scheme, then the resulting global scheme remains well-balanced. We perform a few numerical tests showing on the one hand the well-balanced property of the resulting global numerical model, on the other hand the accuracy and robutness of our superposition approach. Therefore, the definition of the coupling term we present allows to superpose a local 2D model over a 1D main channel model, with non-flat topographies and mix incoming-outgoing lateral fluxes, using independent grids and finite volume solvers. Key words: shallow-water equations, superposition, coupling, finite volume, well-balanced, river, flood plain. ∗Corresponding author: Jerome.Mo[email protected] Preprint submitted to Elsevier 2 February 2010 1 Introduction In river hydraulics, operational models are generally based on the Saint-Venant equations (1D shallow-water). If overflowing, flood plain are represented in the 1D model by storage areas, that are defined by using empirical laws and/or terms to be calibrated, see [6] or e.g. [14]. Obviously, flow dynamic inside the storage areas is not computed; also the empirical laws can be difficult to calibrate. If for any reason, the end-user has to model the flow in the flood plain, a 2D model must be used. Then, the classical approach is to decompose the domain, re-define the mesh, then couple the 1D model (in the non-flooded areas) with 2D models (in flooded areas) at interfaces, see e.g. [15,13]. Coupling conditions have to be imposed at interfaces only. An efficient coupling procedure may be a Schwarz-like algorithm. Nevertheless, this approach presents some drawbacks. It requires to re-define the 1D hydraulic model (mesh, boundaries etc) and very probably the related topography data. The 1D model (which is potentially a complex network) must be segmented (decomposed) in order to combine it with the 2D models. It can be a heavy task. A superposition approach is proposed in [8] and [12]. In such an approach, instead of decomposing the original 1D (network) model, one superposes the 2D model (socalled ”local zoom model”). The superposition approach presents some advantages. The original 1D model remains intact and the 2D local models can be performed with their own dynamics (typically, time steps and mesh grids are much smaller for 2D solvers than for 1D solvers). Nevertheless, an accurate definition of the coupling terms between both models is required. At interfaces, incoming characteristics are still good conditions, but along the 1D main channel one must introduce a coupling term in the 1D equations (modelling the loss or gain of mass and momentum). This coupling term has to take into account the outgoing and incoming fluxes if overflowing. From a continuous point of view, the coupling source term can be derived formally from the 3D Navier-Stokes equations, see [12]. Then, next step is to define a stable and well-balanced global scheme . An important difficulty is to discriminate between the 1D-topography graph Zb(˜x) and the 2D-topography graph zb(x, y), since Zbdepends on the curvilinear coordinate ˜x while zbdepends on the cartesian coordinates (x, y). In addition for real cases, data are sparse, uncertain, and the 1D topography and the 2D topography do not have to respect the same hydrological constraints. If in addition, one wants to consider different meshes and schemes for the 1D model and the 2D model, the discretization of the coupling source term must be such that it leads to a consistent, stable and well-balanced scheme. This is the problem we address in the present study, while focusing on explicit finite volume schemes. Let us point out that we present the 1D-2D coupling in the context of river hydraulics, but this could also be apply to any other flows involving 1D and 2D shallow-water equations with non-flat topographies. We consider the possibility of using different finite volume schemes for the 1D model 2 and the 2D model. In the numerical analysis presented in next sections, they can be based on different time-space grids but they must be explicit in time. More precisely, we consider finite volume methods in conservative form with source terms (the topography terms and the coupling term). Then, for both models, we can consider any solver belonging to a whole familly of approximate Riemann solvers. We prove that the resulting global scheme is well-balanced in the sense that it preserves water at rest, with and without overflowing. We present some numerical results for an academic test case with a non-constant topography in which there are outgoing and incoming lateral fluxes. In order to couple the models, we use a Schwarz coupling algorithm (global in time). This could be done also by using an optimal control approach as in [8], [12]. The numerical results show that after convergence, the coupling source term Ψ defined in the present study leads to a global solution as accurate as a full 2D solution in case of matching grids, and leads to a robust and accurate solution if grids are mismatching. The paper is organized as follows. In Section 2, we present the two mathematical models. Their discretization using well-balanced finite volume scheme is presented in Section 3. The discretization of the coupling source term in 1D equations is described in Section 4. We begin with the simplest case (matching grids and 1D linear axis). Then, we consider the case of 1D curvilinear geometry with matching grids. Finally, the most general case (curvilinear and mismatching grids) is considered. We prove in Theorem 1 that the here introduced discrete source term leads to a global wellbalanced scheme, whatever the choice of the well-balanced finite volume method used for the 1D and 2D models. In Section 5 we present some numerical experiments to validate the definition of the discrete source term, and to show the efficiency of the present superposition approach. We recall briefly in Appendix A the derivation of the coupling source term in the 1D equations from the 3D Navier-Stokes equations (we refer to [12] for more details). 2 Mathematical models 2.1 The 1D model with source term The 1D model is based on Saint-Venant equations (1D shallow-water equations). Nevertheless, since our goal is to couple this 1D model to a 2D shallow-water model, we must take into account transfers through the two lateral boundaries of the main channel. If we integrate the 3D Navier-Stokes equations over the vertical wetted area, in the presence of lateral transfer terms, we obtain some source terms in the 1D equations. The derivation of these source terms is presented in Appendix A. 3 The result is the following. Let us denote the channel-following coordinates by: ˜x. We denote the unidimensional variables (i.e. depending on (˜x, t) only) as follows: Sthe wet cross section, Qthe discharge, Hthe water depth. And Zbdenotes the unidimensional topography (depending on ˜xonly). We assume that: the channel width variations are small, uis nearly constant over the cross section, and (u, v) does not depend on zon boundaries b1and b2. Furthermore, for the sake of simplicity, we consider rectangular cross sections only. Hence S=b.H, where bis the channel width, see Fig. A.3. The derivation of the 1D shallow-water equations with source term is presented in Appendix A (replace xby ˜x). Under the assumptions above, the equations are the following:          ∂S ∂t +∂Q ∂˜x=−(qη1+qη2) ∂Q ∂t +∂ ∂˜x Q2 S+P!−g∂b ∂˜x H2 2+gS ∂Zb ∂˜x=−(qη1ut1+qη2ut2) (1) where P=gS H 2is a pressure term. In the right hand side, qηirepresents the discharge normal to the lateral boundary iof the main channel, i= 1,2; utirepresents the tangential velocity at lateral boundary i, see Fig. A.4. (It is the projection of the horizontal velocity vector (u, v) on the lateral boundary direction). System (1) is closed with appropriate boundary conditions and initial conditions. If we set w= [S, Q]Tthe unknown vector, then the 1D model writes as an hyperbolic system with source term: ∂tw+∂˜xf1(b(˜x), w) = g1(˜x, w)∂˜xZb+g2(˜x, w)∂˜xb+ Ψ (2) where f1(b(˜x), w) =   Q Q2 S+gS2 2b   , g1(˜x, w) =   0 −gS   , g2(˜x, w) =   0 g S2 2b2    and Ψ =   −(qη1+qη2) −(qη1ut1+qη2ut2)  (3) is the coupling term in the 1D model i.e. normal discharges and tangential velocities normal to the (two) lateral boundaries of the main channel are provided by the 2D model. 4 2.2 The 2D model The 2D hydraulics model is based on the bidimensional shallow-water equations in their conservative formulation. The unknowns are the water depth hand the local discharge q=hu, where u= (u, v)Tis the depth-averaged velocity vector. In a bi-dimensional domain Ω and for a computational time interval [0, T], equations are:        ∂th+ div( q) = 0 in Ω×] 0, T] ∂tq+ div 1 hq⊗q+1 2g∇h2+gh∇zb+gn2kqk2 h7/3q= 0 in Ω×] 0, T] (4) where gis the magnitude of the gravity, zbthe (bi-dimensional) topography and nthe Manning roughness coefficient. Initial conditions h(0) = h0,q(0) = q0, and boundary conditions are given. In all the sequel, c=√gh, denotes the local wave celerity. Remark 2.1 a) Let us point out that we have to discriminate between the 1D- topography graph Zb(˜x)and the 2D-topography graph zb(x, y). As a matter of fact, Zbdepend on the curvilinear coordinate while zbdepends on the cartesian coordinates. Furthermore, for actual data, the 1D topography and the 2D topography do not have to respect the same hydraulic constraints. b) If the 2D equations are locally coupled to 1D equations (the 2D model is playing the role of a zoom) then appropriate conditions at interfaces 1D-2D can be incoming characteristics, see e.g. [8,12]. If we define the unknown vector W=W(t, ~x) with ~x = (x, y) by W= (h, qx, qy)T, q= (qx, qy), then the system can be written as a 2D hyperbolic system with source terms. By simplicity we will not consider the friction term gn2kqk2 h7/3qin this study so that (4) can be rewritten as follows : ∂tW+∂xF1(W) + ∂yF2(W) = κ1(~x, W)∂xzb+κ2(~x, W)∂yzb(5) where F1(W) =                qx q2 x h+1 2gh2 qxqy h                , F2(W) =                qy qxqy h q2 y h+1 2gh2                5 and κ1(W) =        0 −gh 0        , κ2(W) =        0 0 −gh        3 Well-balanced finite volume methods We present the finite volume schemes considered to discretize the 1D and 2D models. First the numerical schemes for the 2D case are presented, in terms of some 1D numerical flux functions. Then, we present the numerical schemes for the one dimensional system. All the schemes considered here are based on well-balanced finite volume methods. 3.1 Finite volume methods for the 2D system For the two-dimensional system, we consider a partition of the 2D domain Ω in control volumes denoted by Ki. By Ei j we denote the common edges between the control volumes Kiand its neighbour number j;ηi j is the unit normal vector to Ei j outward to the control volume Ki, see Figure A.5. We also denote by |Ki|the area of Ki,|Eij|the length of Eij and Wn i≈1 |Ki|Z Ki W(tn, x, y)dxdy, with tn=tn−1+ ∆t, being ∆tthe time step considered for the 2D finite volume method. The structure of 2D finite volume schemes is as follows: Wn+1 i=Wn i−∆t |Ki|X j∈Ki|Ei j|Φ2D(Wn i, Wn j, ηi j) + ∆tGn i(6) where Gi≈1 |Ki|Z Ki (κ1∂xzb+κ2∂yzb)dxdy. and Φ2D i j = Φ2D(Wn i, Wn j, ηi j) is an approximation of the normal flux related to the edge Ei j. 6 The basic principle of 2D finite volume methods is to neglect tangential variations and to consider a projected 1D Riemann problem. Concretely, one starts from the 2D equations in the form (5). If we denote: F= (F1, F2), κ= (κ1, κ2), if η= (η1, η2) is an unit vector, and η⊥= (−η2, η1), then we can rewrite the previous system as follows: ∂tW+∂η(F ·η) + ∂η⊥(F ·η⊥) = (κ·η)∂ηzb+ (κ·η⊥)∂η⊥zb. Moreover, we can use the invariance rotation property of shallow-water equations (see [16]): (F ·η)(W) = T−1 ηF1(TηW) where Tηis the rotation matrix: Tη=       1 0 0 0η1η2 0−η2η1        . Using this invariance rotation property and multiplying the previous system by Tη we obtain: ∂t(TηW) + ∂ηF1(TηW) + Tη∂η⊥(F ·η⊥) = Tη(κ·η)∂ηzb+Tη(κ·η⊥)∂η⊥zb(7) The simplification of 2D finite volume method consists in neglecting the tangential variations. With this simplification, and by noticing that: Tη(κ·η) = κ1, we obtain the 1D shallow-water equations plus a passive linear transport equation over the direction η: ∂t(TηW) + ∂ηF1(TηW) = κ1(TηW)∂ηzb,(8) where: TηW= (h, qη, qη⊥)T. The third equation of this system writes as follows: ∂tqη⊥+∂ηqηqη⊥ h= 0.(9) Then, in order to approximate the normal flux of the 2D shallow-water equations, (F·η), we consider numerical fluxes of the form: Φ2D(Wi, Wj, ηi j) = T−1 ηi j Φ(Tηi j Wi, Tηi j Wj).(10) 7 where Φ(Vi, Vj) is a numerical flux associated to the following 1D Riemann problem:              ∂tV+∂ξF1(V) = κ1(V)∂ξzb V(ξ, 0) =      V1if ξ < 0 V2if ξ > 0 (11) Let us point out that the 1D numerical flux Φ depends on the source term κ1(W)∂ξzb (i.e. Φ is not a numerical flux function corresponding to the homogeneous problem). In the numerical schemes considered here the expression of Giis the following: Gi=1 |Ki|X j∈Ki |Eij| 2T−1 ηij κ1,i j(¯zb,j −zb,i), where κ1,ij =κ1(Vi+Vj 2). The definition of zb,j takes into account if flood is produced or not, since its definition includes a treatment of wet/dry fronts. Following [3], we propose to define it as follows: zb,j =     zb,j if hj6= 0 zb,i +hiif hj= 0 (12) Different approximate Riemann solvers. We denote by Φ the numerical flux function, corresponding to the 1D system (11). We consider the following family of numerical solvers (see e.g. [5]): Φij = Φ(Vi, Vj) = F1(Vj) + F1(Vi) 2−1 2Pi j{Ai,j (Vj−Vi)−κ1,i j(¯zb,j −zb,i)}(13) By Ai,j we denote the Roe matrix associated to F1, such that: F1(Vj)−F1(Vi) = Aij(Vj−Vi). We suppose that matrix Ai j can be diagonalized, that is Ai j =XijΛijX−1 ij , where Xij is a matrix whose columns are a basis of eigenvectors and Λij is the diagonal matrix defined by the eigenvalues. We denote the eigenvalues of Ai,j by: λ1,ij < λ2,ij < λ3,ij. Depending on the definition of Pi,j, we obtain different methods : 8 •Roe method corresponds to Pij =Xijsgn(Λij)X−1 ij where: sgn(Λij) = diag(sgn(λ1,ij),sgn(λ2,ij),sgn(λ3,ij)).(14) •HLL method corresponds to Pij =S1I+S2A−1 ij , where: S1=Sr+Sl Sr−Sl , S2=−2SrSl Sr−Sl , Sl= min{λ1,ij,0}, Sr= max{λ3,ij,0}.(15) •Rusanov method corresponds to: Pij = max(l=1,2,3){|λl,ij|}A−1 ij . •HLLC scheme. System (11) has a linearly degenerated field associated to the third equation (see (8), (9)), which corresponds to a passive linear transport. To reduce the numerical diffusion on the approximation of the contact discontinuites, related to the linearly degenerated field, we can consider the HLLC scheme, see [16]. It is based on the definition of the physical flux function that verifies: [F1]3= [F1]1qη⊥/h Then, the third component of the numerical flux function is not defined by (13), but as a function of the first component of Φ and an upwind approximation of the passive scalar ϕ∗. We denote the resulting numerical flux by Φc. We consider the following extension of the HLLC method: Φc i j =       [Φij]1 [Φij]2 ([Φij]1)ϕ∗ ij        (16) where ϕ∗ ij =     qη⊥i/hiif S∗<0 qη⊥j/hjif S∗>0 and [Φij]1and [Φij]2are the two first components of the 1D flux defined by (13). Although in HLLC method these two first components are given by the HLL method, any other solver defined by (13), as Roe or Rusanov, can be also used. A consistent definition of the intermediate wave speed S∗for not flat topography has been introduced in [7]. •Some other choices are possible; for example, flux limiter functions or Lax-Wendroff scheme (see [5]). 9 On the other hand, we remember that the definition of the coupling term is an approximation of:   −[F1(TηW)]1 [F1(TηW)]3   =  −[Tη⊥T−1 ηF1(TηW)]1 −[Tη⊥T−1 ηF1(TηW)]2   . The idea is to change η⊥by another vector β. We have: [−TβT−1 ηF1(TηW)]1=−qηand [−TβT−1 ηF1(TηW)]2=−qηuβ+g 2h2β·η. The vector βis introduced in order to approximate utby uβ. If β=±η⊥then [−TβT−1 ηF1(TηW)]3=qηuβ; otherwise, it remains the term (g 2h2β·η). Hence, a correction must be introduced into the definition of Ψiin order to rid of the term (g 2h2β·η). We have: g 2h2β·η=−h 2[TβT−1 ηκ1(W)]2and [TβT−1 ηκ1(W)]1= 0. So, we propose the following definition of the coupling term: Ψi=b Ψ(βi j1, ηi j1) + b Ψ(βi j2, ηi j2), with [b Ψ(β, ηi j)]l=−TβΦ2D(Wi, Wj, ηi j)−1 2T−1 ηi j κ1,i j(zb,j −zb,i) + hi 2T−1 ηκ1,il ,(24) for l= 1,2 and κ1,i =κ1(Wi). Of course, if β=±η⊥then b Ψ equals to Ψ. Remark 4.2 Let us recall that 1D numerical flux is denoted by φ, 2D numerical flux by Φ2D, and Φ2D=T−1 ηΦor Φ2D=T−1 ηΦc.Φand Φcare 1D numerical fluxes defined by (13) and (16) respectively. We point out that definition (24) of the discrete coupling term does not imply any particular relationship between the 1D and 2D solvers involved. Therefore, it is possible to use for example low computational cost scheme for 2D model (e.g. Rusanov method) with a more accurate scheme for 1D model (e.g. Roe solver). Let us state the well-balanced properties of the resulting global scheme with the proposed discrete coupling term (24). We have Theorem 4.1 The scheme defined by (17)-(18) and the coupling term (24) exactly preserves water at rest with or without overflowing. This is true for any choice of the upwinding matrix Pi+1/2for the 1D solver (see (18)) and for any choice of the 2D numerical solver. 16 Proof.- As the 1D method is well-balanced (see [2]), it is enough to prove that in case of water at rest the coupling term vanish, independently if overflowing or not. Let us define: C=F1(Wj)−F1(Wi) 2−1 2κ1,i j(zbj −zbi) (25) For water at rest, we have: C= 0. The second term in the expression of Cwrites: g 2 hj+hi 2(hj+zbj −hi−zi). If hj= 0, in vertu of the definition of zbj to treat dry/wet areas, see (12), we have: zbj =zi+hi. Hence this term equals to zero. Otherwise, the stationary solution of water at rest verifies: h+zb=constant. Hence, again, this equals to zero. Using the property of Roe matrix, the term that multiplies Pi j in the definition of Φi j, is A−1 i j C. Hence, it is equal to zero. So, Φi+1/2=F1(Wi) + F1(Wj) 2. and [b Ψ(β, ηi j))]l=−[TβT−1 ηi j F1(Wi) + F1(Wj) 2−1 2κ1,i j(zbj −zbi) + hi 2κ1,i]l, for l= 1,2. Finally, the term that multiplies TβT−1 ηi j in the previous expression, coincides with C (defined by (25)); therefore it is equal to zero. Then, b Ψ(β, ηi j) = 0 and Ψi= 0. 2 Remark 4.3 The present 2D finite volume methods do not preserve all stationary solutions of the system. As a matter of fact, the 2D finite volume methods neglect the tangential variations of the unknown at interfaces. Nevertheless, for the family of methods presented in Section 3, it is possible to prove that they preserve all regular 1D stationary solutions, up to second order in ∆x(see [4]), if the finite volume mesh is defined by rectangles oriented in the same sense than the solution variations. Furthermore, we can prove the same property for the present coupled system. It consists to prove that Ψi=O(∆x2)for any stationary solution. 17 4.2.4 Subgrid and mismatching grids cases In this section we consider the case of subgrid meshes or mismatching 1D/2D grids, in the sense described in subsection 4.2.1. In both cases, the generalization of the source term Ψ is straightforward. Let us define the projection operator Rthat computes average values over ˜ Kiand δt : Ri(tm, v) = 1 δtk˜ Kik tm+δt Z tmZ ˜ Ki v dΩdt, where ˜ Kiis a 2D volume defined as a subset matching with the 1D mesh. Remember that δt is the time step corresponding to the discretization of the 1D model. We consider a piecewise constant function, FΨi(x, t) defined by the value {Ψij}jat the volumes Kj. Thus, we define: Ψn i=Ri(tn, FΨi). The main difference between the subgrid case and mismatching one, is the computation of the previous integral. For the subgrid case, ˜ Kiis the union of several volumes of the 2D mesh, then the integral is easily computed. In the mismatching case, ˜ Ki is not exactly defined by the union of several volumes Kj, then it is necessary to estimate the corresponding intersection areas. The integration in time in the denition Ri(tm, v) is due to the fact that the time steps may be different for the 1D and the 2D solvers. The resulting scheme verifies the same well-balanced properties than for the previous cases. The proof is straightforward since each integrand vanishes for water at rest (with or without overflowing), as it has been seen in proof of Theorem 4.1. Then Ψn i= 0, and the well-balanced property is obtained. 5 Numerical results We perform two types of numerical tests. Test 1 is a numerical verification of the well-balanced property when water is at rest. Test 2 concerns a steady-state flow (in a non-constant topography) involving incoming and outgoing lateral fluxes. The 1D solver is the HLL scheme, see (15), while the 2D solver is the corresponding HLLC scheme, see (16). These solvers are implemented into our software DassFlow [11,10]. The coupling algorithm used is a global in time Schwarz-like algorithm, see Fig. A.9. Inputs of the coupling algorithm are: initial conditions for both 1D and 2D models, boundary conditions of the 1D model and a first guess of boundary conditions for the 2D model. 18 5.1 Test 1. Water at rest The first test case considered corresponds to a water at rest solution with and without overflowing. All tests done are successfull since water stays at rest (norm of velocities are about 10−8). The used mesh is presented in Fig. A.10 (2274 cells). 5.2 Test 2. Steady state solution with incoming - outgoing flow 5.2.1 Validation procedure The definition of a non trivial steady state analytical solution for the coupled model is a difficult task. We decided to use the following approach to validate our coupling process: •First, we compute a ”reference” steady flow using the validated 2D HLLC solver for the entire computational domain on a fine grid. •Then, we use our coupling algorithm and compare the solution with the previous generated ”reference” solution. For the 1D model, we use either the same fine grid (matching grids case) or a coarse grid (subgrids case). 5.2.2 Test case description This test is presented in Fig. A.10. The channel length is 200m, the main channel is 2m wide and the circle diameter 100m. The 2D mesh is made of triangular cells, except in the main channel where cells are rectangular. The bathymetry is the following: zb(x, y) =      0.5 if 80 ≤x≤120 0 otherwise Boundary conditions are the following: •Inflow (West) : discharge imposed qin = 2m3/s, •Outflow (East) : ∂h ∂n=∂q ∂n= 0, •Walls boundary conditions are imposed on the other boundaries. This test case includes a non-constant bathymetry, incoming lateral fluxes and outgoing ones. Thus, it is suitable to show the robustness of our definition of the coupling term Ψ. 19 Using a time step ∆t= 0.1 (CFL number is about 0.1), we obtain the steady state presented in Fig. A.10 b). For the coupling test, the simulation time is T= 50s 5.2.3 Matching grids case The 2D ”reference mesh” (see Fig. A.10) is made of 2274 cells. Both 1D and ”2D zoom” meshes are restriction of the ”reference mesh” on the corresponding area. Validation of the source term. We validate the expression of the coupling source term (in the 1D equations) and its implementation as follows. •Initial conditions for both the 1D and 2D models are extracted from the 2D ”reference” solution. •Boundary conditions of the 2D model are extracted from the 2D reference solution. Therefore, the 2D zoom model must reproduce exactly the solution of the 2D reference model (in one iteration of Schwarz algorithm). Furthermore if Ψ, the coupling source term in 1D equations, is accurate, then the 1D model solution should fit perfectly with the reference solution. The numerical results show that both the 1D solution and the 2D zoom solution match perfectly with the reference solution (differences in percent are about 10−6). Thus, this test case shows the accuracy of the coupling source term Ψ and validates its implementation. Coupling with Schwarz algorithm. Now, we consider the coupling algorithm of Schwarz described in Fig. A.9. The convergence threshold used is ǫSchwarz = 10−11. We have a-priori no guarantee that the coupling algorithm will converge to the 2D reference solution. It depends partially on the first guess for the 2D incoming characteristics (thus on the initial condition). Also, in the present test case, the flow computed is steady-state and the iterates are unsteady. In our numerical tests, starting from a ”reasonable initial condition”, the algorithm converged with 4 iterates at maximum. We present in Fig. A.11 the three solutions: the reference one, the 1D solution and the 2D local zoom solution. The latters match perfectly with the ”reference” solution. These last two numerical tests show that in case of matching grids, a convergent coupling procedure based on the coupling source term Ψ leads to a global solution as accurate as a full 2D solution. Concerning the CPU times, let us remark that in the present simple configuration 20 the faster approach is obviously to compute a full 2D solution (since the 1D part is negligible) and not the present superposition strategy. Roughly, the coupling algorithm which can require up to p= 4 iterates (depending on the first guess) is more expansive (about ptimes) than a full 2D simulation. In an operational context, difficulties and time costs are different. As mentionned in the introduction, the 1D model can be a complex network (with empirical laws calibrated by hand etc), its decomposition can be complex (1D data are not the same than 2D ones), the human time to decompose the 1D model can be important etc. Thus, in some operational contexts, the present superposition strategy has few advantages, including in a CPU time point of view. 5.2.4 Subgrids case Now, the reference mesh is much finer, it has 10 298 cells. The 2D zoom mesh is a restriction of this reference mesh. For the 1D model, we consider two coarser meshes: submesh either of a ratio Rspac =∆x1D ∆x2D= 2 or of a ratio Rspac = 10. Concerning time discretization, we keep the same time step for both models: Rtemp =∆t1D ∆t2D= 1. Error due to the mismatching grids. As previously, as a first step, we proceed as follows: •Initial conditions for both 1D and 2D models are extracted from the 2D ”reference” solution, •The boundary conditions of the 2D model is extracted from the 2D ”reference” solution. Therefore, we perform one iteration only of the Schwarz algorithm, the 2D zoom model reproduces the reference solution while the coupling source term in 1D equations include errors (with respect to the reference solution) due to the mismatching grids. In other words, errors are due to the integration of the source term over the 1D cells. Results obtained are presented in Fig. A.12. Coupling with Schwarz algorithm. Starting from a ”reasonable” initial condition, we deduce a first guess for the 2D incoming characteristics and we iterate the algorithm of Schwarz described in Fig. A.9. After convergence, the errors obtained are very similar than the latter, they are errors due to the grid mismatch. Summary of the numerical results (Test 2). These numerical results show that after convergence, the coupling source term Ψ defined in the present article leads to a global solution as accurate as a full 2D solution in case of matching grids, and leads to an accurate solution if grids are mismatching (for a sake of simplicity, only the subgrid has been implemented). This feature remains true with a more complex 21 topography (see the toy test case presented in [12]). 5.3 A practical point of view: superposition of a 2D local zoom model vs full 2D In river hydraulics, the 1D areas (main channels of the river branches considered) can be complex networks. Furthermore, elaborating an operational 1D numerical model is a long and heavy task since data are not dense and precise (topography for example). Thus, in case of flooding, the present principle consisting to superpose locally a 2D model seems to be a good alternative to a re-definition of the whole ”model” (that means to define a new full 2D model). In a discrete point of view, this principle of superposition becomes possible only if coupling term leading to an accurate and well-balanced global scheme can be defined. In this last section, we illustrate both approaches: superposition of a 2D local zoom model vs full 2D with the following constraint: the 1D main channel (mesh and topography) is given and cannot be changed. Then, we compare the following two computations: •”Full 2D”: we keep intact the 1D main channel (mesh and topography), and we extend the mesh inside the flood plain, see Fig. A.13 a). The 2D solver is applied all over the domain. •”Superposition”: Over the 1D model (mesh, topography, solver), we superpose the 2D model, using a finer mesh (those with Rspac = 10) but with the constraint of meshing the lateral boundaries of the channel. The results obtained with the full 2D solver are presented in Fig. A.13 b) and c) (legend ”2D coarse-fine”). A comparison with Fig. A.12 shows that, as expected, we obtain a more precise solution using the principle of superposition. This result is obvious since in the superposition procedure, we have defined a 2D mesh which is finer next to the main channel (Rspac = 10) and similar next to the circle boundary. Nevertheless, this test case illustrates the following practical approach. If for any reason we must keep intact the 1D mesh, the present computations show that the superposition strategy leads to a more accurate solution compared to the full 2D solution based on the ”existing” 1D mesh. 6 Conclusion In this paper, we address the difficulty to couple and superpose numerically 1D and 2D shallow-water equations with non-flat topographies. We focus on the discrete definition of the coupling source term in the 1D equation, in order to obtain a wellbalanced and consistent approximation. To do so, we begin by writing a general form of the discrete the 1D and 2D problems by using a family of well-balanced 22 finite volume solvers. We propose a discretization of the coupling source term that allows to use different finite volume schemes (explicit in time) for the 1D and the 2D problems. This allows, for example, to use a more accurate solver for 2D equations, and a less CPU time-consuming solver for the 1D equations, or the contrary. We prove that the proposed discretization is well-balanced independently of the choice of each solver. We present the problem for curvilinear 1D channel and mismatching grids. Our discretization technique introduces naturally an upwind definition of the tangential velocity at boundaries of the 1D channel. Thus, the global scheme remain stable whatever if lateral fluxes are locally incoming or outgoing. Numerical test cases show the efficiency and robustness of the discretization done, and show that the superposing approach is relevant. 23 A Derivation of the 2D source term in 1D equations We recall the derivation of the coupling source term in the 1D equations from the 3D Navier-Stokes equations. We refer to [12] for more details. We consider the incompressible Navier Stokes equations with no-slip boundary conditions at the bottom, and the kinematic equation at the free surface. We consider a 1D symmetric channel, centered in x-axis. By b(x) we denote the channel width, x∈[0, L]. Lateral boundaries of the channel in the (x, y)-plan are defined by: y=−b(x)/2 and y=b(x)/2 (see Fig. A.2). We define the 1D wet cross section Sas follows: S=Rb2 b1Rzb+h zbdzdy. For the sake of simplicity, we denote: b1(x) = −b(x)/2 and b2(x) = b(x)/2. Mass conservation. Using standard notations, the mass conservation equation writes: ∂u ∂x +∂v ∂y +∂w ∂z = 0. By integration over Swe get: b2 Z b1 zb+h Z zb ∂u ∂x dz dy |{z } A + b2 Z b1 zb+h Z zb ∂v ∂y dz dy |{z } B + b2 Z b1 zb+h Z zb ∂w ∂z dz dy |{z } C = 0 (A.1) We apply the Leibnitz’s integral rule and we use the no slip condition at bottom; this gives: A= b2 Z b1   ∂ ∂x zb+h Z zb u dz −uS ∂(zb+h) ∂x   dy =∂ ∂x b2 Z b1 zb+h Z zb u dz dy −   zb+h Z zb u dz  y=b2 ∂b2 ∂x +   zb+h Z zb u dz  y=b1 ∂b1 ∂x − b2 Z b1 uS ∂(zb+h) ∂x dy |{z } CL1 B= b2 Z b1   ∂ ∂y zb+h Z zb v dz −vS ∂(zb+h) ∂y   dy =   zb+h Z zb v dz  y=b2 −   zb+h Z zb v dz  y=b1 − b2 Z b1 vS ∂(zb+h) ∂y dy |{z } CL2 C= b2 Z b1 wSdy |{z } CL3 where (uS, vS, wS)Tdenotes the surface velocity. We set: Q=Rb2 b1Rzb+h zbu dz dy. Since the free surface boundary condition gives: CL1 + CL2 + CL3 = Rb2 b1 ∂(zb+h) ∂t dy = 24 ∂tS, we obtain: (A.1) ⇐⇒ ∂S ∂t +∂Q ∂x =K1 with: K1=−     zb+h Z zb v dz  b2 −   zb+h Z zb u dz  b2 ∂b2 ∂x −   zb+h Z zb v dz  b1 +   zb+h Z zb u dz  b1 ∂b1 ∂x    Without overflowing, these two terms vanish since h= 0 at lateral boundaries (standard 1D shallow-water equations). If 2D flow information is available at boundaries, hRzb+h zbu dzibk and hRzb+h zbv dzibk , k= 1,2, represent lineic discharges at boundaries in x-direction and y-direction respectively. We denote by η1and η2the unit external normal vectors to the boundaries y= b1(x) = −b(x)/2 and y=b2(x) = b(x)/2 respectively (see Fig. A.2), that is: η1=1 δb (∂xb1,−1)Tand η2=1 δb (−∂xb2,1)T,with δb=v u u t1 + ∂xb 2!2 . If the channel width variation is small then: δb≈1. We set: δb= 1. We set: qηi=   zb+h Z zb u dz  bi (ηi)1+   zb+h Z zb v dz  bi (ηi)2, i = 1,2.(A.2) Then we obtain: K1=−(qη1+qη2) and the mass conservation equation is: ∂S ∂t +∂Q ∂x =−(qη1+qη2) (A.3) The values qηi, i = 1,2,must be provided by a 2D model. Momentum conservation. As previously, we integrate the 3D momentum equations over the 1D wet cross section Sand we use Leibniz integration rule. We consider the following hydrostatic pressure: p=−ρg(z−(H+Zb)) where (H+Zb) is the mean transverse water elevation (Zbis the mean transverse topography i.e. it is the 25 x Ki Kj1 ηi j1 ηi j1 Kj2 i j2 η ηi j2 Fig. A.6. Simplest case and matching grids. j2 Real domain Analytical 1D channel Ki Kj1 K Fig. A.7. General domain with matching grids. Fig. A.8. Information transfer. 1D main channel Ω1, 2D flood plain Ω2, overflowing boundaries Γ1,Γ2 32 2D incoming charac. Coupling var. : 2D charac interpolated from 1D charac IF > εSchwarz o Averaged 2D charac. 1D charac. at interface 2D Model 0 −−> T 0 −−> T 1D Model 1D source term Initial conditions Discrepancy between: Convergence test (not converged) Boundary conditions Fig. A.9. Coupling algorithm based on a Schwarz method. a) b) Fig. A.10. Test case: a) Bathymetry and mesh of 2274 cells, b) ”Reference” steady state solution (velocity field). 33 0.25 0.3 0.35 0.4 0.45 0.5 0.55 0.6 0.65 0 20 40 60 80 100 120 140 160 180 200 reference 2D 1D -3e-06 -2e-06 -1e-06 0 1e-06 2e-06 3e-06 4e-06 5e-06 0 20 40 60 80 100 120 140 160 180 200 2D 1D a) b) Fig. A.11. Matching grids case. Comparison of velocity values (u) in the 1D main channel (common area) and after Schwarz algorithm convergence. a) 2D reference solution and values computed by the 1D model and by the local 2D zoom model; b) Differences in percent. 34 0.25 0.3 0.35 0.4 0.45 0.5 0.55 0.6 0 20 40 60 80 100 120 140 160 180 200 reference 2D 1D 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 0 20 40 60 80 100 120 140 160 180 200 2D 1D a) b) 0.25 0.3 0.35 0.4 0.45 0.5 0.55 0.6 0 20 40 60 80 100 120 140 160 180 200 reference 2D 1D 0 2 4 6 8 10 12 0 20 40 60 80 100 120 140 160 180 200 2D 1D c) d) Fig. A.12. Subgrids case. Comparison of velocity values (u) in the 1D main channel (common area) and after Schwarz algorithm convergence. a) Rspac = 2. 2D reference solution and values computed by the 1D model and by the local 2D zoom model; b) Differences in percent. c) Rspac = 10. 2D reference solution and values computed by the 1D model and by the local 2D zoom model; d) Differences in percent. 35 a) 0.25 0.3 0.35 0.4 0.45 0.5 0.55 0.6 0 20 40 60 80 100 120 140 160 180 200 reference 2D coarse-fine 0 5 10 15 20 25 0 20 40 60 80 100 120 140 160 180 200 2D coarse-fine b) c) Fig. A.13. Superposition vs full 2D. a) The full 2D mesh is defined from the 1D mesh in the channel (R= 1). b) Velocity uin the main channel: ”full 2D” solution (R= 1, legend ”2D coarse-fine”) and the 2D coupled solution with Rspac = 10 (legend ”reference”). c) Differences in percent. 36