M Ciallella
email: mirco.ciallella@inria.fr

L Micalizzi
email: lorenzo.micalizzi@math.uzh.ch

P Öffner
email: poeffner@uni-mainz.de

D Torlo

An Arbitrary High Order and Positivity Preserving Method for the Shallow Water Equations

Keywords: positivity preserving, well-balanced, WENO, modified Patankar, shallow water, deferred correction

In this paper, we develop and present an arbitrary high order well-balanced finite volume WENO method combined with the modified Patankar Deferred Correction (mPDeC) time integration method for the shallow water equations. Due to the positivity-preserving property of mPDeC, the resulting scheme is unconditionally positivity preserving for the water height. To apply the mPDeC approach, we have to interpret the spatial semi-discretization in terms of production-destruction systems. Only small modifications inside the classical WENO implementation are necessary and we explain how it can be done. In numerical simulations, focusing on a fifth order method, we demonstrate the good performance of the new method and verify the theoretical properties.

Introduction

In the last years, the development of structure preserving high-order methods for hyperbolic conservation/balance laws have been an active field of research [START_REF] Abgrall | Reinterpretation and extension of entropy correction terms for residual distribution and discontinuous Galerkin schemes[END_REF][START_REF] Abgrall | High order asymptotic preserving deferred correction implicit-explicit schemes for kinetic models[END_REF][START_REF] Chen | Review of entropy stable discontinuous Galerkin methods for systems of conservation laws on unstructured simplex meshes[END_REF][START_REF] Kuzmin | Entropy conservation property and entropy stabilization of high-order continuous Galerkin approximations to scalar conservation laws[END_REF][START_REF] Gaburro | A unified framework for the solution of hyperbolic PDE systems using high order direct Arbitrary-Lagrangian-Eulerian schemes on moving unstructured meshes with topology change[END_REF][START_REF] Thomann | An all speed second order well-balanced IMEX relaxation scheme for the Euler equations with gravity[END_REF][START_REF] Veiga | An arbitrary high-order spectral difference method for the induction equation[END_REF]. In the context of the shallow water equations, one is mainly interested in maintaining positive levels of water height, in conserving the equilibrium/stationary states and in entropy conservation/dissipation methods. There exists various ways to obtain these desired results, e.g. the applications of limiters for the positivity is only one example, cf. [START_REF] Berberich | High order well-balanced finite volume methods for multi-dimensional systems of hyperbolic balance laws[END_REF][START_REF] Ricchiuto | On the C-property and generalized C-property of residual distribution for the shallow water equations[END_REF][START_REF] Mantri | Well-balanced discontinuous Galerkin scheme for 2× 2 hyperbolic balance law[END_REF][START_REF] Noelle | High-order well-balanced finite volume WENO schemes for shallow water equation with moving water[END_REF][START_REF] Ranocha | Shallow water equations: Split-form, entropy stable, well-balanced, and positivity preserving numerical methods[END_REF][START_REF] Cheng | A new approach for designing moving-water equilibria preserving schemes for the shallow water equations[END_REF][START_REF] Xing | A survey of high order schemes for the shallow water equations[END_REF] and references therein. In this paper, we also deal with these issues and we present a new high-order, well-balanced, positivity preserving method for the shallow water equation starting from a classical WENO scheme. In order to obtain well-balanced (WB) solutions, we subtract the residual of the a priori known stationary solution from our numerical scheme as shown in [START_REF] Berberich | High order well-balanced finite volume methods for multi-dimensional systems of hyperbolic balance laws[END_REF]. Then, to ensure the positivity of the water height, the modified Patankar Deferred Correction (mPDeC) method is used for the time-integration of this variable. Even if (modified) Patankar (mP) methods have been already used inside a numerical method for fluid simulations [START_REF] Huang | Positivity-preserving time discretizations for production-destruction equations with applications to non-equilibrium flows[END_REF][START_REF] Huang | A third-order unconditionally positivity-preserving scheme for productiondestruction equations with applications to non-equilibrium flows[END_REF][START_REF] Meister | A positivity preserving and well-balanced DG scheme using finite volume subcells in almost dry regions[END_REF], those mP schemes have been based on extensions of classical RK methods and they are of maximum order three 1 . By applying the 1 modified Patankar trick inside the Deferred Correction (DeC) framework, the authors of [START_REF] Öffner | Arbitrary high-order, conservative and positivity preserving Patankar-type deferred correction schemes[END_REF] were able to construct a conservative, arbitrarily high-order and positivity preserving method for production-destruction systems (PDS) of ordinary differential equations. To obtain a positive WENO spatial reconstruction, a positive limiter must be used [START_REF] Zhang | On positivity-preserving high order discontinuous Galerkin schemes for compressible Euler equations on rectangular meshes[END_REF][START_REF] Perthame | On positivity preserving finite volume schemes for Euler equations[END_REF]. In this work the mPDeC is applied for the first time to a PDE problem, with a finite volume WENO spatial discretization. This lead to a high order accurate method enjoying all the previously cited properties (WB, positivity preservation). Finally, in order to apply mPDeC on the semi-discretized problem, the finite volume method must be rewritten into a PDS and we explain in details how this can be done. It should be stressed out that only small modifications inside the classical finite volume implementation are necessary as it can be seen in the reproducibility repository [START_REF] Ciallella | Modified Patankar Deferred Correction WENO Code for Shallow Water Equations[END_REF]. To our opinion the modifications can be adapted to most WENO codes in a straightforward manner and the approach is a good alternative to already existing methods. The paper is structured as follows: In Section 2, we introduce the considered model, the classical shallow water equations, and we repeat the basic properties focusing on steady state solutions and the positivity of the water height. Next, in Section 3 we describe the used classical finite volume WENO approach from [START_REF] Shu | Essentially non-oscillatory and weighted essentially non-oscillatory schemes for hyperbolic conservation laws[END_REF][START_REF] Shu | Efficient implementation of essentially non-oscillatory shock-capturing schemes[END_REF] and its well-balanced modification from [START_REF] Berberich | High order well-balanced finite volume methods for multi-dimensional systems of hyperbolic balance laws[END_REF]. We focus on the fifth order WENO method. However, all the ingredients to go to arbitrary high order can be found in the related repository [START_REF] Ciallella | Modified Patankar Deferred Correction WENO Code for Shallow Water Equations[END_REF]. In Section 4, the time-integration is considered focusing on Deferred Correction (DeC), the modified Patankar approach and its combination to the modified Patankar DeC (mPDeC) method developed in [START_REF] Öffner | Arbitrary high-order, conservative and positivity preserving Patankar-type deferred correction schemes[END_REF]. In Section 5, we describe how mPDeC can be combined with the WENO approach. It is important to interpret the semi-discretization in terms of production-destruction systems. Details of the implementation are given with additional algorithms. Then, in Section 6, we verify the theoretical properties of the scheme with numerical simulations with WENO5 focusing on the high-order accuracy, the well-balanced and the positivity preserving properties. In addition, we demonstrate also the excellent performance for more challenging test cases. Finally, in Section 7 we summarize the obtained results and perspectives for future works. In Appendixes A and B we describe in details the WENO reconstruction and apply it for WENO5 with 4-points Gaussian quadrature rule, which, up to our knowledge, is not available in literature.

Shallow Water Equations 2.1 Model

The shallow water equations (SWE) model the behaviour of shallow free surface flows under the action of gravity. They are used to simulate the flows in rivers and coastal areas, and can be applied to predict tides, storm surge levels and coastline changes from hurricanes and ocean currents. They are also used in atmospheric flows, debris flows, and certain hydraulic structures like open channels and sedimentation tanks. SWEs take the form of non-homogeneous hyperbolic conservation laws with source terms modeling the effects of bathymetry and viscous friction. In this paper, we will consider the effect of the bathymetry as the only source term. If the bottom topography is assumed to be constant with respect to time, the SWEs can be recast in balance law form as:

∂u ∂t + ∇ • F(u) = S(u, x, y) on Ω T = Ω × [0, T ] ⊂ R 2 × R + (1)
with conserved variables, flux and source terms given by

u =   h hu hv   , F(u) = F G =   hu hv hu 2 + g h 2
where h represents the relative water height, u = (u, v) are the flow speed components, g is the gravity acceleration and b(x, y) is the local bathymetry. The source term helps modeling the effects induced on the flow caused by the bathymetry changes in space. Finally, it is also convenient to introduce the free surface water level η := h + b. All the aforementioned variables can be better interpreted by looking at Figure 1.

Properties of the model

As it is described inter alia in [START_REF] Bollermann | A well-balanced reconstruction of wet/dry fronts for the shallow water equations[END_REF][START_REF] Berthon | A fully well-balanced, positive and entropy-satisfying Godunov-type method for the shallow-water equations[END_REF][START_REF] Lukáčová-Medvidová | Well-balanced finite volume evolution Galerkin methods for the shallow water equations[END_REF][START_REF] Ricchiuto | An explicit residual based approach for shallow water flows[END_REF][START_REF] Xing | A survey of high order schemes for the shallow water equations[END_REF], the construction and development of effective and accurate numerical methods for the shallow water equations have received much interest in the last decades and it is still ongoing. In particular, one is interested in schemes that preserve physical quantities or structures from the continuous level. In this paper, we are targeting two types of difficulties which are often encountered when working with SWE: the preservation of steady state solutions and water height positivity2 .

Steady state solutions

The SWE system (2) is known to admit some steady state solutions whose form depends on the equilibrium between the source terms S and the remaining terms of the equations. The numerical simulations should be able to capture these behaviors even on coarse grids. Without additional techniques, many methods fail at balancing the source terms and the flux, resulting in small pertubations of the steady state. The perturbations could be then amplified by the method causing a bad approximation of the exact behaviour. This situation is sometimes called numerical storm in such context. To prevent it, one is interested in schemes that are capable of exactly balancing the flux and the source terms to obtain the desired steadystate solution. Numerical schemes enjoying this property are called well-balanced schemes. The still water surface is often the first equilibrium taken in consideration. It is given by u = v = 0; η(x, y, t) = h(x, y, t) + b(x, y)

≡ η 0 ∈ R + 0 , ∀ (x, y) ∈ Ω, t ∈ [0, T ]. (3)
It represents a steady-state solution, and is referred to as lake at rest. However, that is only a special case of the moving water equilibrium. Provided we verify the compatibility condition for the bathymetry which is (-v, u) • ∇b = 0, the steady state solutions are characterized by the invariants

h(x, y, t)u(x, y, t) = Const. and E(x, y, t) = 1 2 ||u(x, y, t)|| 2 + g(h(x, y, t) + b(x, y, t)) = Const. (4)
Here, E is the specific total energy (moving water equilibrium variable), cf. [START_REF] Ricchiuto | An explicit residual based approach for shallow water flows[END_REF] for more details.

Obviously, the lake at rest (3) is a special case of (4) when the velocity reduces to zero. Because of these reasons, it is desirable to solve (2) with well-balanced schemes and also our described method has this property as we shall describe in Section 3.2.

