C Chalons

A Del Grosso

A second-order well-balanced Lagrange-projection scheme for Shallow Water Exner equations in 1D and 2D

The present work is devoted to the numerical approximation of the Shallow water Exner system in both one and two dimensions, where the Exner equation expresses the evolution in time of the bed sediment. In particular the Grass formula is taken into account to model the solid transport discharge contributions. The numerical scheme is based on the Lagrange-projection formalism which consists in splitting the mathematical model into the acoustic and transport system. In particular, the Exner equation is taken into account only at the transport level; both a decoupled and weakly coupled formulation are proposed. The method is designed in such a way to satisfy the well-balanced property as well. Details to reach the second-order of accuracy are given; numerical results are shown to validate the numerical schemes.

Introduction and governing equations

This work deals with the design and implementation of a second-order well-balanced Lagrange-projection scheme applied to the 1D and 2D shallow water system with bed sediment not-constant in time. Lagrange-projection approach consists in splitting the acoustic and transport waves of the model, leading to the possibility of approximating the two resulting systems with di erent approaches. This reveals itself to be useful for instance in subsonic regimes, where the acoustic waves are the reason of the restrictive CFL condition one has to employ in order to have a stable numerical scheme. Indeed, the Lagrange-projection decomposition makes possible to implicitly approximate only the acoustic system and, thus, to circumnavigate the problem of restrictive time-steps.

Nowadays, the Lagrange-projection approach have been studied in order to satisfy di erent properties and as applied to several models. Giving few examples, we refer for instance to the work [START_REF] Chalons | An All-Regime Lagrange-Projection Like Scheme for the Gas Dynamics Equations on Unstructured Meshes[END_REF], where all-regime rst-order explicit and semi-implicit Lagrange-projection schemes have been applied to the gas dynamics model in several dimensions, or to [START_REF] Chalons | An all-regime Lagrange-Projection like scheme for 2D homogeneous models for two-phase ows on unstructured meshes[END_REF], where the scheme was extended to the 2D two-phase ows model. Another possible reference is [START_REF] Bouchut | An entropy satisfying two-speed relaxation system for the barotropic Euler equations[END_REF] with the numerical approximation of low Mach number ows of the barotropic Euler equations where the asymptotic-preserving property is satis ed as well. On the other hand, when it comes to the modeling of the shallow water system in the Lagrange-projection formalism, we can refer to [START_REF] Chalons | A large time-step and well-balanced Lagrange-Projection type scheme for the Shallow-water equations[END_REF] for an implicit well-balanced rst-order scheme, to [START_REF] Manuel | A fully well-balanced Lagrange-Projection type scheme for the Shallow-water equations[END_REF] for a fully well-balanced rstorder explicit method, and nally to [START_REF] Morales | High order fully well-balanced Lagrange-Projection scheme for Shallow-water[END_REF] for high-order fully well-balanced schemes. Last but not least, and without being exhaustive, we refer for instance to [START_REF] Duboc | High-order dimensionally split Lagrangeremap schemes for compressible hydrodynamics[END_REF][START_REF] Maire | A Cell-Centered Lagrangian Scheme for Two-Dimensional Compressible Flow Problems[END_REF][START_REF] Raphaël | Une Méthode Particulaire Lagrangienne de type Galerkin Discontinu Application à la Mécanique des Fluides et l'Interaction Laser/Plasma[END_REF][START_REF] Billaud Friess | A second order antidi usive Lagrange-remap scheme for two-component ows[END_REF][START_REF] Bürger | Antidi usive and Random-Sampling Lagrangian-Remap Schemes for the Multiclass Lighthill-Whitham-Richards Tra c Model[END_REF] for other interesting studies in this framework.

While details for the Lagrange-projection splitting are given in the following section, let us present now the mathematical model we are interested in. It is composed of two di erent ones, the hydrodynamic and the morphodynamic model. The former is simply given by the well-known shallow water system, which is derived from the Navier-Stokes system under the hypothesis that the vertical scale is much smaller than the horizontal dimension. As such, it is composed of the continuity and momentum equations, which are expressed as in the following

∂ t h + ∂ x (hu) = 0 ∂ t (hu) + ∂ x (hu 2 + gh 2
2 ) = -gh∂ x z

(1.1) * Laboratoire de Mathématiques de Versailles, UMR 8100, Université de Versailles Saint-Quentin-en-Yvelines, UFR des Sciences, bâtiment Fermat, 45 avenue des Etats-Unis, 78035 Versailles cedex, France. E-mail: christophe.chalons@uvsq.fr † Laboratoire de Mathématiques de Versailles, UMR 8100, Université de Versailles Saint-Quentin-en-Yvelines, UFR des Sciences, bâtiment Fermat, 45 avenue des Etats-Unis, 78035 Versailles cedex, France. E-mail: alessia.del-grosso2@uvsq.fr 1 where h(x, t) > 0 is the water depth, u(x, t) the averaged velocity and z(x, t) the bed level. In particular, H = h + z is the free surface elevation. Finally, g is the gravitational acceleration, t > 0 represents the time and x the axial coordinate. Let us brie y recall that this system is strictly hyperbolic with real eigenvalues given by u ± c with the sound speed c = √ gh. Then, we aim to simulate the interaction between the sediments and the ow, thus we consider the topography z not constant in time and we make use of the so-called Exner equation, which reads ∂ t z + ζ∂ x q b = 0.

(1.2)

Here q b = q b (h, u) is the solid transport discharge, ζ = 1 1-ρ0 and ρ 0 is the porosity of the sediment layer. The coupling of (1.1) and (1.2) leads to the nal system,

     ∂ t h + ∂ x (hu) = 0 ∂ t (hu) + ∂ x (hu 2 + gh 2
2 ) + gh∂ x z = 0 ∂ t z + ζ∂ x q b = 0, (1.3) which in compact form reads

∂ t Q + ∂ x F(Q) + A(Q)∂ x Q = O where Q =   h hu z   , F(Q) =   hu hu 2 + p ζq b   , A(Q) =   0 0 0 0 0 gh 0 0 0  
with the pressure term p = gh 2 2 . For details about shallow-water equations with and without not-constant in time bed sediment, we refer for instance to [START_REF] Audusse | Sediment transport modelling : Relaxation schemes for Saint-Venant -Exner and three layer models[END_REF][START_REF] Briganti | An e cient and exible solver for the simulation of the morphodynamics of fast evolving ows on coarse sediment beaches[END_REF][START_REF] Castro Díaz | Sediment transport models in Shallow Water equations and numerical approach by high order nite volume methods[END_REF][START_REF] Berthon | An analytical solution of Shallow Water system coupled to Exner equation[END_REF] and [START_REF] Vreugdenhil | Numerical Methods for Shallow Water Flow[END_REF][START_REF] Audusse | A fast and stable well-balanced scheme with hydrostatic reconstruction for shallow water ows[END_REF][START_REF] Berthon | E cient well-balanced hydrostatic upwind schemes for shallow-water equations[END_REF][START_REF] Manuel | A fully well-balanced Lagrange-Projection type scheme for the Shallow-water equations[END_REF][START_REF] Morales De Luna | On a Shallow Water Model for the Simulation of Turbidity Currents[END_REF].

It is known that there exist di erent formulations to express the solid transport discharge q b , depending on the characteristics of the sediment and the ow, for instance the Froude number, the slope of the bottom or the grain size. One of the easiest and most frequently used formulation is the well-known Grass model, which expresses the instantaneous sediment transport as a power law of the averaged velocity u, namely

q b = A g u|u| mg-1 , 1 ≤ m g ≤ 4, (1.4)
refer to [START_REF] Murillo | An Exner-based coupled model for two-dimensional transient ow over erodible bed[END_REF][START_REF] Audusse | A simple three-wave approximate Riemann solver for the Saint-Venant-Exner equations[END_REF]. This deterministic Grass formulation is suitable to model non-cohesive granular sediment. Here A g ∈ [0, 1] is a dimensional calibration constant which is usually measured experimentally and expresses the kind of interaction between the uid and the sediment, whose strength increases as A g approaches to 1. In particular, the value of the constant A g is related to factors as the grain size and the kinematic viscosity. Finally, we set m g = 3.

In this work we will take into account only the Grass model, but for details about other possible formulations for the solid transport discharge, refer to [START_REF] Murillo | An Exner-based coupled model for two-dimensional transient ow over erodible bed[END_REF]. It is important to stress that, depending on the formulation used for q b , system (1.3) may be hyperbolic or not. In particular, in [START_REF] Castro Díaz | Sediment transport models in Shallow Water equations and numerical approach by high order nite volume methods[END_REF] it has been con rmed that considering the Grass formula leads to a strictly hyperbolic system with all real eigenvalues. Indeed, de ning the quantities

a 1 = -2u, a 2 = u 2 -c 2 (1 + ζ∂ hu q b ) and a 3 = -ζc 2 ∂ h q b ,
one can easily see that the eigenvalues are given by the solution of the following equation

λ 3 + a 1 λ 2 + a 2 λ + a 3 = 0.
Hence, the three eigenvalues read

λ k = 2 √ -p cos( θ + 2kπ 3 ) - a 1 3 with k = 0, 1, 2. (1.5)
where ).

p = 3a 2 -a 2 1 9 , r = 9a 1 a 2 -27a 3 + 2a
We remark that, in order to have real eigenvalues, we need p 3 + r 2 ≤ 0, which can be proved in the case of the Grass model.

As mentioned at the very beginning of this work, we are also looking for a well-balanced numerical scheme, namely able to preserve the smooth stationary solutions of the system, that is to say the steady states which satisfy the ordinary di erential equations

∂ x (hu) = 0, ∂ x (hu 2 + gh 2 2
) + gh∂ x z = 0 and ∂ x q b = 0, and obey to q = hu = constant = q 0 , q 2 0 2h 2 + g(h + z) = constant and q b = constant. Note that if we use the Grass formula for q b , then the only non trivial stationary solutions are given by u = 0 and h + z = constant, (1.6) which is called "lake at rest" equilibrium. For well-balanced schemes for the shallow-water equations, see for instance to [START_REF] Audusse | A fast and stable well-balanced scheme with hydrostatic reconstruction for shallow water ows[END_REF][START_REF] Berthon | E cient well-balanced hydrostatic upwind schemes for shallow-water equations[END_REF][START_REF] Michel-Dansac | A well-balanced scheme for the shallowwater equations with topography[END_REF][START_REF] Michel-Dansac | A well-balanced scheme for the shallowwater equations with topography or Manning friction[END_REF], while for well-balanced methods in the Lagrange-projection formalism we refer to [START_REF] Chalons | A large time-step and well-balanced Lagrange-Projection type scheme for the Shallow-water equations[END_REF][START_REF] Manuel | A fully well-balanced Lagrange-Projection type scheme for the Shallow-water equations[END_REF][START_REF] Morales | High order fully well-balanced Lagrange-Projection scheme for Shallow-water[END_REF]. As far as the evolution in time of the bed sediment is considered as well, for general well-balanced schemes we refer for instance to [START_REF] Audusse | A simple three-wave approximate Riemann solver for the Saint-Venant-Exner equations[END_REF][START_REF] Liu | Well-balanced central-upwind scheme for a fully coupled shallow water system modeling ows over erodible bed[END_REF].

Outline of the paper. For the sake of clarity, let us brie y give the paper structure. In the next section 2, the Lagrangeprojection splitting is presented considering both Eulerian and Lagrangian variables. A brief summary of the approximate Riemann solver for the acoustic system is outlined as well. In section 3, both rst and second-order numerical schemes are described, distinguishing between acoustic and transport step. Then, the two-dimensional extension of the mathematical model and numerical method is illustrated in section 4. Finally, section 5 is exploited to present the numerical results and validate our numerical schemes. Concluding, nal remarks are drawn in section 6.

Operator splitting and Lagrangian coordinates

