Efficient numerical schemes for viscoplastic avalanches. Part 1: The 1D case.

Efficient numerical schemes for viscoplastic avalanches. Part 1: the 1D case E.D. Fern´andez-Nieto 1, J.M. Gallardo 2, P. Vigneaux 3 January 13th, 2014 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 to the yield threshold: finitevolume schemes are used together with duality methods (namely Augmented Lagrangian and Berm´udez-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. We derived such schemes and numerical experiments are presented to show their performances. Contents 1 Introduction 2 2 Model and resolution approaches 4 3 Treating the velocity inequality with two duality methods 7 3.1 StudyoftheBMapproach.............................. 8 3.1.1 TheBMalgorithm .............................. 8 3.1.2 Study of the optimal parameter . . . . . . . . . . . . . . . . . . . . . . . 9 3.1.3 Some remarks on the treatment of wet/dry fronts . . . . . . . . . . . . . 10 4 Well-balanced discretization with general bed and wet/dry fronts 10 4.1 The global coupled scheme . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15 4.2 Well-balanced properties . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17 5 Numerical comparison of both approaches 18 5.1 Well-balancedtests.................................. 19 5.2 Avalanchewithobstacle ............................... 24 1Departamento de Matem´atica Aplicada I, Universidad de Sevilla. E.T.S. Arquitectura. Avda, Reina Mercedes, s/n. 41012 Sevilla, Spain ([email protected]) 2Departamento de An´alisis Matem´atico, Universidad de M´alaga. F. Matem´aticas, Campus Teatinos S/N, Spain ([email protected]) 3Unite de Mathematiques Pures et Appliquees, Ecole Normale Superieure de Lyon. 46, Allee d’Italie, 69364 Lyon Cedex 07, France ([email protected]) 1 6 Conclusions 25 Appendix 34 A Derivation of the AL approach 34 B Derivation of the BM method 36 C Study of the BM optimal parameter 38 D The coupled scheme in the case of the AL 39 E Proof of Theorem 1 40 F A duct flow case 42 1 Introduction The objective of this paper is to present improved numerical schemes for shallow water models for viscoplastic materials on variable topography (or bathymetry). The associated difficulties are twofold. First, we will place ourselves in a context –increasingly in use– where no regularization method is used and thus the variational inequalities traducing the plastic nature of the material are handled directly through duality methods. Second, we will describe new well-balanced schemes in this viscoplastic context to take into account general bottoms and wet/dry fronts. In recent years, an increasing interest has been developed for shallow water models in the context of simulations for the flow of viscoplastic materials down inclined planes. 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. Such flow behaviour can be encountered in many practical situations such as food pastes, heavy oils, lavas and avalanches. As a consequence, the theory of the fluid mechanics of such materials has applications in a wide variety of fields such as chemical industry, energy industry and geophysical fluid dynamics. From the mathematical viewpoint, the non-linearity associated to viscoplastic models (such as the Bingham model, as we will see below), leads to feasible but very expensive computational times for the full 3D equations (see e.g. [28]). Consequently, numerous reduced 2D model using the shallow flow approximation have been derived. In the context of avalanches, we refer to the article of Ancey [2] and references therein for a detailed review on these developments. Recently, in [25], an interesting shallow water model based on a Bingham-like constitutive law together with Coulomb frictional condition on the bottom was derived, in local coordinates for the case of a non-planar topography. But the algorithm, presented to solve associated equations, does not take into account either well-balanced properties or the treatment of wet/dry fronts. Another shallow model based on the Herschel-Bulkley constitutive law (which generalizes the Bingham law) was derived in [1] and a new well-balanced scheme was introduced to take into account both non-linearities of this constitutive law, leading to a scheme which preserves more accurately stationary states. This point is important when it comes to determine the stopping time of the flow, when the material enters in its rigid state. And this kind of property is also linked to the use of duality methods which allow to properly deal with the plasticity. Indeed, a common point of an increasing part of the recently developed numerical methods for viscoplastic flows (see e.g. [33], [37, 36], [24], [35], [25]) is that they use decomposition- 2 coordination methods to solve the variational inequality associated to constitutive laws with a so-called plastic threshold (the most simple and iconic one being the Bingham model). This kind of approach takes its roots in the seminal works of Duvaut-Lions [15] and the series of papers of Glowinski and coworkers (see the recent book [22] for a detailed review), initiated at the end of the seventies and anchored in the Augmented Lagrangian formalism. One of the crucial advantage of these methods over regularized approaches (see e.g. Papanastasiou’s [29] or the so-called bi-viscosity methods [14]) is that they rigorously take into account the plastic threshold. Of note, it is well known that in Augmented Lagrangian (AL) methods, the optimal values of the parameters are not easy to determine in the general case. These parameters ((r, ρ) in Glowinski’s nomenclature) influence the speed of convergence of the iterative process towards the saddle-point, solution of the problem. As a consequence, a study of some sort of optimality for such parameters is of real interest when it comes to improve the computational efficiency. As an alternative duality method, one can use the so-called Berm´udez-Moreno (BM) algorithm. This method, which is built upon some properties of the Yosida regularization of maximal monotone operators, has been extensively used for a wide range of applications (see [19] and the references therein). In order to apply the method, the Yosida regularization of the subdifferential associated to the non-differentiable operator appearing in the formulation of the considered model needs to be determined. As for the AL, the performance of the algorithm strongly depends on the choice of two constant parameters. Fortunately, several ways to overcome this problem have been proposed in the literature ([31, 30, 19]), and they will be considered in this paper. Another difficulty that appears when it comes to couple shallow-water models and viscoplastic constitutive laws, is the adequate coupling between the discretization associated to the duality method and the one associated with the spatial terms, in such a way the global scheme is wellbalanced. For shallow water type equations, finite volume methods have proved their efficiency and we adopt them in the present work. In this context, a careful treatment must be made to design a well-balanced scheme when coupling the finite volume scheme and the duality method. This idea was first introduced in [8], in the context of a Bingham fluid treated with an AL method. We extend here this idea for a Shallow-Water-Bingham model on a general topography and in the presence of wet/dry fronts. The well-balanced properties are related to the stationary solutions of the system. In our case, we seek numerical schemes which preserve exactly two types of stationary solutions. For hyperbolic systems with source terms, a discretization of the source terms compatible with the one of the flux term must be performed. Otherwise, stemming from the numerical diffusion terms, a first order error in space takes place. This error, after time iteration, may yield large errors in wave amplitude and speed. The pioneering work by Roe [32] relates the choice of the approximation of the source term with the property of preserving stationary solutions. Berm´udez and V´azquez-C´end´on introduced in [5] an extension of Roe’s solver, in the context of shallow water equations, which preserves exactly the stationary solution of water at rest. This work originated the so-called well-balanced solvers, in the sense that the discrete source terms balance the discrete flux terms when computed on some (or all) of the steady solutions of the continuous system. Different extensions have been done: see for instance Greenberg-Leroux [23], LeVeque [27], Chac´on et al. [12]. An additional difficulty in the simulation of free surface flows comes from the appearance of dry areas in the computational domain, due to the fluid evolution or to the initial conditions. Standard numerical schemes may compute spurious solutions in the presence of wet/dry fronts, unless appropriate modifications are made. See [7, 34] for a review on some methods appearing in the literature to deal with this problem. Moreover, in the context of shallow water equations, Roe-type schemes lose their well-balanced properties when wet/dry transitions appear. 3 Indeed, they may produce nonphysical negative values of the thickness of the water layer near the wet/dry front. Some ways to modify Roe’s method to fix these problems have been proposed in [10, 11]. The contribution of this paper is twofold. First, by adapting the guidelines in [30], we determine, in a theoretical way, an optimal choice of parameters in the sense that they provide the highest rate of convergence for the BM algorithm. For the AL, we perform some numerical studies of the optimal choice of parameters and we then compare both methods on various problems to give insight on their respective behaviours. To our best knowledge this is the first time that BM is applied to such kind of models and that such a systematic comparative study of the behaviour of the number of iteration in duality methods is done for several very different viscoplastic flows. Second, both for the AL and BM methods which are embodied in a general framework, we design a well-balanced scheme which takes into account wet/dry fronts on general bottoms, for such viscoplastic free surface flows. Again, this is the first time that dry area treatment is proposed for flows with plastic behaviour. This is crucial when it comes to study real applications where there are always a flow with wet/dry front and a rigorous treatment is needed to compare qualitatively numerical simulations and physical experiments. This work is organized as follows. In Section 2, we introduce a typical model for viscoplastic free-surface flows on general bottoms and the general resolution approach. Then in Section 3, we present in detail the two duality methods used to treat the viscoplastic behaviour of the material, namely the AL and the BM methods. In Section 4, we describe the design of the overall well-balanced scheme which takes into account general bottoms and the presence of wet/dry fronts. Numerical tests are finally presented in Section 5 to illustrate the various properties of the scheme and compare both methods. A duct flow test and a stationary test allow to make a convergence analysis thanks to the availability of analytic solutions. Then, several tests of avalanches, academic but very challenging from the numerical viewpoint, are performed to show the robustness of the scheme. Conclusions are drawn in Section 6. 2 Model and resolution approaches As a model problem for viscoplastic shallow flows, we use the Bingham shallow-water model derived in [8], but we add here the fact that the bottom is more general than a plane slope. The derivation being very close, we refer to [8] for more details and here we directly present the resulting model. Its physical characteristics and relevance are briefly described in the following. The geometry is as shown on Figure 1. We consider a fluid domain of height 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 4 b(x) H(x) z x α Figure 1: Domain in 2D (left) and 1D (right). 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. [15]) 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) 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) + V·∇xV(Ψ −V)dx +ZΩ βV ·(Ψ −V)dx +ZΩ 2ηHD(V) : D(Ψ −V)dx +ZΩ 2ηHdivxV(divxΨ−divxV)dx +ZΩ τyHp|D(Ψ)|2+ (divxΨ)2−p|D(V)|2+ (divxV)2dx ≥ZΩ H(fΩ+fz∇xb)·(Ψ −V)−ZΩ (H)2 2fz(divxΨ−divxV)dx, (4) where D(U) := 1 2∇xU+ (∇xU)t,(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) Of note, usual Sobolev embeddings and the fact that His bounded allow aforementioned problem to be well-posed. As said previously, this model is based on a so-called shallow water approximation, i.e. the height of the fluid is assumed to be much smaller than the characteristic 5 length of the domain. Furthermore, this model was derived through an asymptotic expansion where the slope is supposed to be small (α1) and the norm of the gradient of b(x) is small (k∇xbk  1). But it is worth noting that this model is also valid for a slope α= 0 (horizontal bottom), which is not generally the case for other models proposed in the literature (see for example [3], [17]). Another interesting feature of the model is that in the case of a plane horizontal slope (α= 0) and with a vanishing yield stress (τy= 0), we recover a viscous shallow water system which has the same structure as the one derived by Gerbeau and Perthame in [20] (note that the hypothesis of friction at the bottom, instead of a no-slip condition is a key point in this degeneracy to [20]). The shallow water formulation (4) is in weak form. It can be rewritten in the strong form to have the expression of the associated (integrated) Bingham constitutive law. Namely, the corresponding formulation is H∂tV+V·∇xV+βV −divx(Hσ) = H(fΩ+fz∇xb)−∇x(H)2 2fz,(8) where          σ= 2η(D(V) + tr(D(V))I) + τy D(V) + tr(D(V))I p|D(V)|2+|tr(D(V))|2if |D(V)| 6= 0 |σ| ≤ τyif |D(V)|= 0, (9) the second invariant |σ|of a tensor σbeing defined here as |σ|:= qΣi,jσ2 ij.(10) 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-a), then the tilted axis (with respect to the horizontal) is x1, i.e. fΩ= (−gsin α , 0), fz=−gcos α. (11) In this paper, we will fulfill our objectives on the 1D version of (2)-(3)-(4). Of note, the ideas presented here can be extended in 2D. This extension being not trivial, it will be presented in another article. In the 1D case, we naturally take x∈Ω = [0, L] and the associated frame of reference becomes (x, z) (cf. Figure 1-b). We have then V(t) = {Ψ∈H1(Ω) /Ψ = 0 on ∂Ω}:= H1 0([0, L]),(12) and the problem (3)-(4) degenerates to: find H∈L2([0, T], L∞([0, L])) and V∈L2([0, T],V(t)), with ∂tV∈L2([0, T], L2([0, L])), such that ∂tH+∂x(HV ) = 0,(13) and ∀Ψ∈ V(t),ZL 0 H∂tV(Ψ −V) + 1 2∂x(V2)(Ψ −V)dx +ZL 0 4ηH∂xV ∂x(Ψ −V)dx +ZL 0 τy√2H|∂xΨ|−|∂xV|dx +ZL 0 βV (Ψ −V)dx ≥ZL 0 H(fΩ+fz∂xb)(Ψ −V)dx −ZL 0 H2 2fz(∂xΨ−∂xV)dx. (14) 6 The gravity becomes (note that, as the velocity, the projection of the force on Ω is now a scalar): fΩ=−gsin α, fz=−gcos α. (15) In this paper, as in many schemes used in the literature for this type of models, we will consider a first order backward semi-discretization in time (explicit Euler’s method). We denote the time step by ∆t. We have: Hn+1 −Hn ∆t+∂x(HnVn) = 0,(16) and ZL 0 HnVn+1 −Vn ∆t(Ψ −Vn+1) + 1 2∂x((Vn)2)(Ψ −Vn+1)dx +ZL 0 βV n+1(Ψ −Vn+1)dx +ZL 0 τy√2Hn|∂xΨ|−|∂xVn+1|dx +ZL 0 4ηHn∂x(Vn+1)∂x(Ψ −Vn+1)dx ≥ZL 0 HnfΩ+fz∂xb(Ψ −Vn+1)dx −ZL 0 (Hn)2 2fz(∂xΨ−∂xVn+1)dx, ∀Ψ.(17) Doing so, we see that problem on the height and problem on the velocity are decoupled. At each time step, supposing that we know (Hn, V n), we need to solve both problems for (Hn+1, V n+1). One of the goals of this article is to compare two duality methods to handle the variational inequality of the problem on the velocity, namely the augmented Lagrangian method and Berm´udez-Moreno method. It is the subject of the next section, where we will also see more clearly the underlying shallow water nature of this system. For this kind of equations, finite volume discretizations are particularly well suited and this is why we want to use them in this context. The other goal of this paper is to show that a careful design of the scheme is needed to obtain a well-balanced property in the case of variable bottom and in the presence of wet/dry fronts. This will be treated in Section 4. 3 Treating the velocity inequality with two duality methods The speed problem of the above shallow Bingham model is a variational inequality and consequently specific methods to solve it are needed. Various types of methods exist in the literature and, in this paper, we choose two of them which have proved to be efficient for such problems (see e.g. [26]), namely the augmented Lagrangian method and the Berm´udez-Moreno method. Their definitions and derivations are different, but interestingly the obtained structure of the algorithms is the same. We will thus make the most of this fact in Section 4, to design an unified scheme for the space discretization. In addition, this common structure in terms of code implementation makes it interesting to compare their efficiency and computational cost. In particular, both algorithms depend on a parameter which influence the speed of convergence to the solution. Consequently, it is interesting to determine if optimal parameters can be derived theoretically. The Augmented Lagrangian algorithm is a method which has been extensively used for a wide range of non linear optimization problems (see [21]). It appears that (17) can be reformulated as 7 a minimization problem and solving for Vn+1 ends up to find the saddle-point of an augmented Lagrangian. This approach allows to elegantly deal with the non-differentiable terms of (17). For sake of brevity, the full derivation of the Augmented Lagrangian algorithm for (17) is given in Appendix A. Of note, the final form of the present Augmented Lagrangian algorithm has a structure which is very similar to the one derived the following subsection and Appendix A is also here for completeness in comparing the two approaches. Concerning the so called augmentation parameter r, it is generally difficult to find a priori an optimal value which leads to the fastest convergence of this iterative method. Some optimal parameters are derived in [18] for very simple model problem and, in the context of seismic reflection tomography, a heuristics was proposed to evolve ralong the iterations of the AL to increase the speed of the algorithm [13]. But overall, a general methodology to find theoretically an optimal value of rseems to remain an open question. This fact was one of the reasons to study in more details another duality method, namely the Berm´udez-Moreno method. 3.1 Study of the BM approach For the sake of brevity, the full derivation of the Berm´udez-Moreno method is given in Appendix B and we will directly give the resulting algorithm in the following. It is important to note that it is as easy to implement as the Augmented Lagrangian algorithm. 3.1.1 The BM algorithm Let us summarize the BM algorithm for the speed problem (17). Berm´udez-Moreno algorithm •Initialization: suppose that Vn,Hnand θnare known. For k= 0, we set Vk=Vnand θk=θn. •Define the parameters λand ω(See Section 3.1.2). •Iterate: –Find Vk+1 ∈ V solution of the following linear problem: Hn ∆t+βVk+1 −∂x((4ηHn+ω)∂xVk+1)−∂x(ω∂xVk+1) =Hn ∆tVn−Hn 2∂x((Vn)2) + 1 2∂x((Hn)2fz) + Hn(fΩ+fz∂xb) + ∂xθk.(18) –Update the so called Berm´udez-Moreno multiplier θk+1 via ξk+1 =∂xVk+1 +λθkand θk+1 =                −ω ξk+1 +τy√2Hn(x) 1−λ ω if ξk+1 > λτy√2Hn(x), ξk+1 λif ξk+1 ∈[−λτy√2Hn(x), λ τy√2Hn(x)], −ω ξk+1 −τy√2Hn(x) 1−λ ω if ξk+1 <−λτy√2Hn(x). (19) Note that this computation is local in space, i.e., it is done at each discretization point. 8 –Check convergence (see below) and update: Vk=Vk+1,θk=θk+1,k7→ k+ 1 and go to the next iteration... •... until convergence is reached: kθk+1 −θkk kθkk≤tol.(20) At convergence, we get the value of Vn+1 by setting Vn+1 =Vk+1 (in the numerical tests presented in this paper, we set tol = 10−5). Up to now, we did not describe the discretization in space. As we said previously, we want to adopt a finite volume approach. Consequently, it is worth realizing that the underlying global problem coupling (16) and (17) involves the following system (we use a slight change of notation which will be useful in the following: Hn+1 is denoted as Hk+1; in spite of this choice, note again that Hk+1 is not involved in the Berm´udez-Moreno algorithm and, so, does not change in this loop): (˜ P)n,k                  Hk+1 −Hn ∆t+∂x(HnVn) = 0, HnVk+1 −Vn ∆t+βV k+1 −∂x4ηHn∂x(Vk+1)−∂x(ω∂xVk+1) =Hn(fΩ+fz∂xb)−Hn 2∂x((Vn)2) + 1 2∂x((Hn)2fz) + ∂xθk. (21) Consequently, even if there is a decoupling of both problems in terms of the time discretization and the Berm´udez-Moreno algorithm, it appears that to obtain a global well-balanced scheme, there must be a coupling between the mass and momentum equations induced by the source terms (involving topography and the multiplier θ). For shallow water type systems with source terms, this has been extensively studied in the literature. In Section 4, we precisely describe the aforementioned coupling for the present problem. 3.1.2 Study of the optimal parameter As it has been documented in the literature, the main drawback of the BM algorithm is that its rate of convergence strongly depends on the choice of parameters λand ω. Several efforts have been made to overcome this problem, allowing the choice of appropriate parameters in different functional frameworks ([30], [31], [19]). We will adapt here the guidelines in [30] in order to deduce an optimal choice of parameters, the main idea being to look for ωminimizing the contractivity constant of the sequence θk−θ. The condition λ ω =1 2 will be assumed throughout the rest of the paper, so only one of the parameters has to be chosen, say ω. The interested reader will find all the details of the derivation in Appendix C. To sum up, a quasi-optimal choice of the parameter ωwould be given by ωopt(Hn max) = Hn max ∆t+βL2 Nπ2+ 4ηHn max,(22) where Hn max =kHnk∞. 9 [Step 3] Update {θk+1 i+1/2}ivia θk+1 i+1/2=                  −ω ξk+1 i+1/2+τy√2Hn i+1/2 1−λ ω if ξk+1 i+1/2> λτy√2Hn i+1/2, ξk+1 i+1/2 λif ξk+1 i+1/2∈[−λτy√2Hn i+1/2, λ τy√2Hn i+1/2], −ω ξk+1 i+1/2−τy√2Hn i+1/2 1−λ ω if ξk+1 i+1/2<−λτy√2Hn i+1/2. (37) [Step 4] Set: for all i,Vk i=Vk+1 i,θk i+1/2=θk+1 i+1/2and return to Step 1. [Step 5] At convergence, when condition (20) is verified, set Vn+1 i=Vk+1 iand θn+1 i+1/2=θk+1 i+1/2∀i. – Resolution of the problem on Hk+1 Compute Hn+1 =Hk+1 with the finite volume method defined by the first component of (28), defined in terms of the most recent Lagrange multiplier {θn+1 i+1/2}iand taking into account the wet/dry treatment presented before where needed. It is worth giving some more details about [Step 1]. Remark that the second component of (28) defines a linear system where the unknowns are {Vk+1 i}i. If we denote Vk+1 the vector whose ith component is Vk+1 i, the aforementioned linear system can be written as AnVk+1 =bn,k,(38) where Anis a matrix defined in terms of {Hn i}i; consequently, Andoes not change during the duality loop (in k). As a matter of fact, Anis a tridiagonal matrix, whose line iis defined by the following entries: An i,i−1=−4η ∆x2Hn i−1/2−δn i−1/2 ∆x2,An i,i+1 =−4η ∆x2Hn i+1/2−δn i+1/2 ∆x2, An i,i =Hn i ∆t+4η ∆x2(Hn i−1/2+Hn i+1/2) + 1 ∆x2(δn i−1/2+δn i+1/2) + β. On the contrary, the right hand side of the linear system (38) changes for each iteration in k. The ith component of bn,k is decomposed as bn,k i=bn,(1) i+bn,k,(2) i+bn,k,(3) i,(39) where bn,(1) i=Hn ifΩ+fz bi+1 −bi−1 2 ∆x,(40) bn,k,(2) i=ζk i+1/2−ζk i−1/2 ∆x,(41) bn,k,(3) i=Hn i [φ(Wn i−1, Wn i,{ζk j+1/2}j=i+1 j=i−1)]2−[φ(Wn i, Wn i+1,{ζk j+1/2}j=i+1 j=i−1)]2 ∆x.(42) Note that all of this is applied only where 1 2(Hn i−1/2+Hn i+1/2)≥Hε. Indeed, on the contrary, we set An i,i−1=An i,i+1 = 0, An i,i = 1, bn,k i= 0. 16 4.2 Well-balanced properties In this section, we study the well-balanced properties of the proposed numerical scheme. Concretely, we are interested in studying stationary solutions defined by a constant free surface or by a free surface parallel to the reference slope with variable bottom, for rigid enough materials (see Figure 2). H(x) b(x) z x α Figure 2: Stationary solution with free surface parallel to the reference plane. The first type of stationary solution corresponds to a material at rest with zero velocity and a constant free surface. It is a stationary solution of the system independently of the rigidity of the material. For the second type, we consider a stationary solution with zero velocity and verifying b(x) + H(x) = C,∀x∈[0, L], where Cis a constant value. Let us study the conditions on such a solution. One should have ZL 0 τy√2H|∂xψ|dx ≥ZL 0H(fΩ+fz∂xb)ψ−fz H2 2∂xψdx, ∀ψ. Taking into account that ∂xb=−∂xHand integrating by parts, this condition is equivalent to ZL 0 τy√2H|∂xψ|dx ≥ZL 0 HfΩψdx =ZL 0−fΩZx 0 H(s)ds +c∂xψ dx. We can set c=fΩRL/2 0H(x)dx. Then, this condition is verified if  fΩZx L/2 H(s)ds≤τy√2H(x)∀x∈[0, L].(43) We can also obtain the analytic value to which (i) µconverges in the case of the augmented Lagrangian method, and (ii) θfor the Berm´udez-Moreno method. For the AL or the BM algorithms, in the case of a stationary solution with zero velocity, we have the following equation: ∂xfΩ 2H2=H(fΩ+fz∂xb) + ∂xζ, with ζ(x) =      H(x)µ(x) in the case of AL, θ(x) in the case of BM. (44) Taking into account that the stationary solution that we consider in this test verifies b+H= C, this equation simplifies to the following one: −HfΩ=∂x(ζ). 17 Then, ζ=−fΩZx L/2 H(x)dx. (45) We have the following result for the proposed algorithms: Theorem 1 The AL and the BM algorithms preserve exactly the following two types of stationary solutions: i) Material at rest with free surface parallel to the reference slope: V= 0, b +H=constant, for any given bottom function b(x), if the material is rigid enough, i.e. if τyis such that the following discrete version of (43) holds:  fΩ∆xi X j=1 H0 j− [M/2] X j=1 H0 j≤τy√2H0 i+1/2∀i. (46) ii) Material at rest with constant free surface: V= 0, fΩx+fz(H+b) = constant, if {ζi+1/2}N i=1 is initialized as follows: ζ1 i+1/2=−∆xi X j=1 H0 j− [M/2] X j=1 H0 jfΩ+fz H0 i+1 +bi+1 −(H0 i+bi) ∆x,(47) where [M/2] is the integer part of M/2, being Mthe number of points of the mesh grid.  The proof is given in Appendix E. 5 Numerical comparison of both approaches In this section, we present numerical tests to illustrate the good properties of the schemes presented above. A first natural test is to check the order of convergence in space on a non trivial stationary problem, in order to test the accuracy of the duality methods for variational inequalities (thus in absence of well-balancing and wet/dry front issues). It appears that the viscoplastic model (16)–(17) degenerates to the well-known Poiseuille-Bingham flow for which an analytic solution is known. All the details of this test are given in Appendix F. The conclusions are that both the augmented Lagrangian and the Berm´udez-Moreno methods are at least of order two (in L2- norm) in space for the velocity of this non-zero stationary solution. In terms of computational cost, for this specific test and a fixed duality parameter (either ror ω), it is shown that the Berm´udez-Moreno is approximately 20% cheaper than the augmented Lagrangian . We then focus on the main novelty of this paper, i.e. the well-balanced properties of the schemes. 18 5.1 Well-balanced tests Test 5.1.a: Analytical stationary solution In this test we study the error and convergence of the numerical results for the stationary solution studied in Subsection 4.2. Concretely, we consider a domain of length L= 10 and a solution defined by V= 0, H(x) = 2 −b(x), b(x) = cos(πx), 0 1 2 3 4 5 6 7 8 9 1 0 −0.5 0 0.5 1 1.5 2 2.5 3 3.5 4 Free surface Bottom !−plane Figure 3: Test 5.1.a. Free surface, bottom and Ω-plane. as shown in Figure 3. This is a stationary solution of the system if the material is rigid enough. For this definition of H(x) we can compute condition (43) exactly; we obtain that it is a stationary solution of the system if τyverifies τy≥max x∈[0,10]  gsin α2(5 −x)−sin(πx) π √2(2 −cos(πx)) . For this test we set α= 10◦, so τymust be approximately greater than 17.03. We compare the numerical results with the analytical ones corresponding to H(x) and the multipliers µand θ. By (44) and (45) we have that µ(x) = 2(x−5) −sin(πx)/π 2−cos(πx)gsin α, θ(x) = 2(x−5) −sin(πx) πgsin α. Following Theorem 1, if we initialize the multipliers with (47) then the stationary solution is exactly preserved, up to machine precision. Then, we initialize both multipliers to zero in order to study their convergence to the analytical values. In Figure 4, we present the convergence of µand θto the analytical solution when both quantities are initialized to zero, for 100 computational cells. In Figure 5, we study the numerical order of convergence in space, through a mesh refinement. For µand θfirst order is reached by computing the error in the L∞norm, and second order for the L2error. For the BM algorithm the parameter ωwas set to ωopt (equation (69)), while for the AL method we have set r= 10. In Figure 6, we present a comparison of the number of iterations kend necessary to converge in the fixed point algorithm of the Augmented Lagrangian and the Berm´udez-Moreno methods, for the first iteration in time. We also picture a vertical line corresponding to the optimal parameter for the BM method. For the BM method we have proposed two different ways to define the parameter ω, as a constant value in space (but variable in time), or depending on the thickness of the material 19 ! " # $ % & ' ( ) * "! !#! !"& !"! !& ! & "! "& #! + + ,-./0123./ 4567823./ (a) µ(x) ! " # $ % & ' ( ) * "! !#! !"& !"! !& ! & "! "& #! + + ,-./0123./ 4567823./ (b) θ(x) Figure 4: Test 5.1.a. Convergence of µand θto the analytical values layer: ω=Hn(x)$, being $a constant value in space. Then, we denote by BM(ω) the results corresponding to the choice of ωas a constant parameter, and by BM($) the results corresponding to ω=Hn(x)$. For BM(ω), the optimal parameter wopt is defined by (69), equals to 20.5 approximately for this test. While for BM($) the optimal parameter $opt is defined by (23), equals to 10.3 approximately for this test. Note that both, ωopt and $opt are near to the optimal value obtained numerically, which are respectively 26 and 21. The results presented in Figure 6 correspond to ∆x= 0.1 and {ω, r, $} ∈ [1,300], concretely a partition with subintervals of length equals to 5. Let us remark that for ω=r=$= 1 the BM($) is the one that need a smaller number of iterations, after the BM(ω) and the one that needs a greater number of iterations is the AL. This behavior is also observed for the values of the parameters smaller than the optimal one. When they are greater than the optimal parameter both versions of the BM algorithm present a similar number of iterations. Nevertheless, for the augmented Lagrangian method we observe that the number of iterations decreases when the value of rincreases. This is a phenomena that has been yet observed in the case of augmented Lagrangian method for other applications. Nevertheless, from a practical point of view, the choice of bigger values of r, when applied to other numerical tests, can imply some problems of stability and ill-conditioning problems in the linear system related to the fixed point algorithm (see [18]). Finally, we study the influence of the angle on the numerical solution. Again we initialize the multipliers to zero and we compare the errors for V(x) and H(x) corresponding to α= 10◦ and α= 45◦through a mesh refinement, see Tables 1 and 2. We observe that errors are almost independent of the angle and the grid. Test 5.1.b: Stationary solution on a random bottom In this test, we consider a random bottom and wet/dry fronts. We also consider two initial conditions to set the position of the free surface: horizontal or parallel to tan(α)x(See Figure 7). Concretely, we set the following two definitions of the height of the material layer: H(x) = max(zref,1−b(x),0),or H(x) = max zref,2 cos(α)−b(x)−tan(α)x, 0. where the bottom function has been defined as b(x) = r1(x)(1+r2(x))er3(x), where rj(x)∈[0,1], j= 1,2,3, are three random numbers for each value of x(See Figure 7). zref,1corresponds to 20 10−310−210−1100 10−6 10−5 10−4 10−3 10−2 10−1 dx Err (a) BM Convergence on θ. L∞error. 10−310−210−1100 10−8 10−7 10−6 10−5 10−4 10−3 10−2 dx Err (b) BM Convergence on θ. L2error. 10−310−210−1100 10−5 10−4 10−3 10−2 10−1 100 101 dx Err (c) AL Convergence on µ. L∞error. 10−310−210−1100 10−6 10−5 10−4 10−3 10−2 10−1 100 dx Err (d) AL Convergence on µ. L2error. Figure 5: Test 5.1.a. Order of convergence of θand µfor the BM and AL algorithms. Black continuous lines show first and second order of convergence; blue line with + is the computed error. 0 50 100 150 200 250 300 0 2000 4000 6000 8000 10000 12000 14000 { !, r, " } k B.M. ( ! ) !opt L.A. B.M. ( " ) "opt Figure 6: Test 5.1.a. Number of iterations in terms of r(for AL), ωand $(for BM). 21 Table 1: Test 5.1.a. Errors for V(x) and H(x) for BM. αCells V L∞Error V L2Error H L∞Error H L2Error 10◦20 2.703E-09 4.271E-10 6.471E-10 8.643E-11 40 1.948E-09 2.196E-10 5.950E-10 5.082E-11 80 9.234E-10 7.323E-11 4.465E-10 2.588E-11 160 1.704E-10 9.552E-12 2.482E-10 1.008E-11 320 3.103E-10 1.222E-11 1.254E-10 3.591E-12 640 5.229E-10 1.456E-11 1.857E-11 3.740E-13 1280 1.974E-10 3.884E-12 2.168E-11 3.106E-13 2560 2.449E-10 3.406E-12 4.742E-12 4.786E-14 45◦20 2.917E-09 4.554E-10 7.021E-10 9.259E-11 40 1.898E-09 2.131E-10 6.851E-10 5.823E-11 80 1.175E-09 9.311E-11 4.553E-10 2.637E-11 160 4.501E-10 2.516E-11 3.043E-10 1.236E-11 320 9.862E-11 3.877E-12 1.567E-10 4.487E-12 640 3.981E-10 1.108E-11 5.827E-11 1.178E-12 1280 3.564E-10 7.011E-12 8.732E-12 1.256E-13 2560 6.904E-11 9.600E-13 1.571E-11 1.589E-13 Table 2: Test 5.1.a. Errors for V(x) and H(x) for AL (r= 10). αCells V L∞Error V L2Error H L∞Error H L2Error 10◦20 3.593E-10 5.038E-11 3.300E-09 4.946E-10 40 4.887E-10 4.150E-11 2.184E-09 2.396E-10 80 4.106E-10 2.401E-11 1.268E-09 9.998E-11 160 2.997E-10 1.218E-11 7.103E-10 3.964E-11 320 1.813E-10 5.172E-12 1.792E-10 7.054E-12 640 8.896E-11 1.789E-12 1.815E-10 5.037E-12 1280 2.219E-11 3.147E-13 3.449E-10 6.764E-12 2560 1.273E-11 1.279E-13 2.105E-10 2.919E-12 45◦20 3.226E-10 4.666E-11 3.328E-09 4.975E-10 40 4.637E-10 3.997E-11 2.148E-09 2.352E-10 80 4.064E-10 2.371E-11 1.315E-09 1.032E-10 160 3.047E-10 1.238E-11 7.871E-10 4.389E-11 320 1.946E-10 5.549E-12 2.903E-10 1.142E-11 640 1.038E-10 2.087E-12 8.149E-11 2.259E-12 1280 3.619E-11 5.137E-13 2.924E-10 5.733E-12 2560 6.308E-12 6.349E-14 2.538E-10 3.518E-12 the height of the material on the Ω-plane, while zref,2corresponds to the level of the horizontal free surface. For example, Figure 7 is obtained with zref,1= 2, zref,2= 3. With the purpose to consider the case with and without wet/dry fronts, we set two different values of zref,1and zref,2, concretely: zref,1= 2, zref,1= 5 and zref,2= 3, zref,2= 5.8. 22 The numerical results presented in this test correspond to α= 10◦. Analogously to the previous tests, we obtain similar results if we increase the angle. For this test we set ∆x= 0.1. The multipliers µand θare initialized to zero. In Figure 8 we present the convergence of µand θto (47). This theoretical value has been used in Theorem 1 to prove that, with this initialization of the multipliers, the proposed numerical scheme preserves exactly both types of stationary solutions. In Tables 3 and 4 are presented the errors for V(x) and H(x), when µand θare initialized to zero. ! " # $ % & ' ( ) * "! ! " # $ % & ' + + ,-..+/0-123. 456657 !!892:. ! " # $ % & ' ( ) * "! ! " # $ % & ' + + ,-..+/0-123. 456657 !!892:. Figure 7: Test 5.1.a. Free surface, bottom and Ω-plane. ! " # $ % & ' ( ) * "! !"! !& ! & "! + + ,-./0.12345 678.02345 (a) µ(x) ! " # $ % & ' ( ) * "! !%! !$! !#! !"! ! "! #! $! %! + + ,-./0.12345 678.02345 (b) θ(x) Figure 8: Test 5.1.b. Initial condition: H(x) = zref,1−b(x) with zref,1= 5. Convergence of µ and θto the theoretical value (47). Table 3: Test 5.1.a. Initial condition: H(x) = max(zref,1−b(x),0). Errors for V(x) and H(x). zref,1{ζ1 i+1/2}iVL∞Error V L2Error H L∞Error H L2Error 2 BM(ωopt) 1.164E-11 1.536E-13 1.257E-08 1.772E-10 BM ($opt) 3.672E-09 5.301E-11 1.210E-08 2.780E-10 AL (r= 10) 4.978E-10 1.117E-11 8.800E-10 1.103E-11 5 BM(ωopt) 3.109E-10 2.191E-11 8.731E-10 2.780E-11 BM($opt) 3.414E-10 2.415E-11 8.733E-10 2.811E-11 AL(r= 10) 8.847E-10 6.203E-11 9.557E-10 2.962E-11 23 Table 4: Test 5.1.a. Initial condition: H(x) = max zref,2 cos(α)−b(x)−tan(α)x, 0. Errors for V(x) and H(x). zref,2Method V L∞Error V L2Error H L∞Error H L2Error 3 BM(ωopt) 6.316E-17 1.078E-18 1.662E-11 3.325E-13 BM ($opt) 6.848E-17 1.695E-18 3.317E-10 6.634E-12 AL (r= 10) 3.933E-17 1.430E-18 2.220E-16 4.680E-18 5.8 BM(ωopt) 1.951E-16 1.094E-17 1.776E-15 4.721E-17 BM($opt) 1.700E-16 1.004E-17 1.332E-15 4.569E-17 AL(r= 10) 7.430E-17 3.760E-18 8.882E-16 3.320E-17 5.2 Avalanche with obstacle For this test we consider the case of an avalanche on an inclined plane with an obstacle. Concretely, we set α= 30◦and b(x) = 14e−x2/1.2+ 2e−(x−5)4/0.1+ 4e−(x−10)2/0.8. As initial condition, we set V= 0 and (see Figure 9), H(x) = (4−b(x) if x∈[7,9], 0 otherwise. We study the influence of the rigidity coefficient on the evolution of the avalanche and the final solution at rest. The length of the domain is L= 10. In Figures 10, 11, 12 the evolution of the avalanche for τy∈ {1,4,8,12}is presented at times t= 1, t= 1.5, and t= 2. In Figure 13, we plot the stationary solution reached for each value of τy. In these figures, we only present the results obtained with the Berm´udez-Moreno method. Of note, the results are exactly the same with augmented Lagrangian method, so we do not present them for sake of brevity. We have considered 200 computational cells. On the left column of these figures, we present the evolution of the free surface. Right column corresponds to the velocity. We can remark zones with evidence of rigidity: we can distinguish clearly some zones with constant velocity, that is, zones where the material moves as a block. We can remark two difficulties of this test related with the wet/dry front. First, in the evolution of the avalanche, the obstacle in the middle of the domain splits the avalanche in two parts. Second, the part of the avalanche arriving at the far left of the domain goes up on a high bed which limits its movement. This leads to a back and forth motion that eventually ends to a stationary state when all the material becomes rigid. This back and forth motion goes faster to stationary state when τyincreases. But the associated free surface has a more complex shape, which is also due to the complex, non linear, interaction of the material when it passes over the obstacle inducing the splitting of the material in the two basins. For τy= 1 we can observe in Figure 13(a) that the stationary solution is close to the one of a fluid, that is, an horizontal free surface. The bump in the middle of the domain produces that the solution is divided in two parts and two different levels of the free surface. In Figures 13(c)-13(g) we can observe the influence of the rigidity of the material on the final stationary solution. Let us also remark that the computed velocity at the stationary solution are in all cases of order 10−9. In Figure 14 a comparison of the free surface at rest with a mesh refinement is presented. The results correspond to 200, 400, 800 and 1600 computational cells. For all cases of τywe 24 − 6 −4−2 0 2 4 6 0 2 4 6 8 10 12 14 Free surface Bottom !−plane Figure 9: Test 5.2. Free surface: initial condition. can remark that the position and form of both parts in which is divided the avalanche by the obstacle in the middle of the domain agree with the mesh refinement. When τyincreases we can observe that it is more difficult to capture exactly the shape of the free surface, although the averaged form is well captured in any case. In Figure 15 we present the computational cost until t= 1 (sum of the number of iterations in the duality loops at each time step) with respect to the duality parameter. We consider the four cases, τy∈ {1,4,8,12}, and {ω, r, $} ∈ [0.1,10]. Let us remark that for each value of {ω, r, $} ∈ [0.1,10], we set it invariant for all time steps. Nevertheless, the computation of the theoretical optimal value of ωand $,ωopt and $opt, defined by (69) and (23), respectively, are variable in time. In order to compare the computational cost of BM algorithms with this optimal choice of the parameters, in Figure 15 we mark two horizontal lines at the level of the sum of the number of iterations in the duality loops at each time step obtained with BM(ωopt) and BM($opt). Let us remark that BM(ωopt) is slightly better than BM($opt) for the case τy= 1, although both are close. Nevertheless, for τy∈ {4,8,12}there is a great difference of efficiency, being BM($opt) close to the optimal computational cost in all situations. Finally, we remark that this difference is produced by the effect of the wet/dry fronts. In the case of numerical tests without wet/dry fronts, the behavior of both versions of BM method are very similar. 6 Conclusions In this work, we proposed a discretization of a shallow Bingham model by a well-balanced finite volume method which is combined with duality techniques. Augmented Lagrangian and Berm´udez-Moreno algorithms have been considered to discretize the momentum equations. For the mass conservation equation we proposed a well-balanced correction which depends on the definition of the multiplier associated to the duality technique. This correction includes the use of a limiter that has been specially designed to recover the well-balanced properties of the numerical method. We prove that the proposed methods are able to preserve exactly two types of stationary solutions. A treatment of wet/dry fronts has also been proposed. It takes into account the rigidity of the material. For the case of the B.M. algorithm the definition of the optimal value of the parameter ωhas been deduced. We also present two different versions of B.M. algorithm, by considering that ωcan be variable in space. In the numerical tests section, we have first compared the algorithms for an analytical solution for a simplified model. Second, we have compared with the analytical solution of the multipliers for the case of a stationary solution. For the case of a random bottom a comparison of the multiplier with the theoretical one, which is proposed in Theorem 1, is also presented. Finally, we consider a test corresponding 25 [10] M. Castro, A. Ferreiro Ferreiro, J. Garc´ıa-Rodr´ıguez, J. Gonz´alez-Vida, J. Mac´ıas, C. Par´es, and M. Elena V´azquez-Cend´on. The numerical treatment of wet/dry fronts in shallow flows: Application to one-layer and two-layer systems. Mathematical and Computer Modelling, 42(3-4):419–439, 2005. [11] M. J. Castro, J. M. Gonz´ales-Vida, and C. Par´es. Numerical treatment of wet/dry fronts in shallow flows with a modified roe scheme. Math. Mod. Meth. App. Sci., 16(6):897–931, 2006. [12] T. Chac´on, M. J. Castro, E. D. Fern´andez-Nieto, and C. Par´es. On well-balanced finite volume methods for non-conservative non-homogeneous hyperbolic systems. SIAM J. Sci. Comput., 29(3):1093–1126, 2007. [13] F. Delbos, J. C. Gilbert, R. Glowinski, and D. Sinoquet. Constrained optimization in seismic reflection tomography: a Gauss-Newton augmented Lagrangian approach. Geophysical Journal International, 164(3):670–684, 2006. [14] J. Dent and T. Lang. A biviscous modified Bingham model of snow avalanche motion. Ann. Glaciol., 4:42–46, 1983. [15] G. Duvaut and J.-L. Lions. Inequalities in mechanics and physics. Springer-Verlag, 1976. [16] I. Ekeland and R. T´emam. Convex analysis and variational problems, volume 28 of Classics in Applied Mathematics. Society for Industrial and Applied Mathematics (SIAM), 1999. [17] E. D. Fern´andez-Nieto, P. Noble, and J.-P. Vila. Shallow water equations for non-newtonian fluids. Journal of Non-Newtonian Fluid Mechanics, 165(13–14):712 – 732, 2010. [18] M. Fortin and R. Glowinski. Augmented Lagrangian methods: applications to the numerical solution of boundary-value problems. North-Holland, 1983. [19] J. M. Gallardo, C. Par´es, and M. Castro. A generalized duality method for solving variational inequalities. Applications to some nonlinear Dirichlet problems. Numer. Math., 100(2):259–291, 2005. [20] J.-F. Gerbeau and B. Perthame. Derivation of viscous Saint-Venant system for laminar shallow water; numerical validation. Discrete Contin. Dyn. Syst., Ser. B, 1(1):89–102, 2001. [21] R. Glowinski and P. Le Tallec. Augmented Lagrangian and operator-splitting methods in nonlinear mechanics, volume 9 of SIAM Studies in Applied Mathematics. Society for Industrial and Applied Mathematics (SIAM), 1989. [22] R. Glowinski and A. Wachs. On the numerical simulation of viscoplastic fluid flow. In R. Glowinski and J. Xu, editors, Numerical Methods for Non-Newtonian Fluids, volume 16 of Handbook of Numerical Analysis, pages 483 – 717. Elsevier, 2011. [23] J. M. Greenberg and A.-Y. Le Roux. A well-balanced scheme for the numerical processing of source terms in hyperbolic equations. SIAM J. Numer. Anal., 33(1):1–16, 1996. [24] R. R. Huilgol and Z. You. Application of the augmented Lagrangian method to steady pipe flows of Bingham, Casson and Herschel-Bulkley fluids. Journal of Non-Newtonian Fluid Mechanics, 128(2-3):126 – 143, 2005. [25] I. R. Ionescu. Augmented lagrangian for shallow viscoplastic flow with topography. Journal of Computational Physics, 242:544–560, 2013. 32 [26] K. Ito and K. Kunisch. Lagrange Multiplier Approach to Variational Problems and Applications. Society for Industrial and Applied Mathematics, 2008. [27] R. J. LeVeque. Balancing source terms and flux gradients in high-resolution Godunov methods: the quasi-steady wave-propagation algorithm. J. Comput. Phys., 146(1):346–365, 1998. [28] K. D. Nikitin, M. A. Olshanskii, K. M. Terekhov, and Y. V. Vassilevski. A numerical method for the simulation of free surface flows of viscoplastic fluid in 3D. Journal of Computational Mathematics, 29(6):605–622, 2011. [29] T. Papanastasiou. Flows of materials with yield. J. Rheol., pages 385––404., 1987. [30] C. Par´es, M. Castro, and J. Mac´ıas. On the convergence of the Berm´udez-Moreno algorithm with constant parameters. Numer. Math., 92(1):113–128, 2002. [31] C. Par´es, J. Mac´ıas, and M. Castro. Duality methods with an automatic choice of parameters. Application to shallow water equations in conservative form. Numer. Math., 89(1):161–189, 2001. [32] P. L. Roe. Upwind differencing schemes for hyperbolic conservation laws with source terms. In Nonlinear hyperbolic problems (St. Etienne, 1986). C. Carraso et al. (Eds), volume 1270 of Lecture Notes in Math., pages 41–51. Springer, Berlin, 1987. [33] N. Roquet and P. Saramito. An adaptive finite element method for bingham fluid flows around a cylinder. Computer Methods in Applied Mechanics and Engineering, 192(31- 32):3317–3341, 2003. [34] E. F. Toro. Shock-capturing methods for free-surface shallow flows. John Wiley & Sons, 2001. [35] G. Vinay, A. Wachs, and J.-F. Agassant. Numerical simulation of weakly compressible bingham flows: The restart of pipeline flows of waxy crude oils. Journal of Non-Newtonian Fluid Mechanics, 136(2-3):93–105, 2006. [36] D. Vola, F. Babik, and J.-C. Latch´e. On a numerical strategy to compute gravity currents of non-Newtonian fluids. J. Comput. Phys., 201(2):397–420, 2004. [37] D. Vola, L. Boscardin, and J. C. Latch´e. Laminar unsteady flows of Bingham fluids: a numerical strategy and some benchmark results. J. Comput. Phys., 187(2):441–456, 2003. 33 Appendix A Derivation of the AL approach Let us describe an augmented Lagrangian method designed to solve problem (17) for Vn+1, supposing that (Hn, V n) are known. We refer to Glowinski and coworkers for details on this duality method `a la Uzawa: the book [21] and the article [22] for a recent review in the context of Bingham flows. The variational inequality (17) can be rewritten as a minimization problem: Jn(Vn+1) = min V∈V Jn(V),(48) where Jn(V) = Fn(B(V)) + Gn(V), with V=H1 0([0, L]), H=L2([0, L]), B:V → H, B(V) = ∂xV, Fn:H → R, F(λ) = ZL 0 τy√2Hn|λ|dx, and Gn:V → R, Gn(V) = ZL 0 HnV2/2−VnV ∆t+1 2∂x((Vn)2)Vdx +ZL 0 βV2 2dx +ZL 0 4ηHn1 2(∂xV)2dx −ZL 0 (fΩ+fz∂xb)HnV+ZL 0 fz (Hn)2 2∂xV dx. As Jn(V) is a non-differentiable function, we consider the Lagrangian Ln:V ×H×H → R, Ln(V, q, µ) = Fn(q) + Gn(V) + ZL 0 Hnµ(B(V)−q)dx, and the augmented Lagrangian function, for a given positive value r∈R: Ln r(V, q, µ) = Ln(V, q, µ) + r 2ZL 0 Hn(B(V)−q)2dx. (49) Then, we search for the saddle point of Ln r(V, q, µ) over V × H × H. Indeed, if we denote by (V∗, q∗, µ∗) this saddle point, then V∗is the solution of the minimization problem (48) (cf. [21]). To do so, we consider the algorithm proposed in [21], based on Uzawa’s algorithm, to approximate the saddle point of (49). Augmented Lagrangian algorithm •Initialization: Suppose that Vn,Hnand µnare known. For k= 0, we set Vk=Vnand µk=µn. Initialize r. •Iterate: –Find qk+1 ∈ H solution of Ln r(Vk, qk+1, µk)≤ Ln r(Vk, q, µk),∀q∈ H. 34 In other words, qk+1 ∈ H is the solution of following minimization problem: min q∈H Hnr 2q2+Hnτy√2|q|−Hn(µk+rB(Vk))q.(50) The solution of this problem is (denoting the sign function as “sgn”): qk+1 =   0 if |µk+rB(Vk)|< τy, 1 r(µk+rB(Vk)) −τy√2 sgn(µk+rB(Vk))otherwise. (51) –Find Vk+1 ∈ V solution of Ln r(Vk+1, qk+1, µk)≤ Ln r(V, qk+1, µk),∀V∈ V. Thus, Vk+1 is the solution of a minimization problem, which can be characterized by differentiating Ln r(V, q, µ) with respect to V. From (49), we deduce that Vk+1 is the solution of the following linear problem (whose resolution is detailed in Section 4): HnVk+1 −Vn ∆t+βV k+1 −∂x4ηHn∂x(Vk+1)−∂xrHn∂x(Vk+1) = (fΩ+fz∂xb)Hn+∂xfz (Hn)2 2−Hn 2∂x((Vn)2) + ∂x(Hn(µk−rqk+1)).(52) –Update the Lagrange multiplier via µk+1 =µk+r(∂xVk+1 −qk+1).(53) –Check convergence (see below) and update: Vk=Vk+1,µk=µk+1,k7→ k+ 1 and go to the next iteration... •... until convergence is reached: kµk+1 −µkk kµkk≤tol.(54) At convergence, we get the value of Vn+1 by setting Vn+1 =Vk+1 (in the numerical tests presented in this paper, we set tol = 10−5). It is shown in [21] that this algorithm converges to the saddle point of (49). Of note, we did not describe the discretization in space yet. As we said previously, we want to adopt a finite volume approach. Consequently, it is worth realizing that the underlying global problem coupling (16) and (17) involves the following system (we use a slight change of notation which will be useful in the following: Hn+1 is denoted as Hk+1; in spite of this choice, note again that Hk+1 is not involved in the augmented Lagrangian algorithm and, so, does not change in this loop): (P)n,k                  Hk+1 −Hn ∆t+∂x(HnVn)=0, HnVk+1 −Vn ∆t+βV k+1 −∂x4ηHn∂x(Vk+1)−∂xrHn∂x(Vk+1) = (fΩ+fz∂xb)Hn+∂x(Hn)2fz 2−Hn 2∂x((Vn)2) + ∂x(Hn(µk−rqk+1)). (55) 35 Consequently, even if there is a decoupling of both problems in terms of the time discretization and the augmented Lagrangian algorithm, it appears that to obtain a global well-balanced scheme, there must be a coupling between the mass and momentum equations induced by the source terms (involving topography and the Lagrange multiplier). For shallow water type systems with source terms, this has been extensively studied in the literature. In Section 4, we precisely describe the aforementioned coupling for the present problem. The natural follow-up of this appendix is Section 3.1. B Derivation of the BM method In this section, the solution of the velocity problem (17) is approached by means of the duality algorithm introduced in [4]. We shall focus on the application of the method to our particular case and refer to [4], [30], [31] and [19] for further details. Define V=H1 0([0, L]) and let h·,·i be the duality pairing between Vand its dual space V0=H−1([0, L]). The variational inequality (17) can be rewritten as: Find V∈ V such that hA(V),Ψ−Vi+j(Ψ) −j(V)≥ hL, Ψ−Vi(56) holds for every Ψ ∈ V. Here, A:V → V0denotes the linear operator hA(V),Ψi=ZL 0Hn ∆t+βV−∂x(4ηHn∂xV)Ψdx, Ψ∈ V. Throughout this section it will be assumed that there exists a constant Hn min such that Hn(x)≥ Hn min >0. Then the operator Ais coercive with constant γ= min Hn min ∆t+β, 4ηHn min>0.(57) The functional j:V → Ris defined by j(V) = ZL 0 Φ(x, B(V)(x))dx, where Φ: [0, L]×R→Ris given by Φ(x, z) = τy√2Hn(x)|z|, and B:V → H is the derivative operator B(V) = ∂xV, where H=L2([0, L]) with the usual scalar product (·,·)L2and norm k·kL2. Finally, L∈ V0represents the functional L=HnVn ∆t−Hn 2∂x(Vn)2+Hn(fΩ+fz∂xb) + 1 2∂x((Hn)2fz). Notice that j(V) = T(B(V)), where T:H → Ris given by T(Z) = ZL 0 Φ(x, Z(x))dx. Let ω > 0 be an arbitrary parameter and define Gω=∂T −ωI, where the subdifferential ∂T is the multivalued operator ∂T(Z) = {W∈ H:T(Q)−T(Z)≥(W, Q −Z),∀Q∈ H}, Z ∈ H. 36 As Tis proper, convex and lower semicontinuous, its subdifferential turns out to be a maximal monotone operator. Thus, if λ > 0 is such that λω < 1, the resolvent Jλ= (I+λGω)−1is an univalued operator defined on H. Moreover, if λω ≤1/2 the Yosida regularization of Gω, Gω λ=I−Jω λ λ, is a Lipschitz function with constant 1/λ (see [9]). Observe that, due to the continuity of T, the subdifferentials of jand Tare related as ([16]) ∂j(V) = B∗(∂T(B(V))) , where B∗denotes the dual operator of B. Thus, problem (56) can be reformulated as follows: Find V∈ V and θ∈ H such that (A(V) + B∗(ωB(V)) + B∗(θ) = L, θ=Gω λ(B(V) + λθ).(58) As Gω λis a Lipschitz function, it makes then sense to define the fixed-point Berm´udez-Moreno algorithm as follows: For k≥0, θkbeing known, compute Vkand θk+1 by solving (A(Vk) + B∗(ωB(Vk)) + B∗(θk) = L, θk+1 =Gω λ(B(Vk) + λθk).(59) As it was proved in [4], the sequence Vkconverges to the solution Vdue to the coerciveness of the operator A. For questions regarding the convergence of the multipliers θkwe refer the reader to [30]. The key point in the definition of the BM algorithm (59) is the construction of the Yosida regularization Gω λ, that will be worked out in what follows. First of all, notice that ([16]) ∂T(Z) = {W∈ H:W(x)∈∂Φ(x, Z(x)) a.e. x∈[0, L]}.(60) Thus, a simple computation shows that ∂T(Z) is the set of elements W∈ H that verify W(x)∈       {τy√2Hn(x)}if Z(x)>0, [−τy√2Hn(x), τy√2Hn(x)] if Z(x)=0, {−τy√2Hn(x)}if Z(x)<0, a.e. x∈[0, L]. After some algebra, we deduce the following expression, which is valid a.e. x∈[0, L] for each Z∈ H: Gω λ(Z)(x) =                −ωZ(x) + τy√2Hn(x) 1−λω if Z(x)> λτy√2Hn(x), Z(x) λif Z(x)∈[−λτy√2Hn(x), λτy√2Hn(x)], −ωZ(x)−τy√2Hn(x) 1−λω if Z(x)<−λτy√2Hn(x). (61) The complete algorithm for solving (17) is then given in Section 3.1.1 of the main text. 37 C Study of the BM optimal parameter Given an uniform partition of [0, L] of size h, let Vhbe the finite-dimensional subspace of Vof standard conforming P1finite elements. Consider now the discrete versions of (58) and (59), where the elements Vand Vkare assumed to belong to Vhinstead of V(for the sake of clarity, the dependence on hwill not be explicitly stated). Combining both expressions, we get hA(Vk−V),Ψi+ω(∂xVk−∂xV, ∂xΨ)L2+ (θk−θ, ∂xΨ)L2= 0,∀Ψ∈ Vh,(62) and also, using that Gω λis Lipschitz with constant 1/λ, kθk+1 −θk2 L2=kGω λ(∂xVk+λθk)−Gω λ(∂xV+λθ)k2 L2 ≤1 λ2k∂xVk−∂xVk2 L2+kθk−θk2 L2+2 λ(θk−θ, ∂xVk−∂xV)L2.(63) Taking Ψ = Vk−Vin (62), we deduce hA(Vk−V), V k−Vi+ωk∂xVk−∂xVk2 L2+ (θk−θ, ∂xVk−∂xV)L2= 0, so (63) can be written as kθk+1 −θk2 L2≤ kθk−θk2 L2−4ωhA(Vk−V), V k−Vi,(64) taking into account that λω = 1/2. Assume that γ1and γ2are positive constants (which may depend on h) verifying γ1kΨkL2≤ k∂xΨkL2≤γ2kΨkL2,∀Ψ∈ Vh,(65) and define Hn max =kHnk∞. Then, from (62) and the definition of A, we have (θk−θ, ∂xΨ)L2=−hA(Vk−V),Ψi−ω(∂xVk−∂xV, ∂xΨ)L2 ≤Hn max ∆t+βkVk−VkL2kΨkL2+ (4ηHn max +ω)k∂xVk−∂xVkL2k∂xΨkL2 ≤Hn max ∆t+βγ−1 1+ (4ηHn max +ω)γ2kVk−VkL2k∂xΨkL2,∀Ψ∈ Vh. Taking now Ψ ∈ Vhsuch that ∂xΨ = θk−θ, we obtain from the above inequality kθk−θkL2≤Hn max ∆t+βγ−1 1+ (4ηHn max +ω)γ2kVk−VkL2.(66) On the other hand, using the coerciveness of A, we have that hA(Vk−V), V k−Vi ≥ γkVk−Vk2 V=γ(kVk−Vk2 L2+k∂xVk−∂xVk2 L2) ≥γ(1 + γ2 1)kVk−Vk2 L2,(67) where the coercivity constant γwas given by (57). Finally, combining (64), (66) and (67), we deduce the following inequality: kθk+1 −θkL2≤ L(ω)kθk−θkL2,(68) where L(ω) = 1−4ωγ(1 + γ2 1)Hn max ∆t+βγ−1 1+ (4ηHn max +ω)γ2−21/2 . 38 The optimal choice of the parameter ωopt will be that minimizing L(ω). An easy computation shows that ωopt =Hn max ∆t+β1 γ1γ2 + 4ηHn max. From a practical point of view, it is necessary to have good estimates of the constants γ1and γ2appearing in (65). To this end, we consider the following spectral problem: Find 0 6=vh∈ Vh and µh∈Rsuch that (v0 h, ϕ0 h)L2=µh(vh, ϕh)L2,∀ϕh∈ Vh. It is well-known ([6]) that, for uniform mesh size h, there exists an orthonormal basis of Vh composed by eigenvectors {ϕ(1) h, . . . , ϕ(N) h}associated to the eigenvalues 0 < µ(1) h≤ ··· ≤ µ(N) hof the spectral problem. Indeed, these eigenvalues have the following form: µ(j) h=6 h2 1−cos(jπh/L) 2 + cos(jπh/L), j = 1, . . . , N, where Ndenotes the number of internal nodes. For arbitrary vh=PN j=1 vjϕ(j) hwe have k∂xvhk2 L2= (∂xvh, ∂xvh)L2= N X j=1 v2 j(∂xϕ(j) h, ∂xϕ(j) h)L2= N X j=1 v2 jµ(j) h(ϕ(j) h, ϕ(j) h)L2(≥µ(1) hkvhk2 L2, ≤µ(N) hkvhk2 L2, so we deduce the inequalities qµ(1) hkvhkL2≤ k∂xvhkL2≤qµ(N) hkvhkL2,∀vh∈ Vh. As noticed in [6], the following optimal estimate holds as h→0: |µ(j)−µ(j) h|=O(h2), where µ(j)= (jπ/L)2. Therefore, we consider the approximations γ1=pµ(1) =π/L and γN=pµ(N)=Nπ/L. Summarizing, a quasi-optimal choice of the parameter ωwould be given by ωopt(Hn max) = Hn max ∆t+βL2 Nπ2+ 4ηHn max.(69) In the case of wet/dry fronts further adaptations can be made and are described in Section 3.1.3. D The coupled scheme in the case of the AL For this duality method the structure of the complete algorithm is very close to the BM one. In particular, the linear problem verified by the speed in the duality loop ([Step 2] below) is embedded in the general formulation (38)-(42). Global numerical scheme for (13)-(14) – Augmented Lagrangian method •Initialization at time t= 0 for n= 0: Vn,Hn,µnare given by the initial conditions. •Time loop: For n= 0, ..., nmax. – Resolution of the problem on Vk+1 {Vn i}i,{Hn i}iand {µn i+1/2}iare known. Compute quantities which are invariant in the following loop: Augmented Lagrangian loop: 39 [Step 0] Initialize for k= 0: for all i,Vk i=Vn iand µk i+1/2=µn i+1/2. [Step 1] Update {qk+1 i+1/2}iby computing qk+1 i+1/2=1 rµk i+1/2+rVk i+1 −Vk i ∆x 1−τy√2 |µk i+1/2+rVk i+1−Vk i ∆x| + .(70) In the case of a wet/dry front, if Hi≥Hε,Hi+1 < Hεand bi+Hi< bi+1 then we set qk+1 i+1/2= 0. Or if Hi< Hε,Hi+1 ≥Hεand bi+1 +Hi+1 < bi, we also set qk+1 i+1/2= 0. [Step 2] Update {Vk+1 i}iby solving the linear system defined by the second component of (28). See details at (38)-(42), with δn i+1/2and ζk i+1/2given by (27)-(a) and (26)-(a). [Step 3] Update {µk+1 i+1/2}ivia µk+1 i+1/2=µk i+1/2+r Vk+1 i+1 −Vk+1 i ∆x−qk+1 i+1/2!.(71) [Step 4] Set: for all i,Vk i=Vk+1 i,µk i+1/2=µk+1 i+1/2and return to Step 1. [Step 5] At convergence, when condition (54) is verified, set Vn+1 i=Vk+1 iand µn+1 i+1/2=µk+1 i+1/2∀i. – Resolution of the problem on Hk+1 Compute Hn+1 =Hk+1 with the finite volume method determined by the first component of (28), defined in terms of the most recent Lagrange multiplier {µn+1 i+1/2}iand taking into account the wet/dry treatment presented before where needed. This can be compared to the Berm´udez-Moreno algorithm presented in Section 4.1. E Proof of Theorem 1 Proof It is enough to prove that if for a time t=tnwe have Hn i=H0 iand Vn i= 0, then Hn+1 i=Hn i and Vn+1 i= 0, ∀i= 1, . . . , M. i) First, let us prove that the numerical scheme preserves exactly the stationary solution defined by Vn i= 0, bi+Hn i=constant, where {Hn i}i={H0 i}iverifies (46). i.a) First, let us prove that Vn+1 i= 0, where Vn+1 i=Vkend i, being kend the number of iterations necessary to converge in the fixed point algorithm. If we prove that bn,1= 0 (see equation (39)), as {V1 j}jis the solution of the linear system (38), it implies that V1 i= 0. Moreover, if we prove that ζ2 i+1/2=ζ1 i+1/2then we obtain that kend = 1 and Vn+1 i= 0. Thus, let us prove that bn,1= 0. The components of the vector bn,1= 0 are bn,1 i=bn,(1) i+bn,1,(2) i+bn,1,(3) i, i = 1, . . . , M. where bn,(1) i,bn,1,(2) iand bn,1,(3) iare defined by (40), (41) and (42), respectively. Taking into account that ζ0 i+1/2is defined by (47) and bi+H0 i=constant, we have that bn,1,(2) i=ζk i+1/2−ζk i−1/2 ∆x=−H0 ifΩ. 40 As V0 i= 0 and bn,1,(3) iis defined in terms of the second component of numerical flux function (30), we have bn,1,(3) i=H0 ifz H0 i+1 −H0 i−1 2 ∆x. Then, taking into account the definition of bn,(1) i(40) we obtain that bn,1 i=H0 ifz 2∆x(bi+1 +H0 i+1 −(bi−1+H0 i−1)) = 0. Thus, bn,1= 0 and consequently V1 i= 0, i= 1, . . . , M. Let us now prove that ζ2 i+1/2=ζ1 i+1/2. ζ1 i+1/2is defined by (47), and by (46) we have that |ζ1 i+1/2| ≤ τy√2H0 i+1/2. Then (see Subsection 4.1) for the AL algorithm we have that |µ1 j+1/2| ≤ τy√2. So, by (70) we have that q2 i+1/2= 0 and consequently, since V1 j= 0 ∀j, by (71) we obtain that µ2 i+1/2=µ1 i+1/2. For the BM algorithm we have that θ1 i+1/2≤τy√2H0 j+1/2; then, by (37) we have that θ2 j+1/2=θ1 j+1/2. Therefore, we deduce that ζ2 i+1/2=ζ1 i+1/2. Consequently, we obtain that Vn+1 i=Vn i= 0. i.b) Now, let us prove that Hn+1 i=Hn i. From (28), and taking into account that Vn j= 0 for j=i−1, i, i + 1, we have that Hn+1 i=Hn i+1 2 ∆t ∆xα0,i+1/2(Hn i+1 −Hn i+ [Gn({ζk j+1/2}j=i+1 j=i−1]1)− α0,i−1/2(Hn i−Hn i−1+ [Gn({ζk j−1/2}j=i+1 j=i−1]1). Then, it is enough to prove that Hn i+1 −Hn i+ [Gn({ζk j+1/2}j=i+1 j=i−1]1= 0,∀i, where [Gn({ζk j+1/2}j=i+1 j=i−1]1=1 fz fΩ∆x+fz(bi+1 −bi) + ∆(ζ+δn∂xV)k i+1/2 Hi+1/2!. In this case, as bi+Hn iis constant, the flux limiter used in the definition of ∆(ζ+δn∂xV)k i+1/2 is equal to one. So, we obtain ∆(ζ+δn∂xV)k i+1/2=1 2(ζkend i+3/2−ζkend i−1/2). By using that ζkend j+1/2=ζ1 j+1/2(see Section i.a) of the proof), for j=i−1 and j=i+ 1, and (47), we obtain that ∆(ζ+δn∂xV)k i+1/2=−∆xHn i+1/2fΩ. As a consequence, [Gn({ζk j+1/2}j=i+1 j=i−1]1=bi+1 −biand Hn i+1 −Hn i+ [Gn({ζk j+1/2}j=i+1 j=i−1]1=Hi+1 +bi+1 −(Hi+bi) = 0. Then, Hn+1 i=Hn i, what concludes the proof for case i). ii) For the case of a the stationary solution with a constant free surface we suppose that fΩxi+fz(H0 i+bi) = constant. (72) 41