Efficient numerical schemes for viscoplastic avalanches. Part 2: the 2D case

Efficient numerical schemes for viscoplastic avalanches. Part 2: the 2D caseI Enrique D. Fernández-Nietoa, José M. Gallardob, Paul Vigneauxc,d,∗ aDepartamento de Matemática Aplicada I, Universidad de Sevilla. E.T.S. Arquitectura. Avda, Reina Mercedes, s/n. 41012 Sevilla, Spain bDepartamento de Análisis Matemático, Estadística e I. O. y Matemática Aplicada, Universidad de Málaga. F. Ciencias, Campus Teatinos s/n, 29071 Málaga, Spain cUnitée de Mathématiques Pures et Appliquées CNRS UMR 5669, Ecole Normale Supérieure de Lyon. 46 allée d’Italie. 69364 Lyon Cedex 07. France dUniv. Lyon, INRIA Team NUMED, Ecole Normale Supérieure de Lyon. 46 allée d’Italie. 69364 Lyon Cedex 07. France Abstract This paper deals with the numerical resolution of a shallow water viscoplastic flow model. Viscoplastic materials are characterized by the existence of a yield stress: below a certain critical threshold in the imposed stress, there is no deformation and the material behaves like a rigid solid, but when that yield value is exceeded, the material flows like a fluid. In the context of avalanches, it means that after going down a slope, the material can stop and its free surface has a non-trivial shape, as opposed to the case of water (Newtonian fluid). The model involves variational inequalities associated with the yield threshold: finite volume schemes are used together with duality methods (namely Augmented Lagrangian and Bermúdez-Moreno) to discretize the problem. To be able to accurately simulate the stopping behaviour of the avalanche, new schemes need to be designed, involving the classical notion of well-balancing. In the present context, it needs to be extended to take into account the viscoplastic nature of the material as well as general bottoms with wet/dry fronts which are encountered in geophysical geometries. Here we derive such schemes in 2D as the follow up of the companion paper treating the 1D case. Numerical tests include in particular a generalized 2D benchmark for Bingham codes (the Bingham-Couette flow with two non-zero boundary conditions on the velocity) and a simulation of the avalanche path of Taconnaz in Chamonix - Mont-Blanc to show the usability of these schemes on real topographies from digital elevation models (DEM). Keywords: Viscoplastic, Shallow Water, Finite Volume, Well-Balanced, Variational inequality, Bingham, Taconnaz 1. Introduction This article is the sequel of the companion paper [10]. It deals with the most general framework for the simulation of viscoplastic (Bingham) avalanches with wet/dry fronts on general 2D bottoms, thanks to a shallow water model discretized with a finite volume method. This approach has several motivations: (i) finite volume methods are used in nearly 40% of the discretizations of geophysical avalanches (40% of the rest being done with finite difference methods) as mentioned in the review [25], (ii) we use the variational inequality framework with duality methods which have proved to be the most accurate for the computation of viscoplastic flows, see [21], (iii) enrichment of geophysical shallow models towards viscoplastic behaviour is increasingly in use to take into account the material ability to rigidify [1,25,21], (iv) the prototypical shallow viscoplastic model covered in this paper is well adapted to (wet) dense snow avalanches (by opposition to powder snow avalanches) which are occurring more and more frequently with the global warming of the atmosphere (see [2]). In the present work, we extend in 2D the 1D schemes developed in [10]. We make a careful study of the ability of the schemes to compute the stationary states of an avalanche. This is done thanks to a coupling between the finite volume method used for the discretization in space and the duality method used to solve the non-Newtonian character of the material. This leads to an extended notion of viscoplastic well-balancing. In this case, we say that the method is well-balanced if it preserves exactly two kinds of stationary solutions: (i) material at rest independently of the rigidity of the material and (ii) rigid enough material with free surface parallel to a reference plane. Let us remark that the latter is more relevant for viscoplastic materials, nevertheless it is necessary to preserve also the first one, material at rest, in order to be consistent with a numerical method for a Newtonian fluid when τytends to zero. We study also the numerical cost and (when possible) the a priori estimation of the optimal intrinsic parameter of two duality methods: the augmented Lagrangian (AL) and the Bermúdez-Moreno (BM) methods. Note that the IThis document is an improved version (some misprints corrected) of the published article whose DOI is given below. ∗Corresponding author Email addresses: [email protected] (Enrique D. Fernández-Nieto), [email protected] (José M. Gallardo), [email protected] (Paul Vigneaux) Preprint submitted to J. of Comput. Phys. – schemes proposed here can be extended straightforwardly to power viscoplastic laws (such as Herschel-Bulkley). Several computational tests illustrate the performances of the schemes. First, we consider a 2D Couette geometry for which we propose a generalized analytic solution with two non-zero boundary conditions on the velocity (usual solutions assume that one of the boundary is fixed). This can be a benchmark for testing classic 2D Bingham codes. Second, we study the well-balanced property for rigid materials by considering a complex random bottom on a 30◦slope reference plane: the free surface is parallel to the reference plane and has complex wet/dry fronts. Third, we build a numerically demanding 2D academic dam break test on a complex topography where the final stationary solution exhibits strong gradients of the free surface at the wet/dry front. Finally, we show the ability of the schemes to compute the final stationary state of an avalanche on a real topography (given by the ASTER Digital Elevation Model): we used historical data of (frequent) avalanches at the Taconnaz path in the Mont-Blanc, one of the longest sites in Europe with a path close to 7km. This test involves a large computation domain and long physical times with a rich dynamics of the progressive stopping of the avalanche. In section 2, we recall the model under consideration. In section 3, we derive the 2D versions of the augmentedLagrangian and the Bermúdez-Moreno methods. For the latter, we give a theoretical a priori estimation of duality parameter which leads to smallest computational time of the duality resolution. This is associated to the computation of the velocity field. We then detail (section 4) the construction and properties of the 2D well-balanced finite volume method for the space discretization. Numerical illustrations are finally presented in section 5. 2. Models The model problem for viscoplastic shallow flows is naturally the one presented in the first part of this work [10]. We refer to [4] for more details. The geometry is as shown on Figure 1. We consider a fluid domain of height Figure 1: Sketch of the 2D domain and convention for the local coordinates. x=(x1,x2) Hover a general bottom b. More precisely, let Ω⊂R2be a given domain for the space variable x. The R2plane generated by Ωis supposed to be sloping at an angle αfrom the horizontal plane. We denote by z∈Rthe variable in the orthogonal direction to Ω. The bottom which bounds the fluid by below is defined by b(x),x∈Ω. We denote by D(t) the fluid domain defined as D(t)={(x,z)∈Ω×R/b(x)<z<b(x)+H(t,x)},(1) where His the time-dependent height of the fluid. As usual for shallow water type models, we denote by V=V(t,x)∈R2the vector of the average of the velocity (orthogonal to the z-axis) along the depth of the fluid (i.e. from z=b(x) to z=b(x)+H(t,x)). We take into account the fact that there may be friction on the bottom through a coefficient β. The fluid undergoes a body force denoted as ( fΩ,fz)∈R2×Rin the Ω×zframe of reference. Note that fΩand fzare both assumed to be constant. Since we are considering a Bingham constitutive law, the material is characterized by a viscosity ηand a yield stress τy. The latter is associated to the plastic behaviour of the material and this leads (cf. [8]) to a variational inequality for the momentum conservation relation (see equation 4). On the contrary, the conservation of mass is rather classic for this type of integrated model (see equation 3). Given the space V(t)={Ψ∈H1(Ω)2/Ψ = 0 on ∂Ω}:=H1 0(Ω)2(2) 2 and some initial conditions at t=0, the problem is to find H∈L2([0,T],L∞(Ω)), V∈L2([0,T]; V(t)), with ∂tV∈L2([0,T]; L2(Ω)2), such that ∂tH+divx(HV)=0,(3) and ∀Ψ∈ V(t),ZΩ H∂tV+(V·∇x)V·(Ψ−V)dx+ZΩ βV·(Ψ−V)dx +ZΩ 2ηHD(V) : D(Ψ−V)dx+ZΩ 2ηHdivxV(divxΨ−divxV)dx +ZΩ √2τyHpD(Ψ): D(Ψ)+(divxΨ)2−pD(V): D(V)+(divxV)2dx ≥ZΩ H(fΩ+fz∇xb)·(Ψ−V)dx−ZΩ H2 2fz(divxΨ−divxV)dx,(4) where D(U) :=1 2h∇xU+(∇xU)ti,(5) ∇xU:= ∂Ui ∂xj!i,j ,i=1,2,j=1,2,(6) divxU:=∂U1 ∂x1 +∂U2 ∂x2 ,∀U(t,·) :=(U1,U2)∈ V(t).(7) The shallow water formulation (4) is in weak form. Note that the derivation of (4) comes from the asymptotic analysis of the integrated 3D equations: this explains why we first present this variational form. From (4), we can then find the strong (i.e. non-variational) form (9)-(10) , which is clearer to read. Obviously, we obtain a Bingham constitutive law but it is modified compared to the canonical law. Indeed, the integration of the 3D equations to the 2D form leads to an integrated Bingham law (10) (with corrector terms tr(D(V))I). Let us consider the space X=R2×2with scalar product (p,q) :=(p:q+tr(p) tr(q)) and associated norm k·k. That is, for p∈X, kpk=v u tX i,j p2 i,j+X i pi,i 2 .(8) By using this notation the strong formulation of (4) can be written as follows: H∂tV+(V·∇x)V−divx(Hσ)=H(fΩ+fz∇xb)−∇x H2 2fz!−βV,(9) where                σ=2η(D(V)+tr(D(V))I)+√2τy D(V)+tr(D(V))I kD(V)kif kD(V)k,0 kσk ≤ √2τyif kD(V)k=0. (10) Let us remark that the √2 factor multiplying τyappears because of the use of Frobenius norm. Then, the actual formulation is equivalent to the one considering a Eulerian norm. Note that in the following, the body force will be the influence of gravity, denoted by g. To write this force, we must decide what is the orientation of the plane generated by Ω; by convention we will say that if (x1,x2,z) is the frame of reference (cf. Figure 1), then the tilted axis (with respect to the horizontal) is x1, i.e. fΩ=(−gsin α , 0),fz=−gcos α. (11) Note also that for numerical accuracy, it is often better to consider the simulation of geophysical flows over large domains (like for instance in Section 5.3 for the Taconnaz avalanche path) by rescaling the equations to simulate the flow on a domain of order one length. Namely, we introduce a characteristic horizontal length Lcand vertical height Hc. We then make the following rescaling (denoting =Hc/Lc): x=Lc˜ x,z=Hc˜z,t=Lc √gHc ˜ t,b=Hc˜ b,V=pgHc˜ V, η =HcpgHc˜η, τy=gHc˜τy, β =pgHc˜ β, f.=g˜ f.. (12) 3 In these new variables (and omitting the tildes), (4) reads: ∀Ψ∈ V(t),ZΩ H∂tV+(V·∇x)V·(Ψ−V)dx+ZΩ βV·(Ψ−V)dx +ZΩ 2ηHD(V) : D(Ψ−V)dx+ZΩ 2ηHdivxV(divxΨ−divxV)dx +ZΩ √2τyHkD(Ψ)k−kD(V)kdx ≥ZΩ H(fΩ+fz∇xb)·(Ψ−V)dx−ZΩ H2 2fz(divxΨ−divxV)dx.(13) As most often done in the literature and since the main objective of this work is to treat the viscoplastic discretization difficulty, we consider a first order backward semi-discretization in time. If we denote by ∆tthe time step, we have from (3)-(4) : Hn+1−Hn ∆t+divx(HnVn)=0,(14) and ZΩ HnVn+1−Vn ∆t+(Vn·∇x)Vn·(Ψ−Vn+1)dx+ZΩ βVn+1·(Ψ−Vn+1)dx +ZΩ 2ηHnD(Vn+1) : D(Ψ−Vn+1)dx+ZΩ 2ηHndivxVn+1(divxΨ−divxVn+1)dx +ZΩ √2τyHnkD(Ψ)k−kD(Vn+1)kdx ≥ZΩ Hn(fΩ+fz∇xb)·(Ψ−Vn+1)dx−ZΩ (Hn)2 2fz(divxΨ−divxVn+1)dx,∀Ψ.(15) Doing so, we see that problems on the height and on the velocity are decoupled. At each time step, supposing that we know (Hn,Vn), we need to solve both problems for (Hn+1,Vn+1). As in the companion paper [10], we compare two duality methods to handle the variational inequality of the problem on the velocity, namely the Augmented Lagrangian method and Bermúdez-Moreno method. It is the subject of the next section. 3. Duality methods in 2D 3.1. The AL approach We will extend in 2D the derivation done in [10]. Supposing that (Hn,Vn) are known, the goal is here to solve the problem (15) for Vn+1. Using ad hoc spaces, variational inequality (15) is now equivalent to a minimization problem Jn(Vn+1)=min V∈V Jn(V),(16) where Jn(V)=Fn(B(V)) +Gn(V), with V=(H1 0(Ω))2. Let us also denote H=L2(Ω)2×2, B: V → H V7→ B(V)=D(V)!,Fn: H → R λ7→ Fn(λ)=RΩτyHnkλkdx!, and Gn:V → R, Gn(V)=ZΩ Hn |V|2/2−Vn·V ∆t+(Vn·∇xVn)·V!dx+ZΩ β|V|2 2dx +ZΩ ηHnkB(V)k2dx−ZΩ Hn(fΩ+fz∇xb)·Vdx+ZΩ fz (Hn)2 2divxVdx. In 2D, we define the Lagrangian functional by Ln:V×H ×H → R, Ln(V,q, µ)=Fn(q)+Gn(V)+ZΩ Hn(µ, B(V)−q)dx, 4 and the augmented Lagrangian functional, for a given positive value r∈R, as: Ln r(V,q, µ)=Ln(V,q, µ)+r 2ZΩ HnkB(V)−qk2dx.(17) Again, we determine the saddle point of Ln r(V,q, µ) over V×H×H thanks to an augmented Lagrangian algorithm (cf. [11]). Augmented Lagrangian algorithm (2D) •Initialization: suppose that Vn,Hnand µnare known. For k=0, we set Vk=Vnand µk=µn. •Iterate: –Find qk+1∈ H solution of Ln r(Vk,qk+1, µk)≤ Ln r(Vk,q, µk),∀q∈ H. In other words, qk+1∈ H is the solution of following minimization problem: min q∈H Hnr 2kqk2+Hn√2τykqk− Hn(µk+rB(Vk)): q−Hntr(µk+rB(Vk)) tr(q).(18) And the solution of this problem is computed locally for all x∈Ω: qk+1=         0 if kµk+rB(Vk)k<√2τy, 1 r(µk+rB(Vk)) −√2τy µk+rB(Vk) kµk+rB(Vk)kotherwise.(19) –Find Vk+1∈ V solution of Ln r(Vk+1,qk+1, µk)≤ Ln r(V,qk+1, µk),∀V∈ V. From (17), by differentiating Ln r(V,q, µ) with respect to V, we deduce that Vk+1is the solution of the following linear problem (whose resolution is detailed in Section 4): Hn Vk+1−Vn ∆t!+βVk+1−(2η+r)divx(HnD(Vk+1)) +∇x(Hndivx(Vk+1)) +Hn(Vn·∇xVn)−(fΩ+fz∇xb)Hn−∇x (Hn)2 2fz! −divxHnµk−rqk+1−∇xHntr(µk−rqk+1)=0.(20) –Update the Lagrange multiplier via µk+1=µk+rB(Vk+1)−qk+1.(21) –Check convergence (see below) and update: Vk=Vk+1,µk=µk+1,k=k+1 and go to the next iteration ... •... until convergence is reached: kµk+1−µkk kµkk≤tol.(22) At convergence, we get the value of Vn+1by setting Vn+1=Vk+1(in the numerical tests presented in this paper, we set tol =10−5). It is also shown in [11] that this algorithm converges to the saddle point of (17). 5 3.2. The BM approach and its optimal parameter 3.2.1. The BM algorithm The BM algorithm in the two-dimensional case follows similar guidelines as in [10], once a proper choice of norms is made. In particular, we use the space V=H1 0(Ω)2endowed with the scalar product (V,W)V=ZΩ D(V): D(W)dx+ZΩ divx(V)divx(W)dx,V,W∈ V. It readily follows from the arithmetic-geometric mean property and Korn inequality that the associated norm kVkV=kD(V)k2 L2+kdivx(V)k2 L21/2 verifies C−1 Kk∇VkL2≤ kVkV≤√3k∇VkL2, where CKis a Korn constant. Therefore, the norm k·kVis equivalent to the norm k∇·kL2in H1 0(Ω)2, so Vturns out to be a Hilbert space. In a similar way, the space H=L2(Ω)2×2is also a Hilbert space with the scalar product (Z,W)H=ZΩ Z:Wdx+ZΩ tr(Z) tr(W)dx, as the associated norm is equivalent to k·kL2. Notice also that kVkV=kB(V)kHfor every V∈ V. Consider now the linear operator A:V → V0defined as hA(V),Ψi=ZΩHn ∆t+βV·Ψdx+ZΩ 2ηHnD(V): D(Ψ)dx+ZΩ 2ηHndivx(V)divx(Ψ)dx, which is coercive with constant γ=2ηHn min (where Hn min =min Hn(x)>0): hA(V),Vi ≥ Hn min ∆t+β!kVk2 L2+2ηHn minkVk2 V≥2ηHn minkVk2 V,∀V∈ V. Define also the functional j:V → Rgiven by j(V)=ZΩ √2τyHnp|D(V)|2+divx(V)2dx, and let L∈ V0be hL,Ψi=ZΩ Hn ∆tVn·Ψdx−ZΩ HnVn· ∇xVnΨdx+ZΩ Hn(fΩ+fz∇xb)·Ψdx−ZΩ (Hn)2 2fzdivx(Ψ)dx. Then, the variational inequality (15) can be expressed as: Find V∈ V such that hA(V),Ψ−Vi+j(Ψ)−j(V)≥ hL,Ψ−Vi,∀Ψ∈ V.(23) Let Φ:Ω×X→Rbe the function Φ(x,p)=√2τyHn(x)kpk, and define T:H → Ras T(Z)=ZΩ Φ(x,Z(x))dx. Using that divx(V)=tr(D(V)) we have j(V)=T(B(V)), where B:V → H is given by B(V)=D(V). Now, reasoning as in [10], the variational inequality (23) can be written as: Find V∈ V and θ∈ H such that        A(V)+ωB∗(B(V)) +B∗(θk)=L, θ=Gω λ(B(V)+λθ),(24) where Gω λis the Yosida approximation of Gω=∂T−ωI; the parameters λand ωare arbitrary positive numbers satisfying λω < 1. The BM method for (24) reads then as follows: For k≥0, θkbeing known, compute Vkand θk+1by solving        A(Vk)+ωB∗(B(Vk)) +B∗(θk)=L, θk+1=Gω λ(B(Vk)+λθk).(25) From now on we will assume the condition λω =1/2, which ensures the convergence of the BM algorithm and it is also fundamental in the computation of the optimal parameters in Section 3.2.2 (see [10] and the references therein). 6 Remark 1. The BM method shares some conceptual properties with the iterative method introduced in [5] (see also [7, Sect. 7.3]). This method is a version of the classical Uzawa’s algorithm, which is based on a projection operator on a closed convex set. In the BM algorithm, the projector is substituted by a Yosida approximation, which can be applied in more general contexts. Indeed, when the functional j(v)is the support function of a closed convex set, BM reduces to Uzawa’s method. Recall that ([9]): ∂T(Z)={W∈ H:W(x)∈∂Φ(x,Z(x)) a.e. x∈Ω}.(26) The subdifferential of Tcan thus be computed pointwise in terms of the subdifferential of Φ. To this end, remember that we have defined the space X=R2×2with scalar product (p,q)=(p:q+tr(p) tr(q)) and associated norm k·k. Then, define φ:X→Ras φ(p)=ckpk, where cis an arbitrary constant. For p,0, the function φis Gâteaux differentiable, so the subdifferential ∂φ(p) consists only of the gradient: ∇xφ(p)=cp kpk. On the other hand, one can see that ∂φ(0) ={q∈X:kqk ≤ c}. Finally, the Yosida approximation Gω λcan be computed as follows: for a.e. x∈Ωand Z∈ H, Gω λ(Z)(x)=                   Z(x) λif kZ(x)k ≤ λ√2τyHn(x), √2τyHn(x)−ωkZ(x)k (1 −λω)kZ(x)kZ(x) if kZ(x)k> λ √2τyHn(x). This expression can be regarded as a generalization of the formula obtained in the one-dimensional case. We notice that kZ(x)k=pZ(x): Z(x)+tr(Z(x))2,a.e. x∈Ω, so the following relation holds: kZkH= ZΩkZ(x)k2dx!1/2 . We end this section by giving the explicit form of the linear problem to be solved at each iteration of (25). After integration by parts, it can be written as follows (and compared to (20)): Hn Vk+1−Vn ∆t!+βVk+1−2η(divx(HnD(Vk+1)) +∇x(Hndivx(Vk+1))) −ω(divx(D(Vk+1)) +∇x(divx(Vk+1))) +Hn(Vn·∇xVn)−(fΩ+fz∇xb)Hn−∇x (Hn)2 2fz!−divx(θk)−∇x(tr(θk)) =0.(27) 3.2.2. Study of the optimal parameter The analysis on the optimal choice of parameters performed in [10] can be adapted to the 2D case. First, let Vhbe a finite-dimensional subspace of Vof standard conforming P1finite elements, being hthe mesh size (dependence on hwill be dropped unless necessary). Now, [10, Equation (62)] and [10, Equation (64)] read as hA(Vk−V),Ψi+ω(B(Vk−V),B(Ψ))H+(θk−θ, B(Ψ))H=0,∀Ψ∈ Vh,(28) and kθk+1−θk2 H≤ kθk−θk2 H−4ωhA(Vk−V),Vk−Vi,(29) respectively. Notice that the condition λω =1/2 has been assumed. Let now CPand CKbe, respectively, the constants in the Poincaré and Korn inequalities (i.e., kΨkL2≤ CPk∇ΨkL2and k∇ΨkL2≤CKkD(Ψ)kL2, for every Ψ∈ Vh), and define γ1=C−1 PC−1 K. Let γ2be such that kVkV≤γ2kVkL2for every V∈ Vh. Then, equation (28) implies, for all Ψ∈ Vh, (θk−θ, B(Ψ))H≤Hn max ∆t+βkVk−VkL2kΨkL2+(ω+2ηHn max)kB(Vk−V)kHkB(Ψ)kH≤Γ(ω)kVk−VkL2kB(Ψ)kH, where Hn max =kHnk∞and Γ(ω)=Hn max ∆t+βγ−1 1+(ω+2ηHn max)γ2. 7 Assuming that Ψ∈ Vhis such that B(Ψ)=θk−θ, it follows that kθk−θkH≤Γ(ω)kVk−VkL2. Finally, using the inequality kVkV≥γ1kVkL2and the coerciveness of A, from (29) we deduce that kθk+1−θk2 H≤ L(ω)kθk−θk2 H, where L(ω)=1−4ωγγ2 1Γ(ω)−2. Minimization of L(ω) leads to the following expression for the optimal parameter ωopt: ωopt =γ−1 1γ−1 2Hn max ∆t+β+2ηHn max.(30) It only remains to estimate the constants γ1and γ2. Following [12, Sect. 5.6], the Korn constant can be simply taken as CK=√2. Assuming that the domain Ωis convex, it is known (see [18]) that the optimal choice for CP is d/π, where dis the diameter of Ω. Thus, γ1=π/ √2d. On the other hand, we have that kVkV≤√3k∇VkL2≤ √3˜γ2kVkL2for a certain constant ˜γ2, so γ2=√3˜γ2. Reasoning as in [10, Appendix C], ˜γ2can be taken as √µmax, where µmax is the maximum eigenvalue of the discrete Laplacian problem; in general, this value has to be computed numerically. Remark 2. In the particular case of a rectangle Ω = [x10,x10+Lx1]×[x20,x20+Lx2]with a uniform discretization (∆x1,∆x2),γ1would be γ1=π q2(L2 x1+L2 x2) , while γ2could be taken as γ2=s3π21 ∆x2 1 +1 ∆x2 2. Note that very recently the FISTA method [23] was introduced for the simulation of viscoplastic flows. It is inspired by proximal gradient methods and allows to speed-up computations, compared to the non-optimized augmented Lagrangian. FISTA could be a complementary approach to the present BM method which has an automatic computation of the optimal duality parameter while giving the same quality results of plastic zones, as the long proven AL method [21]. 4. Discretization in space and Well-Balancing in 2D In this section, we define the spatial discretization for the conservation equation (14) and the velocity equation associated to the iterative algorithm of the AL (equation (20)) and the BM (equation (27)). As mentioned in the companion paper [10] in 1D, there is a rather subtle coupling between both equations through the well-balanced property of the global scheme. This need to be carefully extended when going to the 2D framework. 4.1. Definitions Note that both equations (20) and (27) can be written under the same structure. Namely, for a time t=tnand known the velocity at iteration n,Vn, the common system is HnVk+1 ∆t+βVk+1−divx((2ηHn+δn)D(Vk+1)) +∇x((2ηHn+δn)divx(Vk+1))= HnVn ∆t−Vn·∇xVn+(fΩ+fz∇x(b+Hn))+divx(HnΠk)+∇x(tr(HnΠk)),(31) where: •for the AL method: δn=rHn,Πk=µk−rqk; (32) •for the BM method: δn=ωn,Πk=θk/Hn,(33) where ωnis defined by the optimal value (30), in terms of Hn. 8 Let us suppose that the domain Ωis a rectangle, Ω = [x10,x10+Lx1]×[x20,x20+Lx2]. Let us consider a partition of Nx1intervals of length ∆x1=Lx1/Nx1along x1; and another partition with Nx2intervals of length ∆x2=Lx2/Nx2 along x2. The 2D mesh is then defined by the union of control volumes {Ki,j}j=1,...,Nx2 i=1,...,Nx1 ,where Ki,j=[x1,i−1/2,x1,i+1/2]×[x2,j−1/2,x2,j+1/2],i=1,...,Nx1,j=1,...,Nx2, with x1,i+1/2:=x10+i∆x1,i=0,...,Nx1, x2,j+1/2:=x20+j∆x2,j=0,...,Nx2. Figure 2: Notations for the discretization To approximate Hnand Vk, solutions of the semi-discrete system defined by (14)-(31), we consider a finitevolume solver. Then, let us denote at (x1,i,x2,j), center of Ki,j, Hn i,j≈1 Ki,jZKi,j Hn(x)dx,Vk i,j≈1 Ki,jZKi,j Vk(x)dx. The duality multiplier Πkis approximated at the vertices of the partition, then let us denote (see figure 2) Πk i+1/2,j+1/2≈Πk(x1,i+1/2,x2,j+1/2). By the definition of Πk, equations (32)-(33), we denote Πk i+1/2,j+1/2=(µk i+1/2,j+1/2−rqk i+1/2,j+1/2,for AL, θk i+1/2,j+1/2/Hn i+1/2,j+1/2,for BM, and Hn i+1/2,j+1/2=(Hn i,j+Hn i+1,j+Hn i,j+1+Hn i+1,j+1)/4.(34) In order to discretize in space the system (14)-(31), we consider a well-balanced finite volume method defined in terms of a diagonal viscosity matrix. This implies that we can present the discretization of the system equation by equation. Nevertheless, their discretizations are not really decoupled because, in order to obtain a well-balanced property, it is necessary to take into account the definition of the duality multipliers Πkin the approximation of the mass conservation equation. Then, we first present the discretization of equation (31) and, second, the discretization of (14). .Discretization of the velocity equation (31) associated to the iterative algorithm Equation (31) is approximated as follows: Hn i,j ∆t+β!Vk+1 i,j−Dk+1 i,j=Hn i,jVn i,j ∆t−1 ∆x1∆x2∆x2(Fn− i+1/2,j+Fn+ i−1/2,j)+ ∆x1(Fn− i,j+1/2+Fn+ i,j−1/2)+Ek i,j,(35) where Fn± i+1/2,j= (V1)n i,j+(V1)n i+1,j 2(Vn i+1,j−Vn i,j)−fΩ ∆x1 4+fz 2(bi+1,j+Hn i+1,j−bi,j−Hn i,j) 1 0!±Si+1/2,j 2(Vn i+1,j−Vn i,j), Fn± i,j+1/2= (V2)n i,j+(V2)n i,j+1 2(Vn i,j+1−Vn i,j)−fΩ ∆x2 4+fz 2(bi,j+1+Hn i,j+1−bi,j−Hn i,j) 0 1!±Si,j+1/2 2(Vn i,j+1−Vn i,j), 9 !!"# !!!$"# !$!%"# !%!&"# & $ ! ' # ( ) *% + , , -.//,01.234/ 567768 !,!,9:3;/ Figure 3: Left: slice of the bottom and the free surface at x2=0.5 in global coordinates (meaning that x1is obtained via the rotation associated to the angle α, see Fig. 1). Right: bottom b(x) (in black) and initial condition b+H(blue), in local coordinates. If we consider the initialization of the multiplier defined in Theorem 1, the stationary solution is preserved up to machine precision for a value of τyverifying condition (48). For sake of conciseness we do not show the illustrations here. In the present test, we instead initialize the multipliers to zero. The simulation is done from t=0 to 1. Actually, at the first time iteration, the multipliers converge (inside the duality loop) to some function which then remains unchanged along the subsequent time iterations. As in 2D the multiplier is (intricately) not uniquely defined, it is not assured that the iterative algorithm converges to the one defined by (45)-(46), which ensures the exact wellbalanced property. Nevertheless, we can see numerically that both AL & BM schemes still preserve the stationary solution with a good accuracy. In Figure 4, the multiplier to which the algorithm converges is plotted for illustration (this is done with the BM method but the results are the same with the AL algorithm). We obtain that the time averaged L2error for t∈[0,1] is 6.3×10−4for Hand 1.4×10−8for the velocity norm. The averaged difference between H(x,t) and the initial condition, for t∈[0,1], is represented in Figure 5, together with the averaged norm of the velocity (also for t∈[0,1]). Figure 4: Multiplier at t=1. Left: Πk 11, center: Πk 12, right: Πk 22, being kthe last iteration of the duality algorithm. 5.2. Well-balanced test on academic avalanche This test is a dam break simulation where the Ω-plane is sloping at α=20◦and the bottom with two obstacles is defined as follows on [0; 1]2(cf. Figure 6): b(x)=1.5e−[20(x1−0.75)]4+6e−[5(x1−1.25)]4+3e−[10(x1−0.5)]2−[18(x2−0.5)]2+10(x2−0.5)2+11.1e−0.9x1+1.1e−9x1.(53) As initial condition, we set V≡0 and (see Figure 6) H(x)=(10 −b(x) if (x1,x2)∈[0.7; 0.89] ×[0.4; 0.6], 0 otherwise .(54) 16 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 0 0.001 0.002 0.003 0.004 0.005 0.006 0.007 0.008 0.009 0.01 ! ! " "#$ "#% "#& "#' ( " "#$ "#% "#& "#' ( "#$ "#% "#& "#' ( (#$ (#% (#& (#' $ $#$ )!("!* Figure 5: Left: Averaged error, |H(x,t)−H(x,0)|for t∈[0,1]. Right: identically, averaged norm of the velocity for t∈[0,1]. Figure 6: Left: slice of the bottom b(x1,0.5) with the plane of reference inclined with an angle of 20◦, in global coordinates (meaning that x1is obtained via the rotation associated to the angle α, see Fig. 1). Right: bottom b(x) (in black) and initial condition b+H(blue), in local coordinates. We set η=10−3m2.s−1,τy=√2/2m2.s−2,β=10−3m.s−1and g=9.81 m.s−2. Even though bis defined analytically and it is an academic avalanche, this test is very demanding due to the high slope and the strong gradients of the bottom as well as the quantity of material in the initial column (note the strong aspect ratios in Figures 6left and right). With these values of the parameters, the material reaches a stationary state around t=1s, so we made simulations up to t=6sto check the ability of the 2D scheme to preserve this rigid free surface. This test has a rich hydrodynamics, as shown on Figures 8,9and 10. Recall that V=0 on ∂Ω. In the first phase of the collapse of the column, there are reflections (principally in the x1direction) on the wall at x1=1 and on both obstacles (ridge at x1=0.75 and Gaussian at (x1,x2)=(0.5,0.5)) inside the domain (see t=0.05 s). This notably leads to counterwaves which collide on the ridge (see t=0.09 sand t=0.11 s). Then, the material essentially separates in two parts on each side of the ridge and then reaches a steady state. In the upper part the material oscillates for a short time (feeding a bit the other side of the ridge), cf. t=0.21 and 0.26 son Figure 9. The lower part of the material separates and goes around the Gaussian to finally meets at the bottom of the hill and reaches a stationary shape with V1 and Hwith a shape with very high gradients (after t≈1s, see Figures 10 and 11). It is very difficult to compute a stationary solution in this kind of configuration with wet/dry fronts but we see in Figure 11 that our scheme performs well to do so: the velocity is very small (kVk2(Ω)=1.5e−9) and the level lines {x|H(x)=0}are indeed very well superposed between t=5 and 6 s. Note that these results are computed on a mesh with 4502points and can be considered as converged in terms of spatial resolution. Indeed, Fig. 11d shows a mesh convergence study of the level lines {x|H(t=6,x)=0}for various ∆x. It can be seen that the results for the two more refined meshes are very close. 17 We also use this test to compare the AL and the BM methods in terms of numerical cost. Note that following our 1D study for BM [10, Section 3.1.3], we directly present the wet/dry front optimized version, $opt, of the optimal choice of the ωparameter. We proceed as follows. For both methods, we simulate from t=0 to t=1, with 1102mesh points for the Ωsquare, and we study the cost of the duality loops as a function of the duality parameter (we use 14 discrete values between 0.01 and 64). This cost is defined as the sum, along all the iterations in time (tn), of the iterations of each duality loop used to compute Vn. We perform this study for four values of τy=√2/20,√2/2,2√2 and 5 √2m2.s−2, which are representative of 4 different dynamics of this 2D test case. Remark that we limited the number of iterations in a dual loop at 10,000 iterations: in practice, it is not reached except for the BM method at τy=√2/20 m2.s−2. Consequently this does not change the following conclusions but allows to perform this study in a more reasonable CPU time. The results are presented in Figure 7. Recall that the duality parameter is taken as a constant for all the time iterations for the AL and the standard BM. However, when optimal BM is used with the a priori derived $opt =$opt(t), one can not give a meaning to the cost for a given $and there is only one value of the duality cost: this leads to the horizontal lines in Fig. 7. We can see that, for this dam break problem, when τyincreases the duality cost decreases for both AL and BM methods. Remark that the BM curve for τy=√2/20 m2.s−2is far less convex than the other curves because, at some time iterations, it reached the maximum number of duality iterations (=10,000) as mentioned above. But, still, we can observe an optimal value of ω, which is furthermore close to the $opt estimation. The AL seems to always be cheaper than the BM, especially at small τy. However, it can be seen that $opt always leads to a good estimation of the observed optimal cost and that when τyincreases, this cost of the BM is closer to the cost of the AL. As a consequence, the optimal BM method can be very competitive w.r.t. the AL, especially at high τy. Indeed the optimal rfor the AL method is not known a priori and the practitioner can be far from it, leading to significantly higher CPU times. 10-2 10-1 100101102 r , ω 104 105 106 107 108 Sum of all iters in dual. loop AL & BM - Study of the cost for various viscoplasticity τy AL τy = p 2 / 20 AL τy = p 2 / 2 AL τy = 2 p 2 AL τy = 5 p 2 BM τy = p 2 / 20 BM τy = p 2 / 2 BM τy = 2 p 2 BM τy = 5 p 2 Figure 7: Duality numerical cost for four τy. The colored continuous curves are for the augmented Lagrangian while the dashed ones are for the standard Bermúdez-Moreno . The horizontal thick lines correspond to the cost for the BM with the optimal duality parameter $opt : their colors (varying with τy) correspond to the same colored dashed curve of the standard BM to which it principally needs to be compared; namely τy=√2/20,√2/2,2√2,5√2m2.s−2is in blue, green, red, cyan, respectively. 18 (a) t=0.03 (b) kVk2(x), t=0.03 (c) t=0.05 (d) kVk2(x), t=0.05 (e) t=0.09 (f) kVk2(x), t=0.09 Figure 8: Left: Free surface (blue) and bottom (black). Right: contours of kVk2(x). From t=0.03 to 0.09. 19 (a) t=0.11 (b) t=0.11 (c) t=0.21 (d) t=0.21 (e) t=0.26 (f) t=0.26 Figure 9: Left: Free surface (blue) and bottom (black). Right: contours of kVk2(x). From t=0.11 to 0.26. 20 (a) t=1 (b) t=1 Figure 10: Left: Free surface (blue) and bottom (black). Right: contours of kVk2(x). At t=1. (a) Free surface at t=6: rotated view to better see high gradients of Hat the wet/dry front. (b) Contours of kVk2(x) at t=6. Note: kVk2(Ω)=1.5e−9. (c) Another stationary evidence. Square mesh with 4502points. Level line {x|H(x)=0}for different times: dataifor i=1 to 11 stands for tfrom 5 to 6 with a time step of 0.1. (d) Mesh refinement study: level line {x|H(x)=0} at t=6 for different mesh sizes. DX0,1,2 and 3 stands for a square grid discretized with respectively 752,1502,3002and 4502mesh points. Figure 11: Details on the stationary state at t=6 and mesh convergence study. 21 5.3. Taconnaz avalanche path, Chamonix - Mont-Blanc (a) Photo of the upper part of the site (Dome du Gouter on the far left and Gros Bechard on the right), courtesy [16]. Note the significant amount of ice and complexity of the terrain. (b) Topography from ASTER GDEM: 431 ×213 mesh resolution. Dome du Gouter approximately at (6500,1600) and Gros Bechard at (3500,1000). Figure 12: Topography of the Taconnaz avalanche path, Chamonix - Mont-Blanc. In this section, we test the ability of the 2D numerical scheme to simulate an avalanche on a real topography. Namely, we choose the Taconnaz avalanche path in the region of Chamonix, France. This site is one of the longest in Europe with a length close to 7000 m (avalanches can start around 3300 m above sea level and stop around 1000 m a.s.l.), a width between 300 and 400 m and a mean slope of 25◦(with some portions in departure areas of avalanches of mean slopes 30◦). Taconnaz is well known for a significant frequency of avalanches (composed of dense and mixed snow, with speed of 70 m/s in the worst case scenario), with 75 events between years 1900 and 2000 [15]. We obtain the topography of the Taconnaz avalanche path thanks to the ASTER Global Digital Elevation Model (GDEM) v2, whose initial resolution in x1and x2is around 25-30 m [3]. For simulation purposes, we interpolate the topography on a finer grid: we built a uniform square mesh with a 16 m resolution (431×213 points in x1×x2), see Fig. 12. Based on historical observations [15], we put on top of this topography a truncated Gaussian of material for Hwhose maximal height is 9 m and volume is 0.6×106m3(observed volumes are between 0.01×106 and 1.5×106m3) at an altitude of 3700 m a.s.l. on the slopes of Dome du Gouter, see Figures 13a and 13b. Namely, H(t=0,x)=max 0,−2+11e−(4.10−5)(x1−6380)2−(2.10−5)(x2−1050)2.(55) This is used as an initial condition to represent the dense snow composing the avalanche; further V(t=0) ≡0. For the material, we set η=10−1m2.s−1,τy=√2m2.s−2,β=2.10−3m.s−1,g=9.81 m.s−2. Of note, for a given real observed avalanche, it is very difficult to give precise values of these parameters so we put reasonable values which lead to an observed deposit of the avalanche in the field (see [6]); in particular, we do not enforce that physical time scales are relevant, focusing only on the localization of the deposit. The objective is here to show that algorithms derived in this paper are applicable on real avalanches data to compute the stopping state. The fitting of these parameters is out of the scope of this paper and is left for future works. The dynamics of this test, which spans from t=0 to 120,000 scan be decomposed in 3 phases, going to the stationary state. A first "fast" phase on a "short" time scale (t=0 to approx. 3000 s) where the deposit reaches the bottom of Taconnaz path. The front is not stationary but it is not far from its stationary localization, see Figure 14. In a second phase, on a longer time scale (from approx. t=3000 sto 40,000 s), the velocities are decreasing but there is a significant motion of the material in the whole deposit from the mountain top to the bottom: this leads to a progressive advance of the front of the avalanche. In a last "slow" phase on a much longer time (from approx. t=40,000 sto 106,000 s), there is essentially no motion on the top 3/4 of the deposit and most of the material is in the 1/4 bottom part where the slope is still significant (see Hon Fig. 18 (b) Bottom): as a consequence there is a slow but progressive sliding motion of this bottom part (as also shown by the time evolution of maxx∈ΩkVk2(t,x) on Figure 16) and the deposit front is moving a bit (compare Figs. 15 and 17) but finally stops at about t=106,000 s. We show that the front {x|H(x)=0}is completely stationary after t=106,000 sby also showing the solution 22 at t=120,000 s(Fig. 17). Note that the final shape of the deposit is very close to one of the biggest deposits shown in [6] and measured from true avalanches at Taconnaz. Our test thus covers all the topographical difficulties associated to the Taconnaz avalanche path. It can be seen that the stationary state is very well computed: the final velocity is locally of order 10−10, and globally kVk2(Ω)≤4.6×10−9, see Figure 17. The position of the wet/dry front is shown to be stationary with superimposed level lines {x|H(x)=0}after t=106,000 swith a very good accuracy (it does not move up to t=120,000 s), see Figure 19. Note that the stationary state is difficult to capture since the major part of the deposit accumulates in a zone where there is a significant slope of b, see Figure 18. The viscoplastic nature of the material with a bumped surface of the deposit (H) is clearly exhibited in this stationary state. These results show the ability of present well-balanced schemes to perform accurate simulations for Bingham type materials with real topographies from digital elevation models (DEM). (a) topography b(x) (black) and free surface H(brown). (b) filled contours of b(top) and Hat t=0s(bottom). Figure 13: Details on topography and initial condition of the simulation on Taconnaz avalanche path. Note that Fig. 18b gives also btogether with its gradient. 6. Conclusions In this article, we presented 2D numerical schemes in the finite volume framework which allow to compute accurately shallow viscoplastic flows: thanks to a specific design coupling duality methods and well-balancing, they preserve with a good accuracy the stationary solutions (naturally associated to the viscoplasticity) on general 2D shapes of bottom and free surfaces. These schemes deal with true wet/dry fronts and there is no need to add a small quantity of material in all the domain (as sometimes done by other methods). The well-balanced property is shown to be exact on two kinds of stationary solutions (Theorem 1). A careful study of the optimal cost of the two duality methods (Augmented Lagrangian and Bermúdez-Moreno) was performed and showed that the BM method can become competitive at high τydue to the fact that the optimal duality parameter is known a priori. Such studies are quite rare in the 2D framework. We finally give numerical evidence that these numerical methods can be successfully applied to real topographies as shown by the avalanche test case in the Taconnaz path obtained from ASTER GDEM. As a by-product of this study, we also provide a 2D benchmark for classic 2D Bingham codes thanks to an analytic solution with totally non-homogeneous boundary conditions on the velocity. Acknowledgments This research has been partially supported by the Spanish Government and FEDER through the Research projects MTM2012-38383-C02-01, MTM2012-38383-C02-02, MTM2015-70490-C2-1R and MTM2015-70490C2-2R. Part of this work was done while P. V. was visiting E.D. F.-N. and J.M. G., during a stay in 2013 thanks to a grant from the Instituto Universitario de Investigación de Matemáticas de la Universidad de Sevilla (IMUS). 23 Figure 14: First times of the avalanche between t=76 and t=5596: topography b(black) and free surface H(brown). See also Fig. 15 for the corresponding velocities. At t=5596, the red arrow shows the localization of the zoom made on Fig. 18a. P. V. wishes to thank everyone at IMUS for their hospitality. A visit of E.D. F.-N. was supported in 2015 by the LABEX MILYON (ANR-10-LABX-0070) of Université de Lyon, within the program "Investissements d’Avenir" (ANR-11-IDEX-0007) operated by the French National Research Agency (ANR). The support of French ANR Grant ANR-08-JCJC-0104, TELLUS Grant from CNRS INSU-INSMI (2016 Call) and InFIniTi Grant from CNRS (2017 Call) is also gratefully acknowledged. 24 Figure 15: First times of the avalanche between t=76 and t=5596: kVk2(x) (filled contours) and level line {x|H(x)=0}(white thick line). The colorbar is the same on all snapshots so we also give, at the bottom right of the figure, the corresponding maximum values of kVk2(x) as a function of time. Note: at t=5596, maxxkVk2(x)=4.22, see also Fig. 16 for the total history. See also Fig. 14 for the corresponding 3D views of band H. Figure 16: History of Vconverging to the stationary state for the Taconnaz test, in semi-log scale. 25