Remark 2.1 (C-Property). As described in [START_REF] Cea | Unstructured finite volume discretisation of bed friction and convective flux in solute transport models linked to the shallow water equations[END_REF][START_REF] Ricchiuto | On the C-property and generalized C-property of residual distribution for the shallow water equations[END_REF][START_REF] Ricchiuto | An explicit residual based approach for shallow water flows[END_REF], instead of speaking of well-balanced schemes, one could alternatively say that a scheme enjoys the C-property if it preserves exactly the steady state (3). However one still speaks of C-property when referring to other steady states (4). When the conservation of the steady state is not exact but is obtained within error rates below the formal accuracy of the scheme, one often speaks of generalized C-property [START_REF] Ricchiuto | On the C-property and generalized C-property of residual distribution for the shallow water equations[END_REF].

Positivity of the solution

The other major difficulty which pops up in many simulations of the SWE is the appearance of dry regions in many real applications as for example dam break problems, flood waves and run-up phenomena at a coast with tsunamis. Here, the water height (h = 0) will be zero. As a result, all eigenvalues of the Jacobian of the flux coincide, cf. [START_REF] Ricchiuto | An explicit residual based approach for shallow water flows[END_REF] and the SWE model will be not strictly hyperbolic anymore. If, by numerical oscillations h becomes negative, the problem is also not well-posed and the calculations will simply break down. It is thus essential for a good scheme to preserve the positivity of h at any time and any point. Especially, in situations when we have dry and wet areas, the scheme has to be constructed such that it can handle these numerical challenges. In this paper, we shall target the issue of positivity of the solutions by combining a positivity preserving time-integration method together with the WENO approach. More about this will follow in Section 4.3.

Remark 2.2. Apart from the wet-dry areas in which the value of h approaches 0, we may have some regions of the space domain in which the bottom topography b overcomes the water level resulting in areas which are completely dry. In this context, the given definition (3) of the lake at rest steady state is unsatisfactory because in the mentioned regions we have η 0 < b which, if we strictly stick to (3), would give h < 0 which does not correspond to the physics of the phenomen that must be modeled. Therefore in such cases we need to modify the definition of the lake at rest steady state introducing the so-called dry lake at rest

u = v = 0; h(x, y, t) = η 0 -b(x, y), if b ≤ η 0 , 0, else, ∀ (x, y) ∈ Ω, t ∈ [0, T ]. ( 5
)
3 Space discretization: Finite Volume method

The system of PDEs considered herein will be solved by means of the Method Of Lines (MOL). Thus, space and time are going to be treated separately. This section has the goal of presenting a standard high order finite volume framework. In particular, we will consider a Cartesian setting but the method can be easily applied to a general unstructured framework in a straightforward way. The computational domain Ω is covered with N x × N y non-overlapping control volumes

Ω i,j = [x i-1/2 , x i+1/2 ] × [y j-1/2 , y j+1/2 ] (6)
with the dimensions given by ∆x = x i+1/2 -x i-1/2 and ∆y = y j+1/2 -y j-1/2 .

Considering the system of hyperbolic balance laws described by ( 1) and ( 2), for the control volume Ω i,j we can define the cell average at time t:

U i,j (t) := 1 ∆x∆y x i+1/2 x i-1/2 y j+1/2 y j-1/2 u(x, y, t) dxdy. (7)
Integrating [START_REF] Abgrall | High order schemes for hyperbolic problems using globally continuous approximation and avoiding mass matrices[END_REF] over Ω i,j provides the semi-discrete evolution formula with respect to U i,j

dU i,j (t) dt + 1 ∆x (F i+1/2,j (t) -F i-1/2,j (t)) + 1 ∆y (G i,j+1/2 (t) -G i,j-1/2 (t)) = S i,j (t) (8)
where S i,j is the cell average of the source terms over cell Ω i,j at time t S i,j (t) := 1 ∆x∆y

x i+1/2 x i-1/2 y j+1/2 y j-1/2 S(x, y, t) dxdy (9)
and F i+1/2,j and G i,j+1/2 are the cell-averages of the physical fluxes over cell boundaries at time t:

F i+1/2,j (t) = 1 ∆y y j+1/2 y j-1/2 F(u(x i+1/2 , y, t)) dy, (10)
G i,j+1/2 (t) = 1 ∆x x i+1/2 x i-1/2 G(u(x, y j+1/2 , t)) dx. (11)
Equation ( 8), along with all the previous definitions, is so far exact.

From now until the end of this section, by an abuse of notation, we will not explicitly write the dependence on t in all variables. U i,j , F i+1/2,j , G i,j+1/2 and S i,j can be then approximated to the desired order of accuracy using appropriate quadrature formulae and reconstruction techniques. From now on, we shall focus only on F i+1/2,j (G i,j+1/2 is obtained in a similar manner). Many low order and high order reconstruction techniques have been developed in the last decades. Few of them can be found in [START_REF] Godunov | A finite difference method for the computation of discontinuous solutions of the equations of fluid dynamics[END_REF][START_REF] Van Leer | Towards the ultimate conservative difference scheme. V. A second-order sequel to Godunov's method[END_REF][START_REF] Colella | The piecewise parabolic method (PPM) for gas-dynamical simulations[END_REF][START_REF] Shu | Efficient implementation of essentially non-oscillatory shock-capturing schemes[END_REF][START_REF] Gaburro | High order direct Arbitrary-Lagrangian-Eulerian schemes on moving Voronoi meshes with topology changes[END_REF][START_REF] Bassi | A high-order accurate discontinuous finite element method for the numerical solution of the compressible Navier-Stokes equations[END_REF][START_REF] Dumbser | Arbitrary high order pnpm schemes on unstructured meshes for the compressible Navier-Stokes equations[END_REF][START_REF] Núnez-De | Hybrid DG/FV schemes for magnetohydrodynamics and relativistic hydrodynamics[END_REF].

Once the reconstruction has been performed, at each face we have two sets of values of U, corresponding to x L i+1/2 and x R i+1/2 , which will be referred to as the left and right extrapolated values:

u L i+1/2,θ = u(x L i+1/2 , y θ ) , u R i+1/2,θ = u(x R i+1/2 , y θ ). ( 12
)
In our case, the weighted essentially non-oscillatory (WENO) [START_REF] Shu | Essentially non-oscillatory and weighted essentially non-oscillatory schemes for hyperbolic conservation laws[END_REF][START_REF] Shu | Efficient implementation of essentially non-oscillatory shock-capturing schemes[END_REF] reconstruction has been considered to avoid severe oscillations at discontinuities. By applying a consistent quadrature rule, and dropping the time dependence, the flux in the x-direction reads,

F i+1/2,j = 1 ∆y N θ θ=1 w θ F(u(x i+1/2 , y θ )) = 1 ∆y N θ θ=1 w θ F(u L i+1/2,θ , u R i+1/2,θ ). ( 13
)
where the subscript θ = 1, . . . , N θ corresponds to different Gaussian points y θ ∈ [y j-1/2 , y j+1/2 ] and weights w θ ∈ [0, 1]. The last step in the evaluation of the fluxes replaces F(u(x i+1/2 , y θ )) in ( 13) with a monotone and consistent numerical flux F(u L , u R ). The consistency is proved only if F(u, u) = F(u).

In our scheme we employ a Rusanov-type Riemann solver:

F(u L , u R ) = 1 2 F(u R ) + F(u L ) - 1 2 s max u R -u L , (14)
where s max is the maximum eigenvalue of the normal flux-Jacobian of the system (1).

In order to express the dependence of some quantities on several cell averages U i,j , it is useful to collect them all in the vector U.

Weighted Essentially Non-Oscillatory (WENO) method

As it is described in Shu's seminal paper [START_REF] Shu | Essentially non-oscillatory and weighted essentially non-oscillatory schemes for hyperbolic conservation laws[END_REF] the classical WENO (as well as ENO) approach contains three major steps:

1. Use the WENO reconstruction procedure which will be described in the following to obtain the values at the Gaussian points. This step involves two one-dimensional reconstructions in the two directions.

2.

Compute the numerical flux at quadrature points and integrate them at each cell interface.

3. Form a classical semidiscrete FV method and apply a time-integration method to update the cell averaged values.

Scalar reconstruction

The goal of the WENO method is to compute point-wise value of variable of interest u(x, y) at Gaussian quadrature points (x i+1/2 , y θ ), in order to have a conservative and high order accurate procedure. In general, two ways can be followed to obtain the same result: genuine multidimensional reconstruction [START_REF] Shi | A technique of treating negative weights in WENO schemes[END_REF] and dimension-by-dimension reconstruction [START_REF] Titarev | Finite-volume WENO schemes for three-dimensional conservation laws[END_REF][START_REF] Casper | A finite-volume high-order ENO scheme for two-dimensional hyperbolic systems[END_REF]. The latter is a procedure made up by successive onedimensional reconstruction sweeps and it is much simpler and less computationally expensive than the genuine multidimensional one. For this reason, we shall only focus on this one. The high order reconstructed variables we are looking for will be referred to as u L i+1/2,θ and u R i+1/2,θ . For the left values, we need to reconstruct the variable inside the cell Ω i,j , while, for the right values, similar arguments apply on the cell Ω i+1,j . We aim at reconstructing the variables with an accuracy of order p (p odd). So, we define a stencil of p cells, {Ω lx,ly , l x = i -r + 1, . . . , i + r -1, l y = j -r + 1, . . . , j + r -1}, where 2r -1 = p. For instance, WENO5 has accuracy p = 5, with r = 3, and uses a 5-cells stencil from i -2 to i + 2. In the first step of the two-dimensional reconstruction, a one-dimensional WENO reconstruction along the x-direction is performed obtaining the averages at cell interface x i+1/2 with respect to the y-direction for l y = j -r + 1, . . . , j + r -

1 v R ly =
1 ∆y

y ly +1/2 y ly -1/2 u(x R i+1/2 , y) dy, v L ly =
1 ∆y

y ly +1/2 y ly -1/2 u(x L i+1/2 , y) dy. ( 15
)
In the second sweep we perform another one-dimensional reconstruction along the y-direction in the Gaussian integration points on the y-axis (x = x i+1/2 , y = y θ ), with y θ ∈ [y j-1/2 , y j+1/2 ]. The reconstructed values can be, more generally, defined for each WENO sweep as the one-dimensional averages

q i = 1 ∆ξ ξ i+1/2
ξ i-1/2 q(ξ) dξ of a function q(ξ) where ∆ξ = ξ i+1/2 -ξ i-1/2 is the cell size. For each one-dimensional step of the procedure, there are r candidate stencils for reconstruction. For each of these stencils, made up by r cells, there is a correspoding polynomial of degree (r -1) referred as p m (ξ) m = 0, . . . , r -1. The goal of the WENO reconstruction is that of using all information coming from the r stencils employed for the reconstruction. For this reason, the WENO approach defines the reconstructed value as a convex combination of the r values of all polynomials in each quadrature point, weighted with positive nonlinear weights. The weights are chosen in order to achieve (2r -1)th order of accuracy when the solution is smooth and prefer the smoother stencils when discontinuities occur in the field. For a given (quadrature) point ξ the design of weights consists of three steps. Firstly, the optimal linear weights d m are sought so that the combination of all polynomials with these weights produces the polynomial of degree (2r -2) corresponding to the large stencil. Then, the nonlinear weights ω m can be defined as