In this section we brie y explain the decomposition which entails the splitting of system (1.3) into two di erent ones, the so-called acoustic and transport systems. The former takes into account the acoustic e ects of the model and the source term related to the topography, while the latter the transport phenomena. We will see further on that this splitting can be interpreted as a Lagrange-projection one, as we rst formulate the shallow-water system in Lagrangian coordinates (acoustic step) and then we project the solution into Eulerian coordinates (transport step). For more the details about this decomposition, the reader can refer to [START_REF] Manuel | A fully well-balanced Lagrange-Projection type scheme for the Shallow-water equations[END_REF][START_REF] Chalons | A large time-step and well-balanced Lagrange-Projection type scheme for the Shallow-water equations[END_REF][START_REF] Morales | High order fully well-balanced Lagrange-Projection scheme for Shallow-water[END_REF].

Considering rst only the shallow water system, it can be reformulated as

∂ t h + h∂ x u + u∂ x h = 0 ∂ t (hu) + hu∂ x u + u∂ x (hu) + ∂ x ( gh 2 2 ) = -gh∂ x z,
where we used the chain rule for space derivatives. Therefore, the acoustic and transport system are respectively given by

∂ t h + h∂ x u = 0 ∂ t (hu) + hu∂ x u + ∂ x ( gh 2 2 ) = -gh∂ x z, (2.1)
and

∂ t X + u∂ x X = 0 (2.2)
with X = h and X = hu. We also observe that the acoustic system (2.1) can be expressed as

∂ t τ -∂ m u = 0 ∂ t u + ∂ m p = -g τ ∂ m z (2.3)
where we have introduced the unknown τ = 1 h and the mass variable m such that 1 h ∂ x = ∂ m , see [START_REF] Manuel | A fully well-balanced Lagrange-Projection type scheme for the Shallow-water equations[END_REF][START_REF] Chalons | Large Time Step and Asymptotic Preserving Numerical Schemes for the Gas Dynamics Equations with Source Terms[END_REF][START_REF] Chalons | An all-regime Lagrange-Projection like scheme for 2D homogeneous models for two-phase ows on unstructured meshes[END_REF][START_REF] Chalons | An All-Regime Lagrange-Projection Like Scheme for the Gas Dynamics Equations on Unstructured Meshes[END_REF][START_REF] Chalons | A large time-step and well-balanced Lagrange-Projection type scheme for the Shallow-water equations[END_REF]. Moreover, it is easy to nd that the eigenvalues of system (2.3) are ±h √ gh = ±hc.

About the Exner equation. Let us now consider the Exner equation as well, for which one could imagine at least three possibilities for numerical treatment. The rst one would account for it at the acoustic level, the second one directly inside the transport step, and the third one by splitting it inside both steps. However, the issue of coupling the Exner equation and the shallow water system has been vastly studied in the literature. In particular, it is known that a fully decoupled scheme can lead to a numerical method which produces spurious oscillations inside the numerical solutions; this issue has been clearly presented in [START_REF] Cordier | Bedload transport in shallow water models: why splitting (may) fail, how hyperbolicity (can) help[END_REF]. However, it is not even necessary to consider a fully coupled scheme in order to avoid this problem; indeed a weakly coupling of the equations at the numerical level can lead to satisfying results, see for instance [START_REF] Audusse | A simple three-wave approximate Riemann solver for the Saint-Venant-Exner equations[END_REF]. In this work we mainly focus on weakly coupled numerical approach. Moreover, in the following we completely take into account the Exner equation in the projection step.

Finally, in this rst part of this work, the acoustic and transport systems respectively read

     ∂ t τ -∂ m u = 0 ∂ t u + ∂ m p = -g τ ∂ m z ∂ t z = 0 (2.4) and      ∂ t h + u∂ x h = 0 ∂ t (hu) + u∂ x (hu) = 0 ∂ t z + ζ∂ x q b = 0. (2.5)
Hence, the numerical strategy will be composed of two steps:

1. Take into account the acoustic e ects of the model by solving system (2.4);

2. Consider and solve the transport system (2.5);

System (1.1) in Lagrangian coordinates. In order to interpret the strategy as a Lagrange-Projection one, we rst de ne the uid particle ξ and the characteristic curves

∂x ∂t (ξ, t) = u(x(ξ, t), t) x(ξ, 0) = ξ (2.6)
which de ne the trajectory : t → x(ξ, t), of ξ as the time goes on. Therefore, any function : (x, t) → ϕ(x, t) in Eulerian coordinates can be written in Lagrangian coordinates, φ(ξ, t) = ϕ(x(ξ, t), t).

Let us now introduce the volume ratio

L(ξ, t) = ∂x ∂ξ (ξ, t) (2.7) such that ∂L ∂t (ξ, t) = ∂ ξ u(x(ξ, t), t) L(ξ, 0) = 1. (2.8)
Consequently, we note that

∂ t L(ξ, t) = ∂ ξ u(x(ξ, t), t) = ∂ ξ ū(ξ, t),
and thus

∂ ξ φ(ξ, t) = L(ξ, t)∂ x ϕ(x, t) and ∂ t φ(ξ, t) = ∂ t ϕ(x, t) + u(x, t)∂ x ϕ(x, t).
We have now all the ingredients to write (1.1) in Lagrangian coordinates. More precisely, observing that smooth solutions of (1.1) satisfy

L(∂ t h + u∂ x h + h∂ x u) = 0 L(∂ t (hu) + u∂ x (hu) + hu∂ x u + ∂ x p + gh∂ x z) = 0 and L∂ t h + h∂ t L = 0 L∂ t hu + hu∂ t L + ∂ ξ p + g h∂ ξ z = 0 we get ∂ t (Lh) = 0 ∂ t (Lhu) + ∂ ξ p = -g h∂ ξ z. (2.9)
Notice that in the following sections, we shall omit the bar over the Lagrangian functions when there can be no confusion. This new formulation of system (1.1) exploiting the Lagrangian coordinates makes the above numerical strategy based on an acoustic-transport splitting strictly equivalent to a Lagrangian-projection splitting, that can be summarized as in the following,

1. Solve the system (2.9) in Lagrangian coordinates;

2. Project the solution in Eulerian coordinates solving the transport system (2.5).

To conclude this section, we derive from (2.9) an evolution equation for Lu that will be useful in the next sections. Removing the bars and using the discharge equation on Lu, we have in particular

h∂ t (Lu) + Lu∂ t h + ∂ ξ p + gh∂ ξ z = 0 and, since 0 = ∂ t (Lh) = h∂ t L + L∂ t h = h∂ ξ u + L∂ t h, we get h∂ t (Lu) -hu∂ ξ u + ∂ ξ p + gh∂ ξ z = 0 and nally ∂ t (Lu) -∂ ξ u 2 2 = -g∂ ξ (h + z).
(2.10)

Note that the source term now involves the quantity h + z which is constant for stationary solutions.

A well-balanced approximate Riemann solver for the acoustic system

In this section we brie y recall the de nition of a well-balanced approximate Riemann solver proposed in [START_REF] Chalons | A large time-step and well-balanced Lagrange-Projection type scheme for the Shallow-water equations[END_REF] for the acoustic system (2.4) and that will be useful in the following. It is based on the Gallice theory [START_REF] Gallice | Solveurs simples positifs et entropiques pour les systèmes hyperboliques avec terme source[END_REF][START_REF] Gallice | Positive and entropy stable Godunov-type schemes for gas dynamics and MHD equations in Lagrangian or Eulerian coordinates[END_REF] which is an extension to balance laws of the Harten, Lax and van Leer formalism [START_REF] Harten | On upstream di erencing and Godunov-type schemes for hyperbolic conservation laws[END_REF] for conservation laws. For more details about this Riemann solver, the reader can refer to [START_REF] Chalons | A large time-step and well-balanced Lagrange-Projection type scheme for the Shallow-water equations[END_REF]. Then, suppose that we want to solve (2.4) with the following Riemann initial data

U(m, t = 0) = U L if m < 0 U R if m ≥ 0
where we have set

U L =   τ L u L z L   and U R =   τ R u R z R   .
The proposed approximate solution has the following form

Û( m t ; U L , U R ) =          U L if m t < -a U * L if -a < m t < 0 U * R if 0 < m t < a U R if m t > a
where a is a constant and where the intermediate states

U * L =   τ * L u * L z L   , U * R =   τ * R u * R z R   (2.11) are de ned thanks to                        τ * L = τ L + 1 a (u * -u L ) τ * R = τ R - 1 a (u * -u R ) u * L = u * R = u * u * = 1 2 (u L + u R ) - 1 2a (Π R -Π L ) - M 2a Π * = 1 2 (Π L + Π R ) - a 2 (u R -u L ) with M = g 2 ( 1 τ L + 1 τ R )(z R -z L ).
Observe that Π is a new variable introduced to be able to de ne the approximate Riemann solver. In particular, Π can be interpreted as a linearization of the pressure term p and its initial data are well-prepared in the sense that Π = p. Then, it is clear that M = 0 if z L = z R so that the classical form of an approximate Riemann solver for a system of conservation law is recovered, while X * L = X L and X * R = X R with X = τ, u, when the "lake at rest" stationary conditions are satis ed, namely

u L = u R = 0, h L + z L = h R + z R .
In that sense, the proposed approximate Riemann solver is said to be well-balanced. Note that this is true whatever the definition of the constant a is. In practice, we will choose a = max((hc) L , (hc) R ) according the well-known subcharacteristic stability condition. For more details, we refer again to [START_REF] Chalons | A large time-step and well-balanced Lagrange-Projection type scheme for the Shallow-water equations[END_REF].

Numerical methods

Before getting into the heart of the matter, we give few details about the time and space discretizations we use in the following sections. Given a constant time step ∆t, we de ne the intermediate times by t n = n∆t for n ∈ N. Then, the mesh interfaces are x j+1/2 = j∆x for j ∈ Z, where ∆x is the constant space step, while x j is the center of the cell [x j-1/2 , x j+1/2 ). Hence, given a variable ϕ, we denote ϕ n j its constant average approximation on each cell [x j-1/2 , x j+1/2 ) at time t n , namely

ϕ n j ≈ 1 ∆x x j+1/2 x j-1/2 ϕ(x, t n ) dx,
with n ∈ N and j ∈ Z. At last, note that, regarding the mass variable m, we use ∆m j = h n j ∆x for all j. Given the sequence {ϕ n j } j , we now look for its approximation at the next time level t n+1 , namely {ϕ n+1 j } j . At this stage, we are able to present the numerical schemes, starting with the rst-order method and proceeding with the secondorder one. For each of them, we will prove the well-balanced condition as well. Let us recall that the numerical schemes are divided into two di erent steps. First we have the acoustic step, in which we numerically solve system (2.4). Then, we exploit its solution as the initial condition for solving system (2.5). We can sum up this procedure in the following way, 1. Acoustic step: solve system (2.4) in order to update Q n to Q n+1-; 2. Transport step: nd Q n+1 from Q n+1-by the approximation of the solution of system (2.5).

Note that in the rst step, we implicitly use the change of variable U = U(Q) to rst de ne U n from Q n before solving (2.4), and then the change of variables Q = Q(U) to de ne Q n+1-from U n+1-.

First-order scheme

Here we give the details for the rst-order scheme distinguishing between the acoustic and transport steps.

Acoustic step and Lagrangian reformulation

As far as the discretization of (2.4) is concerned, we suggest to use a classical Godunov-type method based on the wellbalanced approximate Riemann solver proposed in section 2.1. As usual, it simply consists in averaging on each cell the juxtaposition of the approximate Riemann solutions set at each interface. Therefore, it follows after easy calculations that the numerical discretization of the acoustic relaxation system (2.4) can be formulated as

τ n+1- j = τ n j + ∆t ∆mj (u * j+ 1 2 -u * j-1 2 ) u n+1- j = u n j -∆t ∆mj (Π * j+ 1 2 -Π * j-1 2 ) -∆t{ g τ ∂ m z} n j (3.1)
where we have set

u * j+ 1 2 = u * j+ 1 2 (Q n j , Q n j+1 ) = 1 2 (u n j+1 + u n j ) - 1 2a n j+ 1 2 (Π n j+1 -Π n j ) - M n j+1/2 2a j+1/2 , Π * j+ 1 2 = Π * j+ 1 2 (Q n j , Q n j+1 ) = 1 2 (Π n j+1 + Π n j ) - a n j+ 1 2 2 (u n j+1 -u n j ), (3.2)
a n j+ 1 2 = max((hc) n j , (hc) n j+1
), while regarding the source term we have

