Thomas Guérin
email: thomas.guerin@univ-lr.fr

Xavier Bertin
email: xavier.bertin@univ-lr.fr

Guillaume Dodet
email: guillaume.dodet@univ-brest.fr

A numerical scheme for coastal morphodynamic modelling on unstructured grids

Keywords: Morphodynamic modelling, unstructured grid, WENO, diffusion, coastal environments, Exner equation

Over the last decade, modelling systems based on unstructured grids have been appearing increasingly attractive to investigate the dynamics of coastal zones.

However, the resolution of the sediment continuity equation to simulate bed evolution is a complex problem which often leads to the development of numerical oscillations. To overcome this problem, addition of artificial diffusion or bathymetric filters are commonly employed methods, although these techniques can potentially over-smooth the bathymetry. This study aims to present a numerical scheme based on the Weighted Essentially Non-Oscillatory (WENO) formalism to solve the bed continuity equation on unstructured grids in a finite volume formulation. The new solution is compared against a classical method, which combines a basic node-centered finite volume method with artificial diffusion, for three idealized test cases. This comparison reveals that a higher accuracy is obtained with our new method while the addition of diffusion appears inappropriate mainly due to the arbitrary choice of the diffusion coefficient. Moreover, the increased computation time associated with the WENO-based method to solve the bed continuity equation is negligible when considering a fully-coupled simulation with tides and waves. Finally, the application of the new method to the pluri-monthly evolution of an idealized inlet subjected to tides and waves

Introduction