ω m = αm r-1 k=0 α k with α m = dm (βm+ ) 2
, where is a small constant introduced to avoid division by zero (we use = 10 -6 in the simulations) and β m are the smoothness indicators

β m = r-1 k=1 ξ i+1/2 ξ i-1/2 d k dx k p m (ξ) 2 ∆ξ 2k-1 dξ , m = 0, . . . , r -1. ( 16
)
If some of d m are negative then a special procedure must be used to tackle the reconstruction problem [START_REF] Shi | A technique of treating negative weights in WENO schemes[END_REF].

The final WENO reconstructed quantity is given by q( ξ) = r-1 k=0 p k ( ξ)ω k . The numerical experiments presented herein have been performed through a piece-wise parabolic WENO5 reconstruction (r = 3), which formally corresponds to fifth order accurate approximation for smooth solutions. However, in order to actually retain the fifth-order accuracy the quadrature formulae must be consistent with the WENO reconstruction. As Titarev and Toro stated in [START_REF] Titarev | Finite-volume WENO schemes for three-dimensional conservation laws[END_REF], the best results in terms of accuracy and computational cost for r = 3 are obtained if the two-point Gaussian quadrature rule is used. However, this leads eventually to a formal fourth order of accuracy. For this reason, we implemented the four-point Gaussian quadrature rule with positive optimal weights. Up to our knowledge, the two-point Gaussian quadrature has already been thoroughly discussed in many references cited above. However, the case with 4-point Gaussian quadrature has not been fully described in literature, hence, we are going to introduce all the coefficients needed to use such formula in Appendix B and a Matlab script to compute the weights and coefficients for all orders is provided in [START_REF] Ciallella | Modified Patankar Deferred Correction WENO Code for Shallow Water Equations[END_REF]. Remark 3.1 (Positivity limiter). We aim at a positive solution and during the reconstruction procedure, it might happen that h(x L i+1/2 ) or h(x R i+1/2 ) become negative. In order to ensure that positive cell averages lead to positive reconstructions at the cell interfaces, we use the positivity limiter introduced by Perthame and Shu [START_REF] Perthame | On positivity preserving finite volume schemes for Euler equations[END_REF] and developed for two dimensional problems in [START_REF] Zhang | On positivity-preserving high order discontinuous Galerkin schemes for compressible Euler equations on rectangular meshes[END_REF]. The limiter is used in the simulation section with a parameter ε = 10 -6 as minimum water height, if not otherwise specified. We refer to [START_REF] Xing | A survey of high order schemes for the shallow water equations[END_REF][START_REF] Zhang | On positivity-preserving high order discontinuous Galerkin schemes for compressible Euler equations on rectangular meshes[END_REF] for details on the implementation. This limiter when used in combination with the forward Euler (FE) method restricts the CFL conditions to CFL FE := w Lobatto 1 the weight of the Gauss-Lobatto quadrature rule of the corresponding space accuracy. For instance, with WENO5, CFL FE = 1/12. The restriction slightly improves for high order SSPRK methods, for example we have CFL SSPRK (5,4) ≈ 1.508•CFL FE . Unfortunately, explicit SSPRK methods are at most fourth order accurate, so for fifth order schemes (as DeC5), there is no warranty that the solution stays nonnegative under any CFL condition. Instead, we highlight that the new presented approach is unconditionally positive and thus not subjected to any CFL restriction.

Well-Balanced modification of the standard Finite Volume method

In order to achieve Well-Balancing with respect to the (eventually dry) lake at rest steady state, in this work we coupled the WENO formulation with a simple modification firstly introduced in [START_REF] Berberich | High order well-balanced finite volume methods for multi-dimensional systems of hyperbolic balance laws[END_REF]. The modification consists in recasting the original problem into an equivalent one in terms of the deviation of the seeked solution U from the reference solution Ũ which must be preserved. In the particular case in which a steady solution ( ∂ Ũ ∂t = 0) must be preserved, the modification leads to the new problem

d dt U i,j + 1 ∆x (F i+1/2,j (U) -F i-1/2,j (U)) - 1 ∆x (F i+1/2,j ( Ũ) -F i-1/2,j ( Ũ))+ 1 ∆y (G i,j+1/2 (U) -G i,j-1/2 (U)) - 1 ∆y (G i,j+1/2 ( Ũ) -G i,j-1/2 ( Ũ)) = S i,j (U) -S i,j ( Ũ), (17)
which can be interpreted as a classical finite volume formulation with modified fluxes and source:

F i+1/2,j (U) = F i+1/2,j (U) -F i+1/2,j ( Ũ), G i,j+1/2 (U) = G i,j+1/2 (U) -G i,j+1/2 ( Ũ), S i,j (U) = S i,j (U) -S i,j ( Ũ). ( 18
)
This approach is very easy to code and, further, the structures related to the steady reference solution can be computed in advance once and then used for every timestep without affecting the computational time.

It must be underlined that, with this technique, all cell average computations, WENO reconstruction and source terms of the reference solution are performed following the same procedures and quadrature rules carried out for solving the balance law. Hence, all the terms always match when at the equilibrium.

Time discretization

After the description of the space discretization, we will introduce the time discretization in the this section. Due to the MOL approach, it is enough to focus here on the simple ODE case. The time discretization is one of the major point of this paper since we want to apply for the first time the arbitrary high-order, conservative and positivity preserving modified Patankar Deferred Correction method (mPDeC) together with the described WENO approach to the shallow water equations resulting in a high-order, conservative, unconditionally positivity preserving, non-oscillatory and well-balanced scheme. In agreement with [START_REF] Meister | On unconditionally positive implicit time integration for the DG scheme applied to shallow water flows[END_REF], we refer to the numerical methods which are provably positive with respect to the water height without any CFL restriction, i.e., for all time steps ∆t > 0, as unconditionally positivity preserving, in the context of the numerical solution of the shallow water equations. Before describing the combined algorithm, we introduce the Deferred Correction (DeC) method [START_REF] Dutt | Spectral Deferred Correction Methods for Ordinary Differential Equations[END_REF] as described in [START_REF] Abgrall | High order schemes for hyperbolic problems using globally continuous approximation and avoiding mass matrices[END_REF][START_REF] Abgrall | High order asymptotic preserving deferred correction implicit-explicit schemes for kinetic models[END_REF][START_REF] Torlo | Hyperbolic problems: high order methods and model order reduction[END_REF] and we repeat its modification using the Patankar trick from [START_REF] Öffner | Arbitrary high-order, conservative and positivity preserving Patankar-type deferred correction schemes[END_REF][START_REF] Öffner | Approximation and Stability Properties of Numerical Methods for Hyperbolic Conservation Laws[END_REF][START_REF] Torlo | Issues with positivity-preserving Patankar-type schemes[END_REF].

Deferred Correction method

The general DeC approach has been introduced in [START_REF] Dutt | Spectral Deferred Correction Methods for Ordinary Differential Equations[END_REF] and has been further developed and applied in [START_REF] Christlieb | Integral deferred correction methods constructed with high order Runge-Kutta integrators[END_REF][START_REF] Liu | Strong stability preserving property of the deferred correction time discretization[END_REF][START_REF] Minion | Semi-implicit spectral deferred correction methods for ordinary differential equations[END_REF] in different contexts, whereas a simplified version was presented in [START_REF] Abgrall | High order schemes for hyperbolic problems using globally continuous approximation and avoiding mass matrices[END_REF]. In [START_REF] Abgrall | High order schemes for hyperbolic problems using globally continuous approximation and avoiding mass matrices[END_REF], a compact operator notation was also introduced and we shall follow this framework herein. However, even if the notation changed, the main idea of DeC is always the same. DeC is based on the Picard-Lindelöf Theorem in the continuous setting and the classical proof makes use of Picard iterations to minimize the error and to obtain convergence. DeC is constructed to mimic these Picard iterations at the discrete level and decreases the approximation error in several iterative steps. To explain the method, we consider the following time-dependent initial value problem

y (t) = f (y(t)), y(t 0 ) = y 0 , (19)
where y : R → R S and f : R S → R S resulting from our MOL approach.

For our description, two operators are introduced: L 1 and L 2 . Here, the L 1 operator represents a low-order easy-to-solve numerical scheme, e.g. the explicit Euler method, and L 2 is a high order operator that can present difficulties in its practical solution, e.g. an implicit RK scheme. The DeC method can be written as a combination of these two operators. Given a time interval [t n , t n+1 ], we subdivide it into M subintervals

t n = t n,0 = t 0 y 0 t n,1 = t 1 y 1 t n,m = t m y m t n,M = t M = t n+1 y M Figure 2: Time interval divided into subintervals {[t n,m-1 , t n,m ]} M m=1
, where t n,0 = t n and t n,M = t n+1 . Therefore, we mimic for every subinterval [t 0 , t m ] the Picard-Lindelöf Theorem for both operators L 1 and L 2 . With an abuse of notation, we drop the dependency on the timestep n for subtimesteps t n,m and substates y n,m as denoted in Figure 2. Then, the L 2 operator is given by

L 2 (y 0 , . . . , y M ) :=        y M -y 0 -∆t M r=0 θ M r f (y r ) . . . y 1 -y 0 -∆t M r=0 θ 1 r f (y r ) . ( 20
)
where θ m r are the weights of a high order quadrature rule in {t m } M m=0 . The L 2 operator represents a high order numerical scheme if set equal to zero, i. e., L 2 (y 0 , . . . , y M ) = 0. The order depends on the distribution of the subtimesteps, for instance, with M equispaced subtimesteps one obtains (M + 1)th order, while with M Gauss-Lobatto quadrature subtimesteps one has (2M )th order. Unfortunately, the resulting scheme is implicit and, further, the terms f may be nonlinear. The L 1 operator is given by the forward Euler discretization for each state y m in the time interval, i. e.,

L 1 (y 0 , . . . , y M ) :=        y M -y 0 -β M ∆tf (y 0 ) . . . y 1 -y 0 -β 1 ∆tf (y 0 ) (21)
with coefficients β m := t m -t 0 t M -t 0 . To simplify the notation and to describe DeC, we introduce the matrix of states for the variable y at all subtimesteps.

y := (y 0 , . . . , y M ) ∈ R (M +1)×S , such that (22)
L 1 (y) := L 1 (y 0 , . . . , y M ) and L 2 (y) := L 2 (y 0 , . . . , y M ).

The DeC algorithm uses a combination of the L 1 and L 2 operators to provide an iterative procedure. The aim is to recursively approximate y * , the numerical solution of the L 2 (y * ) = 0 scheme, similarly to the Picard iterations in the continuous setting. The successive states of the iteration process will be denoted by the superscript (k), where k is the iteration index, e.g. y (k) ∈ R (M +1)×S . The total number of iterations (also called correction steps in the following) is denoted by K. To describe the procedure, we have to refer to both the m-th subtimestep and the k-th iteration of the DeC algorithm. We will indicate the variable by y m,(k) ∈ R S . Finally, the DeC method can be written as DeC Algorithm