{ g τ ∂ m z} n j = 1 2 ∆m j+1/2 ∆m j { g τ ∂ m z} n j+1/2 + ∆m j-1/2 ∆m j { g τ ∂ m z} n j-1/2 with { g τ ∂ m z} n j+1/2 = M n j+1/2 ∆m j+1/2 ∀j (3.3)
where ∆m j+1/2 = (∆m j + ∆m j+1 )/2, ∆m j = τ n j ∆x and

M j+1/2 = g 2 ( 1 τ n j + 1 τ n j+1 )(z n j+1 -z n j ).
Let us note that z remains constant in this step, thus z n+1- j = z n j for all j.

Lagrangian reformulation of (3.1). Let us observe that (3.1) reveals to be strictly equivalent to

L n+1- j h n+1- j = L n j h n j L n+1- j (hu) n+1- j = L n j (hu) n j -∆t ∆x (Π * j+ 1 2 -Π * j-1 2 ) + ∆ts n j (3.4)
where we have set

L n+1- j = L n j + ∆t ∆x (u * j+ 1 2 -u * j-1 2 ) with L n j = 1 (3.5) and s = -gh∂ x z, s n j = 1 2 s n j+1/2 + s n j-1/2 and s n j+1/2 = - M n j+1/2 ∆x ∀j.
Considering that the Lagrangian variable ξ is discretized using the same mesh step as the one we used for x, namely ∆ξ = ∆x, and ξ j+1/2 = x j+1/2 , ξ j = x j for all j, it is clear that (3.4) and (3.5) respectively approximate (2.9) and (2.8). This Lagrangian reformulation turns out to be crucial in order to derive hereafter a second-order extension of the propose numerical scheme. Note that we still have of course z n+1- j = z n j for all j.

Projection step

As already mentioned, the aim of this step is to turn into Eulerian coordinates the solution obtained at the end of the previous step thanks to (3.4) and given in Lagrangian coordinates. This amounts to solve the transport system (2.5) which also contains the evolution equation for z. Recall indeed that z stayed constant in the rst step.

In order to project X = h, hu on the Eulerian grid, we use the following identity ξ2 ξ1 L(ξ, t)X(ξ, t)dξ =

x(ξ2,t)

x(ξ1,t)

X(x, t)dx

where we recall that the trajectories t → x(ξ, t) and the volume ratio L(ξ, t) are de ned by (2.6) and (2.7). Therefore, it is natural to de ne ξj+1/2 such that x( ξj+1/2 , t n+1 ) = x j+1/2 and x( ξj+1/2 , t n ) = ξj+1/2 for all j and to write

X n+1 j = 1 ∆x x j+ 1 2 x j-1 2 X(x, t n+1 )dx = 1 ∆x x( ξj+ 1 2 ,t n+1 ) x( ξj-1 2 ,t n+1 ) X(x, t n+1 )dx = 1 ∆x ξj+ 1 2 ξj-1 2 L(ξ, t n+1-)X(ξ, t n+1-)dξ. (3.6)
Splitting the last integral into three parts, namely

X n+1 j = 1 ∆x ξ j-1 2 ξj-1 2 L(ξ, t n+1-)X(ξ, t n+1-)dξ+ + 1 ∆x ξ j+ 1 2 ξ j-1 2 L(ξ, t n+1-)X(ξ, t n+1-)dξ + 1 ∆x ξj+ 1 2 ξ j+ 1 2 L(ξ, t n+1-)X(ξ, t n+1-)dξ, (3.7)
and approximating ξj+1/2 by

x j+1/2 = x( ξj+1/2 , t n+1 ) x( ξj+1/2 , t n ) + ∆t∂ t x( ξj+1/2 , t n ) ξj+1/2 + ∆tu * j+1/2 ,
it is natural to set, using rst-order approximations of the integrals,

X n+1 j = ξ j-1 2 -ξj-1 2 ∆x (LX) n+1 j-1/2 + (LX) n+1- j + ξj+ 1 2 -ξ j+ 1 2 ∆x (LX) n+1 j+1/2 (3.8)
where for all j

(LX) n+1 j+1/2 = (LX) n+1- j if u * j+1/2 ≥ 0 (LX) n+1- j+1 if u * j+1/2 < 0. After easy manipulations, (3.8) is equivalent to X n+1 j = (LX) n+1- j - ∆t ∆x u * j+ 1 2 (LX) n+1- j+ 1 2 -u * j-1 2 (LX) n+1- j-1 2 . (3.9)
Therefore, taking X = h and X = hu concludes the projection on the Eulerian grid.

Let us now consider the Exner equation for which we propose two di erent strategies. On one hand, we simply update the topography as in the following

z n+1 j = z n j -ζ ∆t ∆x u * j+ 1 2 q b u n j+ 1 2 -u * j-1 2 q b u n j-1 2 , (3.10)
with

( q b u ) n j+1/2 =        ( q b u )(u n j+1 ) if u * j+ 1 2 ≤ 0 ( q b u )(u n j ) if u * j+ 1 2 > 0,
where we see that the numerical uxes are evaluated simply using the solution at time t n , and not the one obtained at the end of the Lagrangian step. This means in some sense that the evolution of the topography is not coupled with the one hydrodynamic model from a numerical point of view and this is the reason why this strategy is said to be decoupled. At last, note that q b u is well-de ned and actually depends on u only, as we suppose the solid transport discharge to be given by the Grass formula with m g = 3. On the other hand, we propose a weakly coupled strategy in which we exploit the solution obtained at time t n+1-, by setting

z n+1 j = z n j -ζ ∆t ∆x u * j+ 1 2 q b u n+1- j+ 1 2 -u * j-1 2 q b u n+1- j-1 2 . ( 3.11)
At this stage, notice that it would be tempting to de ne

q b u n+1- j+ 1 2
using the velocity u n+1-simply de ned by u n+1-= Lhu n+1- Lh n+1-. However, even if this option results to be natural, dividing by Lh rises di culties when considering the secondorder extension. For this reason, we prefer to set

( q b u ) n+1- j+1/2 =        ( q b u ) (Lu) n+1- j+1 if u * j+ 1 2 ≤ 0 ( q b u ) (Lu) n+1- j if u * j+ 1 2 > 0,
where a possible discretization of the evolution equation (2.10) for Lu reads

(Lu) n+1- j = (Lu) n j + ∆t 2∆x ((u * j+ 1 2 ) 2 -(u * j-1 2 ) 2 ) + ∆t ŝj+ 1 2 + ŝj-1 2 2 (3.12)
with

ŝj+ 1 2 = -g (h + z) n j+1 -(h + z) n j ∆x .
As a last remark we observe that, since we are also interested in a 2D formulation and (2.10) could not be extended in two dimensions as it is (see hereafter), we will propose an alternative approximation of Lu which reads

(Lu) n+1- j = (Lu) n j + u n j + u n j+1 2 ∆t ∆x (u * j+ 1 2 -u * j-1 2 ) + ∆t ŝj+ 1 2 + ŝj-1 2 2 .
Both formulations turn out to give the same results in 1D.

Overall scheme and well-balanced property

Next, we give an equivalent formulation of our rst-order scheme which takes into account both the acoustic and transport steps. This formulation is interesting in the sense that it clearly shows that the scheme is indeed conservative when there is no source term. More precisely, considering together (3.4) and (3.9) we easily get

     h n+1 j = h n j - ∆t ∆x u * j+ 1 2 (Lh) n+1- j+ 1 2 -u * j-1 2 (Lh) n+1- j-1 2 (hu) n+1 j = (hu) n j - ∆t ∆x u * j+ 1 2 (Lhu) n+1- j+ 1 2 + Π * j+ 1 2 -u * j-1 2 (Lhu) n+1- j-1 2 -Π * j-1 2 + ∆ts n j (3.13)
while the evolution equations (3.10) and (3.11) for z are clearly conservative. Let us now prove the well-balanced property.

Theorem 1. The rst-order numerical scheme with updating formula (3.13) and (3.10) or (3.11) preserves the "lake at rest" stationary solution (1.6).

Proof. Assuming to be under the "lake at rest" condition, that is to say u n j = 0 and h n j + z n j = constant, for all j, it is straightforward to demonstrate that this stationary solution is preserved. Indeed, few algebraic computations show that u * j+ 1 2 = 0 as u n j = u n j+1 = 0 and

Π n j+1 -Π n j = g 2 ((h n j+1 ) 2 -(h n j ) 2 ) = g 2 (h n j + h n j+1 )(h n j+1 -h n j ) = - g 2 (h n j + h n j+1 )(z n j+1 -z n j ).
Similarly, it can be proved that 1 ∆x (Π *

j+ 1 2 -Π * j-1 2
) = s n j and thus (Lhu) n+1- j = (Lhu) n j = (hu) n j = 0. We can also note that (Lu) n+1- j = 0 as (h + z) n j = (h + z) n j+1 for all j. Finally, it is easily seen that under the property u * j+1/2 = 0 for all j, then the transport step gives h n+1 j = h n j , (hu) n+1 j = 0, and z n+1 j = z n j which concludes the proof.

Second-order scheme

We now explain how to reach the second order of accuracy in both space and time. While increasing the order of accuarcy is a standard process, the key issue is to preserve the well-balanced property.

In order to construct a high-order approximation in space, we will make use of classical rst-order polynomial reconstructions, but applied to the so-called uctuations which is non-standard and has been introduced in [START_REF] Morales | High order fully well-balanced Lagrange-Projection scheme for Shallow-water[END_REF] to combine both the well-balanced property and the higher order property.

Regarding the second-order discretization in time, we simply consider Runge-Kutta TVD scheme at second order [START_REF] Gottlieb | Total variation diminishing RUNGE-KUTTA schemes[END_REF]. In particular, we apply it to the overall scheme (Lagrangian and remap step together) in order to avoid di usion due to the splitting.

Lagrangian step

In order to reach the second order of accuracy in space, we begin by de ning at time t n and for each cell j a stationary solution denoted by x → Q n,e j (x) and de ned for all x by

h n,e j (x) = h n j + z n j -z n (x), u n,e j (x) = u n j and z n,e j (x) = z n (x), (3.14)
where x → z n (x) is nothing but the piecewise constant approximation of z at time t n , namely such that z n (x) = z n j for all x in [x j-1/2 , x j+1/2 ). Such a reconstructed solution satis es the in-cell conservativity property 1 ∆x

x j+ 1 2 x j-1 2 Q n,e j (x)dx = Q n j .
Next, we follow [START_REF] Morales | High order fully well-balanced Lagrange-Projection scheme for Shallow-water[END_REF] and introduce the so-called j-uctuations de ned as

D n k,j = Q n k - 1 ∆x x k+1/2 x k-1/2
Q n,e j (x)dx, for all k. Observe that D n j,j = 0 by construction, while D n k,j = 0 for all k if the approximate solution at time t n satis es the "lake at rest" condition (1.6).

At last, for each cell I j we make use of a reconstructed polynomial vector P n j (x) de ned by

P n j (x) = Q n j + ∆ n j (x -x j ), where ∆ n j = ∆ n j (D n j-1,j , D n j,j , D n j+1,j
) is the ENO [START_REF] Toro | Riemann Solvers and Numerical Methods for Fluid Dynamics[END_REF] or the MINMOD [START_REF] Toro | PRICE: Primitive Centred Schemes for Hyperbolic Systems[END_REF] slope applied to the uctuations. The numerical uxes u * j+ 1 2 and Π * j+ 1 2 are then de ned in a very classical way using the interfaces values

Q n j+ 1 2 L = P n j (x j+ 1 2 ) and Q n j+ 1 2 R = P n j+1 (x j+ 1
2 ), and formula (3.2), namely

u * j+ 1 2 = u * j+ 1 2 (Q n j+ 1 2 L , Q n j+ 1 2 R ) and Π * j+ 1 2 = Π * j+ 1 2 (Q n j+ 1 2 L , Q n j+ 1 2 R ), (3.15) with a n j+ 1 2 = max h n j+ 1 2 ,L c n j+ 1 2 ,L , h n j+ 1 2 ,R c n j+ 1 2 ,R .
Regarding the source term, once again we exploit formulas (3.3). Let us note that, thanks to formula (3.14), z j+ 1 2 L = z j and z j+ 1 2 R = z j+1 as the uctuations related to the topography are null. Finally, the discretization of the Lagrangian system (2.9) reads as in the rst-order step, namely

