Some Geophysical Applications with Finite Volume Solvers of Two-Layer and Two-Phase Systems

Some Geophysical Applications with Finite Volume Solvers of Two-Layer and Two-Phase Systems E.D. Fernández-Nieto Abstract There exists in the literature a huge range of geophysical applications that have been modeled trough two-layer or two-phase models. In this work first some averaged two-layer and two-phase models are presented. We focus on applications to submarine avalanches, debris flows and sediment transport in rivers. Secondly, their numerical approximation by a finite volume method is discussed and a numerical test is presented. Keywords Saint-Venant model · Two-layer · Two-phase · Finite volume · Well-balanced MSC (2010) 35L40 · 35L65 · 65Z05 · 76T25 1 Introduction In these notes firstly we describe several depth-averaged models which have been presented in the bibliography to study submarine avalanches, debris flows and sediment transport problems. Secondly, we describe a well-balanced finite volume solver to approximate the solution of the models. There are several difficulties related to the discretization of these systems, which can be written under the structure of a hyperbolic system with a conservative term, a non-conservative product and source terms. Some of the difficulties to discretize these systems are the following: (i) in some models the flux function not only depend on the vector of unknowns, but it can also depend on a given function. (ii) the coupling term between the layers or the phases is usually written as a non-conservative product. Then, it is not well defined as a distribution and the choice of a family of paths is These notes are part of several joint works with F. Bouchut, M.J. Castro Díaz, A. Mangeney, T. Morales de Luna and G. Narbona-Reina. E.D. Fernández-Nieto (B) University of Seville, 2 Av. Reina Mercedes, 41012 Seville, Spain e-mail: [email protected] necessary (see [3]). (iii) the source terms are defined in terms of the fixed topography. (iv) all the models considered in this work include a source term corresponding to a Coulomb friction law. One of the difficulties of this term is that it is multievaluated for the case of a material at rest. The notes are organized as what follows: In Sect. 2we describe the model proposed in [4] to study submarine avalanches. In Sect. 3a model to study debris flows is presented. Section4is devoted to detail the influence of the Coulomb friction term in sediment transport models. In Sect.5a finite volume method that can be used to discretize these models is presented. In Sect.6we present a numerical test. And in Sect.7some conclusions and perspectives. 2 A Depth-Averaged Model for Submarine Avalanches Some submarine avalanches can be produced by a collapse of the sediment layer. This can be produced in areas near the coast with a high bathymetry gradient and with a high amount of sedimentation. That can be produced for example near a river mouth. A high concentration of sediment at the bed can be also destabilized with a small earthquake. In some cases the evolution of the sediment layer can degenerate to a submarine avalanche which produces finally a tsunami. In some other cases an aerial avalanche near the coast or a lake can produce a tsunami. In [4] it is proposed a two-layer SWE system where the submarine avalanche and the eventual generated tsunami can be studied. The first layer corresponds to the fluid and the second one to the sediment layer. Heinrich et al. proposed also a two-layer SWE [8] but without taking into account the effects of the fluid on the landslide dynamics, the sea-bottom deformation induced by the landslide is used as input data in the tsunami model. For the sediment layer a Savage–Hutter type model is considered. The pionering work of Savage–Hutter [15] derives a model to describe granular flows over a slopping plane based on Mohr–Coulomb considerations: a Coulomb friction law is assumed to reflect the avalanche/bottom interaction and the normal stress tensor is defined by a constitutive law relating the longitudinal and the normal stresses through a proportionality factor. One of the characteristics of the model proposed in [4] is that the definition of the Coulomb friction term takes into account buoyancy effects, because we are studying submarine avalanches. Another characteristic is that, depending on the aspect ratio between the water density and the sediment density, the movement of the sediment avalanche can be more or less influenced by the presence of the fluid. The submarine avalanches produce a movement on the fluid layer, whose consequence can be a tsunami. With index 1 is denoted the upper layer, composed of a homogeneous inviscid fluid of constant density ρ1, and with index 2 the grain layer of density ρ2=(1−ψ0)ρs+ψ0ρ1, Water surface θ b x z X=(x,z) Z X 1 h 2 h Fig. 1 Two-layer: local coordinates where ρsis the grain density and ψ0is the porosity of the layer. Let us remark that this definition corresponds to the case in which the grain layer is fully immersed, where we suppose that the pores in the grain layer are filled with the fluid of the upper layer. Otherwise, if we consider a test of an aerial avalanche that enters into a fluid, the value of ρ2is variable in space. Its definition is discontinuous at the interface between the aerial and the submarine avalanche. The model is described in local coordinates over a non-erodible bottom (See Fig.1). Then, h1and h2denote the thickness of the fluid and grain layers, respectively, measured orthogonally to the bottom. Ui,i=1,2, is the velocity parallel to the bottom and Wi,i=1,2, is the velocity perpendicular to the bottom, with ireferring to layers 1 and 2. Then, the model is defined by the following set of equations: ⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ ∂th1+∂X(h1U1)=0, ∂t(h1¯ U1)+∂Xh1U2 1+gh2 1 2cos θ=−gh1dXb+gsin θdXθh2 1 2− −gh1∂X(cos θh2)+1 ρ1 fric(U1,U2), ∂th2+∂X(h2U2)=0, ∂t(h2U2)+∂Xh2U2 2+gh2 2 2cos θ(rλ2+K(1−rλ2))=−gh2dXb− −rgh2(λ1+K(1−λ1))∂X(h1cos θ)−1 ρ2 fric(U1,U2)+gh2 2 2sin θdXθ+T, (1) the term fric(U1,U2)models the friction between the two layers. It is proportional to the relative velocity between the layers and the harmonic average of ρ1h1and ρ2h2 (see [16] and references therein). The term Tis defined by a Coulomb friction law (see [15]). We observe that this term must be understood as (see [9]): If |T|≥σc⇒T=−(g(1−r)h2cos θ+h2U2 2dXθ) U2 |U2|tanδ0,(2) If |T|<σ c⇒U2=0,(3) where σc=g(1−r)h2cos θtanδ0. The density aspect ratio is denoted by r=ρ1/ρ2, where ρ1is the density of the fluid and ρ2is the density of the mixture sediment layer. Note that the definition of rcan be variable in space and time in the case of tsunamis produced by aerial avalanches near the coast. 3 Depth-Averaged Models for Debris Flows One of the simplifications considered in the model presented in previous sections is done on the mixture sediment layer. Instead to consider a two-phase model it is simulated as a one-layer where the solid and fluid phases move at the same velocity. Moreover the porosity is considered to be constant. An improvement of the previous model is to couple the fluid layer with a two-phase model for the mixture sediment layer. Two-phase models are also considered to study debris flows. A two-phase depth averaged model was proposed by Pitman and Le in [14], and reformulated by Pelanti et al. in [13]. One of the difficulties of two-phase models is the definition of a closure relation. In these model the closure was related to the pressure at the free surface of the mixture layer. Being the problem that this closure does not imply that the model has a dissipative energy. In [1] a two-phase model has been proposed with a dissipative energy balance. The model proposed in [1] can be rewritten in terms of the solid pressure at the free surface, the fluid pressure at the free surface, the solid pressure at the bed, or the fluid pressure at the bed. This extra unknown can be seen as a Lagrangian multiplier associated to the closure relation, that is the incompressibility of the solid phase. Although all the reformulations are equivalent, from the numerical point of view we propose to rewrite the model in terms of the solid pressure at the free surface. Because, by considering a projection method at the first step the hyperbolic model has better properties. In order to work in a framework well suited to avalanches problems it is interesting to write the model in local coordinates. In previous section we have described the model in local coordinates over a variable non-erodible bottom. Another possibility is the one considered in [1]. The model is written in an inclined plane frame. Thus, we consider a fixed slope with constant angle θwith respect to the horizontal, −π/2< θ<π/2, and the coordinates (x,z)are respectively tangent to and normal to this slope, the xaxis being along the steepest direction, and the yaxis being horizontal (see Fig.2). Fig. 2 Local coordinates on a plane of reference and bathymetry z x b h b ~ Then, for a given bathymetry, firstly a reference plane is defined, by taking into account the mean slope of the domain. We denote this plane by ˜ b(x). Secondly, by a projection of the bathymetry on this reference plane we define b(x). We consider a bottom topography b(x), and a thin layer of material over it with thickness h(t,x). The material thus occupies the domain b(x)<z<b(x)+h(t,x). By denoting by ψthe solid pressure at the free surface divided by the solid density, the two-phase system can be written as follows: ∂t(hϕ) +∂x(hϕv)=0,(4a) ∂t(h(1−ϕ)) +∂x(h(1−ϕ)u)=0,(4b) ∂t(hϕv)+∂x(hϕv2)=−h(1−ϕ)∂xψ−ϕgh cos θ∂ x(˜ b+b+h) −1 2(1−rm)gh2cos θ∂ xϕ −ϕgh sin θ+βh(u−v), +T,(4c) ∂t(h(1−ϕ)u)+∂x(h(1−ϕ)u2)=h r(1−ϕ)∂xψ−(1−ϕ)gh cos θ∂ x(˜ b+b+h) −(1−ϕ)gh sin θ−1 rm βh(u−v), (4d) ∂x(h(1−ϕ)(u−v)) =0.(4e) As in the model presented in previous section the term Tdenotes the contribution in the model of the Coulomb friction law. In this case it is defined as follows: If |T|≥σc,m⇒T=−gcos θ(1−rm)hϕtan δ0 v |v|,(5) If |T|<σ c,m⇒v=0,(6) where σc,m=gcos θ(1−rm)hϕtan δ0. The velocities are denoted by vfor the solid phase and ufor the fluid phase. Moreover, βis the friction coefficient between the phases and rmis the density aspect ratio in the mixture layer, rm=ρf/ρs,beginρfthe fluid density. Let us remark that by coupling this model with the one presented in previous section we have ρf=ρ1. Moreover the weight of the upper fluid layer on the mixture layer can be directly imposed in the definition of ψ. A depth-averaged two-phase model including dilatancy effects has been proposed in [2]. It can be seen as an extension with slope aligned variable dependency of the model proposed in [10](seealso[7]). 4 Sediment Transport: Influence of the Coulomb Friction Law Finally, let us present a simplified version of one of the models deduced in [6]. It has been deduced through an asymptotic analysis, and a coupling between a Shallow Water system for the fluid layer with a Reynolds equation to model the evolution of the sediment layer. Interestingly, this model can be deduced under very similar assumptions than that of the model presented in Sect.2for submarine avalanches. The main difference is the characteristic time at which moves the sediment layer. For the case of sediment transport in rivers or coastal areas the time scale of the movement of the sediment bed is very different from the one of the fluid layer. The second main part in the deduction of this model is the definition of the friction term between the fluid and the sediment layer. The friction coefficient is defined in terms of the ratio between the Coulomb friction angle and the critical Shield parameter. The main difference with the submarine avalanche model is that in this case we do not consider local coordinates, it is enough to write the model in cartesian coordinates. While for the case of avalanches model it is crucial to consider some kind of local coordinates, as the ones presented in previous sections. For this model index 1 denotes again the fluid layer and index 2 the sediment layer. For the case of quasi-uniform flows the model can be written under the form of a SWE system as follows (for more details see [6]): ⎧ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎩ ∂th1+∂x(h1u1)=0, ∂tq1+∂x(h1u2 1)+1 2g∂xh2 1+gh1∂x(b+h2)+ghm rP=0, ∂th2+∂xh2vb(1/r−1)gds=0, (7) with P=∇ x(rh1+h2+b)+(1−r)sgn(u2)tan δ0.(8) Instead to have a momentum equation for the sediment layer, the velocity of the sediment layer is deduced from the asymptotic analysis, it is defined by vb(1/r−1)gds. The main part in the deduction of the non-dimensional bedload velocity vbis the friction coefficient that defines the friction law between the fluid layer and sediment layer. It has to depend on tan δ0,δ0being the same friction angle that defines the Coulomb friction law of the models presented in previous sections for the case of submarine avalanches and debris flows. For the case of a linear friction law between the fluid and the sediment vbis defined as follows: vb=1 (1/r−1)gds u1−ϑ 1−rP,(9) where ϑ=θc/tan δ0. Note that the sign of the velocity of the sediment layer, sgn(u2), has still to be defined. Observe that this coefficient comes from the contribution of the Coulomb friction law at the interface between moving and static sediment particles (see (8)). In order to specify the sign of u2, it is necessary to note that Coulomb friction force has the same sign as the net force acting on the sediment. In fact, note that the definition of Tin the models presented in the two previous sections can be reinterpreted in this way. That is, Tis defined in such a way that the contribution of the Coulomb friction term does not change the sign of the velocity of the layer without this contribution. From the definition of Twe obtain the definition of sgn(u2)for this model corresponding to sediment transport. Concretely, for this model we have sgn ϑhm 1−r(1−r)sgn(u2)tan δ=sgnhmu1 (1/r−1)gds −ϑhm 1−r∂x(rh1+h2+b). Then, using that ϑ=θc/tan δ0, sgn(u2)=sgn u1 (1/r−1)gds −ϑ 1−r∂x(rh1+h2+b).(10) 5 Well-Balanced Finite Volume Methods In this section we present a finite volume method that can be considered to discretize the models presented in previous sections, that is for submarine avalanches, debris flows and sediment transport. The common parts of all these models are the presence of nonconservative products, geometric source terms, and a Coulomb friction law. Concretely, all the models presented previously can be written under the form of a hyperbolic system with a conservative product, a non-conservative term and source terms with the following structure: ∂tW+∂xF(θ, W)=G1(x,W)∂xb+G2(x,W)∂xθ+B(W)∂xW+T.(11) For the case of the model presented in Sect.3this corresponds to the first step in the projection method that we can consider to approximate the Lagrangian multiplier associated to the restriction of incompressibility of the solid phase. The only model for which the flux function depends on a function θ=θ(x)is the case of submarine avalanches. Because it is written in local coordinates on a given non-erodible bathymetry. The source terms modeling the friction between the two layers or the two phases are discretized semi-implicitly (see [13]). Then, we do not detail it in the finite volume discretization. For the discretization of the system, computing cells Ii=[xi−1/2,xi+1/2]are considered. For simplicity, we suppose that these cells have constant size Δx. Let us define xi+1 2=iΔxand by xi=(i−1/2)Δx, the center of the cell Ii.LetΔtbe the constant time step and define tn=nΔt. We denote by Wn ithe approximation of the cell averages of the exact solution provided by the numerical scheme. The discretization of B(W)∂xWfirstly requires to interpret this term as a Borel measure (see [3]), depending on the choice of a family of paths linking given states. Here the family of segments are considered as in [12]. The dependence of the flux function on θ(x), makes it difficult to obtain the desired exact well-balanced property for water at rest. That is, to preserve exactly the stationary solution of the system with zero velocity. For the case of the submarine avalanche model it is helpful to define the flux function F(θ, W)as a function of cos θand cos2θ(see [4] for more details). Finally, as mentioned before, the discretization of the source term T(W)corresponding to the Coulomb friction term is crucial to simulate properly the landslides for the case of submarine avalanches and debris flows. For the sediment transport model, its discretization implies that the sediment moves only when the friction is greater than a threshold. We propose a two-step numerical scheme to treat the Coulomb friction term. Let us suppose that the values Wn iare known. In order to advance in time we proceed as follows: •First Step. Let us denote the unknown approximation as W∗ i=[H∗ 1,iQ∗ 1,iH∗ 2,i Q∗ 2,i]T. For the case of the sediment transport model we do not have the unknown Q2. W∗ i=Wn i−Δt ΔxDFn,+ i−1/2+DFn,− i+1/2,(12) where DFn,± i+1/2=DF± i+1/2(Wn i,Wn i+1)are the generalized Roe flux difference computed using a family of segments (see [11,12]). DF± i+1/2=1 2F(θi+1/2,Wi+1)−F(θi+1/2,Wi)+∂θF(θi+1/2,Wi+1/2)(θi+1−θi) +G1,i+1/2(bi+1−bi)+G2,i+1/2(θi+1−θi) ±Pi+1/2[Ai+1/2(Wi+1−Wi)+∂θF(θi+1/2,Wi+1/2)(θi+1−θi) +G1,i+1/2(bi+1−bi)+G2,i+1/2(θi+1−θi)−Ti+1/2Δx]}, where G1,i+1/2,G2,i+1/2and Ti+1/2are approximations at the interface i+1/2ofthe corresponding source terms. And matrix Pi+1/2is an approximation of the sign of the Roe matrix Ai+1/2(see [11,12]). For the three models presented in this note we consider the same kind of numerical diffusion, which is based on the IFCP method (Intermediate Field Capturing Parabola method) introduced in [5]. Matrix Pi+1/2is defined as follows: Pi+1/2=α0Ci+1/2+α1I+α2Ai+1/2, where, Iis the identity matrix and Ci+1/2is an approximation of the inverse of Ai+1/2. The well-balanced properties of the scheme depend on this approximation. For example, if we are only interested into preserving exactly stationary solutions at rest it is enough to defined Ci+1/2as the inverse matrix of Ai+1/2evaluated with zero velocity. The coefficients α0,α1and α2are the solutions of the following linear system ⎛ ⎝ 1λ− ext (λ− ext )2 1λ+ ext (λ+ ext )2 1χint χ2 int ⎞ ⎠⎛ ⎝ α0 α1 α2⎞ ⎠=⎛ ⎝ |λ− ext | |λ+ ext | |χint|⎞ ⎠, where χint =sgn(λ− ext +λ+ ext )max(|λ− int|,|λ+ int|). For the models corresponding to the submarine avalanches model and for debris flows we have four eigenvalues associated to the Roe matrix, two externals and two internals, verifying λ− ext <λ − int <λ + int <λ + ext . For the case of the sediment transport model we have only three eigenvalues, then we apply previous definition by identifying λ− int =λ+ int, the intermediate eigenvalues. •Second step. For the case of submarine avalanches model and debris flows it is still necessary to introduce the approximation of the Coulomb friction term. For the case of the sediment transport model it is not necessary because it has been directly taken into