y 0,(k) := y(t n ), k = 0, . . . , K, y m,(0) := y(t n ), m = 1, . . . , M, L 1 (y (k) ) = L 1 (y (k-1) ) -L 2 (y (k-1) ) with k = 1, . . . , K, (24)
where K is the number of iterations that we want to compute.

Using the procedure (24), we need, in particular, as many iterations as the desired order of accuracy p, i. e., K = p. This means that we choose the number of subtimesteps in a way that the order of the L 2 operator is itself equal to p. In practice, for each correction and each subtimestep, L 1,m (y

(k) ) = L 1,m (y (k-1) ) - L 2,m (y (k) ) reduces to solve y m,(k) α -y 0 α -∆t M r=0 θ m r f α (y r,(k-1) ) = 0, ∀α = 1, . . . , I. (25)
For more information and properties of the DeC approach, we refer to [START_REF] Abgrall | Relaxation deferred correction methods and their applications to residual distribution schemes[END_REF][START_REF] Veiga | DeC and ADER: Similarities, Differences and a Unified Framework[END_REF] and references therein. In the following, we explain how to adapt the DeC approach to obtain a conservative and positivity preserving time integration scheme.

Patankar method for production-destruction systems

Many problems [START_REF] Colella | The piecewise parabolic method (PPM) for gas-dynamical simulations[END_REF] in nature can be written as a production destruction system (PDS) for the unknown

y ∈ R S f α (y) = S β=1 (p α,β (y) -d α,β (y)), (26)
where p α,β , d α,β ≥ 0 are the production and destruction terms, respectively. The production and destruction terms are conveniently written as matrices. Applications for PDS are for example the biological and/or chemical reactions such as algal bloom [START_REF] Burchard | A high-order conservative Patankar-type discretisation for stiff systems of production-destruction equations[END_REF]. Also parts (or all) of the semi discretization of hyperbolic conservation/balance laws can be interpreted in such PDS system as described in [START_REF] Huang | Positivity-preserving time discretizations for production-destruction equations with applications to non-equilibrium flows[END_REF][START_REF] Huang | A third-order unconditionally positivity-preserving scheme for productiondestruction equations with applications to non-equilibrium flows[END_REF][START_REF] Meister | On unconditionally positive implicit time integration for the DG scheme applied to shallow water flows[END_REF] and also later in this work. The calculated solutions are often describing physical quantities that enjoy some properties, for instance concentrations of chemicals or water height in the context of SWE should be nonnegative. The following definition may be introduced for ODE systems:

Definition 4.1. An ODE (19) is called positive, if positive initial data y 0 > 0 result in positive solutions y(t) > 0, ∀t.
Here, inequalities for vectors are interpreted componentwise, i.e., y(t) > 0 means ∀α :

y α (t) > 0. A PDS (26) is conservative, if p α,β (y) = d β,α (y), ∀ α, β, y.
These properties should be preserved by the numerical scheme as well. Thus, we introduce the following discrete counterpart.

Definition 4.2. A numerical method computing y n+1 ≈ y(t n+1 ) given y n ≈ y(t n ) is called conservative, if α y n+1 α = α y n α . It is called unconditionally positive, if y n > 0 implies y n+1 > 0 for any time step ∆t > 0.
From literature [START_REF] Burchard | A high-order conservative Patankar-type discretisation for stiff systems of production-destruction equations[END_REF], it is well-known that the implicit Euler method is conservative and unconditionally positive preserving whereas the explicit Euler method is only conservative (it might be positive under time step restrictions). To avoid solving a fully nonlinear system of equations, the so-called Patankar modifications have been applied to the explicit Euler method. To build an unconditionally positive numerical scheme, Patankar had the idea [START_REF] Patankar | Numerical heat transfer and fluid flow[END_REF] of firstly weighting the destruction term in the original explicit Euler method with a coefficient as follows

y n+1 α = y n α + ∆t   S β=1 p α,β (y n ) - S β=1 d α,β (y n ) y n+1 α y n α   , α = 1, . . . , S. (27)
Indeed, the resulting scheme ( 27) is unconditionally positive and the implicit terms can be collected on the left hand side, but the conservation relation is violated. Burchard et al. had the idea [START_REF] Burchard | A high-order conservative Patankar-type discretisation for stiff systems of production-destruction equations[END_REF] not only to weight the destruction term but also the production term:

y n+1 α = y n α + ∆t   S β=1 p α,β (y n ) y n+1 β y n β - S β=1 d α,β (y n ) y n+1 α y n α   , α = 1, . . . , S. (28)
They called their constructed scheme (28) modified Patankar scheme and proved that it is is unconditionally positive and conservative. The resulting scheme is linearly implicit, meaning that collecting all the implicit terms on the left hand side, we obtain a linear system at each time iteration. Based on this technique, extensions to second and third order modified Patankar Runge-Kutta (MPRK) methods have been made by several researchers in such context, cf. [START_REF] Huang | Positivity-preserving time discretizations for production-destruction equations with applications to non-equilibrium flows[END_REF][START_REF] Huang | A third-order unconditionally positivity-preserving scheme for productiondestruction equations with applications to non-equilibrium flows[END_REF][START_REF] Kopecz | Unconditionally positive and conservative third order modified Patankar-Runge-Kutta discretizations of production-destruction systems[END_REF][START_REF] Kopecz | On the existence of three-stage third-order modified Patankar-Runge-Kutta schemes[END_REF]. Also the semi implicit RK methods proposed in [START_REF] Chertock | Steady state and sign preserving semi-implicit Runge-Kutta methods for ODEs with stiff damping term[END_REF] can be interpreted as Patankar methods as they weight only the destruction terms [START_REF] Torlo | Issues with positivity-preserving Patankar-type schemes[END_REF]. Finally, in [START_REF] Öffner | Arbitrary high-order, conservative and positivity preserving Patankar-type deferred correction schemes[END_REF] an arbitrarily high-order, conservative and positivity preserving scheme based on the DeC framework has been constructed. Herein we describe the main idea. For the details of the properties and proofs, we refer again to [START_REF] Öffner | Arbitrary high-order, conservative and positivity preserving Patankar-type deferred correction schemes[END_REF].

Modified Patankar Deferred Correction method

The modified Patankar Deferred Correction (mPDeC) is based on the DeC algorithm [START_REF] Gaburro | High order direct Arbitrary-Lagrangian-Eulerian schemes on moving Voronoi meshes with topology changes[END_REF] and it consists in a modification of the L 2 operator through the modified Patankar trick. This amounts to weight the production-destruction terms with respect to the intermediate approximations.

Using the fact that initial states y 0,(k) α are identical for any correction k, the mPDeC correction steps can be rewritten [START_REF] Öffner | Arbitrary high-order, conservative and positivity preserving Patankar-type deferred correction schemes[END_REF] for k = 1, . . . , K, m = 1, . . . , M and ∀α = 1, . . . , S as

y m,(k) α -y 0 α - M r=0 θ m r ∆t S β=1   p α,β (y r,(k-1) ) y m,(k) γ(β,α,θ m r ) y m,(k-1) γ(β,α,θ m r ) -d α,β (y r,(k-1) ) y m,(k) γ(α,β,θ m r ) y m,(k-1) γ(α,β,θ m r )   = 0, (29)
where θ m r are the DeC quadrature weights in time and

γ(α, β, θ) := α if θ ≥ 0, β if θ < 0.
Finally, the new numerical solution is y n+1 = y M,(K) . As in the classical DeC framework, the choice of the distribution, the number of subtimesteps M and the number of iterations K determines the order of accuracy of the scheme. To reach order p, classically M = p -1 equispaced subintervals and K = p corrections should be used. As proven in [START_REF] Öffner | Arbitrary high-order, conservative and positivity preserving Patankar-type deferred correction schemes[END_REF], the scheme is conservative, positivity preserving and can reach arbitrary high order.

A description of the assembly of the mass matrix of the system (29) is described in Algorithm 1 and one timestep of the mPDeC algorithm is sketched in Algorithm [START_REF] Gaburro | A unified framework for the solution of hyperbolic PDE systems using high order direct Arbitrary-Lagrangian-Eulerian schemes on moving unstructured meshes with topology change[END_REF] where the evolution formula is given by Algorithm 2.

Remark 4.3 (Subtimestep distribution).

In our numerical simulations, we apply Gauss-Lobatto nodes in every timestep. They have the advantage of requiring less subtimesteps to reach pth order of accuracy. In the following we will use M = 3 Gauss-Lobatto subtimesteps, which guarantee 6th order of accuracy for the operator L 2 and K = 5 iterations aiming at a 5th order scheme to match the spatial discretization accuracy of WENO5.

Remark 4.4 (Solution of the linear system). At each subtimestep m and iteration (k) we need to solve the linear system given by [START_REF] Veiga | DeC and ADER: Similarities, Differences and a Unified Framework[END_REF]. The mass matrix obtained has the following form:

M(y m,(k-1) ) α,β =            1 + ∆t M r=0 S β=1 θ m r y m,(k-1) α d α,β (y r,(k-1) )1 {θ m r >0} -p α,β (y r,(k-1) )1 {θ m r <0} , for α = β, -∆t M r=0 θ m r y m,(k-1) β p α,β (y r,(k-1) )1 {θ m r >0} -d α,β (y r,(k-1) )1 {θ m r <0} , for α = β, ( 30
)
where 1 is the indicator function. The mass matrix assembly algorithm is described in Algorithm 1. The linear system will then read

M(y m,(k-1) )y m,(k) = y(t n ). ( 31
)
Remark 4.5 (Division on almost wet areas). When the water height is low, we might encounter troubles in computing the divisions in [START_REF] Huang | Positivity-preserving time discretizations for production-destruction equations with applications to non-equilibrium flows[END_REF] as the denominator might be very small. They hypothesis behind the production and destruction system that says that as h α → 0 also d α → 0, can be difficult to be obtain at a numerical level. Hence, to be sure that those divisions do not lead to extremely high values when they should go to 0, we slightly modify the way we implement the division as suggested in [START_REF] Meister | A positivity preserving and well-balanced DG scheme using finite volume subcells in almost dry regions[END_REF]. Given any numerator n and denominator d of (30), we approximate the division by

n d ≈ 0 d < ε, 2d•n d 2 +max{d 2 ,ε} d ≥ ε, (32)
with ε a small tolerance value. Along the computations, if not specified, we will use ε := 10 -6 . This formulation allow to smoothly pass from n d to 0 as d → 0. Moreover, when d 2 ≥ ε the division will be exact. In the following, we apply the mPDeC time marching algorithm [START_REF] Veiga | DeC and ADER: Similarities, Differences and a Unified Framework[END_REF] to the WENO finite volume semidiscretization to solve the shallow water equations. Below we describe the actual implementation procedure.

Implementation of well-balanced mPDeC-WENO Scheme

In the following part, we will describe how the semi-discretization, using the WENO approach, can be written and interpreted as a PDS in order to apply the mPDeC scheme. We will see that only small modifications in the WENO code are necessary. Furthermore, this interpretation is not only restricted to the WENO procedure but can also applied to other high-order FV/FD discretizations. Our WENO approach is working only as a generic example. Finally, the complete algorithm will be presented.

WENO and production-destruction systems