L n+1- j h n+1- j = L n j h n j L n+1- j (hu) n+1- j = L n j (hu) n j -∆t ∆x (Π * j+ 1 2 -Π * j-1 2
) + ∆ts n j .

2D shallow water system reads

∂ t h + ∇ • (hu) = 0 ∂ t (hu) + ∇ • (hu ⊗ u) + ∇p = -gh∇z, (4.1)
while the Exner equation is given by ∂ t z + ζ∂ x q b,x + ζ∂ y q b,y = 0 where q b,x and q b,y are the solid transport discharges in the x and y direction respectively. Exploiting once again the Grass model, their formula are the following,

q b,x = A g u(u 2 + v 2 ) and q b,y = A g v(u 2 + v 2 ).
Notice that the "lake at rest" stationary solutions now satisfy

u = 0, v = 0 and ∇(h + z) = 0.
It is still possible to consider a Lagrangian formulation of these equations. More precisely, let us introduce the Lagrangian coordinates (ξ 1 , ξ 2 ) and consider the map : (ξ 1 , ξ 2 ) → (x, y), with x = x(ξ 1 , ξ 2 , t), y = y(ξ 1 , ξ 2 , t) and such that ∂x ∂t = u(x, y, t) and ∂y ∂t = v(x, y, t),

x(ξ 1 , ξ 2 , 0) = ξ 1 , y(ξ 1 , ξ 2 , 0) = ξ 2 .
We also assume that for each t > 0, this map is invertible and its Jacobian (determinant of the Jacobian matrix) is given by

L(ξ 1 , ξ 2 , t) = ∂ ξ1 x ∂ ξ2 x ∂ ξ1 y ∂ ξ2 y
with L(ξ 1 , ξ 2 , 0) = 1 and, after easy calculations,

∂L(ξ 1 , ξ 2 , t) ∂t = L∇ • u = L∂ x u + L∂ y v. (4.2)
Then, it can be shown that the Lagrangian formulation of the system writes

∂ t (Lh) = 0 ∂ t (Lhu) + L∇p = -ghL∇z, (4.3)
where the gradient is still taken with respect to the Eulerian variables (x, y), while on the other hand we also have

∂ t (Lu) -uL∇ • u + gL∇(h + z) = 0,
which will be useful hereafter. We refer the reader to [START_REF] Raphaël | Une Méthode Particulaire Lagrangienne de type Galerkin Discontinu Application à la Mécanique des Fluides et l'Interaction Laser/Plasma[END_REF][START_REF] Maire | A Cell-Centered Lagrangian Scheme for Two-Dimensional Compressible Flow Problems[END_REF] for more details about Lagrangian coordinates in 2-dimensions.

We now give the basic formulas based on a dimensional splitting and introduce some notations. First of all, the computational domain Ω ⊂ R is divided into M x × M y rectangular cells with constant space steps ∆x and ∆y in the x and y directions respectively. Then, the mesh interfaces are given by x i+1/2 for i ∈ {0, . . . , M x } and y j+1/2 for j ∈ {0, . . . , M y }. Thus, ϕ n i,j denotes the piecewise constant approximation of the variable ϕ in the cell

[x i-1/2 , x i+1/2 ) × [y j-1/2 , y j+1/2 ) at time t n , namely ϕ n i,j ≈ 1 ∆x x i+1/2 x i-1/2 y j+1/2 y j-1/2
ϕ(x, y, t n ) dx dy.

Acoustic step

Following the 1D scheme, the numerical approximation of (4.3) is taken to be

         (Lh) n+1- i,j = (Lh) n i,j (Lhu) n+1- i,j = (Lhu) n i,j - ∆t ∆x (Π * i+ 1 2 ,j -Π * i-1 2 ,j ) + ∆ts 1,n i,j (Lhv) n+1- i,j = (Lhv) n i,j - ∆t ∆y (Π * i,j+ 1 2 -Π * i,j-1 2 ) + ∆ts 2,n i,j (4.4)
where we have set

L n+1- i,j = L n i,j + ∆t ∆x (u * i+ 1 2 ,j -u * i-1 2 ,j ) + ∆t ∆y (v * i,j+ 1 2 -v * i,j-1 2 ) (4.5) with u * i+ 1 2 ,j = 1 2 (u n i+1,j + u n i,j ) - 1 2a n i+ 1 2 ,j (Π n i+1,j -Π n i,j )- M n i+1/2,j 2a i+1/2,j , v * i,j+ 1 2 = 1 2 (v n i,j+1 + v n i,j ) - 1 2a n i,j+ 1 2 (Π n i,j+1 -Π n i,j )- M n i,j+1/2 2a i,j+1/2 , Π * i+ 1 2 ,j = 1 2 (Π n i+1,j + Π n i,j ) - a n i+ 1 2 ,j 2 (u n i+1,j -u n i,j ), Π * i,j+ 1 2 = 1 2 (Π n i,j+1 + Π n i,j ) - a n i,j+ 1 2 2 (v n i,j+1 -v n i,j ), (4.6)
with a n i+1/2,j = max((hc) n i,j , (hc) n i+1,j ), a n i,j+1/2 = max((hc) n i,j , (hc) n i,j+1 ), while regarding the source term we have for all j s 1,n i,j =

1 2 s n i+1/2,j + s n i-1/2,j with s n i+1/2,j = - M n i+1/2,j ∆x s 2,n i,j = 1 2 s n i,j+1/2 + s n i,j-1/2 with s n i,j+1/2 = - M n i,j+1/2 ∆y with M n i+1/2,j = g 2 ( 1 τ n i,j + 1 τ n i+1,j )(z n i+1,j -z n i,j ), M n i,j+1/2 = g 2 ( 1 τ n i,j + 1 τ n i,j+1 )(z n i,j+1 -z n i,j ).
It is clear that the numerical scheme (4.6) is a natural extension of the one used for the one-dimensional system.

Projection step

As before, the second step of the Lagrange-projection scheme consists in projecting the solution obtained at the end of the acoustic step onto the Eulerian grid, that is to say in approximating the transport system

∂ t ϕ + u • ∇ϕ = 0 or equivalently ∂ t ϕ + ∇ • (ϕu) -ϕ∇ • u = 0,
where we assumed ϕ = h, hu, hv. Here and analogously to the 1D formulation (3.9), we simply set

ϕ n+1 i,j =(Lϕ) n+1- i,j - ∆t ∆x u * i+ 1 2 ,j (Lϕ) n+1- i+ 1 2 ,j -u * i-1 2 ,j (Lϕ) n+1- i-1 2 ,j - ∆t ∆y v * i,j+ 1 2 (Lϕ) n+1- i,j+ 1 2 -v * i,j-1 2 (Lϕ) n+1- i,j-1 2 , (4.7)
where

(Lϕ) n+1- i-1 2 ,j = (Lϕ) n+1- i-1,j if u * i-1 2 ,j > 0 (Lϕ) n+1- i,j if u * i-1 2 ,j ≤ 0, and (Lϕ) n+1- i,j-1 2 = (Lϕ) n+1- i,j-1 if v * i,j-1 2 > 0 (Lϕ) n+1- i,j if v * i,j-1 2 ≤ 0.
The Exner equation. As a natural extension of (3.11), we set

z n+1 i,j = z n i,j -ζ ∆t ∆x u * i+ 1 2 ,j q b,x u n+1- i+ 1 2 ,j -u * i-1 2 ,j q b,x u n+1- i-1 2 ,j -ζ ∆t ∆y v * i,j+ 1 2 q b,y v n+1- i,j+ 1 2 -v * i,j-1 2 q b,y v n+1- i,j-1 2 with ( q b,x u ) n+1- i+1/2,j =        ( q b,x u ) (Lu) n+1- i+1,j if u * i+1/2,j ≤ 0 ( q b,x u ) (Lu) n+1- i,j if u * i+1/2,j > 0,
and

( q b,y v ) n+1- i,j+1/2 =        ( q b,y v ) (Lu) n+1- i,j+1 if v * i,j+1/2 ≤ 0 ( q b,y v ) (Lu) n+1- i,j if v * i,j+1/2 > 0,
where a possible discretization of the evolution equations for Lu and Lv read

(Lu) n+1- i,j = (Lu) n i,j + ∆t u i+1,j + u i,j 2
1 ∆x (u * i+ 1 2 ,j -u * i-1 2 ,j ) + 1 ∆y (v * i,j+ 1 2 -v * i,j-1 2 ) -∆t ŝi+ 1 2 ,j + ŝi-1 2 ,j 2
and

(Lv) n+1- i,j = (Lv) n i,j + ∆t v i,j+1 + v i,j 2
1 ∆x (u * i+ 1 2 ,j -u * i-1 2 ,j ) + 1 ∆y (v * i,j+ 1 2 -v * i,j-1 2 ) -∆t ŝi,j+ 1 2 + ŝi,j-1 2 2
where ŝi+ 1 2 ,j = g((h + z) i+1,j -(h + z) i,j )/∆x and ŝi,j+

1 2 = g((h + z) i,j+1 -(h + z) i,j )/∆y.

2D extension of the second-order scheme

We now brie y discuss the extension of the second-order scheme, distinguishing among the Exner equation and the Lagrangian and projection steps for the shallow water system. Once again we reach the second order of accuracy in time exploiting the Runge-Kutta procedure, which is applied to the Lagrangian and projection step together. As expected, the overall strategy is analogous to what we have done for the 1D case.

Regarding the Lagrangian step, we proceed as usual and compute the numerical uxes u *

i+ 1 2 ,j , v * i,j+ 1 2 , Π * i+ 1 2 ,j and Π * i,j+ 1 2
, but also the speeds a i+1/2,j and a i,j+1/2 using left and right interfaces values de ned by means of reconstructed polynomials as in the 1D case, namely

V n i+1/2L,j = V n i,j + ∆ x,t i,j ∆x 2 , V n i+1/2R,j = V n i+1,j -∆ x,t i+1,j ∆x 2 (4.8)
in the x direction, and

V n i,j+1/2L = V n i,j + ∆ y,t i,j ∆y 2 , V n i,j+1/2R = V n i,j+1 -∆ y,t i,j+1 ∆y 2 (4.9)
along the y axis, where we have set V = (h, hu, hv) T . At this stage, we only need to de ne the slopes ∆ x,t i,j , ∆ y,t i,j which are computed exactly as in the 1D case thanks to the de nition of uctuations in the x, respectively y, direction and considering that the y = y j , resp. x = x i , is xed and using reconstructed stationary solutions direction by direction. The details are left to the reader. In particular, such a strategy guarantees the well-balanced property of the numerical scheme, since the slopes turn out to be null under the "lake at rest" conditions.

Regarding the transport step, we consider a direct 2D extension of (3.17) namely

X n+1 i,j =(LX) n+1- i,j - ∆t ∆x u * i+ 1 2 ,j (LX) n+1- i+ 1 2 ,j ξ 1,i+ 1 2 + ξ1,i+ 1 2 2 -u * i-1 2 ,j (LX) n+1- i-1 2 ,j ξ 1,i-1 2 + ξ1,i-1 2 2 - ∆t ∆y v * i,j+ 1 2 (LX) n+1- i,j+ 1 2 ξ 2,j+ 1 2 + ξ2,j+ 1 2 2 -v * i,j-1 2 (LX) n+1- i,j-1 2 ξ 2,j-1 2 + ξ2,j-1 2 2 ,
where we have used clear notations which are based on classical rst-order polynomial reconstructions of the Lagrangian unknowns (LX) in each direction as in the 1D case, with X = h, hu, hv.

Finally, exploiting also the reconstructed values for Lu and Lv and using again classical notations, we suggest a direct 2D extension of (3.11) namely

z n+1 i,j = z n i,j -ζ ∆t ∆x u * i+1/2,j q b,x u n+1- i+1/2,j -u * i-1/2,j q b,x u n+1- i-1/2,j -ζ ∆t ∆y v * i,j+1/2, q b,y v n+1- i,j+1/2 -u * i,j-1/2 q b,y v n+1- i,j-1/2
with a natural de nition for the numerical uxes q b,x u and q b,y v . Again, the details are left to the reader since there is no ambiguity.

