On a shallow water model for the simulation of turbidity currents

On a shallow water model for the simulation of turbidity currents T. Morales de Luna1, M. J. Castro D´ıaz2, C. Par´es Madro˜nal3, E.D. Fern´andez Nieto4 1 , 2 , 3An´alisis Matem´atico Facultad de Ciencias - Univ. M´alaga Campus de Teatinos, s/n 29071-M´alaga, Spain 4Dpto. Matem´atica Aplicada I ETS Arquitectura - Univ. Sevilla Avda. Reina Mercedes N. 2 41012-Sevilla, Spain E-mail: 1[email protected], 2[email protected], 3[email protected], 4[email protected] Abstract We present a model for hyperpycnal plumes or turbidity currents that takes into account the interaction between the turbidity current and the bottom, considering deposition and erosion effects as well as solid transport of particles at the bed load due to the current. Water entrainment from the ambient water in which the turbidity current plunges is also considered. Motion of ambient water is neglected and the rigid lid assumption is considered. The model is obtained as a depth-average system of equations under the shallow water hypothesis describing the balance of fluid mass, sediment mass and mean flow. The character of the system is analyzed and numerical simulations are carried out using finite volume schemes and path-conservative Roe schemes. Key words: Turbidity currents, Hyperbolic systems, Finite Volume Methods, Path-conservative schemes, Numerical modeling AMS subject classifications: 74S10, 35L60, 35L65, 74G15 1 Introduction When a river contains an elevated concentration of suspended sediment, to the extent that the river density is greater than that of the receiving water body, the river can plunge and create a hyperpycnal plume or turbidity current. This hyperpycnal plume can travel significant distances until it loses its identity by entraining surrounding ambient water and dropping its sediment load. A sketch of a turbid underflow is presented in Figure 1. There is great interest in turbidity currents because of their profound impact on the morphology of the continental shelves and ocean basins of the world. It is commonly accepted that they are one of the potential processes through which sediments can be transferred to the deep sea environments. These bottom currents influence the sea bed morphology by depositing, eroding and dispersing large quantities of sediment particles. The resultant deposit often form porous layer of rocks which are potential sources of hydrocarbon. Therefore, understanding and predicting the geometry of these deposits is crucial for effectively exploring and exploiting these reservoirs. An additional concern is the destructive effect that turbidity currents have on underwater structures, such as cables, pipelines and foundations. Large-scale hyperpycnal flow or turbidity currents in the natural environment are difficult to monitor because of the unpredictable nature of the events. As a result, most of our knowledge about these flows is derived from small scale laboratory experiments like the ones described in [19], [11],[1], [17] and [16]. In this paper, we present a numerical model of hyperpycnal flow generated by the plunging of a river. It incorporates the interaction between the turbidity current and bottom, considering eroding and deposit effects as well as solid transport due to the velocity of the current. Here, the along-shore 1 H(t, x) h(t, x) Hyperpycnal flow Water entrainement b(t, x) and deposition ρ=Pns j=0 ρjcj ρw Ambient fluid Erodible bed H0−h(t, x)−b(t, x) Free surface Sediment erosion u(t, x) Figure 1: Sketch of hyperpycnal flow current has been neglected. The equations are obtained as a layer-averaged system from the twodimensional equations under the shallow water hypothesis, namely the hydrostatic distribution of pressure. A number of studies have attempted to predict the distribution of turbidities, mainly based on smallscale tank experiments of particle-driven density currents (see [8], [23], [15], [17], [16]). Numerical models have been developed as well, and some of them applied to the tank experiments (e.g. [3], [5], [14]). The organization of this paper is as follows: In Section 2 we describe the model presented here and in particular the empirical relations and physical constants taken into account. In Section 3 the model is obtained from the two dimensional equations under the shallow water hypothesis and hydrostatic distribution of pressure. The character of the system as well as its eigenstructure is studied in Section 4. An entropy inequality is also obtain under some hypothesis. The theory presented in [22] and [21] is adapted in Section 5 in order to define a Roe scheme for the model considered here. Finally, numerical results and comparisons to experimental data are shown in Section 7. 2 Model description We introduce a layer-averaged model in which equations for conservation of fluid mass, sediment mass and fluid momentum are solved simultaneously. Vast applicability, low computational cost and conceptual simplicity are the main reasons that justify the option for layer averaged model in a onedimensional framework. This type of model has been widely used in the simulation of turbidity currents (e.g., [8], [23], [5], [17]) and some layer-averaged models have already been proposed in [23], [16], [17], [5] and [14]. In the framework proposed by these authors, ns≥1 species of sediments are considered with constant density ρj,for j= 1, . . . , ns. These sediment species are transported by a river with freshwater of constant density ρ0. The river plunges into an ambient fluid (in general the ocean) of density ρwgenerating a turbidity current. The models presented by the cited authors may be written under the following general formulation: 2                ∂th+∂x(hu) = φη, ∂t(hu) + ∂xhu2+gRc h2 2=g(R0+Rc)h∂xH+τ, ∂t(hcj) + ∂x(hu cj) = φj b,for j= 1, . . . , ns ∂tH=ξφb. (2.1) where his the current thickness; uis the depth-averaged velocity; cjfor j= 1, . . . , nsrepresents the vertically averaged volume concentration of the jth sediment; and His the bottom depth from a fixed level. Finally, Rj=ρj−ρ0 ρ0 ,for j= 1, . . . , ns; and Rc= ns X j=1 Rjcj.(2.2) The source term φηrepresents the amount of ocean water entrained and mixed by turbulence to the plume. The sediment flux at the bed of the jth specie is determined from the rates of deposition (Fj d) and erosion (Fj e), φj b=Fj e−Fj d, φb= ns X j=1 φj b.(2.3) We denote ξ= 1/(1 −γ),being γthe porosity of the sediment layer. Finally, τrepresents the friction term. While there are some discrepancies in the physical laws that describe the source terms φηand φb, we shall use here the expressions that are the most extended. In particular, following [23], [17] and [14], the amount of ocean water entrained and mixed by turbulence to the plume, φη, may be described as Ewu, where the water entrainment coefficient, Ew, is determined from the empirical relationship Ew=0.00153 0.0204 + Ri ,(2.4) where the Richardson number, Ri is defined as Ri =Rcgh u2.(2.5) The rate of deposition is described as the product of the settling velocity of sediment, vsj, and the fractional concentration of suspension near by the bed, cbj. For multiple grain sizes, the sum of this product for each size population is used as the net rate of deposition Fj d=vsjcbj.(2.6) This definition is found for example in [23], [5] and [14] while other formulation is presented in [16]. The near bed concentration of sediment, cbj, can be related to the layer averaged concentration, cj, by a factor 2, as it is the case in [12] and [2], or by the more extended expression used in [14] and [5]. cbj cj = 0.4Dj Dsg 1.64 + 1.64,(2.7) where Djis the characteristic grain size and Dsg denotes the geometric mean size of the suspended sediment mixture. The rate of sediment entrainment from the bed is described as Fj e=vsjpjEsj,(2.8) where the sediment entrainment coefficient, Esj, is given by the expression developed by Garcia and Parker [11], Esj=1.3·10−7Z5 j 1+4.3·10−7Z5 j ,(2.9) 3 where Zj=α1 √cD|u| vsj Rα2 pj,(2.10) and Rpj=pRjgDjDj ν(2.11) is the particle Reynolds number, with νthe kinematic viscosity. The parameters (α1, α2) take respective values (1,0.6) for Rpj>2.36 and (0.586,1.23) for Rpj≤2.36. We remark again that this is not the only expression used in models describing the physics of the problem. Some variants can be found in [5], [17] and [16]. In (2.8), pjrepresents the volume fraction of the jth sediment in the bed. It could be described as a function of time as it is the case in [5]. The friction term τis described as a function of the velocity, τ=−(1 + α)cD|u|u. (2.12) cDis the bed drag coefficient which ranges from 0.002 to 0.05 depending on the flow type and αis the ratio of the drag force at the upper flow surface to that at the bed. Finally, the fall velocity vsjis calculated using the empirical relationship introduced by Dietrich [10]. We remark that there are some difficulties with system (2.1). First, let us denote by c0the freshwater concentration present in the plume and by cthe summation of all of the sediment fractions. The relation c0= 1 −cis satisfied and from (2.1) one gets ∂t(c0h) + ∂x(c0hu) = φη−φb.(2.13) This means that erosion and deposition effects modify the freshwater mass while it is clear that in the case that there is no water entrainment from above, φη= 0, freshwater mass should be preserved. Another point is that deposition/erosion and water entrainment are not taken into account for the momentum equation. We will see in Section 3 that source terms related to deposition/erosion and water entrainment appear naturally in the momentum equation. We would also like to take into account the ambient fluid considered and in particular its density ρw, so that we could model a plume that plunges into the ocean as well as a current flowing over a given bottom. Finally, we would like to include into the model the sediment transport of bed-load particles due to the velocity of the current. Taking into account these considerations, we are going to introduce the following model for hyperpycnal plumes:                ∂th+∂x(hu) = φη+φb, ∂t(hu) + ∂xhu2+g(R0+Rc)h2 2=g(R0+Rc)h∂xH+uφη+u 2φb+τ, ∂t(hcj) + ∂x(hu cj) = φj b,for j= 1, . . . , ns ∂tH−ξ∂xqb=ξφb. (2.14) Here, R0=ρ0−ρw ρ0 .(2.15) Remark 2.1.In the general case we consider a river with a mixture of sediments that plunges into the ocean: ρw(saltwater density) is slightly greater than ρ0(freshwater density) and R0<0.The case where ambient fluid is air can also be considered: in that case ρw≈0and R0≈1. Remark 2.2.The system (2.17) presents some major differences with respect to the model presented in [23],[5], and [14]. 4 •From (2.17) one gets ∂t(c0h) + ∂x(c0hu) = φη,(2.16) with c0= 1 −cthe freshwater concentration. Thus, if there is no ocean water entrainment, the total freshwater mass is preserved. •The pressure term gR0h2/2in the second equation allows to consider a turbid current that plunges into the ocean (R0≈0) as well as currents flowing over a given bottom (R0≈1). •The effects of deposition/erosion and water entrainment on the momentum equation are taken into account. •Solid transport of particles at the bottom is allowed via a solid transport flux qb. Taking into account the empirical relations described before, (2.14) reads                                  ∂th+∂x(hu) = Ewu+ ns X j=1 vsj(pjEsj−cbj), ∂t(hu) + ∂xhu2+g(R0+Rc)h2 2=g(R0+Rc)h∂xH +Ewu2+u 2 ns X j=1 vsj(pjEsj−cbj)−(1 + α)cD|u|u, ∂t(hcj) + ∂x(hu cj) = vsj(pjEsj−cbj), ∂tH−ξ∂xqb=ξ ns X j=1 vsj(pjEsj−cbj). (2.17) Sediment transport of bed-load particles due to the velocity of the current is represented by ξqbwhere ξ= 1/(1 −γ),being γthe porosity of the sediment layer and qbthe solid transport flux, that depends on the fluid velocity u,qb=qb(h, hu). Several formulae for solid transport flux have been proposed, see [13], [18], [20], etc. Here, we shall consider the simplest one, which is the Grass model. According to it, the solid transport fluid equation is qb=Agu|u|mg−1,1≤mg≤4,(2.18) where Agis a constant that represents the effects due to grain size and kinematic viscosity and is usually determined by experimental data. Usually, the constant mgis set to mg= 3. One could also consider other expression for the solid transport flux and most of what is going to be said remains valid. The election of Grass model will only play an essential role in Lemma 4.3, where the explicit expresion of the corresponding eigenvalues can be given, as well as in the explicit definition of a Roe matrix in Section 5. 3 Derivation of the layer-averaged equations 3.1 Governing equations for a dilute suspension We consider the situation described in Section 1 and, in particular, in Figure 1. We shall assume the rigid lid hypothesis for the ambient fluid in order to simulate just the turbidity current. The width of the plume is assumed to be small compared to the typical length phenomena in the horizontal direction so that the shallow water hypothesis can be made. The equations of motion for a suspension are considered. The suspension is assumed to be sufficiently dilute to justify the use of the Boussinesq approximation, and the assumption of a kinematic viscosity νequal to the value for clear water. The equation of momentum conservation is ∂~u ∂t + div(~u ⊗~u) = −1 ρ0∇p+1 ρ0∇·T+ (1 + Rc)~g, (3.1) 5 where ~u(~x, t) represents the instantaneous velocity of the mixture, p(~x, t) denotes the instantaneous pressure and ~g = (0,0,−g) is the gravity and T=µ(∇~u + (∇~u)t) (3.2) is the stress tensor. We denote c=Pns j=1 cj,R0= 1 −ρw/ρ0,and Rc=Pns j=1 Rjcj,where Rj= ρj/ρ0−1 and cj(x, t),for j= 1, . . . , nsdenotes the instantaneous volumetric concentration of the jth suspended sediment. The fluid is assumed to be incompressible and the continuity equation is also considered: div~u = 0,(3.3) ∂tρ+ div(ρ~u)=0,(3.4) where ρis the density of the turbidity current. Following the notation of Section 2, we may rewrite the density of the current as ρ= ns X j=0 (ρjcj) = ρ01 + ns X j=1 Rjcj(3.5) where c0(~x, t) denotes the portion of water in the turbidity current. Now, as Pns j=0 cj= 1, from equations (3.3) and (3.4) we get 0 = ∂t(c0ρ0+ ns X j=1 ρjcj) + div((c0ρ0+ ns X j=1 cjρj)~u) = ∂t( ns X j=1 cj(ρj−ρ0)) + div( ns X j=1 cj(ρj−ρ0)~u). (3.6) In order to satisfy (3.6), we shall assume that the volume fraction of each sediment specie satisfies the equation ∂cj ∂t + div(cj~u)=0,for j= 1, . . . , ns.(3.7) As a consequence, this relation is also satisfied by the volume fraction of fresh water and the equations of mass conservation are thus ∂cj ∂t + div(cj~u)=0,for j= 0,1, . . . , ns.(3.8) 3.2 Vertically integrated equations Henceforth ~x = (x, y, z) and ~u = (u, v, w). Variations in the horizontal ydirection are neglected and we suppose v= 0. It is assumed that the thickness of the turbidity current is small compared to the typical length phenomena in the horizontal direction so that the shallow water hypothesis can be made. As the current is essentially two-dimensional, the domain considered is (t, x, z)∈R+×R2/ b(x, t)< z < b(x, t) + h(x, t),(3.9) where b(x, t) is the bottom elevation from a defined reference level and h(x, t) is the length of the turbidity layer. As it has been said, perturbations of the free surface due to the turbidity current are neglected, so that the water surface is considered to be constant and equal to z=H0. It will be useful to define the bottom depth from the surface H(x, t) = H0−b(x, t),(3.10) 6 and the interface between the turbidity current and the ambient water η(x, t) = b(x, t) + h(x, t).(3.11) Some kinematic conditions should be imposed at the bottom and at the interface. The turbidity current lies between an upper non-material interface, Fs=z−η(x, t) = 0, and a basal interface, Fb=b(x, t)−z= 0. These curves are defined so that their unit normals ns=∇Fs/|∇Fs|and nb=∇Fb/|∇Fb|point outwards from the turbidity current. Assuming that these interfaces have velocities ~vsand ~vb,the kinematic boundary conditions are Fs(x, z, t) = 0 : ∂tFs+~vs·∇Fs= 0,(3.12) Fb(x, z, t) = 0 : ∂tFb+~vb·∇Fb= 0.(3.13) The upper and basal interfaces are modified by entrainment of water from the above and deposition/erosion of sediment at the bottom. Let the normal entrainment rate, ds, be the equivalent volume of ambient water entrained to the turbid current at the upper interface per unit area per unit time. Then the velocity of the interface is ~vs=~us+ds~ ns. Similarly, if the normal deposition/erosion rate, db, is the equivalent volume of granular material deposited/eroded at the bottom per unit area per unit time, the velocity of the basal interface is ~vb=~ub+db~nb. It follows that the kinematic conditions (3.12)-(3.13) are z=η(x, t) : ∂tη+us∂xη−ws= (1 + (∂xη)2)1/2ds,(3.14) z=b(x, t) : ∂tb+ub∂xb−wb=−(1 + (∂xb)2)1/2db.(3.15) Thus, we shall consider the kinematic equations ∂tη+u|z=η∂xη−w|z=η=e φη,(3.16) ∂tb+u|z=b∂xb−w|z=b=−e φb,(3.17) where ~u = (u, w),e φη= (1 + (∂xη)2)1/2ds,and e φb= (1 + (∂xb)2)1/2db. Under these hypothesis, mass conservation equations yield ∂tZη b cjdz +∂xZη b cjudz =cj|z=ηe φη+cj|z=be φb,for j= 0,1, . . . , ns.(3.18) We shall only retain the vertical components of the stress tensor and assume that the other terms may be neglected. Let us denote by eτ=µ ρ0∂u ∂z +∂w ∂x (3.19) the (1,2) component of the stress tensor T.Let us also suppose that the pressure preduces to hydrostatic effects, p(x, z, t) = Zη z ns X j=0 (ρjcj)gdz0+gρw(H0−h−b).(3.20) By integrating the momentum equation, we obtain: ∂tZη b udz +∂xZη b u2dz +g ρ0Zη b ∂xpdz −gρw ρ0 h∂xη=eτ|z=η−eτ|z=b+u|z=ηe φη+u|z=be φb.(3.21) Definition 3.1. For any function f(x, z, t)we define the mean of the function as f=1 hZη b fdz, (3.22) The fluctuating part of frelative to the mean is defined by b f=f−f. (3.23) 7 Using this notation, it is clear now that for any functions f1, f2one has hf1f2+Zη bb f1b f2dz =Zη b f1f2dz. (3.24) Thus, equation (3.18) writes as follows: ∂t(hcj) + ∂x(hcju) + ∂xZη bbcjbujdz=cj|z=ηe φη+cj|z=be φb,for j= 0,1, . . . , ns.(3.25) From equation (3.20) we get p(x, z, t) = g ns X j=0 (ρjcj(η−z)) + gρw(H0−h−b) + gZη z ns X j=0 (ρjbcj)dz0,(3.26) and Zη b 1 ρ0 ∂xpdz −gZη b ∂x  Zη z ns X j=0 (ρjbcj)dz0  dz =g ρ0Zη b ∂x  ns X j=0 (ρjcj)(η−z) dz −gρw ρ0 h∂x(h+b) =g ρ0 ∂x Zη b ns X j=0 (ρjcj(η−z))dz +g ρ0 ns X j=0 ρjcjh∂xb−gρw ρ0 h∂x(h+b) =g∂x R0+ ns X j=1 Rjcj h2 2 +g R0+ ns X j=1 Rjcj h∂xb. (3.27) (3.21) can be thus written as follows: ∂t(hu) + ∂x hu2+g R0+ ns X j=1 Rjcj h2 2 +g R0+ ns X j=1 Rjcj h∂xb +Zη bbu2dz +gZη b ∂xZη z ns X j=0 (ρjbcj)dz0dz =eτ|z=η−eτ|z=b+u|z=ηe φη+u|z=be φb. (3.28) 3.3 Similarity assumptions In this point, it is necessary to make some assumptions in order to write the integrated equations in terms of mean variables. We will suppose that the variables cjmaintain approximately similar profiles in the z-direction so that the terms bcjmay be neglected. The term Rη bbu2dz corresponds to the kinematic Reynolds stress. We shall assume that this term and the tensor stress may be expressed in terms of uso that we write eτ|z=η−eτ|z=b−Zη bbu2dz =τ(u).(3.29) Let us consider the approximations φj η≃cj|z=ηe φη, φj b≃cj|z=be φb, j = 0,1, . . . , ns,(3.30) and φη= ns X j=0 φj η≃e φη, φb= ns X j=0 φj b≃e φb.(3.31) 8 In particular, as there is no water entrainment at the bottom, we shall assume φ0 b= 0,and φj bfor j= 1, . . . , nswill be given by the empirical relations described in Section 2. In the same way, there is no sediment transfer at the interface z=ηand we shall suppose φj η= 0,for j= 1, . . . , nsand φ0 η≡φηdescribed in Section 2. Following [4], we shall assume u|z=η=uand u|z=b=1 2u. Thus, using the fact that c0= 1 −cand ∂xb=−∂xHwe arrive to the final model      ∂th+∂x(hu) = φη+φb, ∂t(hc) + ∂x(hu c) = φb, ∂t(hu) + ∂xhu2+g(R0+Rc)h2 2=g(R0+Rc)h∂xH+uφη+u 2φb+τ. (3.32) The second equation will be replaced by one equation for each sediment specie, ∂t(hcj) + ∂x(hu cj) = φj b,for j= 1, . . . , ns,(3.33) where Pns j=1 φj b=φb. 3.4 Evolution of the topography Sediment and erosion of the bottom as well as sediment transport due to the velocity of the current are considered here. The bed-sediment conservation equation has the form ∂tb+ξ∂xqb=−ξφb,(3.34) where qb=qb(h, hu) denotes the solid transport flow and ξ=1 1−γ, where γis the porosity of the bed. 4 Reformulation and properties of the model The system may be rewritten in the form: ∂tW+∂xe F(W) = B(W)∂xW+S(W),(4.1) where W= (h, hu, hc1, . . . , hcns, H)t,e F(W)=(hu, hu2+g/2(R0+Rc)h2, hc1u, . . . , hcnsu, −ξqb)t, B(W) = g(R0+Rc)he2,ns+3 and S(W) =         Ewu+Pns j=1 vsj(pjEsj−cbj) Ewu2+u 2Pns j=1 vsj(pjEsj−cbj)−(1 + α)cD|u|u vs1(p1Es1−cb1) . . . vsns(pnsEsns−cbns) ξPns j=1 vsj(pjEsj−cbj)         (4.2) Here er,s = (δirδjs)ns+3,ns+3 i=1,j=1 is the canonical basis of square matrices of order ns+ 3. We will also use the notation U= (h, hu, hc1, . . . , hcns)t, F(U)=(hu, hu2+g/2(R0+Rc)h2, hc1u, . . . , hcnsu)t. The problem may also be written in the more general form ∂tW+A(W)∂xW=S(W),(4.3) where A=∇e F(W)−B(W). More explicitly, A(W) =           0 1 0 0 . . . 0 0 gh(R0+1 2Rc)−u22ug 2R1hg 2R2h . . . g 2Rnsh−g(R0+Rc)h −c1u c1u0. . . 0 0 −c2u c20u . . . 0 0 . . . . . . . . . . . . . . . . . . . . . −cnsu cns0 0 . . . u . . . −ξ∂qb ∂h −ξ∂qb ∂(hu)0 0 . . . 0 0           (4.4) 9 Dal Masso, LeFloch, and Murat [9], a rigorous definition of weak solutions can be performed using a family of paths. We refer to [9], [22] and [21] for further details. In [24] a generalization of Roe methods to systems of the form (5.4) was introduced. These methods are based on the following general definition of a Roe linearization: Definition 5.1. Given a family of paths Ψ, a matrix function AΨ: Ω ×Ω→ MN(R)is called a Roe linearization if it satisfies: •for any WL, WR∈Ω,AΨ(WL, WR)has Nreal eigenvalues; • AΨ(W, W) = A(W), for all W∈Ω; •for any WL, WR∈Ω: AΨ(WL, WR)·(WR−WL) = Z1 0AΨ(s;WL, WR)∂Ψ ∂s (s;WL, WR)ds. (5.6) A family of paths is a Lipschitz continuous function Ψ : [0,1] ×Ω×Ω→Ω which must satisfy some regularity conditions, (see [9] for the details) and Ψ(0, Wl, Wr) = Wl,Ψ(1, Wl, Wr) = Wr), Ψ(s, W, W) = W. Once a Roe linearization is given, Roe schemes can be generalized for systems of the form (5.4). We refer to [22] for further details. The scheme can be written in the form Wn+1 i=Wn i−∆t ∆xA+ i−1/2·(Wn i−Wn i−1) + A− i+1/2·(Wn i+1 −Wn i),(5.7) where L± i+1/2=    (λi+1/2 1)±0 ... 0 (λi+1/2 N)±    ,A± i+1/2=Ki+1/2L± i+1/2K−1 i+1/2,(5.8) λi+1/2 1, . . . , λi+1/2 Nbeing the eigenvalues of Ai+1/2and Ki+1/2a matrix whose columns form a basis of associated eigenvectors. We consider here the canonical choice of the family of paths Ψ(s;WL, WR) = WL+s(WR−WL).(5.9) In order to build a Roe linearization based on the family of paths (5.9), it is natural to look for matrices with the same structure of A(W): Ai+1/2=Ji+1/2−bi+1/2 −ξqi+1/20.(5.10) For Roe matrices with this structure, (5.6) can be rewritten as follows: Ji+1/2(Ui+1 −Ui)−(Hi+1 −Hi)bi+1/2=F(Ui+1)−F(Ui)−bΨ(Wi, Wi+1) (5.11) qi+1/2(Ui+1 −Ui) = qb(Ui+1)−qb(Ui),(5.12) where bΨ(Wi, Wi+1) = (0,Z1 0 g R0Ψ1(s, Wi, Wi+1) + ns X j=1 RjΨj+2(s, Wi, Wi+1) ∂Ψ2 ∂s ds, 0,...,0)t. (5.13) 16 In the particular case that qbis given by Grass model, properties (5.11)-(5.12) are satisfied if we define: Ji+1/2=       0 1 0 . . . 0 g(R0hi+1/2+1 2γi+1/2)−(ui+1/2)22ui+1/2g 2R1hi+1/2. . . g 2Rnshi+1/2 −(c1)i+1/2ui+1/2(c1)i+1/2ui+1/2 ...0 . . . . . . . . . . . . . . . −(cns)i+1/2(cns)i+1/20. . . ui+1/2        ,(5.14) −ξqi+1/2=ui+1/2di+1/2−di+1/20. . . 0,(5.15) bi+1/2=      0 g(R0hi+1/2+γi+1/2) 0 . . . 0      ,(5.16) where hi+1/2=hi+hi+1 2, ui+1/2=√hiui+phi+1ui+1 √hi+phi+1 ,(5.17) (cj)i+1/2=√hi(cj)i+phi+1(cj)i+1 √hi+phi+1 , γi+1/2=Pns j=1 Rj((hcj)i+ (hcj)i+1) 2,(5.18) di+1/2=Agξ(√hi+phi+1) √hihi+1 +phi+1hi mg−1 X k=0 (ui+1)k(ui)mg−(k+1).(5.19) The eigenstructure of matrices (5.10) is given by lemma 4.2, and thus Ai+1/2is a Roe matrix. Remark 5.1.If we neglect the solid transport flow (qb= 0), the system reduces to a system of conservation laws with a source term. Thus, following [22], we have A± i+1/2=J± i+1/2−J± i+1/2J−1 i+1/2bi+1/2 0 0 (5.20) and the scheme (5.7) may be rewritten under the form Un+1 i=Un i−4t 4xhP+ i−1/2Ji−1/2(Ui−Ui−1)−bi−1/2(Hi−Hi−1) +P− i+1/2Ji+1/2(Ui+1 −Ui)−bi+1/2(Hi+1 −Hi)i, (5.21) where P± i+1/2=1 2Id ±|Ji+1/2|J−1 i+1/2.(5.22) 5.2 The non-hyperbolic case From Lemma 4.2, if R0hi+1/2+1 2γi+1/2+1 2Pns j=1(Rjcj)i+1/2hi+1/2>0 and di+1/2is sufficiently small, Ai+1/2has ns+ 3 real eigenvalues that can be approached by λ1≃ui+1/2−rg(R0hi+1/2+γi+1/2 2+1 2(Rc)i+1/2hi+1/2),(5.23) λ2≃ui+1/2+rg(R0hi+1/2+γi+1/2 2+1 2(Rc)i+1/2hi+1/2),(5.24) λ2+j=ui+1/2, j = 1, . . . , ns(multiplicity ns),(5.25) λns+3 ≃0.(5.26) 17 As R0<0, λ1and λ2may become complex when 1 2γi+1/2+1 2(Rc)i+1/2hi+1/2is small. Indeed, this situation may arise near the plume front or when most of the sediment has been deposited. In this situations the density of the turbidity current is lower than that of the ambient water an thus it should go upwards. Therefore, the model is no longer valid in such situations. Nevertheless, one would like to compute, for example, numerical simulations of dam-break type, where eigenvalues may become complex near the front but the model makes sense far from the front. One way to do so is to set ρ0=ρwif Ai+1/2has complex eigenvalues and to neglect the pressure terms. Thus, whenever Ai+1/2given by (5.10), (5.14)-(5.16) has complex eigenvalues, we replace locally the scheme by an up-winding method for the flux e Fwithout pressure terms. More explicitly, if the Roe matrix has complex eigenvalues at the inter-cell xi0+1/2,then the terms A± i0+1/2(Wi0+1 −Wi0) in the scheme (5.7) are replaced by D± i0+1/2where D− i0+1/2= 0, D+ i0+1/2=e F∗(Wi0+1)−e F∗(Wi0) if ui0+1/2>0,(5.27) D− i0+1/2=e F∗(Wi0+1)−e F∗(Wi0), D+ i0+1/2= 0 if ui0+1/2<0,(5.28) and e F∗(W) = (hu, hu2, hc1u, . . . , hcnsu, −ξqb)t.(5.29) 5.3 Addition of a non-erodible bed Usually, the bottom is formed by a thin layer of sediments that may be eroded or transported by the plume that lays over a rigid bottom which is not erodible. In that case, an initial function M0(x) must be given that establishes the depth for the non-erodible bottom (see Figure 2). Moreover, in some cases a fraction of the deposited sediment may solidify and become part of the non-erodible bottom so that M0could depend on time. In this cases, the source terms (φj b)n=vsj(pn jEn sj−cn bj) have to be limited so that the amount of eroded sediment do not exceed the total mass of the erodible layer. Thus, for each cell Iiwhenever (φj b)n>0,this source term is replaced by the modified expression (b φj b)n=((φj b)n,if pn j(M0−Hn)−ξ4t(φj b)n≥0, pn j(M0−Hn)(ξ4t)−1,otherwise. (5.30) The volume fraction of sediments at the bottom, pn+1 iis given by pn+1 j(M0−Hn+1) = pn j(M0−Hn)−ξ4t(b φj b)n≥0,(5.31) for each cell Ii. Note that (5.30) assures the positivity of pn+1 j. 6 Extension to 2D The model (2.17) can be easily extended to the 2D case. Following the same procedure of Section 3 adapted to the 2D case, one can obtain                        ∂th+∂x(hu) + ∂y(hv) = φη+φb, ∂t(hu) + ∂xhu2+g(R0+Rc)h2 2+∂y(huv) = g(R0+Rc)h∂xH+uφη+u 2φb+τu, ∂t(hv) + ∂x(huv) + ∂yhv2+g(R0+Rc)h2 2=g(R0+Rc)h∂yH+vφη+v 2φb+τv, ∂t(hcj) + ∂x(hu cj) + ∂y(hv cj) = φj b,for j= 1, . . . , ns ∂tH−ξ∂xqx b−ξ∂yqy b=ξφb, (6.1) 18 M0(x) Ambient fluid Free surface Erodible bed Non-erodible bed Hyperpycnal flow H(t, x) h(t, x) Figure 2: Sketch of hyperpycnal flow over a non erodible bed where his the current thickness; ~u = (u, v) is the depth-averaged velocity; cjfor j= 1, . . . , ns represents the vertically averaged volume concentration of the jth sediment; and His the bottom depth from a fixed level. We shall adopt the same notation as in Section 2 adapting the source terms to the 2D framework. So, the expression of φηis replaced by φη=Ewk~uk,(6.2) while the definition of φj bremains essentially the same but replacing |u|by k~ukin (2.10). We shall denote τu=−(1 + α)cDuk~uk, τv=−(1 + α)cDvk~uk.(6.3) (qx b, qy b) is the corresponding 2D solid transport flux which, in the case of Grass model, corresponds to qx b=Aguk~ukmg−1, qy b=Agvk~ukmg−1,1≤mg≤4,(6.4) The system may be rewritten in the form: ∂tW+∂xe F1(W) + ∂ye F2(W) = B1(W)∂xW+B2(W)∂yW+S(W),(6.5) where W= (h, hu, hv, hc1, . . . , hcns, H)tand e F1(W)=(hu, hu2+g/2(R0+Rc)h2, huv, hc1u, . . . , hcnsu, −ξqx b)t,(6.6) e F2(W)=(hv, huv, hv2+g/2(R0+Rc)h2, hc1v, . . . , hcnsv, −ξqy b)t,(6.7) B1(W) = g(R0+Rc)he2,ns+4,(6.8) B2(W) = g(R0+Rc)he3,ns+4 (6.9) and S(W) =           Ewk~uk+Pns j=1 vsj(pjEsj−cbj) Ewuk~uk+u 2Pns j=1 vsj(pjEsj−cbj)−(1 + α)cDuk~uk Ewvk~uk+v 2Pns j=1 vsj(pjEsj−cbj)−(1 + α)cDvk~uk vs1(p1Es1−cb1) . . . vsns(pnsEsns−cbns) ξPns j=1 vsj(pjEsj−cbj)           (6.10) Here er,s = (δirδjs)ns+4,ns+4 i=1,j=1 is the canonical basis of square matrices of order ns+ 4. The problem may also be written in the more general form 19 ∂tW+A1(W)∂xW+A2(W)∂yW=S(W),(6.11) where Ak=∇e Fk(W)−Bk(W), k = 1,2. Some easy calculations allow us to prove the following property: Lemma 6.1. System (6.1) is invariant under rotations. More explicitly, suppose any unitary vector ~n = (nx, ny)and define R~n =    1 0 0 0 0nxny0 0−nynx0 0 0 0 Id    .(6.12) Then the following relations hold e F~n(W) := nxe F1(W) + nye F2(W) = R−1 ~n e F1(R~nW),(6.13) B~n(W) := nxB1(W) + nyB2(W) = R−1 ~n B1(R~nW)R~n,(6.14) A~n(W) := nxA1+nyA2=R−1 ~n A1(R~nW)R~n.(6.15) In order to solve (6.11), the computational domain is divided into discretization cells or finite volumes, Viand we proceed as it was done in Section 5 by solving first ∂tW+∂xe F1(W) + ∂ye F2(W) = B1(W)∂xH+B2(W)∂yH, W(x, y, t =tn) = Wn ifor (x, y)∈Vi.(6.16) Then, let Wn+1/2 ibe the approximations of the cell averages of the solution of this problem. We define Wn+1 i=Wn+1/2 i+4tS(Wn i).(6.17) The following notation is considered: given a finite volume Vi,Niis the set of indexes jsuch that Vj is a neighbor of Vi;Eij is the common edge to two neighbor cells Viand Vj, and |Eij|represents its length; ~nij = (nij,x, nij,y) is the normal unit vector of the edge Eij pointing towards the cell Vj. Following [7], [6], system (6.16) will be solved using a path-conservative Roe scheme which can be written under the form Wn+1 i=Wn i−4t |Vi|X j∈Ni|Eij|A− ij(Wn j−Wn i),(6.18) where Ai,j =A(Wi, Wj,~nij) in now a Roe linearization of nxA1(W) + nyA2(W). We refer to [7] and [6] for further details. Thanks to Lemma 6.1, the definition of Roe matrix is an easy task by using the expression given in the 1D case. The eigenstructure of Roe matrix is also straightforward from Lemma 4.2. 7 Numerical simulations 7.1 Comparison with laboratory experiments First, we compare the model with the data presented in [16] corresponding to a laboratory experiment using a 10-m-long, 0.2-m-wide and 0.5-m-deep flume with a 0.5-m-long gate box at the upstream end. Lock exchange experiments where studied with different bottom topographies and different sediment concentration at the gate box. For the experiment shown here, the initial water depth is 0.2 m at the gate box. The topography in the flume is given by Figure 3. Particle-driven density currents were generated by releasing suspensions, which were composed by siliceous, non-cohesive, sand to silt sized particles (density 2650 kg/m3). The initial volume fraction of the particles was 2% and the solid transport of the bottom is neglected (qb= 0). The sediments 20 00000000000 00000000000 00000000000 00000000000 00000000000 00000000000 11111111111 11111111111 11111111111 11111111111 11111111111 11111111111 0000000000000 0000000000000 0000000000000 1111111111111 1111111111111 1111111111111 0000000000000 0000000000000 0000000000000 1111111111111 1111111111111 1111111111111 000000000000 000000000000 000000000000 111111111111 111111111111 111111111111 length gate ramp length hump length ramp height hump height Figure 3: Comparison with laboratory experiments: Hump topography are represented by a mixture of 20% of fine sand (vs= 0.01m/s,φ= 3.0), 50% very fine sand (vs = 0.005m/s,φ= 3.5) and 30% of coarse silt (vs = 0.003m/s,φ= 4.0). We show the simulations corresponding to experiments B5, B9, B11, C5, C7 and C13 in [16]. The corresponding parameters for each experiment are given in Table 1. Experiment Ramp height Ramp length Number of humps Hump height Hump length B5 0.1 1 0 - - B9 0.1 2 0 - - B11 0.05 1 0 - - C5 0.1 1 3 0.036 1 C7 0 1 3 0.036 1 C13 0.1 1 3 0.036 2 Table 1: Experiment description In general, the model predictions are comparable to the experimental results. In some cases, we remark that some disagreements arise, specially in the proximal area, where the model predicts less deposition than found in experimental results. This behaviour was also observed in the simulations shown in [16]. This is more likely due to the fact that the motion of the upper ambient layer has been neglected. When a turbidity current intrudes shallow water, an opposing current is generated in the ambient fluid, which reduces its forward momentum. In Figures 4, 5 and 6 we see that deposits decrease from a maximum at the beginning of the flume with a slight increase downstream of the slope break. This increase is sensible to the slope of the ramp. In the case of Figures 7, 8 and 9, the topography influence can be seen as an increase on upslope of the humps. 0 1 2 3 4 5 6 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 Distance from gate (m) Deposit mass (g/cm2) 0 1 2 3 4 5 60 0.05 0.1 Bottom elevation Deposit mass (Simulation) Deposit mass (Experimental) Initial topography Figure 4: Comparison with laboratory experiments: Deposit distribution for experiment B5 21 0 1 2 3 4 5 6 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 Distance from gate (m) Deposit mass (g/cm2) 0 1 2 3 4 5 60 0.05 0.1 Bottom elevation Deposit mass (Simulation) Deposit mass (Experimental) Initial topography Figure 5: Comparison with laboratory experiments: Deposit distribution for experiment B9 0 1 2 3 4 5 6 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 Distance from gate (m) Deposit mass (g/cm2) 0 1 2 3 4 5 60 0.05 Bottom elevation Deposit mass (Simulation) Deposit mass (Experimental) Initial topography Figure 6: Comparison with laboratory experiments: Deposit distribution for experiment B11 7.2 Successive turbidity currents over an initial flat bottom Inspired by the experiments shown in [17], we intend now to simulate spillover turbidity currents on submarine channel levees, that is a succession over a long time period of several turbidity currents. Following this idea, we generate successive turbidity currents over an initial 25mlong flat bottom. Each turbidity current is generated by supplying a suspension of sediments at the upstream end of the flow field. The mixture is composed by three sediment species with the same characteristics as the ones used on experiments given by Table 1. The height is fixed to 0.5mand the flux is 0.001m2/s. The supply duration is set at 4sand the turbidity is let to evolve for 100sbefore the next supply arrives. The volume of deposits is then converted to bed thickness so that the final bed after each supply is considered as the initial topography for the following and is supposed to solidify and to be non-erodible when the next plume arrives. In Figure 10 we represent the bottom evolution for this successive flow events. Each layer represents the resulting bottom after 500 flow events. We remark the appearance of wavy structures by repeated deposition of turbidities over a flat bottom as it would be expected in physical situations. 22 0 1 2 3 4 5 6 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 Distance from gate (m) Deposit mass (g/cm2) 0 1 2 3 4 5 60 0.05 0.1 Bottom elevation Deposit mass (Simulation) Deposit mass (Experimental) Initial topography Figure 7: Comparison with laboratory experiments: Deposit distribution for experiment C5 0 1 2 3 4 5 6 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 Distance from gate (m) Deposit mass (g/cm2) 0 1 2 3 4 5 60 0.01 0.02 0.03 0.04 Bottom elevation Deposit mass (Simulation) Deposit mass (Experimental) Initial topography Figure 8: Comparison with laboratory experiments: Deposit distribution for experiment C7 7.3 Solid transport of topography We shall now show the influence of solid transport on the topography. We consider the case of a river (assume R0= 0) over an initial topography given by H(x, t = 0) = (0.1−0.04 sin2(π(x−3)/2),for 2 ≤x≤5, 0.1,elsewhere. (7.1) This initial topography is assumed to be erodible up to the depth level 0.102.The initial height is set to h(x, t = 0) = H(x, 0) and the initial velocity is such that hu(x, t =0)=0.01.We suppose that there is a suspension of siliceous particles (density 2650 kg/m3and vs = 0.002m/s) with initial volume fraction of 1%. At the upstream end, we set a supply with concentration 1% for the first 100s. Solid transport is considered using Grass model with Ag= 0.005. First, assume that there is no erosion of the topography and no deposition of particles, so that only solid transport is considered. Figure 11 shows the topography evolution. Now, let us consider the same situation but adding erosion and deposition effects. The results are shown on Figure 12. 23 0 1 2 3 4 5 6 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 Distance from gate (m) Deposit mass (g/cm2) 0 1 2 3 4 5 60 0.05 0.1 Bottom elevation Deposit mass (Simulation) Deposit mass (Experimental) Initial topography Figure 9: Comparison with laboratory experiments: Deposit distribution for experiment C13 0 2 4 6 8 10 12 14 16 18 20 −2 −1.5 −1 −0.5 0 0.5 Bottom topography Figure 10: Successive turbidity currents over an initial flat bottom: Topography evolution. Each layer corresponds to 500 turbidity beds 7.4 A 2D example As it has been said in Section 6, the model and numerical schemes can be easily extended to 2D turbidity currents. Here we shall show a test case related to experiment C5 in Section 7.1. We consider the domain given by Figure 13. A 0.5-m-long by 1-m-wide gate box has been placed at the upstream end with a mixture of sediments with same characteristics as the ones used in Section 7.1. The initial height of the mixture is set to 0.2m and we consider a hump topography like the one on Figure 14 with ramp and hump height and length corresponding to experiment C5 in Section 7.1. Figures 16-17 show the deposit density along the longitudinal and transversal sections on Figure 13. 24 2 2.5 3 3.5 4 4.5 5 5.5 6 6.5 7 −0.1 −0.095 −0.09 −0.085 −0.08 −0.075 −0.07 −0.065 −0.06 Initial Topography Topography evolution Figure 11: Solid transport of topography: Bottom evolution without erosion nor deposition of particles. Each layer represents 400s 2 2.5 3 3.5 4 4.5 5 5.5 6 6.5 7 −0.1 −0.095 −0.09 −0.085 −0.08 −0.075 −0.07 −0.065 −0.06 Initial topography Topography evolution Figure 12: Solid transport of topography: Bottom evolution with erosion and deposition of particles. Each layer represents 400s 25