First of all we must underline the fact that in order to preserve the positivity of the water height h, the mPDeC scheme is going to be applied only to the first equation of system (1) as h is the only variable that must stay nonnegative. Thus, a simple DeC approach is going to be used to evolve the momentum equations. So from now on, we shall only talk about the modifications introduced for the first equation to turn it into a production-destruction system. Given the foundations of finite volume schemes, each control volume has fluxes entering and exiting its boundary and, for each boundary face, the flux going from element α = [i, j] to element β = [l, r] is going to be equal in module and opposite in sign to the flux from element β to element 1: Compute the mass matrix M(y m,(k-1) ) ←Mass(p α,β (•), d α,β (•),∆t, y (k-1) , m) using Algorithm 1 2: Compute y m,(k) solving the linear system M(y m,(k-1) )y m,(k) = y n given by [START_REF] Veiga | DeC and ADER: Similarities, Differences and a Unified Framework[END_REF] with Jacobi Algorithm 3 3: return y m,(k) Algorithm 3 Jacobi iterative method Require: D diagonal of the matrix, L off-diagonal terms of the matrix, r right hand side of the system, tol tolerance.

1: err ← 2 • tol, k ← 0, x k ← r 2: while err > tol do 3: k ← k + 1 4: x k+1 ← D -1 (r -Lx k ) 5: err ← ||x k -x k-1 || 6: end while 7: return x k+1 Ω i-1,j Ω i,j Ω i+1,j Ω i,j+1 Ω i,j-1 F i-1/2,j (U) F i+1/2,j (U) G i,j-1/2 (U) G i,j+1/2 (U) Ω i-1,j Ω i,j Ω i+1,j Ω i,j+1 Ω i,j-1 p [i,j],[i-1,j] d [i,j],[i-1,j] p [i,j],[i+1,j] d [i,j],[i+1,j] p [i,j],[i,j-1] d [i,j],[i,j-1] d [i,j],[i,j+1] p [i,j],[i,j+1]
Figure 3: Cell Ωi,j, its four neighbors and its production and destruction terms.

α. This is the key feature for turning a finite volume schemes into a PDS. Note that the source is zero on h component. Therefore, we can define the production and destruction terms for a general α = [i, j] and the neighboring β = [l, r] as

p [i,j],[i-1,j] (U) = + 1 ∆x F i-1/2,j (U) + , d [i,j],[i-1,j] (U) = - 1 ∆x F i-1/2,j (U) -, p [i,j],[i+1,j] (U) = - 1 ∆x F i+1/2,j (U) -, d [i,j],[i+1,j] (U) = + 1 ∆x F i+1/2,j (U) + , p [i,j],[i,j-1] (U) = + 1 ∆y G i,j-1/2 (U) + , d [i,j],[i,j-1] (U) = - 1 ∆y G i,j-1/2 (U) -, p [i,j],[i,j+1] (U) = - 1 ∆y G i,j+1/2 (U) -, d [i,j],[i,j+1] (U) = + 1 ∆y G i,j+1/2 (U) + , (33)
where with the superscript + and -we denote the positive and the negative part respectively. All the other p α,β and d α,β not defined here are set to 0. Clearly, this define a conservative and positive PDS, as the properties in Definition 4.1 are verified. The visualization of the production and destruction terms in Figure 3 may help the reader. We clearly observe that the matrices (p α,β ) and (d α,β ) are S × S sparse matrices, with, at most, 4 nonzero entries per row and S = N x • N y .

Once the production and destruction matrices have been assembled, the next step consist in running the mPDeC algorithm in Eq. [START_REF] Veiga | DeC and ADER: Similarities, Differences and a Unified Framework[END_REF]. Note that the matrix built in ( 30) is sparse as well with at most 5 nonzero entries for each row (4 nonzero entries of the production/destruction terms and the diagonal term). Hence, in the numerical computations we will use the classical Jacobi iterative method, see Algorithm 3, to obtain the solution of system (29) at each iteration. Indeed, it is provable [START_REF] Öffner | Arbitrary high-order, conservative and positivity preserving Patankar-type deferred correction schemes[END_REF][START_REF] Burchard | A high-order conservative Patankar-type discretisation for stiff systems of production-destruction equations[END_REF] that the Jacobi iteration algorithm converges when applied on the matrix defined in [START_REF] Huang | Positivity-preserving time discretizations for production-destruction equations with applications to non-equilibrium flows[END_REF]. In all calculations, we set the tolerance to machine precision and the algorithm converge towards the solution in few iterations. Experimentally, we have seen that usually 10-20 iterations suffice, in the worst cases 40 iterations are needed and the number of iterations do not depend much on the mesh size. Overall, considering the assembly the mass matrix, the inversion of the system with Jacobi iterations, the whole mPDeC procedure increases the computational costs of around 18% with respect to the original DeC algorithm using the same CFL. The code has not been construct with the goal of optimizing all the procedures, so it might well be that this extra computational cost can be decreased with better implementations. The greatest advantage is anyway that no CFL restrictions are required in order to guarantee the positivity of the solution. So, a CFL= 1 suffices to guarantee the stability of the method. In the numerical test section, we will also detail the number of Jacobi iterations and the overhead computational cost that the method brings into the system with respect to the explicit method.

Remark 5.1 (Efficiency and properties of Jacobi iterative method). Though being a very simple algorithm, Jacobi iterative method is particularly effective for this application. Indeed, the mass matrix in R S×S is very sparse, i.e., only 5S non-zero elements, hence, at each iterations, only 5S multiplications are computed. Moreover, experimentally, we have seen that the overhead computational cost of the Jacobi solver is around 10% of the whole computational cost and that the iterations needed are usually below 20, see section 6. Moreover, Jacobi guarantees the positivity of the solution at every iteration of the procedure.

Clearly, there are other many iterative methods to solve linear system [START_REF] Saad | Iterative methods for sparse linear systems[END_REF], e.g. Krylov preconditioning methods. All of these methods have larger complexity per iteration, but may converge faster to the solution of the system. Nevertheless, for most of these methods it is not possible to guarantee the positivity of the solution of the iterative solver along the process. That is why, we will use the Jacobi iterative method.

Full Algorithm

After describing how the WENO (FV) procedure can be interpreted as a PDS, we give a more precise description of the full algorithm which is used to calculate the numerical solution.

In our version of this approach, we have to adapt the steps taking into account the re-interpretation of the WENO approach as a production destruction system and the well balancing approach and this is described in Algorithm 4. Finally, we can combine all the ingredients described above in a full algorithm as in Algorithm 5. There, we simply use the mPDeC to evolve in time and the production and destruction functions are given by the WENO description from above. U n+1 =DeC(U n , ∆t, WENOPDS) with DeC Algorithm [START_REF] Gaburro | A unified framework for the solution of hyperbolic PDE systems using high order direct Arbitrary-Lagrangian-Eulerian schemes on moving unstructured meshes with topology change[END_REF] where update formula ( 29) is used for h and ( 25) is used for hu and hv and the WENO PSD function are given by Algorithm 4 5: t = t + ∆t 6: end while Conclusion 5.2. The described method is high order accurate, positivity preserving, conservative, nonoscillatory and well-balanced for the shallow water equation.

Remark 5.3 (Difference with respect to classical WENO).

We want to highlight the differences between a classical WENO and the proposed algorithm are minimal. Indeed, once the spatial discretization is performed with a simple WENO step we need to apply two easy modifications. The first one consists in subtracting the flux related to the steady state variables. The second one consists in the definition of the production and destruction terms. Then the mPDeC can be applied as a simple time integration scheme. The code and these modifications are available at the reproducibility repository [START_REF] Ciallella | Modified Patankar Deferred Correction WENO Code for Shallow Water Equations[END_REF]. Remark 5.4 (Advantages of mPDeC). We shall remark that the presented method does not require any CFL constraint to obtain positive solutions for h, while the classical positivity limiter for WENO5 requires a CFL number of 1/12 ≈ 0.083. Clearly, a CFL number for a classical explicit method must be anyway used (between 1 and 1.5), but this allow to run the simulation with much less time steps than a classical explicit WENO scheme, with the extra computational cost of the Jacobi iterative method, which is negligible with respect the cost of decreasing the CFL number of factor of 12. Moreover, the procedure allows to have a provably positive method with arbitrarily high order of accuracy, while with positive limiters applied to WENO schemes only SSPRK methods guarantee the positivity of the solutions and they exist only up to order 4 [START_REF] Gottlieb | Strong stability preserving Runge-Kutta and multistep time discretizations[END_REF].

Remark 5.5 (Stability and accuracy of combined time-integration methods). We combined two timeintegration methods, the mPDeC and the DeC approach. How can we assess the stability and accuracy properties the combined method. For the accuracy of the combined scheme, we can exclude a loss of accuracy as long as the two methods are consistent up to the same order of accuracy, even if applied separately to each equation. For instance, one could perform a simple Taylor expansion to prove the error behavior for each component and also in the combination step-after-step. In practice, as done in [START_REF] Öffner | Arbitrary high-order, conservative and positivity preserving Patankar-type deferred correction schemes[END_REF], it is sufficient to develop the Taylor expansion on the modified Patankar weighting coefficients to obtain a high order approximation of the DeC. This is also verified in our numerical tests. Stability is more problematic since it is not clear how to analyze it even for modified Patankar schemes in the ODE case, where only few preliminary works are available [START_REF] Izgin | On Lyapunov stability of positive and conservative time integrators and application to second order modified Patankar-Runge-Kutta schemes[END_REF][START_REF] Torlo | Issues with positivity-preserving Patankar-type schemes[END_REF]. This is also due to the nonlinear character of the scheme itself. Therefore, a stability investigation for PDEs is far beyond the goals of this paper. However, in general explicit methods have a more restrictive stability region than implicit methods. Therefore, we will simple assume that our stability region is determined by the underlying explicit method (DeC).

Numerical Simulations

The goal of this section is to present the results obtained with the fifth order positivity-preserving mPDeC scheme, compared to that given by the classical fifth order DeC time integration method. The first test case consists in assessing the convergence properties of the spatial and temporal discretization on an unsteady vortex-type solution [START_REF] Ricchiuto | Analytical travelling vortex solutions of hyperbolic equations for validating very high order schemes[END_REF]. Afterwards, we focus on testing the well-balanced implementation for the lake at rest solution by showing its impact on perturbation analysis. Finally, three challenging simulations are performed to prove its capabilities to cope with wet-dry fronts. It is important to underline that in all simulations, especially those involving shocks, when there is no comparison with DeC5 (dam break problems), the mPDeC5 method obtains consistent results with respect to classical methods. For all the simulations carried out herein periodic boundary conditions have been considered together with the local Lax-Friedrichs (Rusanov) numerical flux.

Unsteady vortex

In order to verify the order of accuracy we consider a moving smooth vortex. The computational domain is the square [0, 3] × [0, 3]. The initial condition is given by some perturbations δ applied on a homogeneous background field (h 0 , u 0 , v 0 ) = (1, 2, 3). Hence, the perturbation for the depth variable h is