To conclude this 2D section, let us mention that both the schemes described here preserve the "lake at rest" stationary solutions. The proof is analogous to the one seen in 1D.

Numerical results

This section is devoted to the presentation of the simulations and outputs of the numerical schemes we described so far. Regarding the 1D time step value, at each time t n we compute two di erent time steps, one for the acoustic system and the other for the transport part, which respectively read

∆t ≤ CFL l ∆x max j {max(τ n j , τ n j+1 )a j+ 1 2 }
, and

∆t ≤ CFL t ∆x max j {u + j-1 2 -u - j+ 1 2 } ,
where CFL l and CFL t are respectively the CFL number for the Lagrangian and the transport system, and

u + j-1 2 = max(u * j-1 2
, 0) and u -

j+ 1 2 = min(u * j+ 1 2 , 0).
Then, the nal time step is taken as the minimum between the two. If working in 2D, the acoustic time step is automatically extended, it is enough to consider both directions, while the transport time step reads

∆t ≤ CFL t min i,j u + i-1 2 ,j -u - i+ 1 2 ,j ∆x + v + i,j-1 2 -v - i,j+ 1 2 ∆y -1
.

Finally, if not otherwise speci ed, we take ζ = 1, q b = A g u 3 with A g = 0.005 for the Exner equation, and transmissive boundary conditions. For the CFL number we use CFL l = 0.45 and CFL l = 0.25 for the rst and second order schemes respectively, while CFL t = 1. When the 1D reference solution is inserted, it is computed exploiting the second-order scheme with decoupled approximation (3.10) for the Exner equation. Then, M = 2000 cells are used, where ∆x = L M with L the length of the channel.

Test of order of accuracy

Here we test the order of accuracy of the numerical schemes described previously. Let us consider a channel of length L = 20m, A g = 0.3, m = 3. The initial condition is given by null velocity and

z IC = 0.1 -0.01e -(x-10) 2 h IC = 2 -0.1e -(x-10) 2 .
We refer to paper [START_REF] Castro Díaz | Sediment transport models in Shallow Water equations and numerical approach by high order nite volume methods[END_REF] for this test case. The reference solution is computed using M = 4096 cells and second order decoupled method. In table 1 we inserted the error in norm L 1 and the empirical order of accuracy (EOA) for the water height h, the discharge q and the topography z of both the weakly coupled and decoupled approach. We can see that both schemes reach the second order of accuracy as expected.

1D Riemann problem: dam break on movable bottom

For this Riemann problem we refer to [START_REF] Audusse | Sediment transport modelling : Relaxation schemes for Saint-Venant -Exner and three layer models[END_REF]. The length of the channel is L = 10m and the dam is placed in the middle. The ending time is t end = 1s. The initial condition is given by null velocity, at topography and water height h L = 2m if x < L/2, h R = 0.125m if x > L/2. In gure 1 we insert the numerical results given by the rst-order weakly coupled scheme exploiting two di erent meshes, in particular M = 200 and M = 2000 cells. In the second case, we observe that the rst-order solution converges towards the reference one, while for M = 200 the solution is less accurate in the plateau zone when considering the topography z. On the other hand, no spurious oscillations appear.

Next, in gure 2, we insert the results for the decoupled and weakly coupled second order numerical methods for di erent mesh sizes, M = 100, M = 200 and M = 500 cells . These two schemes give similar results and, in the topography outputs of both of them, we note some oscillations which decrease as we re ne the mesh. We also observe that these oscillations are more accentuated in the decoupled scheme outputs. This is indeed expected as in the decoupled scheme the topography is not coupled with the hydrodynamic model.

1D sub-critical and supercritical regions

For the following two numerical tests we refer to paper [START_REF] Cordier | Bedload transport in shallow water models: why splitting (may) fail, how hyperbolicity (can) help[END_REF]. As initial condition we consider the sub-critical steady state      hu(x, t = 0) = 0.5 z(x, t = 0) = 0.1(1 + e -(x-5) 2 )

u 2
2 + g(h + z) = 6.386, while the length of the channel is L = 10.0m. In gure 3 we insert the results of the rst and second-order numerical schemes computed with M = 200. We consider both A g = 0.005 (top) and A g = 0.07 (bottom). In both cases, the results are satisfying. In work [START_REF] Cordier | Bedload transport in shallow water models: why splitting (may) fail, how hyperbolicity (can) help[END_REF], the authors noted as a splitting numerical scheme could produce oscillations in the solution of this numerical test, but to reduce the CFL number could remedy the problem. By splitting numerical scheme, they mean a method which solves before the shallow water system for a xed topography, and then update the bed level according to the Exner equation. For our method it is not necessary to further decrease the CFL due to the Lagrange-projection splitting whose numerical di usion is su cient as coming from both steps. Then, in the same paper [START_REF] Cordier | Bedload transport in shallow water models: why splitting (may) fail, how hyperbolicity (can) help[END_REF], the authors presented another test case in which the oscillations of the numerical solution of the splitting method could not be removed, even reducing the CFL number. In gure 4 we show that our solution does not present any oscillation. The decoupled method outputs are not reported as they are very close to the weakly coupled ones. For this last test case, the coe cient A g is kept null until time t = 15s is reached, then the value A g = 0.0005 is used. As IC we considered

     hu(x, t = 0) = 0.6 z(x, t = 0) = 0.1(1 + e -(x-5) 2 ) h(x, t = 0) + z(x, t = 0) = 0.4.
Note that in this test case we have both sub-critical and supercritical regions.

1D "lake at rest" solution and perturbation

This numerical test is useful to check the well-balanced property of the numerical scheme as we start considering a stationary solution and then we insert some perturbations. For the numerical tests of this section, see [START_REF] Chalons | A large time-step and well-balanced Lagrange-Projection type scheme for the Shallow-water equations[END_REF]. Thus, let us consider as initial condition null velocity, h(x, t = 0) + z(x, t = 0) = 3m and z(x, t = 0) = 2 + 0.25(cos(10π(x -0.5)) + 1) if 1.4 < x < 1.6 2 otherwise.

The length of the channel is L = 2.0m. Both rst and second-order schemes maintain the steady state up to an error of order 10 -15 . As a second step, let us introduce small perturbations, namely we impose

h(x, t = 0) = 3 -z(x, t = 0) + 0.001 if 1.1 < x < 1.2 3 -z(x, t = 0) otherwise.
In gure 5 we compare the results of rst and second order numerical schemes against the reference solution. We observe that the outputs are satisfying and in agreement with the ones showed in work [START_REF] Chalons | A large time-step and well-balanced Lagrange-Projection type scheme for the Shallow-water equations[END_REF], no spurious oscillations appear. Clearly, the second order results are less di usive than the rst-order one.

Circular dam break on wet bed

Let us now consider test problems in 2 dimensions. In this rst simple test, the Exner equation is not taken into account, refer to [START_REF] Seaïd | Non-oscillatory relaxation methods for the shallow-water equations in one and two space dimensions[END_REF]. The domain is a L × L square with L = 50m. Here as initial condition we consider a at topography, null velocities in both the x and y directions and water height h(x, y, t = 0) = 10 if r ≤ 11m 1 if r > 11m, with r = (x -25) 2 + (y -25) 2 . Thus, we are considering a cylindrical dam that instantaneously breaks at the beginning time t = 0s. The ending time t End = 0.69s. Satisfying results of the rst order scheme are reported in gure 6.

Water drop in a basin

For this numerical test, see [START_REF] Audusse | A 2d Well-balanced Positivity Preserving Second Order Scheme for Shallow Water Flows on Unstructured Meshes[END_REF]. Here we simulate a water drop in a basin and consequently re ective boundary conditions are used. A L-side square domain with L = 20m is considered. The topography is still taken at and constant in time. At initial time we assume u = (0, 0) t and h(x, y, t = 0) = 2.4(1 + e -0.25((x-10.05) 2 +(y-10.05) 2 ) )

The outputs are shown in gure 7 at time t = 1s, t = 2s, t = 3s and t = 4s respectively. The results agree with the ones reported in [START_REF] Audusse | A 2d Well-balanced Positivity Preserving Second Order Scheme for Shallow Water Flows on Unstructured Meshes[END_REF]. In particular, in the same picture 7, we compare the results at time t = 4s obtained using the 2D extensions of the rst-order and second-order schemes. We can clearly see that the latter scheme gives less di usive solutions under the same mesh M = 100 cells.

2D squared dam break

Once again we assume at topography constant in time. The domain is a square of side L = 200m. The dam position is represented in gure 8. We note a breach of length 75m which is instantaneously opened at time t = 0. At initial time we also have null velocities and h(x, y, t = 0) = 10 if x ≤ 100m 5 if x > 100m.

For more details refer to paper [START_REF] Seaïd | Non-oscillatory relaxation methods for the shallow-water equations in one and two space dimensions[END_REF]. The results for the water height and the velocity eld are shown in gure 9 and they appear to be in agreement with the reference outputs given in [START_REF] Seaïd | Non-oscillatory relaxation methods for the shallow-water equations in one and two space dimensions[END_REF].

2D ow over a smooth bump

The following test problem is useful to check the well-balanced property of the scheme, see [START_REF] Seaïd | Non-oscillatory relaxation methods for the shallow-water equations in one and two space dimensions[END_REF]. The domain is a square of side L = 1m and we consider the Grass formulation for the Exner equation with A g = 1 and ζ = 1. At initial time we Thus, the initial solution satis es the "lake at rest" condition. Indeed, our 2D numerical schemes are able to preserve this kind of stationary solutions up to an error of order 10 -15 .

Conical dune of sand

This test case has been vastly used to validate numerical schemes for shallow water Exner system, here we do refer for instance to [START_REF] Audusse | Parallelization of a relaxation scheme modelling the bedload transport of sediments in shallow water ow[END_REF][START_REF] Hudson | Numerical techniques for morphodynamic modelling[END_REF]. When considering the Grass formulation for the sediment discharge, we take porosity ρ 0 = 0.4, where we recall that ζ = 1 1-ρ0 . The domain is a L × L square with L = 1000m. At time t = 0, we impose z(x, y, t = 0) = 0.1 + (sin( π(x-300)

200

)) 2 (sin( π(y-400)

200

)) 2 if 300 ≤ x ≤ 500 and 400 ≤ y ≤ 600 0.1 otherwise, h(x, y, t = 0) = 10 -z(x, y, t = 0), u(x, y, t = 0) = 10 h(x, y, 0) and v(x, y, t = 0) = 0.

Here, as boundary conditions, we assume that at the upstream we impose u(x, y, t) = 10 h(x,y,t) while the other boundaries are usual transmissive conditions.

Then, we consider two di erent cases: in the rst one we take A g = 1 and ending time t End = 500s, thus we are assuming a fast interaction between the ow and the sediments. As second case, we diminish the value of A g , namely we impose A g = 0.1, thus the strength of the interaction decreases. The outputs for these two test cases can be found in pictures 10 and 11 respectively. The result are in agreement with the ones reported in [START_REF] Audusse | Parallelization of a relaxation scheme modelling the bedload transport of sediments in shallow water ow[END_REF] even if more di usive due to a coarser mesh size (M = 100 cells).

Concluding remarks

In this work a second-order well-balanced Lagrange-projection scheme for the shallow water Exner system has been presented in 1D. Two-dimensional extension is proposed as well. Numerical results proved the validity of the scheme. Generally, no oscillations are present in the numerical outputs, with the exception of one test case, in which it is shown that to slightly re ne the mesh solve the mentioned problem. Further work to solve implicitly the acoustic step is required. Other strategies to take into account the Exner equation in the Lagrange-projection formalism are being explored and will be shortly submitted.

Figure 1 :

 1 Figure 1: RP: dam break on movable bottom; free surface (up) and bed level (bottom). Reference solution (red) and rst-order weakly coupled solution (blue). M = 200 cells (up) and M = 2000 cells (bottom).

  (a) 2nd order weakly coupled scheme (b) 2nd order decoupled scheme

Figure 2 :

 2 Figure 2: RP: dam break on movable bottom; free surface (up) and bed level (bottom). "Decoupled" (bottom) and "weakly coupled" (up) solutions. Mesh of size M = 100 (blue line), M = 200 (magenta line) and M = 500 (green line) cells. Red line for reference solution.