Coastal zones often display fast morphological changes, which can lead to socio-economical and environmental issues since a large part of the population lives in these areas. Moreover, sea-level rise and potential increase in storminess are likely to impact strongly these environments (IPCC, 2013). As a consequence, coastal management such as sediment dredging or erosion control plans becomes increasingly challenging. To better address these problems, morphodynamic modelling systems appeared as attractive tools and have experienced significant improvement during the last decades [START_REF] De Vriend | Analysis of horizontally two-dimensional morphological evolutions in shallow water[END_REF][START_REF] De Vriend | Approaches to long-term modelling of coastal morphology: A review[END_REF][START_REF] Cayocca | Long-term morphological modeling of a tidal inlet: the Arcachon Basin, France[END_REF][START_REF] Fortunato | A modeling system for tidally driven longterm morphodynamics[END_REF][START_REF] Bertin | Simulating morphodynamics with unstructured grids: Description and validation of a modeling system for coastal applications[END_REF][START_REF] Zhang | A coupled modeling scheme for longshore sediment transport of wave-dominated coastsa case study from the southern baltic sea[END_REF]. However, a common problem of these models is the development of numerical oscillations, due to both the decoupled way of solving the hydrodynamic and the sediment continuity (or Exner) equations, and the inherently unstable nature of the non-linear coupling between the sediment transport module and the bed evolution module [START_REF] Fortunato | Improving the stability of a morphodynamic modeling system[END_REF][START_REF] Long | A numerical scheme for morphological bed level calculations[END_REF]. In order to overcome this problem, fully coupled approaches, where the Exner-Saint-Venant system is solved simultaneously, have been successfully applied (e.g. Castro [START_REF] Castro Díaz | Twodimensional sediment transport models in shallow water equations. a second order finite volume approach on unstructured meshes[END_REF][START_REF] Soares-Frazão | HLLC scheme with novel wave-speed estimators appropriate for two-dimensional shallow-water flow on erodible bed[END_REF][START_REF] Bouharguane | Minimisation principles for the evolution of a soft sea bed interacting with a shallow sea[END_REF]. Unfortunately, this type of approach requires that the sediment flux only depends on the water depth and the fluid velocity (e.g. as in [START_REF] Meyer-Peter | Formulas for bed-load transport[END_REF] or [START_REF] Grass | Sediment transport by waves and currents[END_REF] formulae), which is not suitable in coastal zones where sediment transport is a much more complex process due to the presence of short waves. For coastal applications, the hydrodynamic and the sediment transport are usually treated separately and the problem of numerical oscillations is rather solved by using bathymetric filters and/or adding artificial diffusion [START_REF] Cayocca | Long-term morphological modeling of a tidal inlet: the Arcachon Basin, France[END_REF][START_REF] Johnson | Controlling spatial oscillations in bed level update schemes[END_REF]). Yet, these methods require the use of arbitrary thresholds or coefficient values, which potentially hides the physical behavior of the bed forms, while the root of the problem remains unsolved. Thus, the development of numerical schemes adapted to morphodynamic modelling has been the concern of extensive research effort during the last decade. [START_REF] Hudson | Numerical approaches for 1D morphodynamic modelling[END_REF] reviewed several methods for 1D morphodynamic systems, and investigated coupled solution of flow and bed-updating equations with Lax-Wendroff and Roe schemes with and without flux-limiting methods. This effort was extended to horizontally two-dimensional (2DH) morphodynamic modelling by [START_REF] Callaghan | Numerical solutions of the sediment conservation law; a review and improved formulation for coastal morphological modelling[END_REF], who applied a non-oscillating centered scheme (NOCS). Latter on, [START_REF] Long | A numerical scheme for morphological bed level calculations[END_REF] compared several numerical schemes to solve the Exner equation and showed that a weighted essentially non-oscillatory (WENO) scheme [START_REF] Liu | Weighted essentially non-oscillatory schemes[END_REF] with an Euler temporal discretization was the best compromise between computational time, accuracy, and numerical stability. However, these efforts concerned finite differences on regular grids whereas a significant tendency for developing unstructured grid (UG) versions of well-established models can been observed over the last years (e.g. SWAN [START_REF] Zijlema | Parallel, unstructured mesh implementation for SWAN[END_REF], DELFT3D [START_REF] Kernkamp | Efficient scheme for the shallow water equations on unstructured grids with application to the Continental Shelf[END_REF], or WaveWatchIII [START_REF] Tolman | User manual and system documentation of WAVEWATCH III version 4[END_REF], and only a few studies concerned morphodynamic modelling on UG (e.g. [START_REF] Kubatko | An unstructured grid morphodynamic model with a discontinuous galerkin method for bed evolution[END_REF][START_REF] Benkhaldoun | Comparison of unstructured finite-volume morphodynamic models in contracting channel flows[END_REF].

Using WENO schemes on UG has been investigated for solving two-dimensional conservation laws (e.g. [START_REF] Friedrich | Weighted essentially non-oscillatory schemes for the interpolation of mean values on unstructured grids[END_REF][START_REF] Hu | Weighted essentially non-oscillatory schemes on triangular meshes[END_REF][START_REF] Wolf | High-order eno and weno schemes for unstructured grids[END_REF], and even in three space dimensions [START_REF] Tsoutsanis | Weno schemes on arbitrary mixed-element unstructured meshes in three space dimensions[END_REF], but applications were restricted to the Euler and Burger equations. In particular, [START_REF] Liu | A robust reconstruction for unstructured weno schemes[END_REF] distinguished two types of finite volume WENO schemes on UG: (1) a first one designed for the purpose of nonlinear stability or to avoid spurious oscillations (being of our interest in the present study), and (2) a second one (more complex) providing higher order of accuracy for equal order of polynomial reconstruction. To our knowledge, the only application of a WENO scheme on UG to morphodynamic modelling was done by [START_REF] Canestrelli | Well-balanced highorder centered schemes on unstructured meshes for shallow water equations with fixed and mobile bed[END_REF], who employed a coupled solution strategy for solving the hydromorphodynamic system. As mentioned above, this approach cannot be applied for simulating morphodynamics in coastal areas because the sediment transport becomes also a function of wave parameters.

Alternatively, this study presents a numerical method for UG morphodynamic modelling based on the WENO formalism in a finite volume framework that is suitable for coastal applications. This method is implemented into the SED2D sediment transport and bed evolution module of [START_REF] Dodet | Morphodynamic modelling of a wave-dominated tidal inlet: the Albufeira Lagoon[END_REF], which was adapted from the sediment transport and bed evolution module SAND2D (Fortunato andOliveira, 2004, 2007), part of the 2DH morphodynamic modelling system MORSYS2D [START_REF] Bertin | Simulating morphodynamics with unstructured grids: Description and validation of a modeling system for coastal applications[END_REF] and the 3D morphodynamic modelling system MORSELFE [START_REF] Pinto | Development and validation of a three-dimensional morphodynamic modelling system for non-cohesive sediments[END_REF]. As in the SAND2D module, the original method for solving the Exner equation in SED2D uses nodecentered control volumes with sediment flux considered as constant inside each element. In the present modelling system, SED2D is coupled with the hydrodynamic model SELFE [START_REF] Zhang | Selfe: a semi-implicit eulerian-lagrangian finite-element model for cross-scale ocean circulation[END_REF], and the spectral wave model WWM-II [START_REF] Roland | A fully coupled 3D wave-current interaction model on unstructured grids[END_REF]. Three test cases are considered to assess the proposed scheme: (1) a migrating sandwave, allowing us to compare numerical and analytical results, (2) a migrating trench, where the robustness of the method in the presence of strong bathymetric gradients is analyzed, and (3) the pluri-monthly evolution of an idealized inlet subjected to tides and waves.

The morphodynamic modelling system

General outline of the modelling system

The core of the system is the Semi-implicit Eulerian-Lagrangian Finite Element (SELFE) modelling system of [START_REF] Zhang | Selfe: a semi-implicit eulerian-lagrangian finite-element model for cross-scale ocean circulation[END_REF], which has now evolved to SCHISM [START_REF] Zhang | Unstructured-grid model for the north sea and baltic sea: Validation against observations[END_REF], and is based on UG. The main feature of the circulation model in SELFE is the combination of an Eulerian-Lagrangian Method with semi-implicit schemes, to treat the advection in the momentum equations while relaxing the numerical stability constraints of the model (i.e. CFL condition can be exceeded). The Wind Wave Model II (WWM-II) of [START_REF] Roland | A fully coupled 3D wave-current interaction model on unstructured grids[END_REF] (third generation, spectral wave model) is coupled to SELFE and simulates gravity waves generation and propagation by solving the wave action equation (WAE) [START_REF] Komen | Dynamics and Modelling of Ocean Waves[END_REF]. WWM-II uses a residual distribution scheme [START_REF] Abgrall | Residual distribution schemes: Current status and future trends[END_REF] to solve the geographic advection in the WAE, which also relaxes CFL constraints and allows using large time step without compromising the numerical stability. The 2DH sediment transport/bottom evolution module SED2D [START_REF] Dodet | Morphodynamic modelling of a wave-dominated tidal inlet: the Albufeira Lagoon[END_REF] computes sediment fluxes (total load, i.e. sum of bed-load and suspended load) with classical semi-empirical formulations based on depth-averaged velocity, water depth, bottom roughness, sediment properties and wave parameters. The bed evolution over the morphological time step is then computed by solving the Exner equation, this part being detailed in the following sections since this is the core of the present study. This modelling system is fully-coupled, parallelized, and the three modules share the same computational grid and domain-decomposition.

Bed evolution equation and finite volume formulation

The bottom evolution module computes the bed change at each grid node by solving the sediment continuity/Exner equation, given by:

∂z b (x, t) ∂t + 1 1 -λ ∇ • Q(x, t) = 0 (1)
where x = (x, y), z b (x, t) is the bed level elevation (positive upwards), λ is the sediment porosity, and Q = (Q x , Q y ) is the depth-integrated sediment transport rate (in m 3 .s -1 .m -1 ) computed at element centres by the sediment transport module.

Considering node-centered control volumes (Fig. 1), the semi-discrete finite volume formulation (continuous in time, discrete in space) of Eq. 1 can be written as:

∂ ∂t Ωi z b dΩ = - 1 1 -λ Γi Q • n dΓ (2)
with Ω i the control volume (or cell) for node i, Γ i the corresponding boundary, and n the outward unit normal to Γ i . Using an Euler explicit time discretization, we have the fully-discrete finite volume form:

Ωi

∆z b dΩ = - ∆t 1 -λ Γi Q • n dΓ (3)
where ∆z b is the bed change during the morphological time step ∆t. where N el is the number of elements neighboring node i, and Ω i,el is the part of Ω i belonging to element el. S(el, nd) is the element linear shape function that equals 1 at node nd = i and 0 at the two other nodes of the element, which 120 gives:

Ω i,el S(el, nd) dΩ = C nd A i,el (5)
where A i,el is the area of element el neighboring node i, and

C nd =      22/108 if nd = i 7/108 if nd = i (6)
Once right-hand side of Eq. 3 is computed (see section 3), a system of N nd equations with N nd unknowns is obtained (N nd is the total number of grid nodes) and eventually solved with a Jacobi conjugate gradient method.

A fourth-order Runge-Kutta (RK) time discretization was also considered in order to increase the morphological time step but this method implies performing four times the WENO scheme described below for spatial discretization, for each time step. Since the subsequent increase in computation time neither balanced the gain in numerical stability nor improved substantially the accuracy, the Euler explicit time discretization was retained. Similarly, it can be noted that [START_REF] Long | A numerical scheme for morphological bed level calculations[END_REF] did not observe any significant quantitative change in results by considering a third-order RK scheme rather than a simple Euler explicit scheme for time discretization, with a WENO scheme for spatial discretization.

The new numerical method

Contrary to the original method implemented in SED2D where the sediment flux is assumed to be constant inside an element, the main feature of the WENO scheme is to compute a reconstruction polynomial P i (x) for each control volume in order to interpolate the sediment flux at the corresponding boundaries.

Spatial discretization

Each control volume Ω i defines a cell which is polygonally bounded, with a finite number of line segments. Therefore, replacing sediment fluxes Q by P i , the integral from Eq. 3 can be decomposed into:

Γi Q • n dΓ = Γi P i • n dΓ = j Γi,j P i • n dΓ (7)
with j the line segment index. Each line integral is then discretized by a q-point Gaussian integration formula:

Γi,j

P i • n dΓ ≈ |Γ i,j | q k=1 ξ k P i (G k ) • n (8)
where G k and ξ k are the Gaussian points and weights. We use q = 2, so with

x 1 and x 2 being the end points of the line segment Γ i,j , the position of

G k are x(G 1 ) = αx 1 + (1 -α)x 2 and x(G 2 ) = αx 2 + (1 -α)x 1 , with α = 1/2 + √ 3/6
and ξ 1 = ξ 2 = 1/2.

Polynomial reconstruction procedure

(a) Following a WENO procedure, we need to select several stencils for each cell Ω i and to compute the corresponding polynomials which interpolate sediment flux over the cell. As we want a numerical method with a relatively low computational cost, each stencil related to Ω i is defined by three elements neighboring node i (Fig. 2), such as a linear polynomial is computed for each stencil, from the values of sediment flux computed at element centers. Only continuous stencils are considered (i.e. for each stencil, there is no gap between the three elements) which avoids interpolation across discontinuities as recommended in case of non-smooth solution [START_REF] Friedrich | Weighted essentially non-oscillatory schemes for the interpolation of mean values on unstructured grids[END_REF]. Consequently, if node i is an interior grid node, the number N of stencils related to Ω i equals the number of elements neighboring node i. Moreover, using these basic stencils facilitates the implementation of the method on parallelized codes since there is no need to reach an element which is not a direct neighbor of node i. (b) For each stencil, the two linear polynomials corresponding to both com-ponents of the sediment flux are computed as:

node i

     p x,m (x) = p x,m (x, y) = a x,m x + b x,m y + c x,m p y,m (x) = p y,m (x, y) = a y,m x + b y,m y + c y,m (9)
where m is the stencil index, and such as for each element ∆ e belonging to stencil m we have:

     p x,m (x c (∆ e )) = Q x (x c (∆ e )) p y,m (x c (∆ e )) = Q y (x c (∆ e )) (10)
where

Q x (x c (∆ e
     p x,m (x) ∆e = p x,m (x c (∆ e )) = Q x (x c (∆ e )) p y,m (x) ∆e = p y,m (x c (∆ e )) = Q y (x c (∆ e )) ( 11
)
where ∆e is the spatial mean operator over ∆ e .

(c) Aiming to measure the smoothness of p m = (p x,m , p y,m ) (i.e. how much p m varies spatially), an oscillating indicator is computed for each stencil based on [START_REF] Friedrich | Weighted essentially non-oscillatory schemes for the interpolation of mean values on unstructured grids[END_REF]:

OI m = OI x,m + OI y,m (12)
For the x-component we have:

OI x,m = Ωi dX -2 ∂p x,m (x, y) ∂x 2 + ∂p x,m (x, y) ∂y 2 dΩ 1/2 (13)
leading in our case to (d) While an Essentially Non-Oscillating (ENO) scheme [START_REF] Harten | Uniformly high-order accurate non-oscillatory schemes, I[END_REF] would only keep the linear polynomial having the lowest OI m value, the WENO scheme considers a weighted combination of the N linear polynomials to compute the reconstruction polynomial. The weights ω m are computed such that their sum is one, following:

OI x,m = |Ω i | dX 2 (a 2 x,m + b 2 x,m ) (14
ω m = ( + OI m ) -r N k=1 ( + OI k ) -r (15)
where is a small value compared to OI m ensuring a non-zero denominator (we take = 10 -10 m.s -1 ), and r is a positive integer. [START_REF] Friedrich | Weighted essentially non-oscillatory schemes for the interpolation of mean values on unstructured grids[END_REF] indicates that the weights should be of magnitude one for stencils in smooth regions while it should be low in discontinuous regions, this condition being fulfilled for any positive r. A sensitivity analysis leads us to take r = 1.

(e) The reconstruction polynomial at node i is finally computed as:

P i (x) = N k=1 ω k p k (x) (16)
with P i (x) = (P x,i , P y,i ) and p k (x) = (p x,k , p y,k ).

Regarding boundary conditions, the two following cases are considered:

1) If the number N el of elements neighboring node i (where i belongs to the grid boundary(ies)) is such that N el ≥ 3, then the number of stencils used to

compute P i is N ≥ 1.
2) If N el < 3, then no stencil is defined, and P i is simply computed such that for the one or two elements ∆ e neighboring node i:

P i (∆ e ) = Q(x c (∆ e )).

Numerical flux

For each line segment Γ i,j of a cell Ω i , the sediment flux at Gaussian points is approximated by the two reconstruction polynomials P i and P l , the latter corresponding to the neighbor cell Ω l (Γ i,j being the shared boundary segment of both cells). This allows to compute the following two values for right-hand side of Eq. 8:

F i,j = |Γ i,j | q k=1 ξ k P i (G k ) • n = |Γ i,j | 1 2 (P i (G 1 ) + P i (G 2 )) • n (17) F l,j = |Γ l,j | 1 2 (P l (G 1 ) + P l (G 2 )) • n (18) with |Γ i,j | = |Γ l,j |.
A flux limiter (FL) is then applied in order to handle the strongest sediment flux gradients, such as:

F F L i,j = F i,j + 1 2 φ(r F L )(F i -F i,j ) (19)
F F L l,j = F l,j + 1 2 φ(r F L )(F l -F l,j ) (20)
with

F i = |Γ i,j |(P i (x i ) • n) and F l = |Γ l,j |(P l (x l ) • n). The FL function of
Chatkravathy and Osher is used [START_REF] Chakravarthy | High resolution applications of the Osher upwind scheme for the Euler equations[END_REF], which reads φ(r F L ) = max(0, min(r F L , β)), with 1 ≤ β ≤ 2. Through the r F L value, the FL function φ(r F L ) quantifies the upwinding which is added to the scheme.

Important care is taken to define r F L , such that it tends to zero for smooth solutions and it increases near discontinuities. Since the sediment flux is a nonlinear function of the water depth h (always positive), we take r F L = |∆h| h with ∆h = h(i) -h(l) and h = 1 2 (h(i) + h(l)). Moreover we take β = 2, allowing a maximum upwinding for the numerical flux. Indeed, we have F F L i,j = F i,j and F F L l,j = F l,j if r F L = 0 (i.e. no effect of the FL on the scheme), whereas we have

F F L i,j = F i and F F L l,j = F l if r F L ≥ β (i.e
. a maximum upwinding is added to the scheme).

Finally, Eq.3 is solved by using an upwind flux formula to compute the final flux at each line segment of cell Ω i :

F f inal i,j =      min(F F L i,j , F F L l,j ) if z b (i) < z b (l) max(F F L i,j , F F L l,j ) if z b (i) ≥ z b (l) (21)
4. Numerical results

Test case 1: Migrating sandwave

We first apply both the original and the new numerical method of SED2D

to the 2DH migration test case of an initially sinusoidal sandwave under unidirectional and stationary flow in a straight channel, similarly to the 1D test case of [START_REF] Hudson | Numerical approaches for 1D morphodynamic modelling[END_REF]. We recall that an uncoupled solution strategy is used in this study, i.e. the hydrodynamic (fluid velocity and surface elevation) is first solved by SELFE, allowing SED2D to compute the sediment transport and to solve the Exner equation. In order to compare the numerical result with the analytical solution, a simple transport rate function is considered, given by:

     Q = (Q x , Q y ) = (au b x , 0) u x = D x (h∆y) -1 (22)
where a and b are constants, u = (u x , 0) is the depth-averaged current velocity (m.s -1 ), D = (D x , 0) is the constant water discharge (m 3 .s -1 ), h = η -z b ≥ 0 is the water depth (with the mean water level η = 0 in the present case), and ∆y = 1.2 m is the channel width.

In order to devise a stringent test for the new method, the bed slope effect on the sediment transport is not considered in this first test case (unlike in the next two test cases), allowing us to obtain the corresponding analytical solution of the Exner equation by using the method of characteristics:

z b (x, t) = z b (x -c z t, 0) (23)
where c z = (c x,z , 0) is the phase velocity of the bedform:

c x,z (z b ) = 1 1 -λ ∂Q x ∂z b = 1 1 -λ abu b x z b (24)
The evolution of the sandwave is simulated using a = 0.001 s 2 .m -1 , b = 3, D x = 1 m 3 .s -1 , and λ = 0.4, which yields a maximum Courant number of about 0.1 if estimated according to [START_REF] Damgaard | Morphodynamic modelling of rip channel growth[END_REF] and [START_REF] Roelvink | Coastal morphodynamic evolution techniques[END_REF] by max|c z |∆t/∆x, with in our case ∆t = 2 s and ∆x = 0.15 m. The Euler-WENO (EW) scheme is compared against the original Euler node-centered finite volume method of SED2D, in which sediment flux is assumed to be constant inside each element. Since this latter scheme is proned to develop numerical oscillations even for Courant numbers below unity, we also include a diffusion-like term in the sediment transport formula which is common practice to stabilize the bed evolution in morphodynamic modelling [START_REF] Rakha | A morphology model for an eroding beach backed by a seawall[END_REF][START_REF] Cayocca | Long-term morphological modeling of a tidal inlet: the Arcachon Basin, France[END_REF][START_REF] Fortunato | Improving the stability of a morphodynamic modeling system[END_REF]. This additional diffusion method consists in replacing the sediment transport rate Q by

Q * = Q -ε(1 -λ)(|Q x | ∂z b ∂x , |Q y | ∂z b ∂y ) ( 25
)
where ε is a dimensionless coefficient, with usually ε ∈ [0, 5]. Fig. 3 (a) shows the bed profiles at time t = 500 s and along y = 0.75 m for the original scheme without and with additional diffusion (ε = 1), and for the EW scheme. While the original scheme without additional diffusion shows the emergence of numerical oscillations at the dune crest, accuracy is well improved with the EW scheme, as confirmed by the associated errors (Fig. 3 (b)). The root-mean-square errors for the original scheme without diffusion and for the EW scheme are 2.8 mm and 0.8 mm, respectively. An over-smoothing of the dune is obtained for the original scheme with additional diffusion, and will be discussed in more details in the next sections. The convergence analysis verifies this increased accuracy obtained with the EW scheme (Fig. 4), especially for dx < 0.08 m where the original scheme becomes highly unstable (for this particular case oscillations are not developing near maximum transport gradients, which would suggest a potential spatial limit for the original scheme).

- In order to test the EW scheme with a more complex sediment transport formula than in test case 1, the formula of van Rijn (2007a,b) is used to compute both bed-load (q b ) and suspended load transport (q s ): The bed slope effect on the sediment transport is considered following the method of [START_REF] Lesser | Development and validation of a three-dimensional morphological model[END_REF], and the Exner equation is finally solved for the total transport q tot = q b + q s . A median diameter of 0.14 mm is used, while the time step is set to 1 s, satisfying the equivalent Courant number stability criterion. The bed profiles at mid-width channel shown on Fig. 5 after 1700 s of simulation confirm the enhanced stability of the EW scheme compared to the original scheme. Unlike the previous test case, the inclusion of artificial diffusion with the same coefficient value (ε = 1) strongly improves the results while no large over-smoothing of the bed profile is observed.

     q b = 0.015uh(d 50 /h) 1.2 M 1.5 e q s = 0.012ud 50 M 2.

Test case 3: Idealized inlet

In order to evaluate the improvement of our new method with a more realistic case, we applied our modelling system to the idealized coastal lagoon of Nahon et al. ( 2012) (Fig. 6) where tides and waves are considered. This test case is more challenging than the previous ones because the combination of waves and tidal forcings yields both a large variability of sediment fluxes and strong gradients over the domain.

The lagoon has an initial depth of 2.5 m relative to Mean Sea Level (MSL)

and is connected to the sea through a 700 m long and 300 m wide shore-normal oriented channel. The beach/shore face profile is alongshore uniform and goes from 2 m above MSL down to 24 m depth, with maximum slopes of 0.014 at the beach berm and 0.004 offshore. The grid resolution ranges from 300 m at the open boundary down to 25 m at the inlet. As for test case 2, bed-load and suspended load transport are computed using van Rijn (2007a,b) formula (see Eq. 26), with:

M e = (max(0, |u| + γU w -u cr ))/ [(s -1)gd 50 ] 0.5 (29)
where U w is the amplitude of the wave orbital velocity and γ = 0.4 for irregular waves. Following van Rijn (2007a), the critical fluid velocity for initiation of sediment motion in the presence of current and waves is:

u cr = βu cr,c + (1 -β)u cr,w (30)
where β = |u|/(|u| + U w ), and u cr,w is the critical wave orbital velocity for initiation of sediment motion computed as:

u cr,w =      0.24 ((s -1)g) 0.
66 (d 50 ) 0.33 (T p ) 0.33 for 0.05 < d 50 < 0.5 mm 0.95 [(s -1)g] 0.57 (d 50 ) 0.43 (T p ) 0.14 for 0.5

< d 50 < 2 mm (31)
where T p is the wave peak period. As in the previous test case, the bed slope effect on the sediment transport is considered following [START_REF] Lesser | Development and validation of a three-dimensional morphological model[END_REF].

A mixed-energy regime is considered for this test case, meaning that the ratio between the yearly-averaged tidal range and the significant wave height is approximately in the range [1, 2] according to [START_REF] Hayes | Barrier island morphology as a function of tidal and wave regime[END_REF]. Oliveira ( 2007) was also added to this original scheme, aiming to eliminate local extrema in the bathymetry after each morphological time step. On the opposite, the EW scheme is applied without any artificial diffusion nor bathymetric filter, as for the previous test cases.

By analyzing the bathymetry simulated with both schemes after 3 and 5 months on Fig. 7 (taking about 20 hours on 24 processors), several differences can be noticed. First, the main channel is found to be about 2 m deeper with the EW scheme than with the original one. Besides, due to the wave-induced southward littoral drift, sediment accretion is observed at the northern (updrift) side of the inlet. This causes a counterclockwise rotation of the main channel axis, in agreement with mixed-energy-straight inlets described in [START_REF] Davis | Morphodynamics of the barrier-inlet system, west-central Florida[END_REF], this evolution being more pronounced with the EW scheme. Moreover, using the EW scheme leads to the development of a secondary flood channel on the updrift side of the ebb-delta, and shore-parallel sandbars on its downdrift side, unlike using the original method (see also Fig. 8 (a), (b)). Finally, we observe the development of shore-oblique sandbars along the adjacent shorelines only with the EW scheme (Fig. 8 (e)). On the other hand, the bathymetry at once for all test cases. Indeed, with ε = 1, the numerical result is oversmoothed for test case 1, correct for test case 2, and oscillating for test case 3

(not shown but leading us to use a higher value in this case). The problem is that this coefficient requires to be arbitrarily user-defined and does not depend on a relevant parameter, such as the local Courant number. This tuning being specific for each test case, the coefficient value will not even suit over the whole computational grid for some test cases, due to the variable bathymetry and hydrodynamic conditions. This implies to choose a relatively high value to overcome the development of numerical oscillations, but with the drawback of over-smoothing some bed features. This behavior is illustrated with the test case of an idealized inlet subjected to tides and waves, where a higher bathymetric complexity is captured when using the EW scheme. It handles relatively strong sediment transport gradients without over-smoothing the bathymetry where these gradients are lower, unlike the additional diffusion method. Moreover, our proposed method constitutes an alternative to the discontinuous Galerkin method of [START_REF] Kubatko | An unstructured grid morphodynamic model with a discontinuous galerkin method for bed evolution[END_REF] which, despite its higher accuracy, may increase the computation time substantially [START_REF] Budgell | Scalar advection schemes for ocean modelling on unstructured triangular grids[END_REF]. As shown on Fig. 9, this is not the case here since using the EW scheme instead of the original one leads to an increase of the SED2D computation time by a factor less than two, which in the end appears negligible when looking at the total computation time (i.e. for a fully-coupled run). This point is of great importance for long-term morphodynamic modelling (as shown in Guérin ( 2016)), and also when multiple sediment classes are considered where the Exner equation is solved for each class.

Implications for real-world applications

Morphological predictions obtained with the EW scheme substantially differ from those obtained with the original method when simulating an idealized inlet subjected to tides and waves. Indeed, after 5 months of simulation, the inlet main channel is about 2 m shallower when using the original method, which can be explained by an over-smoothing effect of the additional diffusion. A detailed analysis also reveals that several bed features only develop with the EW scheme.

First, a secondary flood channel develops on the updrift side of the ebb-delta while this morphological unit is commonly observed at many tidal inlets, such as the Maumusson inlet (Fig. 8 (b) and (c), marker 1). Secondly, ebb-delta sandbars develop on the downdrift side of the inlet and migrate onshore until they eventually weld onto the beach (Fig. 8 (b) and (c), marker 2). This common behaviour of tidal inlets is also well documented while the modeled migration rate of 1.5 to 3 m.day -1 is coherent with some observations (e.g. [START_REF] Pianca | Mobility of meso-scale morphology on a microtidal ebb delta measured using remote sensing[END_REF]. Finally, periodic oblique sandbars develop along the adjacent shorelines only with the EW scheme. As studied by [START_REF] Garnier | Generation and nonlinear evolution of shore-oblique/transverse sand bars[END_REF] with a 2DH morphodynamic model, these bed features can emerge by self-organization of the coupling between waves, currents and morphology via sediment transport.

A wavelength range of about 350 to 500 m is obained in our case, which is consistent with observations (e.g. [START_REF] Castelle | Double bar beach dynamics on the high-energy meso-macrotidal french aquitanian coast: a review[END_REF] measured a range of 360 to 470 m; see Fig. 8 (f) for illustration). Although their physical significance cannot be formally demonstrated from this study, we expect that applications to realistic sites will greatly benefit from our proposed method. Moreover, the mean intertidal cross-shore bed slope obtained with the EW scheme after several months (∼ 0.01) remains close to the initial one, whereas it reaches very large values (∼ 0.1) with the original method while the bathymetry turns unrealistic. Indeed, the increase of cross-shore bed slope reduces the surfzone width, which increases the gradients of wave radiation stress and in turn increases the waveinduced longshore current. As sediment transport is a non-linear function of the current velocity, this problem may cause large errors in longshore transport rates and impact the evolution of the inlet significantly.

Conclusion

In order to improve an existing unstructured grid, 2DH, morphodynamic modelling system, a numerical scheme combining an Euler temporal discretization and a WENO formalism for spatial discretization is used to solve the Exner equation. Through three idealized test cases, this numerical method is compared to the original one of SED2D module, which stability is guaranteed through the inclusion of additional diffusion. The first two test cases demonstrate the enhanced accuracy of the EW scheme over the original one. Indeed, the additional diffusion method is shown to be inappropriate since it remains arbitrary and does not solve the problem locally. The advantages of the new method are also evaluated through the pluri-monthly morphodynamic simulation of an idealized inlet subjected to tides and waves. Non-oscillating and realistic bed evolutions were obtained, as partly attested when confronting the development and evolution of several bedforms (e.g. ebb-delta sandbars, secondary flood channel, or oblique sandbars at adjacent beaches) to related studies and satellite images. Moreover, the additional computation time due to the use of the EW scheme appears negligible when considering the total computation time (i.e. for a fully-coupled run with waves and tidal forcings). Our new method can be implemented in any UG, 2DH, parallelized, morphodynamic modelling system, but also in 3D models where the Exner equation is solved for bedload transport. Future work will be to use the EW scheme in realistic test cases and to compare its advantages with alternative methods, such as the residual distribution schemes [START_REF] Abgrall | Residual distribution schemes: Current status and future trends[END_REF] which proved their efficiency in the wave model WWM-II.

Figure 1 :

 1 Figure 1: Node-centered control volume Ω i and associated variables.

115

  Bed level elevation z b (known at grid nodes) is assumed to vary linearly within each element, allowing us to express left-hand side of Eq. 3 as:

Figure 2 :

 2 Figure 2: Example of a stencil (gray color) defined by three elements neighboring node i.

  ) and Q y (x c (∆ e )) are the sediment flux components computed by the sediment transport module at the centre x c of element ∆ e . Considering these two values as the mean values of each sediment flux component over element ∆ e , they are conserved by (p x,m , p y,m ) since:

  ) with the grid spacing dX = |∆ e | m , |∆ e | being the area of each element belonging to stencil m. OI y,m is computed by replacing (a x,m , b x,m ) by (a y,m , b y,m ) in Eq. 14. Since OI m is function of a 2 m = (a 2 x,m , a 2 y,m ) and b 2 m = (b 2 x,m , b 2 y,m ), it vanishes in areas of constant sediment fluxes whereas it increases in areas of variable fluxes. The stencils corresponding to the lowest values of OI m will then be favored for computing the reconstruction polynomial, through the weighted average procedure described in the following.

Figure 3 :

 3 Figure 3: Comparison of Euler-original without and with additional diffusion, and Euler-WENO scheme results to analytical solution at t = 500 s and y = 0.75 m: bed profiles (a), and associated errors (b).

Figure 4 :

 4 Figure 4: Convergence plot for test case 1: mean order of convergence is 1.22 for the Euler-WENO scheme.

  is the median sediment diameter, and D * = d 50 g(s -1)/ν 2 1/3 is the dimensionless grain diameter, with ν the kinematic fluid viscosity and s = ρ s /ρ the specific sediment density (ρ and ρ s are the density of water and sediment respectively). Following van Rijn (2007a), the mobility parameter M e is computed as: M e = max(0, |u| -u cr,c )/ [(s -1)gd 50 ] d 50 ) 0.1 log(4h/d 90 ) for 0.05 < d 50 < 0.5 mm 8.5(d 50 ) 0.6 log(4h/d 90 ) for 0.5 < d 50 < 2 mm (28)

Figure 5 :

 5 Figure 5: Comparison of Euler-original without and with additional diffusion, and Euler-WENO scheme results for test case 2 at t = 1700 s and y = 0.55 m.

  Figure 6: Computational grid of the idealized inlet test case, with zoom on initial bathymetry of the inlet.

Figure 8 :

 8 Figure 8: Bathymetry of the idealized inlet (t = 4 months) and the updrift coast (t = 7 months) simulated using the original method with diffusion (a), (d), and the EW scheme (b), (e). (c) The mixed-energy inlet of Maumusson (Atlantic coast, Charente-Maritime, France ; Landsat image), exhibiting a secondary flood channel (1) and an emergent ebb-delta sandbar (2). (f) Shore-oblique sandbars near Cap Ferret (Atlantic coast, Gironde, France ; Google Earth, august 2012).

Figure 9 :

 9 Figure9: Computation times for the idealized inlet test case (6-hours evolution), for the original scheme with additional diffusion and non linear filter, and for the EW scheme. The computation times for SED2D module and the fully-coupled modelling system (SELFE-WWM-SED2D) are plotted.

Acknowledgements

The development of the model with respect to sediment transport and morphodynamics was performed in the context of the project DYNAMO, funded by the French National Research Agency (Grant agreement no. ANR-12-JS02-00008-01). The developers of SELFE and WWM are also greatly acknowledged.

unrealistic (Fig. 8 (d)).

Discussion

Improvements compared to alternative methods

The three test cases clearly show that the additional diffusion method appears problematic since no unique value of the diffusion coefficient is suitable