h(r) = h 0 -δh(r) = h 0 -γ e - 1 arctan 3 (1-r 2 ) , if r < 1, 0, else, with r = (x -1.5) 2 + (y -1.5) 2 (34)
and the vortex amplitude is γ = 0.1. The velocity field is affected by the following perturbations where ∂ r h = ∂ r h(r) is a function of the radial distance from the center of the vortex

δu δv = 2 g ∂ r h (y -1.5) -(x -1.5) , (35)
∂ r h(r) = 3 γ e - 1 arctan 3 (1-r 2 ) arctan 4 (r 2 -1)((r 2 -1) 2 + 1) . (36)
It is important to highlight the fact that this solution is C ∞ , which is a fundamental property for testing arbitrarily high order schemes. Many vortex-type solutions can be found available online but most of them can only be used to test lower order schemes. The exact solution of this problem is given by h(x, y, t) = h(x -u 0 t, y -v 0 t, 0), u(x, y, t) = u(x -u 0 t, y -v 0 t, 0), v(x, y, t) = v(x -u 0 t, y -v 0 t, 0). [START_REF] Lukáčová-Medvidová | Well-balanced finite volume evolution Galerkin methods for the shallow water equations[END_REF] For the unsteady vortex, two convergence tests are run for WENO5 coupled with both the time integration schemes DeC5 and mPDeC to corroborate the fact that, for smooth flows, the results should be almost identical. We used CFL=0.7. The convergence tests are run on cartesian meshes of size 25 × 25, 50 × 50, 100 × 100, 200 × 200, 300 × 300, 400 × 400, 500 × 500 and 600 × 600. The computational mesh of the coarsest grid and initial condition for this test case are shown in Figure 4. For these convergence tests we used a tolerance ε = 10 -30 both for the positivity limiter and the mPDeC divisions, since the errors that we obtain are of the order of 10 -8 . The error || h (u)|| is the L 1 norm of the difference between the exact solution and the approximated one. Figure 5 points out the predicted fifth order behavior for both time integration schemes.

In Figure 6 we study the computational costs of the two methods (mPDeC5 and DeC5) with respect to different CFL numbers and mesh refinements. In Figure 6(a) we plot the ratio of computational time of mPDeC5 over the computational time of DeC5 needed to finish the simulation for the same CFL and mesh. We see that for fine meshes, when the computational times are more reliable, all ratios are close to 1.1. This means that the overhead that mPDeC5 requires, with respect to an explicit method, is of about 10%. In Figure 6(b) we plot errors with respect to computational times and we see a very small difference between mPDeC5 and DeC5, while there is a huge difference in computational costs when changing the CFL. The mPDeC5 is guaranteed to run at any CFL FE < 1, while the positivity for other SSPRK methods is guaranteed only for CFL FE < 1/12 with WENO5. Hence, there is huge advantage with mPDeC. Finally, in Figure 6(c) we can observe on average how many iterations are needed to solve the linear system and a confidence interval defined by the average ± 1 2 standard deviation. It is clear that Jacobi iterations are mainly driven by the CFL number, which explicitly appear as a coefficient of the mass matrix minus the identity. Indeed, all factors ∆t d h or ∆t p h are linearly dependent on the CFL as production and destruction terms are proportional to 1 ∆x .

Lake at rest

As already introduced in the theoretical part, this test case is needed to prove the presented scheme is well-balanced. The computational domain is the square [0, 1] × [0, 1] and the steady solution of this problem and bathymetry are briefly summarized below b(x, y) = 0.1 sin (2 π x) cos (2 π y) , h(x, y) = 1 -b(x, y),

u = v = 0. ( 38
)
This benchmark can also be used to test once again the order of accuracy of our discretization. Indeed, we expect the method to converge with a fifth order slope when not well-balanced and we expect machine precision errors for all the well-balanced tests. This simulation has been performed with four different settings: DeC5 and mPDeC5, well-balanced and not well-balanced. For all cases, we employed a fifth order WENO discretization for the spatial derivatives and CFL=0.9. Also for this test, we chose ε = 10 -30 to check the accuracy of the scheme with errors of the order of 10 -10 . As expected, the error computed for the well-balanced simulations is exactly zero therefore we did not plot them along with the other results. The interested reader can run the simulations and test the properties of our method by downloading the code available at the reproducibility repository [START_REF] Ciallella | Modified Patankar Deferred Correction WENO Code for Shallow Water Equations[END_REF]. The Cartesian mesh employed for this convergence test are 4 × 4, 8 × 8, 16 × 16, 32 × 32, 64 × 64 and 128 × 128. The exact solution is presented in Figure 7, along with the 32 × 32 mesh. As can be noticed from Figure 8, mPDeC5 allows a fifth order convergence rate as theoretically proved with results almost identical to those of DeC5.

Wet-dry lake at rest

Now we test the capability of dealing with wet and dry regions of the scheme in a very simple context. We consider a bathymetry given by a bump b(x, y) = e

1-1 1-r 2 , if r 2 < 1, 0,
else , where . The lake at rest solution is