Figure 3 :

 3 Figure 3: Flow over a movable bump; bed level z and free surface z + h. First (blue) and second-order (magenta) weakly coupled solutions with M = 200 cells. Reference solution in red.

Figure 4 : 20 Figure 5 :

 4205 Figure 4: Free surface and bed level, M = 200 cells. Reference solution (red) and rst-order weakly coupled solution (blue).

Figure 6 :

 6 Figure 6: 2D circular dam break on wet bed; water height (left) and its contour plot (right). 2D extension of rst-order scheme solution. M = 100 cells, t End = 0.69s and CFL= 0.45.

  (a) t = 0s, 2D extension of 2nd-order scheme (b) t = 1s, 2D extension of 2nd-order scheme (c) t = 2s, 2D extension of 2nd-order scheme (d) t = 3s, 2D extension of 2nd-order scheme (e) t = 4s, 2D extension of 2nd-order scheme (f) t = 4s, rst order scheme

Figure 7 :

 7 Figure 7: Water drop in a basin; water height at time t = 0s, t = 1s, t = 2s, t = 3s and t = 4s. M = 100 cells, CFL= 0.25 and CFL= 0.45 for the 2D extension of second and rst order scheme respectively.

Figure 8 :

 8 Figure 8: Geometry of the dam in the squared dam break test problem.

  (a) t = 3.6s (b) t = 3.6s (c) t = 7.2s (d) t = 7.2s (e) t = 10.8s (f) t = 10.8s

Figure 9 :

 9 Figure 9: Squared dam break; water height (left) and contour plot and velocity eld (right) at time t = 3.6s, t = 7.2s and t = 10.8s. Dam in black. 2D extension of second order scheme with M = 100 cells and CFL= 0.25.

Figure 10 :

 10 Figure 10: Conical dune of sand; fast interaction Ag = 1. Bed level (left) and contour plot (right) at time t = 500s. 2D extension of second order scheme with M = 100 cells and CFL= 0.25.

Figure 11 :

 11 Figure 11: Conical dune of sand; fast interaction Ag = 1. Bed level (left) and contour plot (right) at time t = 500s. 2D extension of second order scheme with M = 100 cells and CFL= 0.25.

Table 1 :

 1 Errors and empirical convergence rates for norm L 1 . Mesh of size M = (64, 128, 256, 512, 1024, 2048), CF L = 0.25. Second-order decoupled and weakly coupled numerical schemes.

	Method	Mesh M	Variable	err L 1	O(L 1 )	Variable	err L 1	O(L 1 )	Variable	err L 1	O(L 1 )
	Decoupled	64	h	0.0268	-	hu	0.1175	-	z	0.1792 × 10 -3	-
		128		0.0083 1.6953		0.0354 1.7320		0.0544 × 10 -3	1.7190
		256		0.0027 1.6185		0.0115 1.6249		0.0182 × 10 -3	1.5785
		512		0.0007 1.8756		0.0031 1.8782		0.0050 × 10 -3	1.8761
		1024		0.0002 1.9781		0.0008 1.9818		0.0012 × 10 -3	1.9932
	Weakly coupled	64	h	0.0268	-	hu	0.1175	-	z	0.1824 × 10 -3	-
		128		0.0083 1.6955		0.0354 1.7320		0.0550 × 10 -3	1.7296
		256		0.0027 1.6182		0.0115 1.6248		0.0183 × 10 -3	1.5890
		512		0.0007 1.8755		0.0031 1.8781		0.0050 × 10 -3	1.8734
		1024		0.0002 1.9782		0.0008 1.9818		0.0013 × 10 -3	1.9792

Acknowledgments

This work was partially supported by a grant from Région Île-de-France.

Projection step

In this step, we exploit again a reconstructed polynomial LP. However, in order to preserve the second-order of accuracy, it is crucial to reconstruct the Lagrangian variables (LX), namely as

with the slope

) and where the variable X denotes h, hu and u. Then, the updating formula for X = h and X = hu are given by a second-order approximation of the three integrals that appear in (3.7). This is achieved by using the mid-point rule, thus we get

where we use the upwind de nition

As far as the topography is concerned, we consider the weakly coupled scheme (3.11) where we naturally set

and u * j± 1 2 is given by (3.15). Afterwards, for the decoupled scheme, we rst de ne the reconstructed polynomial for the water height h and the ow hu at time t n , P (X) n j (x) = X n j + ∆ n j (x -x j ) with X = h, hu and ∆ n j the slopes (either ENO or Minmod). Then, we use formula (3.10) where we impose

Theorem 2. The second-order numerical scheme described preserves the "lake at rest" stationary solution (1.6).

Proof. Since we already proved the well-balanced property for the rst-order scheme, it is straightforward to show it for the second-order method as well. Indeed, it is enough to observe that the slopes ∆ n j = ∆ n j (D n j-1,j , D n j,j , D n j+1,j ) are null under the hypothesis of the "lake at rest condition" thanks to de nition of the uctuations. Hence, once again we obtain u *

Finally, it is only worth to specify that the Runge-Kutta TVD procedure automatically preserves the stationary solutions.

Two-dimensional extension

In this section, we brie y describe how we extend the proposed approach in two space dimensions using dimensional splitting. Let us rst recall that if we denote (x, y) ∈ R 2 the space variables and u = (u, v) T the velocity vector, then the


Paolo Censi
email: paolo.censi@unipa.it393479662844

Ahmet Sasmaz

Pierpaolo Zuddas

Marianna Cangemi

Daniela Piazzese

Gulsah Ozek

Marco Venturi

Thomas H Darrah
email: darrah.24@osu.edu

Yigal Erel
email: yerel@vms.huji.ac.il

rare earth elements in neutralalkaline waters: Case study of Lake Van hydrothermal system, Turkey

Keywords: Zr/Hf ratio, REE, Alkaline lakes, Turkey Zr/Hf ratio, REE, Alkaline lakes, Turkey

We investigated the distribution of Zr, Hf, and rare earth elements (REE) as the sum of lanthanides plus Y) in the hydrothermal system in the Lake Van area of south-eastern Turkey. This system is characterised by water with variable pH in alkaline conditions resulting from hydrothermal CO 2 upraise and neoformation of calcite minerals in near equilibrium with the interacting waters. Zr, Hf, and REE determinations were carried out for aqueous phases and suspended solids in lake water and surrounding thermal springs. We found that dissolved Hf is partitioned relative to Zr during calcite formation and that such fractionation is a function of the Ca 2+ activity in warm water. The observed Zr-Hf fractionation is explained by coulombic interactions that occur between suspended solid particles and dissolved phases at the calcite-water interface. There, the surfaces of carbonate minerals demonstrated greater reactivity towards aqueous Hf-bearing species relative to Zr-complexes. This evidence involves a coulombic mechanism of reactivity at the calcite-water interface because Hf complexes are negatively charged while Zr compounds are uncharged. Thus, authigenic calcite can behave as a suitable host for dissolved metal ion species to adsorb on crystal surfaces to remediate waste waters from mine drainage.

Introduction

Zirconium, Hf, and REEs are used in specialised technology industries notwithstanding their limited production [START_REF] Chakhmouradian | Rare earth elements-Minerals, mines, magnets (and more[END_REF][START_REF] Jones | Critical mineral resources of the United States-Economic and environmental geology and prospects for future supply[END_REF]. Demand for these elements is high leading an effort by the mineral exploration industry to understand genetic processes and discover new resources. Geochemically, these elements were considered to be 'immobile' in subsurface conditions while some works indicate easy water transport under subsurface hydrothermal conditions. Accordingly, REE enrichment has been observed in deposits formed from carbonisation in alkaline environments where Zr concentrations were enriched by a factor of 3 to 4 [START_REF] Estrade | The Origin of Skarn-Hosted Rare-Metal Mineralization in the Ambohimirahavavy Alkaline Complex, Madagascar[END_REF]. The behaviour of these strategic elements under surface and subsurface conditions is often neglected because of difficulties in accurate analytical determination at the level of pmol/kg in natural water. In a recent survey, the dissolved Zr/Hf ratio was different from that found in the rock minerals [START_REF] Censi | Trace element fractionation through halite crystallisation. Geochemical mechanisms and environmental implications[END_REF][START_REF] Inguaggiato | Geochemistry of REE, Zr and Hf in a wide range of pH and water composition: the Nevado del Ruiz volcano-hydrothermal system (Colombia)[END_REF][START_REF] Inguaggiato | Zirconium-hafnium and rare earth element signatures discriminating the effect of atmospheric fallout from hydrothermal input in volcanic lake water[END_REF][START_REF] Parisi | Rare earths, zirconium and hafnium distribution in coastal areas: The example of Sabella spallanzanii (Gmelin, 1791)[END_REF][START_REF] Tepe | Distribution of rare earth elements and other high field strength elements in glacial meltwaters and sediments from the western Greenland Ice Sheet: evidence for different sources of particles and nanoparticles[END_REF]Zuddas et al., 2017[START_REF] Censi | Zr, Hf and REE distribution in river water under different ionic strength conditions[END_REF] suggesting that the subsurface cycle is important in the accumulation of these elements. In fact, the Zr/Hf molar value in rock minerals was 71.4±5.6 which corresponds to the so-called 'chondritic' signature [START_REF] Jochum | The solar-system abundances of Nb, Ta, and Y, and the relative abundances of refractory lithophile elements in differentiated planetary bodies[END_REF], while natural water had a wider range of Zr/Hf [START_REF] Censi | Zr, Hf and REE distribution in river water under different ionic strength conditions[END_REF]. The supergene geochemical cycle of Zr and Hf participates in the accumulation and stability of these elements through the process of dissolution and precipitation of the minerals. Because the composition of subsurface and surface waters is often near equilibrium with carbonate minerals like calcite, calcite is expected to play a significant role in the supergene behaviour of these elements. Although carbonates are among the most common minerals in equilibrium with natural water, the formation of neogenic calcite also determined the change in REE distribution observed by Moller and De Lucia (2020) in controlled laboratory investigations. To explain Zr-Hf decoupling during calcite crystallisation and the concurrent behaviour of REEs, we investigated the distribution of these elements in a natural hydrothermal system characterised by the formation of calcite under the neutral and alkaline water system found at Lake Van and characterised by low thermal and lake waters with pH ranging from 6 to 10, resulting from a variable CO2 flux [START_REF] Degens | The Geology of Lake Van[END_REF][START_REF] Moller | Rare-earth patterns with positive cerium anomaly in alkaline waters from Lake Van, Turkey[END_REF]. The Lake Van system is characterised by neogenic precipitation in the form of Ca-Mg carbonates [START_REF] Moller | Rare-earth patterns with positive cerium anomaly in alkaline waters from Lake Van, Turkey[END_REF].

Material and Methods

Field sampling

Lake Van is located in south-eastern Turkey (Fig. 1) which formed during the collision between European and Arabian plates. Field study and water sampling were conducted in November 2018 and consisted of a 16 samples collected from lake water and thermal springs. Coordinates and geologic units of the collected samples are given in Table 1. Lake water was sampled at a depth of 15-20 cm. Physicochemical parameters (pH, redox potential (Eh), temperature, and electrical conductivity) were measured with an Eh oxytrode Pt probe (Hamilton TM ) using a standard reference solution buffer at 475±5 mV. Accuracy of determinations was ±10 mV. For the major ion determinations, water samples were filtered in the field through 0.45 μm Millipore® cellulose acetate filters into 50 mL polyethylene bottles and acidified with ultrapure HNO3. Water samples for trace and rare element analyses were collected in 2,000 mL Nalgene bottles at each sampling site.

Water analysis and thermodynamic modelling

In the laboratory, samples were filtered (Millipore® manifold filter, 47 mm diameter with a 0.45 µm pore size cellulose nitrate membrane), then 5% ultrapure HNO3 acid solution was added to reach pH ≈ 2. Zirconium, Hf, and REEs in aqueous samples were analysed according to the methods of [START_REF] Raso | Simultaneous determinations of zirconium, hafnium, yttrium and lanthanides in seawater according to a co-precipitation technique onto iron-hydroxide[END_REF]. Briefly, an excess of FeCl3 (1%) solution was added to each sample (1 L), and a suitable volume of NH4OH (25%) solution was added to attain a pH of 8 to induce precipitation of solid Fe(OH)3. During this process, Zr, Hf, and REEs were held onto the surface of the crystallising solid. To ensure that crystallisation of Fe(OH)3 was complete, the solution was left in a closed flask for 48-72 h in a stirrer. Next, the iron concentration was measured to assess recovery and was consistently greater than 95%. Precipitated Fe(OH)3, together with Zr, Hf, and REEs, was collected on a membrane filter (Millipore® manifold filter, 47 mm diameter with 0.45 µm pore size). The solid filtrate was then dissolved in 3 M HCl, and the obtained solution was diluted to 1 M HCl prior to analysis using a Quadrupole-ICP-MS (Agilent 7800 series) with an external calibration procedure.

The overall strategy and a detailed evaluation of the analytical errors were reported in [START_REF] Raso | Simultaneous determinations of zirconium, hafnium, yttrium and lanthanides in seawater according to a co-precipitation technique onto iron-hydroxide[END_REF]. The overall enrichment factor measured for Zr, Hf, and REE concentrations corresponded to 66.7.

Lanthanum, Ce, and Eu anomalies were calculated using eqations 1-3 [START_REF] Alibo | Rare earth elements in seawater: Particle association, shalenormalization, and Ce oxidation[END_REF]:

𝐿𝑎 𝐿𝑎 * = [𝐿𝑎] 𝑛 3 * [𝑃𝑟] 𝑛 +2 * [𝑁𝑑] 𝑛 , ( 1
)
𝐶𝑒 𝐶𝑒 * = 2 * [𝐶𝑒] 𝑛 [𝐿𝑎] 𝑛 +[𝑃𝑟] 𝑛

, and

(2)

𝐸𝑢 𝐸𝑢 * = 2 * [𝐸𝑢] 𝑛 [𝑆𝑚] 𝑛 +[𝐺𝑑] 𝑛 , (3)
where the subscript n indicates shale-normalised elemental concentrations calculated relative to the Post Archean Australian Shale (PAAS, [START_REF] Taylor | The geochemical evolution of the continental crust[END_REF].

All chemicals used during laboratory manipulations were ultrapure grade. Ultrapure water (resistivity of 18.2 MΩ cm or better) was obtained from an Arium® mini system (Sartorius, Italy).

Nitric acid 65% (w/w), ammonia solution, and HCl acid were purchased from J.T. Baker chemicals.

Working standard solutions of the elements studied were prepared daily by stepwise dilution of standard multi-element stock solutions from DBH, Merck, or CPI International (1,000±5 mg/L) in a 1 M HCl medium. All labware was polyethene, polypropylene, or Teflon® and the calibration of all volumetric equipment was verified. Assessment of the analytical precision of Zr, Hf, and REEs was performed on five aliquots (500 mL each) of NASS-6 (seawater certified reference material for trace metals) distributed by the National Research Council of Canada. These certified reference materials were treated as water samples according to the procedure mentioned above and the concentrations obtained were compared with those previously reported by [START_REF] Raso | Simultaneous determinations of zirconium, hafnium, yttrium and lanthanides in seawater according to a co-precipitation technique onto iron-hydroxide[END_REF] and [START_REF] Lemaitre | Trace element behaviour at cold seeps and the potential export of dissolved iron to the ocean[END_REF]. Analysed values are reported in Supplementary Data 1.

Thermodynamic modelling of the chemical composition of the water performed was using PHREEQC software and adding additional Lawrence Livermore National Laboratory (LLNL) database (version 3.0.6, Parkhurst and Appelo, 2010).

Suspended solids

X-ray diffraction investigations were carried out with a Philips PW14 1373 X-ray spectrometer (Cu-Kα radiation 2Θ range 10-40°, step size 0.02°, and a 1-min step time). The X-ray spectra were analysed by the Rietveld method (Program DiffracPlus TOPAS®, Version 4.0, Bruker AXS Inc., Karlsruhe, Germany) using parameters for the Rietveld refinement method obtained from the Inorganic Crystal Structure Database (ICSD) program. This method consists of fitting the experimental XRD spectrum to the theoretical spectrum (Young, 1993). Scanning electron microscopy (SEM) of suspended particles was carried out on materials collected after filtration onto cellulose nitrate membranes. Dried solids were mounted on an aluminium stub and gold coated. The SEM analyses were carried out using a LEO 440 SEM equipped with an EDS system OXFORD ISIS Link and Si (Li) PENTAFET detector at the SIDERCEM SRL laboratory (Caltanissetta, Italy).

Chemical analyses were carried out by digesting 100 mg of each sample in 10 mL of a 1:1 HNO3:H2O2 mixture in a sealed Teflon® TFM bomb using a microwave oven (CEM MARS 5 device). After digestion, the TFM bomb was rinsed, and the solution was filtered onto previously acid-cleaned 0.45 μm Millipore® filters to remove silicate residue. The filtered solution was diluted to 50 mL and stored for chemical analyses.

Results and Discussion

System description

The concentration of major elements, pH, Eh, electric conductivity, and temperature are reported in the Supplementary Data 2. Water samples had ionic strength ranging from 0.024 to 0.222 M in thermal water and between 0.004 and 0.443 M in lake water. Thermal water had pH between 5.8 and 7.1 and Eh between 0.013 and 0.172 V, while lake water had pH between 8 and 10 and Eh from -0.288 and 0.012 V. Figure 2a shows that lake water preferentially falls in the field of alkaline elements (only two samples have higher Ca 2+ and Mg 2+ concentrations) while thermal water fell between alkaline and Ca 2+ -rich field (few samples fell in the Mg 2+ corner).

The SEM images of the suspended particulates revealed that suspended particulate matter in lake water consisted of authigenic carbonates associated with biogenic detritus and lithic fragments (Fig. 3a-f), while suspended solids in thermal waters were authigenic carbonates and Fe-rich particles and lithogenic fragments (Fig. 3g-i). Iron-bearing particulates in the samples were often mixed with authigenic calcite and dolomite (see elemental maps in Fig. 3i). The X-ray diffraction analyses confirmed the presence of carbonate minerals identified as calcite and low-Mg calcite representing up to 80% of the authigenic minerals (Table 2). Other phases collected were Feoxyhydroxides and detritic silicate particles reaching 10-15% of the total mass.

Using the chemical composition of the collected samples, we estimated the saturation state with respect to the mineral in equilibrium and found both lake water and thermal water were saturated or oversaturated with respect to calcite or low-Mg calcite. Some water was also saturated with respect to dolomite and goethite in agreement with the results of the X-ray determinations of the collected suspended matter.

Zr and Hf behaviour

Zr and Hf concentrations in both aqueous and suspended particulate matter are reported in Supplementary Data 3 and 4, respectively. In the aqueous phase, the concentration of Zr and Hf was not homogenous. Lake waters had Zr concentrations between 79 and 4,958 pmol/L while Hf concentrations fell between 2 and 53 pmol/L, and thermal waters had Zr concentrations between 32 and 1.425 pmol/L; whereas Hf concentrations fell between 0.5 and 14 pmol/L. The concentrations of Zr and Hf in thermal waters were 1 to 3 orders of magnitude lower than those found in lake water;

however, the concentration of Zr as a function of Hf was linearly correlated with an intercept equal to zero and parametrising a single slope of 90.4±1.7 (Fig. 4a). This value was higher than the Zr/Hf value of 71±5.6 for crustal materials and chondrites [START_REF] Jochum | The solar-system abundances of Nb, Ta, and Y, and the relative abundances of refractory lithophile elements in differentiated planetary bodies[END_REF] indicating that the calcite precipitation fractionates the Zr/Hf ratio in water. In the suspended solids from lake water, Zr concentration ranged from 99.5 to 2,533.2 nmol/kg, whereas Hf fell between 2.6 and 77.1 nmol/kg.

Suspended solids from thermal water had Zr concentrations between 20.7 and 346.7 nmol/kg and Hf between 0.6 and 6.8 nmol/kg. Similarly, Zr and Hf in solids collected from thermal water were lower compared to suspended solids found in lake water. Figure 4B shows the concentration of Zr as a function of Hf in neogenic calcite from lake water and describes a linear slope of 33.5 similar to calcite in equilibrium with thermal water where the slope was 37.

These results show that the Zr/Hf ratio in water is different from that in crustal rocks and that

Hf is preferentially enriched in authigenic carbonate minerals. Preferential enrichment of Hf relative to Zr in newly-formed suspended carbonates was unexpected and may result from processes occurring at the surface-water interface. The calcite-water interface is generally represented by an electric double-layer of weak mineral solubility buffered by CO3 2-ions with the  potential of the {104} carbonate surface independent of pH [START_REF] Herberling | Structure and reactivity of the calcite-water interface[END_REF]. Laboratory measurements indicate the positive surface potential (up to 24 mV) when the Ca 2+ concentration in solution is between 10 -2 and 10 -4 M (Al-Mahrouqui et al., 2017). [START_REF] Byrne | Inorganic speciation of dissolved elements in seawater: The influence of pH on concentration ratios[END_REF] made a critical review of stability constants for Zr and Hf hydroxyl complexes based on the work of [START_REF] Sillen | Stability Constants of Metal-Ion Complexes[END_REF], [START_REF] Baes | The Hydrolysis of Cations[END_REF], [START_REF] Ekberg | Studies on the hydrolytic behavior of zirconium (IV)[END_REF], and [START_REF] Brown | Chemical Thermodynamics of Zirconium[END_REF] and estimated that the most abundant Zr and Hf complexes in aqueous phase were [Zr(OH)4] 0 and [Hf(OH)5] -, respectively [START_REF] Koschinsky | Uptake of elements from seawater by ferromanganese crusts: Solid-phase associations and seawater speciation[END_REF]. Reporting Zr/Hf ratios as a function of pCa (pCa=-log[Ca 2+ ]) (Fig. 5), we found two distinct correlation lines, one for thermal water and one for lake water, suggesting that negatively-charged [Hf(OH)5] -complexes were preferentially scavenged onto {104} positive carbonate surfaces and were, in turn, removed from the aqueous phase. This potential mechanism explains the Zr/Hf fractionation observed during calcite precipitation under pH conditions between 7

and 10 and should correspond to the higher binding energy between Hf and {104} calcite surface relative to Zr as estimated by theoretical ab initio calculations [START_REF] Pinto | Adsorption and migration of single metal atoms on the calcite (10.4) surface[END_REF]. Our results support the Zr/Hf fractionation observed by [START_REF] Firdaus | Strong elemental fractionation of Zr-Hf and Nb-Ta across the Pacific Ocean[END_REF][START_REF] Firdaus | Distribution of zirconium, hafnium, niobium and tantalum in the North Atlantic Ocean, north-eastern Indian Ocean and its adjacent seas[END_REF] during Fe-oxyhydroxide formation. Fe-oxyhydroxide was in fact positively charged under weak alkaline conditions [START_REF] Koschinsky | Uptake of elements from seawater by ferromanganese crusts: Solid-phase associations and seawater speciation[END_REF]. On the contrary, Zr-Hf fractionation did not occur onto surfaces of Mn-oxyhydroxides that were negatively charged at the pH up to acidic values [START_REF] Tripathy | The electrical double layer at hydrous manganese dioxide/electrolyte interface[END_REF].

The Zr-Hf fractionation here observed during calcite formation and interpreted as water-interface interaction is consistent with the Zr-Hf fractionation found in stromatolites [START_REF] Censi | The behaviour of rare-earth elements, Zr and Hf during biologically-mediated deposition of silica-stromatolites and carbonate-rich microbial mats[END_REF], halite (Censi et al., 2017a[START_REF] Censi | Trace element fractionation through halite crystallisation. Geochemical mechanisms and environmental implications[END_REF], and other groundwater environments (Censi et al., 2017b;[START_REF] Inguaggiato | Geochemistry of REE, Zr and Hf in a wide range of pH and water composition: the Nevado del Ruiz volcano-hydrothermal system (Colombia)[END_REF][START_REF] Inguaggiato | Zirconium-hafnium and rare earth element signatures discriminating the effect of atmospheric fallout from hydrothermal input in volcanic lake water[END_REF].

REE behaviour

Dissolved REE concentrations in lake water and thermal water in the Lake Van system ranged between 7,135 and 76,140 pmol/L and from 168 to 17,426 pmol/L, respectively (Supplementary Data 3). In lake waters the shale normalised pattern (Fig. 6A) reveals the presence of La positive anomaly as well an 'ascending' feature leading to enrichment in heavy REE (HREE, from Ho to Lu) while thermal waters have a single slight medium REE (MREE, from Sm to Dy) enrichment (Fig. 6B). In both types of water, Eu anomalies were 0.3 and 4.8 without evidence of a Ce anomaly (0.1 ≤ Ce/Ce * ≤ 1). The progressive increase of REEs observed in the normalised patters in these waters may result from the values of the stability constants of di-carbonate and carbonate REE-complex ions [START_REF] Lee | Complexation of trivalent rare earth elements (Ce, Eu, Gd, Tb, Yb) by carbonate ions[END_REF][START_REF] Liu | Comprehensive investigation of yttrium and rare earth element complexation by carbonate ions using ICP-mass spectrometry[END_REF] potentially producing [REE(CO3)2] - complexes (Johannesson et al., 1994;Johannesson and Lyon, 1994;[START_REF] Kerrich | Trace-element composition of cherts from alkaline lakes in the east African rift: a probe for ancient counterparts[END_REF][START_REF] Moller | Rare-earth patterns with positive cerium anomaly in alkaline waters from Lake Van, Turkey[END_REF]. Our data showed a higher abundance of La and therefore inconsistent with earlier REE measurements from Van Lake waters [START_REF] Moller | Rare-earth patterns with positive cerium anomaly in alkaline waters from Lake Van, Turkey[END_REF]. The water level of Lake Van has dropped due to climate change over the last 30 years while the input of lithogenic solids increased.

The main difference between the REE distribution recognised in lake Van waters and that reported by [START_REF] Moller | Rare-earth patterns with positive cerium anomaly in alkaline waters from Lake Van, Turkey[END_REF] is the lack of positive Ce anomaly in our samples. This may be caused by a larger dissolution of lithogenic particles typical of continental and coastal waters. The dissolution of lithogenic solids usually provides La enrichment in the aqueous phase that can conceal the concurrent Ce enrichment completely or partially [START_REF] Censi | The behaviour of REEs in Thailand's Mae Klong estuary: Suggestions from the Y/Ho ratios and lanthanide tetrad effects[END_REF][START_REF] Greaves | Aeolian sources of rare earth elements to the Western Pacific Ocean[END_REF][START_REF] Grenier | Differentiating lithogenic supplies, water mass transport, and biological processes On and Off the Kerguelen Plateau using rare earth element concentrations and neodymium isotopic compositions[END_REF].

Suspended particulate material collected from the lake and thermal waters had a total concentration of REEs ranging between 508.2 and 835.2 nmol/kg and from 102.3 to 278.5 nmol/kg, respectively (Supplementary Data 4). These values are 3 to 4 orders of magnitude higher than the concentrations found in related waters. The shale-normalised REE patterns of suspended calcites collected in lake waters (Fig. 6c) and thermal waters (Fig. 6d) display positive Ce anomalies ranging from 1.1 to 3.1, an enrichment in MREEs with positive Eu anomalies (0.8 < Eu/Eu* < 4.2), and slight progressive depletion of HREEs. Suspended solids from thermal waters showed 'bulged' REE patterns (Fig. 6D) while solids from lake waters had flatter patterns (Fig. 6C). The REE distribution in shale-normalised patterns usually show MREE enrichment mainly ruled by crystal-chemical reasons, as demonstrated by the review of distribution coefficient values reported by Moller and De Lucia (2020). In suspended solids from thermal water, the MREE enrichment was more evident than in materials collected from lake water. This evidence is probably justified by the concurrent presence of a detrital fraction in the 'bulge' of suspended solids from lake water with a 'flatter' shale-like trend in the normalised concentration [START_REF] Taylor | The significance of the rare earths in geochemistry and cosmochemistry[END_REF].

The behaviour of Y-Ho twin is different in lake and thermal waters. In lakes, Y and Ho concentrations ranged from 3,845 to 62,279 pmol/L and 36 to 686 pmol/L respectively, whereas in thermal springs, Y ranges up to 3,577 pmol/L, whereas Ho fell between 2 and 68 pmol/L. The related Y/Ho ratios spanned from 51 to 111 in lake water and between 53 and 162 in thermal water.

These ratios fell between chondritic and super-chondritic values. On the contrary, in carbonates from suspended solids, the Y/Ho ratio ranged from 32.4 to 49.6 in lake water and from 28.6 and 47.8 in thermal water. Then Y/Ho values in suspended solids spanned between slight subchondritic to chondritic values (Y/Ho = 52±5.1, [START_REF] Jochum | The solar-system abundances of Nb, Ta, and Y, and the relative abundances of refractory lithophile elements in differentiated planetary bodies[END_REF].

The role played by the formation of calcite in the distribution of REEs in alkaline waters was highlighted by the super-chondritic value of the Y/Ho ratio. Our study revealed that Ho was 2 to 3 orders of magnitude higher than Y during calcite formation causing enrichment. The preferential enrichment of Ho in calcite compared to Y observed in our study may reflect the larger value of the Misono softness parameter of Ho 3+ relative to Y 3+ [START_REF] Thompson | Fractionation of yttrium and holmium during basaltic soil weathering[END_REF] due to the different external electronic configuration of these metal ions [START_REF] Qu | Experimental study on the fractionation of yttrium from holmium during the coprecipitation with calcium carbonates in seawater solutions[END_REF].

The different fractionation of these elements found during the formation of calcite also reflects other processes. In thermal water, the intensity of the water-rock interaction may cause MREE enrichment (Fig. 6c). Dissolution of lithogenic rocks may be the source of MREE [START_REF] Bau | Scavenging of dissolved yttrium and rare earths by precipitating iron oxyhydroxide: Experimental evidence for Ce oxidation, Y-Ho fractionation, and lanthanide tetrad effect[END_REF][START_REF] Censi | Zr/Hf ratio and REE behaviour: A coupled indication of lithogenic input in marginal basins and deep-sea brines[END_REF][START_REF] Greaves | Aeolian sources of rare earth elements to the Western Pacific Ocean[END_REF][START_REF] Haley | Rare earth elements in pore waters of marine sediments[END_REF][START_REF] Inguaggiato | Zirconium-hafnium and rare earth element signatures discriminating the effect of atmospheric fallout from hydrothermal input in volcanic lake water[END_REF]Lin et al. 2019) and explain MREE enrichment observed in authigenic hydrothermal carbonates [START_REF] Jakubowicz | Tracing the composition and origin of fluids at an ancient hydrocarbon seep (Hollard Mound, Middle Devonian, Morocco): A Nd, REE and stable isotope study[END_REF][START_REF] Phan | Rare earth elements and radiogenic strontium isotopes in carbonate minerals reveal diagenetic influence in shales and limestones in the Appalachian Basin[END_REF]Wang et al., 2020), phosphates [START_REF] Censi | The behaviour of REEs in Thailand's Mae Klong estuary: Suggestions from the Y/Ho ratios and lanthanide tetrad effects[END_REF][START_REF] Hannigan | The development of middle rare earth element enrichments in freshwaters: Weathering of phosphate minerals[END_REF][START_REF] Reynard | Crystal-chemical controls on rare-earth element concentrations in fossil biogenic apatites and implications for paleoenvironmental reconstructions[END_REF]Zhang et al., 2016), sulphates [START_REF] Censi | Authigenic phase formation and microbial activity control Zr, Hf, and rare earth element distributions in deep-sea brine sediments[END_REF]2018;[START_REF] Kagi | Tetrad effect in lanthanide partitioning between calcium sulphate crystal and its saturated solution[END_REF][START_REF] Playà | Non-marine evaporites with both inherited marine and continental signatures: The Gulf of Carpentaria, Australia, at similar to 70 ka[END_REF][START_REF] Toulkeridis | Tracing the source of gypsum in New Caledonian soils by REE contents and S-Sr isotopic compositions[END_REF], and Fe-oxyhydroxides [START_REF] Bau | Scavenging of dissolved yttrium and rare earths by precipitating iron oxyhydroxide: Experimental evidence for Ce oxidation, Y-Ho fractionation, and lanthanide tetrad effect[END_REF][START_REF] Davranche | Rare earth element patterns: a tool for identifying trace metal sources during wetland soil reduction[END_REF].

Implications and Conclusions

The Based on previous researches addressed to the larger Hf affinity relative to Zr towards crystal surfaces of halite, potash salts and Fe-oxyhydroxides, this study contributes to demonstrate that the larger surface reactivity of Hf relative to Zr is a widespread phenomenon through the sedimentary processes. The evidence of this phenomenon is detected measuring the Zr/Hf ratio in natural waters and in authigenic minerals.

  results of our investigation of water oversaturated with calcite indicate that calcite plays a crucial role in constraining Zr-Hf fractionation and the geochemical behaviour of REEs in alkaline lake waters. The strong reactivity of positively charged {104} surfaces of calcite coupled with the aqueous Zr and Hf speciation lead to preferential accumulation of Hf in calcite relative to Zr. This indication agrees with the larger affinity of Hf towards solid surfaces relative to Zr already evidenced onto the surfaces of Fe-oxyhydroxides. But the Hf removal onto the carbonate surfaces is ruled by pCa, whereas onto the Fe-oxyhydroxide surfaces is ruled by pH. Therefore, changes of the Zr/Hf ratio in natural waters relative to the chondritic signature indicate that these elements are involved in processes occurring at the interface with crystallising authigenic minerals. Differences in the external electronic configuration of Zr and Hf is similar to that between Y and Ho due to differences in softness/hardness between Zr and Hf and between Y and Ho that, in turn, influence hydrolysis during formation of Zr and Hf hydroxyl complexes. This process could be a crucial aspect of the fractionation of these metal ions during calcite formation. Differences in softness between Y and Ho lead increase stability of Ho-O relative to Y-O bonds on solid surfaces and rule Y-Ho fractionation at the solid-liquid interface.

FIGURE CAPTIONS Figure 1 -

 CAPTIONS1 FIGURE CAPTIONS

Figure 2 -

 2 Figure 2 -Major element concentrations analysed in studied waters.

Figure 3 -

 3 Figure 3 -SEM images of suspended solids in studied water samples filtered with membrane

Figure 4

 4 Figure 4  a: Zr and Hf concentration measured in studied water samples and the linear trend

Figure 5 -

 5 Figure 5-Zr/Hf values measured in the dissolved phase of the lake and thermal waters reported towards Ca 2+ concentration (mol/L) expressed as pCa (see equation 1).

Figure 6 

 6 Figure 6  Shale-normalised REE pattern calculated in the dissolved phase of the lake (a) and

TABLE CAPTIONS Table 1 -

 CAPTIONS1 Sampling sites and geographical localisation of studied waters. Table 2 Mineralogical composition of studied suspended particulates in lake waters. Values are given in weight %.

ACKNOWLEDGMENT This research was partially funded by the contracts CORI 2017 (University of Palermo) and CON-0037 funded by the SIDERCEM S.R.L. -University of Palermo Agreement and The Science and Technological Research Council of Turkey (TUBITAK 118Y319).

Zhang, L., Algeo, T.J., Cao, L., Zhao, L., Chen, Z.-Q., Li, Z., 2016. Diagenetic uptake of rare earth elements by conodont apatite. Palaeogeogr. Palaeoclimatol. Palaeoecol. 458, Zuddas, P., [START_REF] Censi | Zr, Hf and REE distribution in river water under different ionic strength conditions[END_REF] The behaviour of zirconium and hafnium during water-rock interaction. Appl. Geochem. 94, 46-52. Zuddas, P., Inguaggiato, C., Censi, P., D'Alessandro, W., 2017. Zr-Hf fractionation during waterrock interaction. Procedia Earth and Planetary Science 17, 670-673.

thermal waters (e). These features are also reported in terms of MREE/HREE and LREE/MREE 463 ratios (f). 464 Lake waters

Thermal waters

Prof. Paolo Censi

Conflict of Interest

Supplementary Material

Click here to access/download Supplementary Material Supplementary data.pdf