r 2 = x 2 + y 2 , ( 39
h(x, y) = max{0.7 -b(x, y), 0}, u = v = 0. ( 40
)
The maximum of the bathymetry is 1, hence, there is a dry island at the center of the domain. For practical purposed, we set the initial conditions to be

h 0 (x, y) = max{0.7 -b(x, y), ε}, u = v = 0, (41)
with ε = 10 -6 . We consider final time T = 1. First, we test the non-well-balanced schemes, to asses the capability of preserving the water height positivity and the accuracy of such methods. The positivity of the classical schemes is not preserved, even with the positivity limiter of remark 3.1. In the dry region, the SW model (2) does not hold and for the time integration schemes it is hard to verify hypotheses that guarantee the positivity of the solution. Hence, for the classical schemes, we force the positivity of the water height every time we need to compute the flux or to convert the variables from conservative to primitive and vice versa. On the other hand, the mPDeC scheme always preserve the positivity of the solution and none of these tricks is required.

In figure 9 we observe that for similar computational times, the mPDeC5-WENO5 gives the same accuracy of all classical schemes. For all schemes the accuracy is 2, as the solution is not C 1 everywhere. The difference is that mPDeC5-WENO5 can be run up to CFL=1, without any problem, while for CFL=1 the Adding the well-balanced technique we obtain machine precision errors for all schemes.

Almost dry lake at rest

Now we modify the previous test, in order to have a smooth solution and to be able to obtain a fifth order accuracy in the schemes. We consider again the bathymetry (39) on the domain [-5, 5] × [-2, 2]. The lake at rest solution is defined, this time, as h(x, y) = max{0.999 -b(x, y), 0},

u = v = 0. ( 42
)
Notice that the bathymetry has a peak with value 1, but, depending on the mesh refinement, the peak will be lower. In most of the simulations, this test will be completely wet and C ∞ . If we discretize the mesh with more than 600 × 180 elements, the water level will be below the threshold ε = 10 -6 . As before, we initialize the water heigth at least equal to ε = 10 -6 and we let the schemes evolve up to final time T = 1. We compare the mPDeC5-WENO5 and the SSPRK(6,4)-WENO5 schemes. The test is steady, hence, we expect only the spatial discretization to be the only responsible of the order of accuracy. For the SSPRK(6,4)-WENO5 to be run, we need to introduce some extra checks on the water height and computations of the flux, so that it does not become negative, while for the mPDeC5 time integration we do not need this type of extra corrections. In both cases we expect to have a small perturbation of the solution while computing the L 1 error, indeed, the initial condition set with minimum level at 10 -6 should introduce an error, in very few cells, of this order.

In table 1 there is the error analysis for the mPDeC5-WENO5 method, while in table 2 there is the one referred to the SSPRK(6,4)-WENO5 method. Despite expecting the error of the initialization to become evident around error of 10 -6 , for the mPDeC5 simulation, this small perturbation confined to very few cells does not propagate much and lead to very accurate results also for errors ≈ 5 • 10 -7 . Even the correction in remark 4.5 does not seem to affect the accuracy of the solution, probably because the water height never reaches values much lower than 10 -6 . So, the order of accuracy stays very close to five, the expected one, even for almost dry solutions. On the other side, in the SSPRK simulation the need of extra corrections in the flux every time the solution is below 10 -6 adds further errors that are visible at level of L 1 error around 10 -6 . Indeed, it seems that the error of the water height starts plateauing close to that value. And there it loses the expected fifth order of accuracy.

Perturbation of the lake at rest

In order to better highlight the improvements one gets with the well-balanced implementation, a perturbation analysis is run on a problem with both wet and dry areas. This test case is run on the rectangular domain [-5, 5] × [-2, 2]. The bathymetry is given by [START_REF] Meister | On unconditionally positive implicit time integration for the DG scheme applied to shallow water flows[END_REF]. The lake at rest solution that we are going to conserve with the well-balanced implementation is

h(x, y) = max{0.7 -b(x, y), ε}, u = v = 0, (43)
with ε = 10 -6 , whereas the perturbation shape that we want to study is

h(x, y) = h(x, y) + 0.05 e 1- 1 (1-ρ 2 ) 2 , if ρ 2 < 1, 0, else , where ρ 2 = 9((x + 2) 2 + (x -0.5) 2 ). (44)
Two simulations have been performed: one with the well-balanced correction and one without. Both simulations use WENO5 as spatial discretization and mPDeC5 for integrating the ODEs coming from the semi-discrete system. From this test on the tolerances of the positivity limiter and of the mPDeC5 divisions is set to ε = 10 -6 . The computational mesh and bathymetry plot are shown in Figure 10.

The results obtained for the two implementations are displayed in Figure 11 where only the isolines of the water height h are shown. Four snapshots are presented at different times of the simulation, t = 0, 0.25, 0.5, 1.

The results on the right-hand side of Figure 11 are those computed without the well-balancing correction. As can be noticed, in this case, the numerical error propagates from the nonconstant bathymetry area around the island placed in the middle of the domain. This error propagates and interacts with the perturbation, making the perturbation waves indistinguishable from the noise. This test case allows to better assess the well-balanced implementation already tested in the previous test case. Indeed, for this case, we have a dry area which involves a jump in the derivative of the water height causing a reduction of the order of accuracy given by the WENO method, whose limiters work with the high order derivatives of the solution. On the other side, the simulation runs with the well-balanced correction allows to exactly preserve the lake at rest solution over which the perturbation analysis is carried out. This leads to a much better capturing of the perturbation, whose evolution is not influenced by the spurious disturbances coming from the wet-dry area.

Circular dry dam break problem

We simulate the break of a circular dam separating two basins with water heights h 1 = 2.5 and h 2 = ε = 10 -6 , meaning that the water in the first basin is falling over a dry area which is all around it. The radius of the discontinuity is r = 7. A sketch of the initial condition, along with the computational mesh, is given in Figure 12. The computational domain is the square [0, 40] × [0, 40] discretized with 100 × 100 cells and the simulation is run until a final time t end = 0.9. This is a somewhat challenging test for the mPDeC5-WENO5 method that has to face both the capture of a sharp discontinuity and the progressing wetting of a dry area while always maintaining its appealing properties. The results are printed for different times, t = 0, 0.3, 0.6, 0.9, in Figure 13 showing the evolution of the water height. The advantages of the mPDeC5 have been clearly proven by running this challenging test case with different CFL conditions. As a matter of fact, we managed to run this test case with a CFL up until 1.5 while ensuring positivity. In particular, the water height equation does not cause any CFL restriction due to the nature of the method used to solve it, which is unconditionally positive. However, since the momentum equations are solved by means of an explicit DeC scheme, the CFL must be bounded. It should be noticed that, in order to retain positivity in the spatial reconstruction, a positivity limiter has been implemented. For explicit SSPRK methods, this limiter has been proven to cause a huge restriction in the CFL condition, which now has to be less than 1 12 CFL SSPRK . Nevertheless, since arbitrary high order DeC cannot be recast as a convex combination of explicit Euler method, there is no proof that the solution would stay positive also under that strict condition. On the other side, in the modified Patankar DeC framework, the limiter is only imposed on the water height equation which is unconditionally positive for any CFL by definition. This means that even when the limiter plays a role in the simulation, like in this case, the scheme stays positive with much higher CFL numbers. Furthermore, it should be underlined the fact that the system is linear implicit, which only requires a simple linear solver, e.g. Jacobi method. Even though implicit methods are much slower than their explicit counterpart, solving a linear problem in this case only yields to a 18% increase in the computational time with respect to the fully explicit DeC5, which is nothing compared to the CFL restriction imposed by the positivity limiter.

Circular wet dam break problem

Next, we consider a wet dam break problem with a setup similar to the previous one. In this case the water height of the two basins are h 1 = 10 and h 2 = 0.5, like it was done in [START_REF] Ricchiuto | Stabilized residual distribution for shallow water simulations[END_REF]. In this case, the computational domain is the square [0, 50] × [0, 50] discretized by a 200 × 200 mesh. The simulation is run until t end = 0.8 with CFL=1 and the solutions is displayed only for one quarter of the domain. Figure 14 displays the final snapshot of the solution at time t end = 0.8 plotted both in 2D and 3D displaying the correct evolution of the water height and the time evolution of the water depth extracted along the diagonal of the portion of the domain studied. The method shows good properties such as discontinuities sharply captured, no oscillations and high order approximation of smooth features.

Wave over dry island

For the last test case, we simulate a wave crashing over a dry island showing the robustness of our method when facing more realistic simulations. The computational domain is the rectangle [-5, 5]×[-2, 2] discretized by a 400 × 120 mesh. The simulation can be reproduced by taking Eq. ( 39) for the bathymetry b(x, y) and the following initial conditions h(x, y) = 0.7 -b(x, y) + 0.5 e

1- 1 (1-ρ 2 ) 2 , if ρ 2 < 1, 0,
else, where in case h(x, y) ≤ ε = 10 -6 , we set h(x, y) = ε and u = 0. The simulation has been run with the WENO5-mPDeC5 scheme until a final time t = 1 with CFL=0.9 and the results are displayed in Figure 16 for different times t = 0, 0.25, 0.5, 0.75, 1. In this case, the variable η = h + b has been chosen to be plotted since it better represents the underlying physics. It should be noticed that the top of the island, which is dry at the beginning, gets wet and dry several times during the simulation while never giving rise to problems due to negative water height. This is instead something that we cannot ensure for the DeC5 case. The simulation starts with a background state moving with speed u = 1 which helps the wave traveling towards the island and immediately starts wetting the island from the left and drying it on the right. The wave breaks into two smaller waves respectively approaching and moving away from the island following the eigenvalues of the flux Jacobian. At time t = 0.25 a run up is happening where the traveling wave is trying to submerge the top of the island while, on the other side, a section of the island is drying. In consecutive times, the wave overtakes the island causing different interacting shock fronts and two symmetrical minimum points highlighted in dark blue at time t = 1. Several structures could been observed in this simulations and the repeated wetting/drying procedures never results in troubles for the mPDeC method. In Figure 15 we depict the average number of Jacobi iteration needed for every linear system. The plot compares different CFL numbers and mesh refinements which keep the aspect ratio. The average is plotted with a confidence interval given by ± 1 2 standard deviation. We observe that the CFL number is very incisive in determining how many Jacobi iterations are needed. For low CFL it seems that larger matrices needs more iterations, but this is not uniform with all the tests, as shown in the unsteady vortex test.

ρ 2 = (x + 2) 2 , (u, v) = (1, 0), (45)

Summary and Outlook

We have presented a new well-balanced, positivity preserving high order numerical method for solving the shallow water equations. By re-writing the WENO semi-discretization in terms of productions-destructions terms, we were able to use the modified Patankar Deferred Correction methods of arbitrarily high-order to ensure the unconditionally positivity of water height. The restriction on the CFL number comes only from the explicit DeC-solver for the momentum equations. The used CFL numbers are of the order of 1 and are much larger than the ones used in explicit SSPRK WENO methods in combination with positivity preserving limiters [START_REF] Xing | A survey of high order schemes for the shallow water equations[END_REF], where the CFL must be lowered to 1/6 or 1/12. One can relax further the CFL constraint by using implicit DeC or RK methods, though introducing more difficulties. Even if we have only presented a fifth order method in this manuscript, the approach is actually of arbitrary high-order. With classical explicit approaches, one must use SSPRK to guarantee positivity and, as known from literature [START_REF] Gottlieb | Strong stability preserving Runge-Kutta and multistep time discretizations[END_REF], explicit (implicit) SSPRK methods exist only up to order four (sixth) for general cases. By applying mPDeC, we avoid those issues and the price to pay is that of solving a (very sparse) linear system for the water height. However, as mentioned before, this increase in computational costs is around 18% percent in our numerical simulations, but the procedure can still be optimized. In the future, a detailed performance test in terms of accuracy and run-time is planned. Additionally, in this work we have only considered positivity preservation and well-balanced properties. In the next step, we would like to extend our investigation to entropy conservative/stable methods. Here, various approaches exist as described inter alia in [START_REF] Abgrall | A general framework to construct schemes satisfying additional conservation relations. application to entropy conservative and entropy dissipative schemes[END_REF][START_REF] Abgrall | Analysis of the SBP-SAT stabilization for finite element methods part II: entropy stability[END_REF][START_REF] Chen | Review of entropy stable discontinuous Galerkin methods for systems of conservation laws on unstructured simplex meshes[END_REF][START_REF] Fisher | Discretely conservative finitedifference formulations for nonlinear conservation laws in split form: Theory and boundary conditions[END_REF][START_REF] Fjordholm | Arbitrarily high-order accurate entropy stable essentially nonoscillatory schemes for systems of conservation laws[END_REF][START_REF] Ranocha | Shallow water equations: Split-form, entropy stable, well-balanced, and positivity preserving numerical methods[END_REF], but from our perspective the convex limiting strategies seems the most promising one [START_REF] Hajduk | Monolithic convex limiting in discontinuous Galerkin discretizations of hyperbolic conservation laws[END_REF][START_REF] Kuzmin | Entropy conservation property and entropy stabilization of high-order continuous Galerkin approximations to scalar conservation laws[END_REF]. In order to do so, the basic stability properties of Patankar methods as ODE solvers have to be first fully understood [START_REF] Izgin | On Lyapunov stability of positive and conservative time integrators and application to second order modified Patankar-Runge-Kutta schemes[END_REF][START_REF] Torlo | Issues with positivity-preserving Patankar-type schemes[END_REF].

Finally, herein we have focused on the shallow water system. However, the method can be easily adapted to more complex models, e.g. the shallow water equations together with biochemical processes like algae bloom in oceans, seas and open water cancels or the Euler equations of gas dynamics, where special treatments for the pressure positivity are necessary [START_REF] Zhang | Positivity-preserving high order finite difference WENO schemes for compressible euler equations[END_REF][START_REF] Wang | Robust high order discontinuous Galerkin schemes for twodimensional gaseous detonations[END_REF]. Those extensions can and will be also considered in future works.

A Reconstruction of the primitive variable

The goal of this section is that of deepening the procedure to compute the coeffiecients of the polynomials and linear weights needed to compute the WENO procedure in Section 3. The results and the actual coefficients for WENO5 with 4 Gaussian quadrature points are written in Section B.

In the following assemble the system that allow to find the coefficients for the WENO polynomial at any quadrature point. The symbolic Matlab script coded for this purpose is also available at [START_REF] Ciallella | Modified Patankar Deferred Correction WENO Code for Shallow Water Equations[END_REF]. With few adjustments, the script can be used to compute all the ingredients needed for a WENO reconstruction of arbitrary high order with arbitrary high order quadrature formulae. The reconstruction of the primitive variables is needed in a high order finite volume method, as only cell averages are available from the previous time step. For the sake of simplicity we are going to work in a simpler one-dimensional scalar framework since the reconstruction is done dimension-by-dimension. We consider a scalar function u(x) whose cell averages u i are known. We aim at reconstructing this variable as a polynomial v(x), where the polynomial may vary in different points. In particular, it will be useful to use a primitive of v(x) which we denote by P(x) =

x x0 v(s)ds. Indeed, we can impose that for every cell average, we have

u i = 1 ∆x x i+1/2 x i-1/2 u(x) dx ! = x i+1/2 x i-1/2 v(x)dx = P(x i+1/2 ) -P(x i-1/2 ) ∆x . (46)
The stencil considered for this example is the one used for the WENO5 reconstruction, which is made up by five cell averages as denoted in Figure 17. The next step consists in using interpolating polynomials ϕ j-1/2 , e.g. Lagrange polynomials at cell interfaces, to approximate the primitive P(x), i.e.,

x i-5/2 x i-3/2 x i-1/2 x i+1/2 x i+3/2 x i+5/2 u i-2 u i-1 u i u i+1 u i+2
P(x) = k2+1 j=k1 a j-1/2 ϕ j-1/2 (x) where ϕ j-1/2 (x) = k2+1 =k1,i =j x -x -1/2 x j-1/2 -x -1/2 (47)
where k 1 and k 2 are the extreme indexes of the considered stencil. So, given then k 2 -k 1 +1 cell averages, the degree of the polynomial P will be k 2 -k 1 + 2. For instance, when working with WENO5, the reconstruction is composed by a linear combination of three lower order polynomials of degree 3, so with P ∈ P 4 , obtained through three cell averages. Then, the 3 polynomials, which then depend on the whole 5 cells stencil, will be be combined with weights which depends on the quadrature points into a fifth order accurate reconstruction.

In order to have this result, we need to compute the aforementioned three lower order polynomials and a high order polynomial made up by information coming from the whole stencil. Let us begin with the procedure to compute the high order polynomial with all available cell averages, i.e., k 1 = -2 and k 2 = 2. The lower order polynomials can be easily computed following the same approach explained hereafter. In this case P(x) is a sixth order polynomial and v(x) = P (x) is a fifth order polynomial which gives the right accuracy order. Using [START_REF] Öffner | Arbitrary high-order, conservative and positivity preserving Patankar-type deferred correction schemes[END_REF] for all cell averages in the stencil with the definition of P given in [START_REF] Patankar | Numerical heat transfer and fluid flow[END_REF] and using the property of Lagrangian polynomials ϕ j-1/2 (x -1/2 ) = δ j, , we obtain a system of equations for the coefficients a j-1/2 with solution

a j-1/2 = 0, if j = k 1 , -k1+j-1 =-k1 u i+ , j = k 1 + 1, . . . , k 2 + 1. (48)
Finally, the expression for the high order, ho, approximation polynomial v ho (x) can be written as

v ho (x) = k2+1 j=k1 a j-1/2 ϕ j-1/2 (x) = k2 =-k1 c ho (x) u i+ = 2 =-2 c ho (x) u i+ , (49)
where c ho are obtained collecting all the coefficients and basis functions related to u i+ . In this way, for any quadrature point ξ ∈ [x i-1/2 , x i+1/2 ] we can evaluate this high order polynomial v ho ( ξ). Following the same procedure for three lower order polynomials, lo, associated to the 3-cells stencils S 0 = {u i , u i+1 , u i+2 }, S 1 = {u i-1 , u i , u i+1 } and S 2 = {u i-2 , u i-1 , u i }, we obtain an expression for these low order polynomials

v lo j (x) = 2-j =-j c lo j (x) u i+ . (50)
The last step concerns the computation of the ideal linear weights. These are the weights that allow to recover the high order reconstruction from a linear combination of the lower order ones for a given quadrature point ξ. Therefore, we need to find the linear weights d j such that

v ho ( ξ) = 2 j=0 d j v lo j ( ξ) ⇐⇒ 2 =-2 c ho ( ξ) u i+ = 2 j=0 2-j =-j d j c lo j ( ξ) u i+ , ∀u i+ . (51)
Since (51) must hold for any quintuplet {u i+ } 2 =-2 , we can write a system of five equations in the 3 linear weights d j for each quadrature point ξ, i.e.,

2 j=0 d j c lo j ( ξ) = c ho ( ξ), ∀ = -2, . . . , 2. (52)
This is an overdetermined system with five equations and only three unknowns d j that can be easily solved by means of a least squares method. Moreover, the solutions found verify exactly all the equations, implying that some of the equations are linearly dependent.

B WENO recontruction (r = 3) with four-point Gaussian quadrature rule

The goal of this section is to present the fifth order WENO reconstruction with four-point Gaussian quadrature rule. Up to our knowledge, there is no reference in literature that explicitly define linear weights and polynomial coefficients needed for this WENO reconstruction. Reference [START_REF] Titarev | Finite-volume WENO schemes for three-dimensional conservation laws[END_REF] well described the fifth-order WENO5 with two-point Gaussian quadrature rule, which unfortunately does not allow to go beyond fourth order. Instead, with the four-point Gaussian quadrature rule, one could reach even eighth order.

Let us consider a one-dimensional cell [ξ i-1/2 , ξ i+1/2 ], we hereafter provide the expressions for q(ξ - i+1/2 ) , q(ξ + i-1/2 ) , q   ξ i ± ∆ξ 2

which are used for the first sweep (first two terms corresponding to the two boundaries) and the second sweep (the last two terms corresponding to the 4 quadrature points). For r = 3 we have only three candidate stencil for the reconstruction S 0 = (i, i + 1, i + 2) , S 1 = (i -1, i, i + 1) , S 2 = (i -2, i -1, i).

The corresponding smoothness indicators are given by: β 0 = 13 12 (q i -2q i+1 + q i+2 ) 2 + 1 4 (3q i -4q i+1 + q i+2 ) 2 , β 1 = 13 12 (q i-1 -2q i + q i+1 ) 2 + 1 4 (q i-1 -q i+1 ) 2 , β 2 = 13 12 (q i-2 -2q i-1

+ q i ) 2 + 1 4 (q i-2 -4q i-1 + 3q i ) 2 .
The optimal weights d m for the left boundary extrapolated value q - i+1/2 at x i+1/2 are

d 0 = 3 10 , d 1 = 3 5 , d 2 = 1 10 (55)
and q - i+1/2 is given by q - i+1/2 = 1 6 ω 0 (-q i+2 + 5q i+1 + 2q i ) + 1 6 ω 1 (-q i-1 + 5q i + 2q i+1 ) + 1 6 ω 2 (2q i-2 -7q i-1 + 11q i ) .

The optimal weights d m for the right boundary extrapolated value q + i-1/2 at x i-1/2 are

d 0 = 1 10 , d 1 = 3 5 , d 2 = 3 10 (57)
and q + i-1/2 is given by

q + i-1/2 = 1 6
ω 0 (2q i+2 -7q i+1 + 11q i ) + 1 6 ω 1 (-q i+1 + 5q i + 2q i-1 ) + 1 6 ω 2 (-q i-2 + 5q i-1 + 2q i ) .

For the first Gaussian quadrature point ξ q 1 = ξ i -∆ξ , the optimal weights are:

d 0 = 269 √ 42
The reconstructed value can be computed from the three polynomials associated to each stencil: For the second Gaussian quadrature point ξ q 2 = ξ i -∆ξ , the optimal weights are:

d 0 = 411 √ 30
The reconstructed value can be then computed from the three polynomials associated to the three stencil: For the third Gaussian quadrature point ξ q 3 = ξ i + ∆ξ , the optimal weights are:

p 0 (ξ q 2 ) =     3
d 0 = 1751 √ 35
The reconstructed value can be then computed from the three polynomials associated to the three stencil: For the fourth Gaussian quadrature point ξ q 4 = ξ i + ∆ξ , the optimal weights are:

d 0 =
The reconstructed value can be then computed from the three polynomials associated to the three stencil: For all Gaussian quadrature points the solution in ξ can be easily built by assembling the three polynomials.

2 huv huv hv 2 + g h 2 2 

 222 , S(u, x, y) = -gh

Figure 1 :

 1 Figure 1: Shallow Water Equations: definition of the variables.

Algorithm 1 5 :θ m r ≥ 0 then 6 : 7 : 9 :

 15679 MassRequire: Production-destruction functions p α,β (•), d α,β (•), ∆t, previous correction variables y(k-1) , current subtimestep m. 1: M ← I 2: for α = 1 to S do 3:for β = 1 to S do 4: for r = 0 to M do if M α,β ← M α,β -∆tθ m r p α,β (y r,(k-1) M α,α ← M α,α + ∆tθ m r d α,β (y r,(k-1) M α,β ← M α,β + ∆tθ m r d α,β (y r,(k-1) ) y m,(k-1) β 10: M α,α ← M α,α -∆tθ mr p α,β (y r,(k-1) ) end for 15: return M Algorithm 2 mPDeC Update formula Require: y (k-1) , ∆t, production-destruction functions p α,β (•), d α,β (•), m.

Algorithm 4 1 : 2 : 3 : 4 : 5 :Algorithm 5 3 :

 41234553 WENO FV with PDS structure Require: U i,j , well balanced fluxes Reconstruct on the quadrature points on cell interfaces the variable U in a high order fashion Compute the numerical fluxes at quadrature points on cell interfaces F(u L , u R ) Subtract the correction for well balanced problems and obtain F and G Integrate over the cell interface to obtain the numerical fluxes F i+1/2,j , F i-1/2,j and G i,j-1/2 , G i,j+1/2 Compute p α,β (U), d α,β (U) as in[START_REF] Kopecz | Unconditionally positive and conservative third order modified Patankar-Runge-Kutta discretizations of production-destruction systems[END_REF] 6: return p α,β (U), d α,β (U) Full algorithm Require: U 0 i,j , T 1: t = 0 2: while t < T do Compute ∆t by CFL restrictions 4:

Figure 4 :

 4 Figure 4: Unsteady vortex: test case setting.

10 - 2 Figure 5 :

 25 Figure 5: Unsteady vortex: convergence tests.

Figure 6 :

 6 Figure 6: Unsteady vortex test: computational time and Jacobi iterations

Figure 7 :

 7 Figure 7: Lake at rest: test case setting.

Figure 8 :10 0 10 1 10 2 10 3 10 4 10 Figure 9 :

 8109 Figure 8: Lake at rest: convergence tests without preserving the exact solution.

Figure 10 :

 10 Figure 10: Perturbation analysis over a steady solution: test case setting

  non WB, t = 1

Figure 11 :Figure 12 :

 1112 Figure 11: Perturbation analysis over a steady solution: water height h isolines for well-balanced (left-hand side) and non well-balanced (right-hand side) results.

Figure 13 :

 13 Figure 13: Circular dry dam break: water height h isocontours.

Figure 14 :

 14 Figure 14: Circular wet dam break

Figure 15 :

 15 Figure 15: Wave over dry island: Jacobi iterations and confidence interval.

Figure 16 :

 16 Figure 16: Wave over dry island: η = h + b isocontours at different times.

Figure 17 :

 17 Figure 17: Stencil of five cell averages for WENO5 reconstruction.

Table 1 :

 1 Almost dry lake at rest: mPDeC5-WENO5 WENO5 scheme, even with the checks on the positivity, can have problems and might have exploding velocities and water heights. Since the reconstruction does not guarantee the positivity for such high CFLs, it is not safe to run such simulations.

	N x	Error h	Order h	Error u	Order u	Error v	Order v
	600	8.346e-05	3.899	7.759e-04	3.057	6.911e-04	3.091
	700	4.166e-05	4.508	4.526e-04	3.498	3.648e-04	4.144
	800	2.134e-05	5.007	2.533e-04	4.346	1.886e-04	4.941
	1000 6.185e-06	5.551	7.406e-05	5.511	5.225e-05	5.753
	1200 2.219e-06	5.621	2.478e-05	6.006	1.774e-05	5.924
	1400 1.120e-06	4.438	1.034e-05	5.669	7.561e-06	5.532
	1600 5.896e-07	4.804	5.200e-06	5.148	3.799e-06	5.155
		Table 2: Almost dry lake at rest: SSPRK(6,4)-WENO5
	N x	Error h	Order h	Error u	Order u	Error v	Order v
	600	8.378e-05	3.891	7.765e-04	3.057	6.907e-04	3.091
	700	4.191e-05	4.493	4.532e-04	3.493	3.644e-04	4.147
	800	2.160e-05	4.963	2.538e-04	4.342	1.884e-04	4.942
	1000 6.426e-06	5.434	7.432e-05	5.504	5.228e-05	5.745
	1200 2.569e-06	5.029	2.502e-05	5.972	1.797e-05	5.857
	1400 1.376e-06	4.051	1.058e-05	5.586	7.805e-06	5.410
	1600 9.062e-07	3.127	5.496e-06	4.901	4.096e-06	4.829
	SSPRK(6,4)-

We do not focus on entropy/energy preservation, see[START_REF] Berthon | A fully well-balanced, positive and entropy-satisfying Godunov-type method for the shallow-water equations[END_REF][START_REF] Ranocha | Shallow water equations: Split-form, entropy stable, well-balanced, and positivity preserving numerical methods[END_REF] and reference therein for this topic.

Acknowledgements

M. Ciallella is funded by an Inria PhD fellowship. L. Micalizzi is supported by SNF, project number 200020 175784. D. Torlo is funded by an Inria Postdoc. P. Öffner gratefully acknowledge the support of the Gutenberg Research College and also wants to thank Mario Ricchiuto for his invitation to Inria Bordeaux. All authors would like to thank Jonatan Núñez for sharing his high-order FV-WENO code on his repository [43]. We have started our work by adapting his code.