Duncan Bossion
email: dbossion@ur.rochester.edu

Pengfei Huo
email: pengfei.huo@rochester.edu

General Formulas of the Structure Constants in the su(N) Lie Algebra

We provide the analytic expressions of the totally symmetric and anti-symmetric structure constants in the su(N ) Lie algebra. The derivation is based on a relation linking the index of a generator to the indexes of its non-null elements. The closed formulas obtained to compute the values of the structure constants are simple expressions involving those indexes and can be analytically evaluated without any need of the expression of the generators. We hope that these expressions can be widely used for analytical and computational interest in Physics.

The su(N ) Lie algebra and their corresponding Lie groups are widely used in fundamental physics, particularly in the Standard Model of particle physics [START_REF] Halzen | Quarks and Leptons: An Introductory Course in Modern Particle Physics[END_REF][START_REF] Georgi | Lie Algebras In Particle Physics: from Isospin To Unified Theories[END_REF]. The su(2) Lie algebra is used describe the spin- 1 2 system. Its generators, the spin operators, are Ŝj = 2 σ j with the Pauli matrices

σ 1 = 0 1 1 0 , σ 2 = 0 -i i 0 , σ 3 = 1 0 0 -1 . ( 1
)
The generators of the su(3) Lie algebra are proportional to the Gell-Mann λ matrices [START_REF] Gell-Mann | Symmetries of baryons and mesons[END_REF] as Ŝj = 2 λ j , with

λ 1 =   0 1 0 1 0 0 0 0 0   , λ 2 =   0 -i 0 i 0 0 0 0 0   , λ 3 =   1 0 0 0 -1 0 0 0 0   λ 4 =   0 0 1 0 0 0 1 0 0   , λ 5 =   0 0 -i 0 0 0 -i 0 0   , λ 6 =   0 0 0 0 0 1 0 1 0   λ 7 =   0 0 0 0 0 -i 0 i 0   , λ 8 = 1 √ 3   1 0 0 0 1 0 0 0 -2   . ( 2
)
These matrices are used in quantum chromodynamics as an approximate symmetry of the strong interaction between quarks and gluons [START_REF] Gell-Mann | Symmetries of baryons and mesons[END_REF]. There are different manners for obtaining the generators of an algebra, but the most commonly used one in physics is based on a generalization of the Pauli matrices of su(2) and of the Gell-Mann matrices [START_REF] Gell-Mann | Symmetries of baryons and mesons[END_REF] of su [START_REF] Gell-Mann | Symmetries of baryons and mesons[END_REF], which is what we used in this work.

The generators of the su(N ) Lie algebra can be expressed as follows [START_REF] Kimura | The bloch vector for n-level systems[END_REF][START_REF] Bertlmann | Bloch vectors for qudits[END_REF]. There are a total of N (N -1)/2 symmetric matrices

ŜSnm = 2 |m n| + |n m| , (3)
as well as N (N -1)/2 anti-symmetric matrices

ŜAnm = -i 2 |m n| -|n m| , (4)
and N -1 diagonal matrices (so-called Cartan generators)

ŜDn = 2n(n -1) n-1 k=1 |k k| + (1 -n)|n n| , (5)
The indexes S nm , A nm and D n indicate generators corresponding to the symmetric, anti-symmetric, and diagonal matrices, respectively. The explicit relation between a projection operator |m n| and the generators can be found in Ref. [START_REF] Bertlmann | Bloch vectors for qudits[END_REF]. Note that these generators are traceless Tr[ Ŝi ] = 0, as well as orthonormal Tr[ Ŝi Ŝj ] = 2 2 δ ij . These higher dimensional su(N ) Lie algebra are commonly used in particle physics. For example, the su(6) Lie algebra in the quark model [START_REF] Kokkedee | The quark model[END_REF]; in the search of a Grand Unification Theory (GUT), su(5) has been proposed as the simplest possible version of GUT by Georgi and Glashow [START_REF] Georgi | Unity of All Elementary-Particle Forces[END_REF]. In atomic and optical physics [START_REF] Hioe | N -Level Coherence Vector and Higher Conservation Laws in Quantum Optics and Quantum Mechanics[END_REF], as well as in physical chemistry [START_REF] Meyer | Classical models for electronic degrees of freedom: Derivation via spin analogy and application to F * +H2 → F+H2[END_REF][START_REF] Runeson | Spin-mapping approach for nonadiabatic molecular dynamics[END_REF], spin analogy is used to map the electronic-nuclear dynamics of open quantum systems involving non-adiabaticity [START_REF] Meyer | Classical models for electronic degrees of freedom: Derivation via spin analogy and application to F * +H2 → F+H2[END_REF][START_REF] Runeson | Spin-mapping approach for nonadiabatic molecular dynamics[END_REF][START_REF] Kuratsuji | Path integral in the representation of SU(2) coherent state and classical dynamics in a generalized phase space[END_REF][START_REF] Runeson | Generalized spin mapping for quantum-classical dynamics[END_REF][START_REF] Bossion | Non-adiabatic ring polymer molecular dynamics with spin mapping variables[END_REF]. The su(N ) Lie algebra is also proposed as general mapping between a multi-state Hamiltonian and a classical-like Hamiltonian [START_REF] Runeson | Generalized spin mapping for quantum-classical dynamics[END_REF]. In quantum computing, the su(N ) Lie algebra [START_REF] Bertlmann | Bloch vectors for qudits[END_REF][START_REF] Patera | The pauli matrices in n dimensions and finest gradings of simple lie algebras of type an1[END_REF] is widely used for describing the qudit [START_REF] Wang | Qudits and high-dimensional quantum computing[END_REF], the state of a d-state quantum system.

The crucial quantities to define an algebra are the generators and what govern their commutation and anticommutation relations. These relations give rise to the structure constants, where the totally anti-symmetric structure constant, f ijk , is defined through the commutation relation

[ Ŝi , Ŝj ] = i N 2 -1 k=1 f ijk Ŝk , (6)
which is anti-symmetric under the exchange of any two indexes. The totally symmetric structure constant, d ijk , is determined through the anti-commutation relation

{ Ŝi , Ŝj } = 2 N δ ij Î + N 2 -1 k=1 d ijk Ŝk , (7)
which is symmetric under the exchange of any two indexes. Using the properties in Eq. 6 and Eq. 7, as well as the fact that the generators are traceless and orthogonal to each other, one can derive the well-known results [START_REF] Haber | Useful relations among the generators in the defining and adjoint representations of SU(N )[END_REF]

f ijk = -i 2 3 Tr [ Ŝi , Ŝj ] Ŝk ; d ijk = 2 3 Tr { Ŝi , Ŝj } Ŝk . (8)
Using the above relation as well as general expressions of the generators (Eq. 3-5) one can compute the numerical values of all the structure constants of the su(N ) algebra through Eq. 8, as has been commonly done in the literature. However, this requires laboriously efforts of different combinations of the N 2 -1 generators of su(N ), which remain numerous even when considering the symmetry properties of the structure constants when N is large. Despite the extensive usage and crucial role of these structure constants, to the best of our knowledge [17], we are not aware of any analytic expression (closed formulas) of f ijk and d ijk .

In this letter, we derive the analytic expressions of the structure constants in the su(N ) Lie algebra. The key results are summarized in Eq. 14 for f ijk and Eq. 21 for d ijk . We first determine the relation of the indexes of the generators (see Eq. 3-5) S nm , A nm , and D nm with the label n and m. We note that one can use recursive relations to obtain the generators of su(N + 1) from the generators of su(N ) [18]. This helps to determine the indexes of the generators as

S nm =n 2 + 2(m -n) -1, (9a)
A nm =n 2 + 2(m -n), (9b)
D n =n 2 -1. (9c)
Note that there is no overlap among the indexes as long as the conditions 1 ≤ m < n ≤ N holds. Hence, there is a one to one correspondence between the value i ∈ {S nm , A nm , D n } with indexes {n, m}, which helps us to identify the type of generator it labels. This is the first key step to determine closed formulas of structure constants.

The Totally Anti-symmetric Structure Constants f ijk . The commutation relation between two symmetric generators is

[ ŜSnm , ŜS n ′ m ′ ] = i N 2 -1 k=1 f SnmS n ′ m ′ k Ŝk (10) = 2 4 δ nm ′ (|m n ′ | -|n ′ m|) + δ nn ′ (|m m ′ | -|m ′ m|) + δ mm ′ (|n n ′ | -|n ′ n|) + δ mn ′ (|n m ′ | -|m ′ n|) = i 2 δ nm ′ ŜA n ′ m + δ nn ′ ( ŜA m ′ m -ŜA mm ′ ) + δ mm ′ ( ŜA n ′ n -ŜA nn ′ ) -δ mn ′ ŜA nm ′ .
With constraints from δ ij , we can identify the index of each generator, and hence obtain the non-zero analytic expression of f SnmS n ′ m ′ k , which are summarized in the first line of Eq. 14.

For a symmetric and an anti-symmetric generator, the commutation relation is

[ ŜSnm , ŜA n ′ m ′ ] = i N 2 -1 k=1 f SnmA n ′ m ′ k Ŝk (11) = i 2 4 δ nn ′ (|m m ′ | + |m ′ m|) -δ nm ′ (|m n ′ | + |n ′ m|) + δ mn ′ (|n m ′ | + |m ′ n|) -δ mm ′ (|n n ′ | + |n ′ n|) = i 2 δ nn ′ ( ŜS m ′ m + ŜS mm ′ ) -δ nm ′ ŜS n ′ m + δ mn ′ ŜS nm ′ -δ mm ′ ( ŜS n ′ n + ŜS nn ′ ) + 2 δ mm ′ δ nn ′ (2|m m| -2|n n|) .
The first two lines of the last equality in Eq. 11 give us directly several structure constants. The last line of Eq. 11 contains diagonal elements, hence we know it will be a combination of diagonal generators. In fact, we can prove (see Supplemental Material [19], Sec I) that

i 2 2 (|m m| -|n n|) (12) = i n 2(n -1) ŜDn + n-1 k>m ŜD k 2k(k -1) - m -1 2m ŜDm .
This helps to determine the rest of the structure constant f SnmA n ′ m ′ k with the expressions documented in Eq. 14.

The commutation relations between symmetric and diagonal generators are not required as we already obtained all the non-zero structure constants involving diagonal and symmetric generators (as we know that we cannot obtain a diagonal matrix through the commutator of a symmetric and a diagonal generator).

Between two anti-symmetric generators, the commutation relation is

[ ŜAnm , ŜA n ′ m ′ ] = i N 2 -1 k=1 f AnmA n ′ m ′ k Ŝk (13) = 2 4 δ nm ′ (|n ′ m| -|m n ′ |) + δ nn ′ (|m m ′ | -|m ′ m|) + δ mm ′ (|n n ′ | -|n ′ n|) + δ mn ′ (|m ′ n| -|n m ′ |) = i 2 -δ nm ′ ŜA n ′ m + δ nn ′ ( ŜA m ′ m -ŜA mm ′ ) + δ mm ′ ( ŜA n ′ n -ŜA nn ′ ) + δ mn ′ ŜA nm ′ ,
which helps to determine the structure constants involving all anti-symmetric generators (second line of Eq. 14). The remaining totally anti-symmetric structure constants are computed through the commutator between two diagonal generators, which is [ ŜDn ,

ŜD n ′ ] = i N 2 -1 k=1 f DnD n ′ k Ŝk = 0 (see proof in Supplemental Ma- terial, Sec II), indicating a zero value for all f DnD n ′ k .
This was a known fact, as the diagonal matrices are generators of the Cartan subalgebra of su(N ), and they commute by definition [START_REF] Georgi | Lie Algebras In Particle Physics: from Isospin To Unified Theories[END_REF].

To summarize, all of the non-zero totally antisymmetric structure constants are expressed as follows

f SnmS kn A km = f SnmS nk A km = f SnmS km A kn = 1 2 , (14)
f AnmA km A kn = 1 2 , f SnmAnmDm = - m -1 2m , f SnmAnmDn = n 2(n -1) , f SnmAnmD k = 1 2k(k -1)
, m < k < n.

One of the interesting usages of these expressions is the construction of an adjoint representation of the su(N ) Lie algebra, which is a defining representation of SO(N 2 -1), whose generators Ti are (N 2 -1)×(N 2 -1) matrices. The (jk)-th matrix element of Ti is [ Ti ] jk = -if ijk . Thus, Ti are anti-symmetric, non-diagonal and traceless, and we provide analytic expressions to obtain them through the totally anti-symmetric structure constants expressions.

Totally Symmetric Structure Constants d ijk . The anti-commutation relation between two symmetric generators is

{ ŜSnm , ŜS n ′ m ′ } = 2 N δ SnmS n ′ m ′ Î + N 2 -1 k=1 d SnmS n ′ m ′ k Ŝk = 2 4 δ nm ′ (|m n ′ | + |n ′ m|) + δ nn ′ (|m m ′ | + |m ′ m|) + δ mm ′ (|n n ′ | + |n ′ n|) + δ mn ′ (|n m ′ | + |m ′ n|) = 2 δ nm ′ ŜS n ′ m + δ nn ′ ( ŜS m ′ m + ŜS mm ′ ) + δ mm ′ ( ŜS n ′ n + ŜS nn ′ ) + δ mn ′ ŜS nm ′ + δ nn ′ δ mm ′ 2 (2|m m| + 2|n n|) . ( 15
)
We know that the last line of Eq. 15 only involves diagonal matrices. In fact, we can prove that (see Supplemental Material, Sec III)

2 2 (|m m| + |n n|) = 2 N Î (16) + N k>n 2 k(k -1) ŜD k + 2 -n 2n(n -1) ŜDn + n-1 k=m+1 1 2k(k -1) ŜD k - m -1 2m
ŜDm .

Thus, we can extract all the non-zero d SnmS n ′ m ′ k , with the expressions summarized in Eq. 21.

Between a symmetric and an anti-symmetric generator, the anti-commutation relation reads

{ ŜSnm , ŜA n ′ m ′ } = N 2 -1 k=1 d SnmA n ′ m ′ k Ŝk (17) = i 2 4 δ nm ′ (|n ′ m| -|m n ′ |) + δ nn ′ (|m m ′ | -|m ′ m|) + δ mm ′ (|n ′ n| -|n n ′ |) + δ mn ′ (|n m ′ | -|m ′ n|) = 2 δ nm ′ ŜA n ′ m + δ nn ′ ( ŜA mm ′ -ŜA m ′ m ) + δ mm ′ ( ŜA n ′ n -ŜA nn ′ ) + δ mn ′ ŜA nm ′ ,
from which one can extract d SnmA n ′ m ′ k . Note that based on Eq. 17, there is no diagonal component ŜD , thus all the d SAD = 0. The anti-commutator between symmetric and diagonal generators is not necessary as we already obtained the structure constants involving those generators by permutation [20].

Computing the anti-commutators between two antisymmetric generators gives

{ ŜAnm , ŜA n ′ m ′ } = 2 N δ AnmA n ′ m ′ Î + N 2 -1 k=1 d AnmA n ′ m ′ k Ŝk = 2 4 δ nn ′ (|m m ′ | + |m ′ m|) -δ nm ′ (|m n ′ | + |n ′ m|) + δ mm ′ (|n n ′ | + |n ′ n|) -δ mn ′ (|n m ′ | + |m ′ n|) = 2 δ nn ′ ( ŜS m ′ m + ŜS mm ′ ) -δ nm ′ ŜS n ′ m + δ mm ′ ( ŜS n ′ n + ŜS nn ′ ) -δ mn ′ ŜS nm ′ + δ nn ′ δ mm ′ 2 (2|m m| + 2|n n|) , (18)
where we recognize that the last line of Eq. 18 is identical to the last line of Eq. 15, which can be expressed as generators in Eq. 18. We do not need to compute the anti-commutator between an asymmetric and a diagonal generator as we already have the result by permutation from Eq. 18 (and Eqs. 17 indicates d SAD = 0).

The remaining d ijk values are obtained through the anti-commutator between two diagonal generators

{ ŜDn , ŜD n ′ } = 2 N δ DnD n ′ Î + N 2 -1 k=1 d DnD n ′ k Ŝk (19) = 2 2n(n -1)2n ′ (n ′ -1) n-1 k=1 δ kk ′ (|k k ′ | + |k ′ k|) + δ kn ′ (1 -n ′ )(|k n ′ | + |n ′ k|) + δ nk ′ (1 -n)(|n k ′ | + |k ′ n|) + δ nn ′ (1 -n)(1 -n ′ )(|n n ′ | + |n ′ n|) = 2n(n -1) 2δ kn ′ ŜD n ′ + 2n ′ (n ′ -1) 2δ nk ′ ŜDn + δ nn ′ 2 n(n -1) ( n-1 k=1 |k k| + (1 -n) 2 |n n|).
One can see that only diagonal matrices are involved in the last line of Eq. 19 and there is no off-diagonal element. We can prove that (see Supplemental Material, Sec IV)

2 n(n -1) ( n-1 k=1 |k k| + (1 -n) 2 |n n|) = 2 N Î (20) + N k>n 2 k(k -1) ŜD k + (2 -n) 2 n(n -1)
ŜDn , which helps to determine all d DnD n ′ k . We summarize all the non-zero totally symmetric structure constants as follows

d SnmS kn S km = d SnmA kn A km = d SnmA mk A nk = 1 2 , ( 21
)
d SnmA nk A km = - 1 2 , d SnmSnmDm = d AnmAnmDm = - m -1 2m , d SnmSnmD k = d AnmAnmD k = 1 2k(k -1) , m < k < n, d SnmSnmDn = d AnmAnmDn = 2 -n 2n(n -1) , d SnmSnmD k = d AnmAnmD k = 2 k(k -1) , n < k, d DnD k D k = 2 n(n -1) , k < n, d DnDnDn = (2 -n) 2 n(n -1)
.

In a recent work on deriving the quantum Liouvillian of coupled electronic-nuclear DOFs based upon the su(N ) representation (see Appendix F in Ref. [START_REF] Runeson | Generalized spin mapping for quantum-classical dynamics[END_REF]), these d ijk are explicitly present in the equation of motion. Having the above analytic expressions will facilitate future theoretical developments.

Mapping Hamiltonian using the su(N ) Lie Algebra. The SU(2) representation of the Lie algebra (spin-1 2 analogy) is often used in quantum dynamics to study systems with two states [START_REF] Meyer | Classical models for electronic degrees of freedom: Derivation via spin analogy and application to F * +H2 → F+H2[END_REF][START_REF] Runeson | Spin-mapping approach for nonadiabatic molecular dynamics[END_REF][START_REF] Bossion | Non-adiabatic ring polymer molecular dynamics with spin mapping variables[END_REF]. For a two level system with the Hamiltonian Ĥ = H

0 Î + 1 H • Ŝ = H 0 Î + k 1 (H k • Ŝk ) = H 0 Î + 1 [2R(V 12 ) • Ŝx + 2I(V 12 ) • Ŝy +(V 11 -V 22 )
• Ŝz ], it can be shown (through the Heisenberg equations of motion (EOMs)) that

d dt S i = 1 3 j,k ε ijk H j S k , (22)
where S i = Tr[ρ Ŝi ] is the expectation value of Ŝi , ρ being the density operator, and ε ijk the Levi-Civita tensor which is the two-dimension totally anti-symmetric structure constant f ijk . This is equivalent to the precession of spin of a spin-1 2 system around a magnetic field M = H, which is a well-known result (eg, Page 424 of Ref. [START_REF] Cohen-Tannoudji | Quantum Mechanics[END_REF]). Those EOMs exactly obey the time-dependent Schrödinger equation (TDSE), ċi = -i 2 j=1 V ij c j . More specifically, with an arbitrary state defined as |Ψ = c 1 |1 + c 2 |2 , by using the transformation S

x = Re{c * 1 c 2 }, S y = Im{c * 1 c 2 } and S z = 1 2 (|c 1 | 2 -|c 2 | 2 )
, one can show that Eq. 22 is equivalent to TDSE.

For a system with N states, one can use the su(N ) Lie algebra for the spin analogy [START_REF] Hioe | N -Level Coherence Vector and Higher Conservation Laws in Quantum Optics and Quantum Mechanics[END_REF][START_REF] Meyer | Classical models for electronic degrees of freedom: Derivation via spin analogy and application to F * +H2 → F+H2[END_REF][START_REF] Runeson | Generalized spin mapping for quantum-classical dynamics[END_REF], describing the precession of the N -states system [22]. This su(N ) analogy is based on a reformulation of the Hamiltonian Ĥ = H 0 Î + Ve ( R) with the generators of su(N ) as follows [START_REF] Hioe | N -Level Coherence Vector and Higher Conservation Laws in Quantum Optics and Quantum Mechanics[END_REF][START_REF] Runeson | Generalized spin mapping for quantum-classical dynamics[END_REF] Ĥ = H 0 Î + 1

N 2 -1 k=1 H k Ŝk , (23)
f ijk H j S k . ( 24
)
where f ijk is the totally anti-symmetric structure constant of su(N ). For an arbitrary state defined as |Ψ =

|c k | 2 - n -1 2n |c n | 2 , (26)
and the analytic expressions of the f ijk , one can show that Eq. 24 is equivalent to the TDSE . Hence Eq. 24 now has a closed formula.

Conclusion.

In this letter, we provide the analytic expressions of the totally symmetric and totally antisymmetric structure constants for the su(N ) Lie algebra. We hope that these expressions can be widely used for analytical and computational interest in Physics, as they are valid for any dimension N of su(N ) Lie algebra without the need to explicitly compute the commutation and anti-commutation relations or use any generator. The structure constants bear important information on the algebra they belong to, and the possibility to obtain those constants with simple relations can bring insight into the high dimensional su(N ) Lie algebra, which might be challenging otherwise.

where H k = 1 1 N Î + 1 2 N 2 - 1 k=1

 11221 Tr[ Ĥ Ŝk ]. Similarly, for the density matrix [8] ρ = S k Ŝk where S k = Tr[ρ Ŝk ]. Plugging the su(N ) generator expression of Ĥ and ρ into the quantum Liouville equation i ∂ ρ/∂t = [ Ĥ, ρ], one arrives at the following equation [23] which can be viewed as the generalization of the spin precession [8

N k=1 c k

  |k , using the transformationS Snm = Re{c * m c n }, S Anm = Im{c * m c n },

ACKNOWLEDGMENTS

This work was supported by the National Science Foundation CAREER Award under Grant No. CHE-1845747.


Mirco Ciallella
email: mirco.ciallella@ensam.eu

Lorenzo Micalizzi

Victor Michel-Dansac
email: victor.michel-dansac@inria.fr

Philipp Öffner

Davide Torlo
email: davide.torlo@sissa.it

A high-order, fully well-balanced, unconditionally positivity-preserving finite volume framework for flood simulations

Keywords: well-balancing, moving steady solutions, positivity preservation, high-order accuracy, flood simulations, shallow water, WENO

In this work, we present a high-order finite volume framework for the numerical simulation of shallow water flows. The method is designed to accurately capture complex dynamics inherent in shallow water systems, particularly suited for applications such as tsunami simulations. The arbitrarily high-order framework ensures precise representation of flow behaviors, crucial for simulating phenomena characterized by rapid changes and fine-scale features. Thanks to an ad-hoc reformulation in terms of production-destruction terms, the time integration ensures positivity preservation without any time-step restrictions, a vital attribute for physical consistency, especially in scenarios where negative water depth reconstructions could lead to unrealistic results. In order to introduce the preservation of general steady equilibria dictated by the underlying balance law, the highorder reconstruction and numerical flux are blended in a convex fashion with a well-balanced approximation, which is able to provide exact preservation of both static and moving equilibria. Through numerical experiments, we demonstrate 1 the effectiveness and robustness of the proposed approach in capturing the intricate dynamics of shallow water flows, while preserving key physical properties essential for flood simulations.

Introduction

The Saint-Venant equations, also known as the shallow water (SW) equations, characterize the behavior of hydrostatic free surface waves influenced by gravity. These nonlinear hyperbolic partial differential equations (PDEs) are applicable under the assumption of either very large wavelengths or very shallow depths. They find extensive use across engineering domains, including river and estuarine hydrodynamics, urban flood management, and tsunami risk evaluation. The numerical approximation of the SW equations remains a highly active area of research. Numerous original methods have been developed across various contexts and settings: finite volume [START_REF] Audusse | A fast and stable well-balanced scheme with hydrostatic reconstruction for shallow water flows[END_REF][START_REF] Gallardo | On a well-balanced high-order finite volume scheme for shallow water equations with topography and dry areas[END_REF][START_REF] Noelle | High-order well-balanced finite volume WENO schemes for shallow water equation with moving water[END_REF][START_REF] Díaz | High order exactly well-balanced numerical methods for shallow water systems[END_REF][START_REF] Bollermann | A well-balanced reconstruction of wet/dry fronts for the shallow water equations[END_REF][START_REF] Kurganov | Central-upwind schemes for the Saint-Venant system[END_REF][START_REF] Kurganov | Finite-volume schemes for shallow-water equations[END_REF][START_REF] Cheng | Moving-water equilibria preserving central-upwind schemes for the shallow water equations[END_REF][START_REF] Xing | On the advantage of well-balanced schemes for moving-water equilibria of the shallow water equations[END_REF][START_REF] Berthon | A fully well-balanced, positive and entropysatisfying Godunov-type method for the shallow-water equations[END_REF][START_REF] Michel-Dansac | A well-balanced scheme for the shallow-water equations with topography or manning friction[END_REF][START_REF] Michel-Dansac | A two-dimensional highorder well-balanced scheme for the shallow water equations with topography and manning friction[END_REF][START_REF] Ciallella | An arbitrary high order and positivity preserving method for the shallow water equations[END_REF][START_REF] Ciallella | Arbitrary High Order WENO Finite Volume Scheme with Flux Globalization for Moving Equilibria Preservation[END_REF], continuous and discontinuous finite element [START_REF] Kashiyama | Three-step explicit finite element computation of shallow water flows on a massively parallel computer[END_REF][START_REF] Hauke | A stabilized finite element method for the Saint-Venant equations with application to irrigation[END_REF][START_REF] Yabe | Conservative semi-lagrangian CIP technique for the shallow water equations[END_REF][START_REF] Bunya | Improvements in mass conservation using alternative boundary implementations for a quasi-bubble finite element shallow water model[END_REF][START_REF] Takase | Space-time SUPG finite element computation of shallow-water flows with moving shorelines[END_REF][START_REF] Song | The shifted boundary method for hyperbolic systems: Embedded domain computations of linear waves and shallow water flows[END_REF][START_REF] Behzadi | An exact source-term balancing scheme on the finite element solution of shallow water equations[END_REF][START_REF] Xing | Exactly well-balanced discontinuous Galerkin methods for the shallow water equations with moving water equilibrium[END_REF][START_REF] Mantri | Well-balanced discontinuous Galerkin scheme for 2×2 hyperbolic balance law[END_REF][START_REF] Arpaia | An efficient covariant frame for the spherical shallow water equations: Well balanced DG approximation and application to tsunami and storm surge[END_REF][START_REF] Bender | Entropy-conservative discontinuous Galerkin methods for the shallow water equations with uncertainty[END_REF][START_REF] Wintermeyer | An entropy stable discontinuous Galerkin method for the shallow water equations on curvilinear meshes with wet/dry fronts accelerated by GPUs[END_REF], residual distribution [START_REF] Ricchiuto | Stabilized residual distribution for shallow water simulations[END_REF][START_REF] Ricchiuto | Application of conservative residual distribution schemes to the solution of the shallow water equations on unstructured meshes[END_REF][START_REF] Ricchiuto | An explicit residual based approach for shallow water flows[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] Arpaia | r-adaptation for Shallow Water flows: conservation, well balancedness, efficiency[END_REF][START_REF] Arpaia | Well-balanced residual distibution for the ALE spherical shallow water equations on moving adaptive meshes[END_REF], and so on.

The ultimate goal of these approaches is to provide reliable and physically meaningful simulations for real-world applications, while demanding minimal computational resources. High-order methods are particularly suitable in this context, as they are able to achieve smaller errors within coarser discretizations. Furthermore, an effective strategy for designing numerical methods with reduced errors is through structure-preserving techniques. These techniques aim to replicate additional consistency conditions beyond those explicitly defined by the system of equations themselves. For the SW equations, the focus is on preserving positive water heights, equilibrium or stationary states, and implementing entropy conservation or dissipation methods. First, the SW equations with source terms are known to admit a family of stationary solutions, which are characterized by a balance between flux divergence and source terms. This concept is connected with that of a well-balanced (WB) discretization, typically characterized by its ability to replicate one or more of these equilibria at the discrete level [START_REF] Berthon | A fully well-balanced, positive and entropysatisfying Godunov-type method for the shallow-water equations[END_REF][START_REF] Berthon | A very easy high-order well-balanced reconstruction for hyperbolic systems with source terms[END_REF][START_REF] Castro | Well-balanced high-order finite volume methods for systems of balance laws[END_REF][START_REF] Gómez-Bueno | High-order well-balanced methods for systems of balance laws: a control-based approach[END_REF][START_REF] Gómez-Bueno | Collocation methods for high-order well-balanced methods for systems of balance laws[END_REF][START_REF] Mantri | Fully well-balanced entropy controlled discontinuous Galerkin spectral element method for shallow water flows: global flux quadrature and cell entropy correction[END_REF][START_REF] Berberich | High order well-balanced finite volume methods for multi-dimensional systems of hyperbolic balance laws[END_REF][START_REF] Chertock | Well-balanced schemes for the Euler equations with gravitation: Conservative formulation using global fluxes[END_REF]. This WB property is crucial for complex, time-dependent simulations, as discretization errors due to the non-preservation of stationary regions could accumulate over time. Second, in the context of flood simulations, it is a necessity to have provably positive discretizations, avoiding negative water heights. To obtain a provably positive reconstruction in the context of highorder weighted essentially non-oscillatory [START_REF] Shu | Essentially non-oscillatory and weighted essentially non-oscillatory schemes for hyperbolic conservation laws[END_REF] (WENO) schemes, an effective positive limiter has been introduced and further developed in [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]. As proven in these references, this limiter achieves the preservation of positive reconstruction, but it restricts the CFL condition, for classical SSPRK [START_REF] Gottlieb | Strong stability-preserving high-order time discretization methods[END_REF] schemes, to the weight of the Gauss-Lobatto quadrature rule of the corresponding space accuracy (e.g., 1/12 for fifth-order schemes). To circumvent this issue, unconditionally positivity preserving time-stepping strategies [START_REF] Meister | A positivity preserving and well-balanced DG scheme using finite volume subcells in almost dry regions[END_REF][START_REF] Ciallella | An arbitrary high order and positivity preserving method for the shallow water equations[END_REF] for the SW equations have been proposed, based on a suitable reformulation of the finite volume semi-discretization in terms of production-destruction terms. These approaches are based on the modified Patankar trick [START_REF] Patankar | Numerical heat transfer and fluid flow[END_REF][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 positivitypreserving scheme for production-destruction equations with applications to nonequilibrium flows[END_REF][START_REF] Öffner | Arbitrary high-order, conservative and positivity preserving Patankar-type deferred correction schemes[END_REF]. The linearly implicit nature of this approach allows for a relaxation of the aforementioned time-step constraint at a reasonable computational cost.

In this paper, we deal with the possibility of integrating the additional preservation of general static and moving equilibria into the arbitrary high-order positivity preserving framework introduced in [START_REF] Ciallella | An arbitrary high order and positivity preserving method for the shallow water equations[END_REF]. To achieve this, we suitably modify the spatial discretization relying on ideas presented in [START_REF] Berthon | A very easy high-order well-balanced reconstruction for hyperbolic systems with source terms[END_REF]. In particular, we perform a convex blending between the original discretization and a WB one, able to exactly capture general families of equilibria. The approach is able to tackle challenging flood simulations, proving to be a good candidate for real-life applications.

The paper is structured as follows. We first introduce the multidimensional SW system in Section 2. Then, the high-order WB positive numerical scheme is detailed in Section 3, where the space discretization is discussed, and in Section 4, where we present the production-destruction formulation in combination with high-order modified Patankar time schemes. The results of the numerical validation are reported in Section 5. Finally, Section 6 is left for conclusions and further developments.

Shallow water equations

The two-dimensional SW equations consist in a hyperbolic system of PDEs, extensively used in many applications to describe the behavior of water flows. Their Eulerian formulation on a space domain Ω ⊆ R 2 , assuming no friction and a time-independent bathymetry, reads

∂u ∂t + ∂F ∂x (u) + ∂G ∂y (u) = S(x, y, u), ∀(x, y) ∈ Ω, ∀t ∈ [0, T f ], (2.1)
where conserved variables, fluxes and source term are respectively given by

u =   h hu hv   , S(x, y, u) = -gh   0 ∂b ∂x (x, y) ∂b ∂y (x, y)   , F(u) =   hu hu 2 + g h 2 2 huv   , G(u) =   hv huv hv 2 + g h 2 2   , (2.2)
with h being the water height, u and v the velocity components of the flow along the x and y directions respectively, g the gravitational constant, and b(x, y) the bathymetry. We also introduce the free surface water level η := h + b, and the discharge variables along the two directions x and y, defined as q x := hu and q y := hv respectively.

Notable properties of the SW equations, which have been drawing the interest of the scientific community in recent years and which play a central role in the context of this paper, are the positivity of the water height and the existence of non-trivial steady solutions. In the context of numerical schemes preserving moving equilibria, one is interested in a detailed capturing of steady solutions satisfying

∂u ∂t ≡ 0 ⇔ ∂F ∂x (u) + ∂G ∂y (u) = S(x, y, u), ∀(x, y) ∈ Ω ⊆ R 2 , ∀t ∈ [0, T f ]. (2.
3)

The simplest and most known steady solution is the so-called "lake at rest" given by

u = v = 0, η ≡ η 0 ∈ R + 0 , ∀ (x, y) ∈ Ω, ∀t ∈ [0, T f ]. (2.4)
Generally speaking, steady solutions are not known in closed-form and they are characterized by the analytical balance (2.3). The smooth steady solutions tackled in this work are the pseudo-monodimensional states in the form

       ∂ ∂s q s = 0, ∂ ∂s q 2 s 2h 2 + g(h + b) = 0, (2.5)
where s is a general handle for the x or the y variable. For more information on these steady solutions, the reader is referred for instance to [START_REF] Michel-Dansac | A well-balanced scheme for the shallow-water equations with topography[END_REF]. For what follows, it is useful to define the so-called equilibrium variables

E s (x, y, u) = q s q 2 s 2h 2 + g(h + b) . (2.6)
After (2.5), steady solutions are characterized by E s (x, y, u) being constant in space.

The system of PDEs under consideration is discretized using the Method of Lines (MOL), a numerical approach that treats space and time independently. In particular, space and time discretizations are the main focus of the next two sections.

Well-balanced space discretization

This section is dedicated to the space discretization. First, in Section 3.1, we describe our classical, non-well-balanced high-order discretization. Then, Section 3.2 is devoted to the generalization of a strategy to achieve a high-order well-balanced (WB) discretization, which was introduced in a one-dimensional setting in [START_REF] Berthon | A very easy high-order well-balanced reconstruction for hyperbolic systems with source terms[END_REF]. Here, we generalize this technique for a two-dimensional WENO framework, applying the basic idea dimension by dimension. The underlying principle consists in a simple blending between a high-order discretization and a WB discretization to be used where a steady state is detected. The main strengths of this approach are its low cost (no nonlinear equations need to be solved) and its ease of use (it consists in multiplying the reconstruction by a suitable coefficient). We emphasize that the resulting scheme will be able to capture and preserve all the moving 1D steady solutions given by (2.5), and not just the so-called lake at rest solution, where velocity vanishes.

Basic high-order discretization

The computational domain Ω is discretized in a Cartesian fashion via N x × N y nonoverlapping control volumes

Ω i,j = [x i-1/2 , x i+1/2 ] × [y j-1/2 , y j+1/2 ],
with uniform spatial steps ∆x = x i+1/2 -x i-1/2 and ∆y = y j+1/2 -y j-1/2 .

Finite volume methods are based on deriving a system of ordinary differential equations (ODEs) for the cell averages of the solution in each control volume Ω i,j 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.
The first step to obtain such a system is to integrate (2.1) over Ω i,j , thus getting

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), (3.1
) where S i,j is the source term average 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, u) dxdy, and F i+1/2,j and G i,j+1/2 are the averages of the fluxes over the cell boundaries

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, 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.

So far, Equation (3.1) has been exactly derived from (2.1). However, in order to obtain the numerical scheme, we need to discretize the fluxes and the source averages.

To that end, we rely on the following ingredients: a high-order reconstruction of the conservative variables in each control volume (WENO [START_REF] Shu | Essentially non-oscillatory and weighted essentially non-oscillatory schemes for hyperbolic conservation laws[END_REF] in our case), consistent quadrature formulas to discretize all integrals (Gauss-Legendre with Q points in our case), and suitable numerical fluxes to compute the fluxes in the boundary integrals (described later on). In the remainder of this section, we drop the time dependency to shorten notation.

Let us first focus on the discretization of the fluxes averages and, more in detail, on F i+1/2,j , as G i,j+1/2 is obtained similarly. Once the reconstruction in each control volume has been performed, at each quadrature point y q ∈ [y j-1/2 , y j+1/2 ] of each edge x i+1/2 we have two high-order reconstructed values for u, corresponding to x L i+1/2 and

x R i+1/2 , which will be referred to as the left and right high-order extrapolated values u L i+1/2,q = u HO (x L i+1/2 , y q ) and u R i+1/2,q = u HO (x R i+1/2 , y q ).

By applying a consistent quadrature rule, the flux in the x-direction reads

F i+1/2,j ≈ Q q=1 w q F(u L i+1/2,q , u R i+1/2,q ),
where F is a consistent numerical flux, and w q is the normalized quadrature weight associated to the quadrature node y q . The choice of F is discussed in Section 3.2.

The high-order source term averages are computed as

S i,j ≈ Q q=1 Q p=1
w q w p S(x q , y p , u HO (x q , y p )), with a surface quadrature obtained as the tensor product of the classical 1D quadrature used for the edges and u HO being the local reconstruction of the solution in the cell. Despite its robustness in capturing discontinuities, while minimizing the oscillations, the WENO reconstruction may provide some negative reconstructed values for the water height, especially close to dry regions. Such negative water heights are not physically admissible, and in fact will immediately lead to the simulation crashing. In order to avoid such an issue, we adopt for the water height reconstruction the positivity limiter introduced in [START_REF] Perthame | On positivity preserving finite volume schemes for Euler equations[END_REF] and further discussed in [START_REF] Zhang | On positivity-preserving high order discontinuous Galerkin schemes for compressible Euler equations on rectangular meshes[END_REF].

As shown in [START_REF] Xing | A survey of high order schemes for the shallow water equations[END_REF], provable positivity preservation for the water height, in the context of this framework, is subjected to severe CFL constraints, when adopting standard time integration techniques. In particular, assuming a simple forward Euler time-stepping and a Lax-Friedrichs numerical flux, the limit CFL guaranteeing positivity preservation is CFL FE := w Lobatto 1 , where w Lobatto 1 is the first weight of the adopted high-order Gauss-Lobatto quadrature rule. This corresponds to CFL FE = 1/12 for a quadrature of order 5. The restriction gets even worse as the order of accuracy increases, e.g., we have CFL FE = 1/20 for a quadrature of order 7. The adoption of high-order SSPRK methods slightly relaxes the constraint, but not significantly. Indeed, for instance, using the SSPRK [START_REF] Audusse | A fast and stable well-balanced scheme with hydrostatic reconstruction for shallow water flows[END_REF][START_REF] Arpaia | An efficient covariant frame for the spherical shallow water equations: Well balanced DG approximation and application to tsunami and storm surge[END_REF] discretization relaxes the condition to CFL SSPRK (5,4) = 1.508 CFL FE . The adopted time discretization, described in Section 4, allows us to drop such limitations and to run simulations at any CFL without violating the positivity constraint on the water height. Due to the explicit nature of the time scheme used for the discharge equations, however, the (far less restrictive) stability constraint CFL ≤ 1 of explicit schemes applies.

Well-balanced blending

We now describe the WB strategy, which makes possible the capture of steady states characterized by constant equilibrium variables (2.6). The key idea comes from the following remark: for the simulation of a steady solution, a well-balanced scheme is exact, and therefore has a better accuracy than any high-order scheme. For unsteady simulations, high-order schemes are more accurate, and should be used whenever the solution is not steady. To achieve a seamless switch between high-order and wellbalanced schemes, we propose a simple blending between the two. This blending is performed according to a suitable steady solution indicator, defined below.

For simplicity, we only derive the reconstruction along the x-direction. The extension to the y-direction is easily performed following a dimension by dimension approach. We replace the reconstructed variables at the interfaces by the convex combination between the high-order extrapolated values and the cell averages

ũL i+1/2,q = (1 -θ i+1/2,j )U i,j + θ i+1/2,j u L i+1/2,q , ũR i+1/2,q = (1 -θ i+1/2,j )U i+1,j + θ i+1/2,j u R i+1/2,q , (3.2)
where θ i+1/2,j is a steady state indicator. On the one hand, it should vanish when the equilibrium variables (2.6) are constant in space; in this case, the modified reconstructed values ũL i+1/2,q are equal to the cell averages U i,j . On the other hand, when far from any equilibrium, ũL i+1/2,q should be an approximation of order P , where P is the order of the discretization (herein, P = 5).

Following [START_REF] Berthon | A very easy high-order well-balanced reconstruction for hyperbolic systems with source terms[END_REF], we define θ i+1/2,j by

θ i+1/2,j = ε i+1/2,j ε i+1/2,j + ∆x C i+1/2,j P , with ε i+1/2,j := ∥E x (x i+1 , y j , θ i+1/2,j ) -E x (x i , y j , θ i+1/2,j )∥,
where C i+1/2,j is a quantity independent of ∆x, which is here chosen, at a given time iteration, as the time residual difference at the previous iteration

C i+1/2,j := 1 2 U n i+1,j -U n-1 i+1,j ∆t + U n i,j -U n-1 i,j
∆t .

We remark that, at equilibrium C i+1/2,j → 0, hence θ i+1/2,j → 0 as well resulting in the low-order WB reconstruction. Similarly, the source term discretization is defined as

Si,j = 1 2 θ i-1/2,j + θ i+1/2,j S i,j + 1 2 1 -θ i-1/2,j S WB i-1/2,j + 1 -θ i+1/2,j S WB i+1/2,j ,
where S WB i-1/2,j and S WB i+1/2,j represent WB discretizations of the source term at the interfaces described in [START_REF] Michel-Dansac | A well-balanced scheme for the shallow-water equations with topography[END_REF].

In order to ensure stability in the context of unsteady wet-dry simulations, we found experimentally useful to introduce, with respect to the classical approach [START_REF] Berthon | A very easy high-order well-balanced reconstruction for hyperbolic systems with source terms[END_REF][START_REF] Michel-Dansac | A well-balanced scheme for the shallow-water equations with topography[END_REF], a similar convex combination in the flux definition

F(u L i+1/2,q , u R i+1/2,q ) = (1 -θ i+1/2,j ) FWB (u L i+1/2,q , u R i+1/2,q ) + θ i+1/2,j FLF (u L i+1/2,q , u R i+1/2,q ),
where FWB represents the WB approximate Riemann solver presented in [START_REF] Michel-Dansac | A well-balanced scheme for the shallow-water equations with topography[END_REF], while FLF is a robust local Lax-Friedrichs numerical flux reading

FLF (u L , u R ) = 1 2 F(u R ) + F(u L ) - 1 2 s max u R + u L ,
where s max is the spectral radius of the normal flux Jacobian of system (2.1). Remark 1. The discretized terms S WB i-1/2,j , S WB i+1/2,j and FWB are designed in such a way to guarantee an exact equilibrium with respect to steady states in the form (2.6), when taking in input the cell averages. The reader can easily verify that, when a steady state of this type is considered, then all θ i+1/2,j are equal to 0 and the scheme reduces to the WB version. Indeed, the modified reconstruction (3.2) degenerates to the cell averages. This means that, despite guaranteeing an exact capturing of the steady states, the basic WB discretization is directly based on cell averages without any reconstruction, and it is, therefore, only first order accurate in general [START_REF] Berthon | A very easy high-order well-balanced reconstruction for hyperbolic systems with source terms[END_REF]. Remark 2. We emphasize an important property of the proposed strategy: the steady solution indicator is defined such that the nonlinear system (2.5) never has to be solved. Instead, it merely relies on evaluating the equilibrium variables (2.6) at the cell interfaces. Remark 3. It should be noticed that the proposed scheme is well-balanced when θ i+1/2,j goes to 0. Numerically speaking, this consists in defining a low enough threshold (10 -10 in the numerical experiments) to set θ i+1/2,j = 0.

Unconditionally positive time discretization

In this section, we describe the time-stepping strategy, which consists in a slight modification of arbitrary high-order deferred correction (DeC) methods for ODEs [START_REF] Abgrall | High order schemes for hyperbolic problems using globally continuous approximation and avoiding mass matrices[END_REF]. In particular, the water height update is reinterpreted as a Production-Destruction System (PDS) and then the modified Patankar trick is applied in order to achieve unconditional preservation of its positivity as in [START_REF] Ciallella | An arbitrary high order and positivity preserving method for the shallow water equations[END_REF]. Both DeC methods and Patankar trick have a long history. In particular, for more information on DeC the interested reader is referred to [START_REF] Dutt | Spectral Deferred Correction Methods for Ordinary Differential Equations[END_REF][START_REF] Micalizzi | A new efficient explicit deferred correction framework: analysis and applications to hyperbolic PDEs and adaptivity[END_REF][START_REF] Micalizzi | Efficient iterative arbitrary high-order methods: an adaptive bridge between low and high order[END_REF][START_REF] Veiga | On improving the efficiency of ADER methods[END_REF][START_REF] Veiga | DeC and ADER: Similarities, Differences and a Unified Framework[END_REF], while Patankar (and modified Patankar) tricks are detailed in [START_REF] Patankar | Numerical heat transfer and fluid flow[END_REF][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][START_REF] Burchard | Application of modified Patankar schemes to stiff biogeochemical models for the water column[END_REF][START_REF] Izgin | A study of the local dynamics of modified Patankar DeC and higher order modified Patankar-RK methods[END_REF][START_REF] Kopecz | On order conditions for modified Patankar-Runge-Kutta schemes[END_REF][START_REF] Kopecz | Unconditionally positive and conservative third order modified Patankar-Runge-Kutta discretizations of production-destruction systems[END_REF][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 positivitypreserving scheme for production-destruction equations with applications to nonequilibrium flows[END_REF][START_REF] Meister | On unconditionally positive implicit time integration for the DG scheme applied to shallow water flows[END_REF][START_REF] Ortleb | Patankar-type Runge-Kutta schemes for linear PDEs[END_REF].

Deferred Correction method

To introduce the DeC method, let us consider the Cauchy problem

d dt c(t) = H(t, c(t)), t ∈ [0, T f ], c(0) = c 0 , (4.1)
where c : [0, T f ] → R Nc is the unknown solution, with N c components, and H : [0, T f ]×R Nc → R Nc is a given function satisfying the classical smoothness assumptions, which guarantee the existence of a unique solution to the Cauchy problem (4.1). As is customary in the context of one-step methods, we focus on a generic interval [t n , t n+1 ] of size ∆t := t n+1 -t n and, given c n ≈ c(t n ), we seek an approximation c n+1 of c(t n+1 ).

Following [START_REF] Abgrall | High order schemes for hyperbolic problems using globally continuous approximation and avoiding mass matrices[END_REF][START_REF] Micalizzi | A new efficient explicit deferred correction framework: analysis and applications to hyperbolic PDEs and adaptivity[END_REF], we introduce M + 1 subtimenodes t m in the interval [t n , t n+1 ], which are such that

t n = t 0 < t 1 < • • • < t M = t n+1 .
The DeC method under consideration consists in an explicit fixed point iterative procedure to compute the approximation of c at all subtimenodes simultaneously. The update formula is given by

c m,(p) := c 0 + ∆t M ℓ=0 θ m ℓ H(t ℓ , c ℓ,(p-1) ), m = 1, . . . , M, p ≥ 1, (4.2)
where c m,(p) is the approximation of the solution in the subtimenode t m obtained at the p th iteration and, for each m, the coefficients (θ m ℓ ) ℓ∈{0,...,M } are the normalized weights of the high-order quadrature formula over [t 0 , t m ] associated to the subtimenodes. In particular, in the previous update formula, we set c m,(p) = c 0 := c n whenever m = 0 or p = 0. One can show that, for small enough ∆t, the iterative process converges. Furthermore, the order of accuracy of c M,(p) with respect to c(t n+1 ) is min (p, R), i.e., each iteration corresponds to an increase in the order of accuracy by one, until a saturation value R, which depends on the number and on the distribution of the adopted subtimenodes. For example, evenly spaced subtimenodes lead to R = M +1, while Gauss-Lobatto subtimenodes yield R = 2M . In this paper, we use Gauss-Lobatto subtimenodes. Therefore, the optimal way of reaching order P is to perform P fixed-point iterations with M + 1 subtimenodes, where M = ⌈ P 2 ⌉. Hence, the arbitrarily high-order time integration method presented in this section, combined with the space discretization described in Section 3, defines an arbitrarily high-order, fully well-balanced framework for the numerical solution of the SW equations (2.1) -(2.2). However, at this level, nothing can be said, in general, about the positivity of the water height. In the next subsection, we present the modification to be performed in the time integration of the water height, guaranteeing unconditional positivity.

Modified Patankar DeC method

In this section, we first focus, in Section 4.2.1, on the unconditionally positive time integration of a specific class of ODEs, namely Production-Destruction Systems (PDSs). Then, we describe in Section 4.2.2 how to apply these notions to the SW equations. A PDS is said to be positive if, starting by a positive initial condition, we get a positive evolution of all the components, i.e.,

Unconditionally positive time integration of PDSs

c(0) > 0 =⇒ c(t) > 0, ∀t ∈ [0, T f ], (4.4)
where the comparison operator, applied to vectors, is meant to be applied to each scalar component.

Conservative and positive PDSs arise in many applications and many numerical methods have been developed to preserve such properties. A successful approach, in this context, is given by the (modified) Patankar trick [START_REF] Patankar | Numerical heat transfer and fluid flow[END_REF][START_REF] Burchard | A high-order conservative Patankar-type discretisation for stiff systems of production-destruction equations[END_REF], which is based on the introduction of some weights on the production and destruction terms. In particular, the application of the modified Patankar trick to the DeC scheme (mPDeC) [START_REF] Öffner | Arbitrary high-order, conservative and positivity preserving Patankar-type deferred correction schemes[END_REF] is characterized by replacing (4.2) with the following update

c m,(p) α = c 0 α + ∆t M ℓ=0 θ m ℓ   Nc β=1 p α,β (c ℓ,(p-1) ) c m,(p) γ(β,α,θ m ℓ ) c m,(p-1) γ(β,α,θ m ℓ ) - Nc β=1 d α,β (c ℓ,(p-1) ) c m,(p) γ(α,β,θ m ℓ ) c m,(p-1) γ(α,β,θ m ℓ )   ,
where c m,(p) = c 0 := c n whenever m = 0 or p = 0, and γ is a switch function defined as

γ(α, β, θ) := α, if θ ≥ 0, β, if θ < 0.
For guidelines concerning the number P of iterations to be performed, and the associated accuracy, the reader is referred to the discussion regarding the standard DeC scheme at the end of Section 4.1. The mPDeC method is positive and conservative, i.e., it satisfies Moreover, the method is linearly implicit and can be recast in compact form as

Mc m,(p) = c n ,
where the matrix M is defined as

M(c (p-1) , m) α,β =            1 + ∆t M ℓ=0 Nc k=1 k̸ =α θ m ℓ c m,(p-1) α d α,k (c ℓ,(p-1) )χ {θ m ℓ ≥0} -p i,k (c ℓ,(p-1) )χ {θ m ℓ <0} , for β = α, -∆t M ℓ=0 θ m l c m,(p-1) β p α,β (c ℓ,(p-1) )χ {θ m ℓ ≥0} -d α,β (c ℓ,(p-1) )χ {θ m ℓ <0} , for β ̸ = α, (4.5)
with χ {•} the indicator function, i.e., a switch with value equal to 1 if the argument condition is true, 0 otherwise. One can prove that the matrix is column diagonally dominant and hence invertible. Furthermore, it is possible to show that, for any b > 0, the solution to Mc = b is such that c > 0. At the implementation level, the system is solved though the Jacobi method, which is provably convergent due to the fact that M is column diagonally dominant. Moreover, in order to avoid divisions by zero, as in [START_REF] Meister | A positivity preserving and well-balanced DG scheme using finite volume subcells in almost dry regions[END_REF][START_REF] Öffner | Arbitrary high-order, conservative and positivity preserving Patankar-type deferred correction schemes[END_REF], the following mollification of the ratios in the matrix (4.5) is considered

n d ≈ 0, if d < 10 -8 , 2d•n d 2 +max{d 2 ,10 -8 } , if d ≥ 10 -8 .
Further details are omitted to avoid lengthening the paper. However, they are thoroughly discussed in [START_REF] Ciallella | An arbitrary high order and positivity preserving method for the shallow water equations[END_REF].

We now explain how the presented notions can be applied to the finite volume semi-discretization of the SW equations.

Application to the Shallow Water equations

The key idea is to reinterpret the water height semi-discretization as a PDS and to apply the modified Patankar trick to the water height DeC update, while performing a standard DeC time-stepping on the updates of the discharge in the x-and y-directions. From (3.1), we note that each cell communicates with the neighboring cells, sharing common edges, via numerical fluxes. Thus, in such a context, the components c α are given by the water height averages h i,j , with indices α identified as couples [i, j], while the production and destruction terms are given by the associated numerical fluxes. Let us recall that the water height equation has no source term contribution.

Considering all the neighbors to the cell [i, j] (i.e., cells sharing an edge with cell [i, j]), one can define the following production and destruction terms

p [i,j],[i-1,j] (U) = + 1 ∆x F(1) i-1/2,j (U) + , d [i,j],[i-1,j] (U) = - 1 ∆x F(1) i-1/2,j (U) -, p [i,j],[i+1,j] (U) = - 1 ∆x F(1) i+1/2,j (U) -, d [i,j],[i+1,j] (U) = + 1 ∆x F(1) i+1/2,j (U) + , p [i,j],[i,j-1] (U) = + 1 ∆y Ĝ(1) i,j-1/2 (U) + , d [i,j],[i,j-1] (U) = - 1 ∆y Ĝ(1) i,j-1/2 (U) -, p [i,j],[i,j+1] (U) = - 1 ∆y Ĝ(1) i,j+1/2 (U) -, d [i,j],[i,j+1] (U) = + 1 ∆y Ĝ(1) i,j+1/2 (U) + , (4.6)
where the superscripts + and -respectively represent the positive and the negative part, while the superscript (1) represents the first component of the numerical fluxes. These production and destruction terms, as well as their relationships (4.6) with the numerical fluxes, are sketched in Figure 1.

Ω i-1,j Ω i,j Ω i+1,j Ω i,j+1 Ω i,j-1 F(1) i-1/2,j (U) F(1) i+1/2,j (U) Ĝ(1) i,j-1/2 (U) Ĝ(1) 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]
Fig. 1 Sketch of the PDS structure for the control volume Ω i,j .

In light of the previous discussion, it is therefore easy to apply the modified Patankar trick to the DeC update of the water height. Let us remark that the strategy provides unconditional positivity of the water height with respect to the time step ∆t. This results in great computational advantages with respect to standard explicit time integration techniques subjected to the typical positivity-preserving CFL constraints. Further details, including a detailed description of a possible implementation, can be found in [START_REF] Ciallella | An arbitrary high order and positivity preserving method for the shallow water equations[END_REF].

Numerical results

In this section, we report the results of several numerical experiments demonstrating the good properties of the scheme, including its robustness. In particular, the tests are meant to verify the high-order accuracy in Section 5.1, the WB property for both stationary and moving equilibria (in Sections 5.2 and 5.3), and the ability to deal with tough flood simulations involving dry areas in Section 5.4. We assume g = 9.81 unless otherwise specified. Let us remark that the basic ingredients of the scheme allow us to reach arbitrarily high-orders of accuracy. Here, we focus on the fifth order version.

Unsteady vortex

Through this test [START_REF] Ricchiuto | Analytical travelling vortex solutions of hyperbolic equations for validating very high order schemes[END_REF], we verify the high-order accuracy of the space and time discretizations, without considering the source term for the moment. To that end, we set b ≡ 0. Therefore, this test is meant to verify the high-order accuracy of the flux discretization; the high-order accuracy of the source term discretization will be checked in a later test.

The considered computational domain is the square Ω := [0, 3] × [0, 3], and the vortex is given by a perturbation δ of a homogeneous background field (h 0 , u 0 , v 0 ) := (1, 2, 3). Let us define the variable r(x, y, t) := (x -x c (t)) 2 + (y -y c (t)) 2 , expressing the distance between (x, y) and vortex center (x c (t), y c (t)) := (1.5, 1.5) + (u 0 t, v 0 t).

The water height is then given by h(r) := h 0 + δh(r), with

δh(r) := -γ      exp - 1 arctan 3 (1 -r 2 ) , if r < 1, 0, otherwise,
where γ := 0.1 is the vortex amplitude. The velocity field, defined by (u, v) := (u 0 , v 0 )+ (δu, δv), is characterized by the following perturbation

δu δv = g r ∂h ∂r y -y c -(x -x c ) ,
where ∂h ∂r is the derivative of h with respect to r, which depends only on the radial distance from the center of the vortex

∂h ∂r (r) =        6 γ r exp(- 1 arctan 3 (1-r 2 ) ) arctan 4 (r 2 -1)(1 + (r 2 -1) 2 ) , if r < 1, 0, otherwise.
We assume periodic boundary conditions and a final time T f := 0.1. It is important to highlight the fact that this solution is C ∞ , which is a fundamental property for testing arbitrarily high-order schemes [START_REF] Ricchiuto | Analytical travelling vortex solutions of hyperbolic equations for validating very high order schemes[END_REF].

The convergence test is run on Cartesian meshes of sizes 25 2 , 50 2 , 100 2 , 200 2 , 300 2 , and 400 2 . The error, denoted by ||ϵ h (u)||, is computed as the L 1 norm of the difference between the approximated solution and the exact one. Figure 2 shows the initial water height for this test case (left panel) and the retrieved fifth order convergence trend expected from theory (right panel). These results are also reported in Table 1, where fifth order accuracy is shown to be achieved.

4.41

Table 1 Errors and orders of accuracy, with respect to the number Nx of cells, for the traveling vortex from Section 5.1. We indeed observe fifth-order accuracy.

Lake at rest

We now focus on showing the capability of the proposed scheme to exactly preserve the lake at rest steady state, governed by (2.4). We first tackle the exact capture of the steady state in Section 5.2.1, and we then perform a perturbation analysis in Section 5.2.2.

Exact capturing

In this section, we demonstrate that the proposed scheme is able to exactly capture the lake at rest steady solution in two situations: a fully wet case, where the water height never vanishes, and a wet-dry case, where the water height may vanish. mesh size ∥ϵ h (u)∥ h q x q y (b) WB scheme Fig. 3 Wet lake at rest from Section 5.2: convergence test. Left panel: non-WB scheme; we observe fifth order accuracy. Right panel: WB scheme; we observe an accuracy up to machine precision for each mesh size.

Wet lake at rest

We first consider the lake at rest steady state given by b(x, y) = 0.1 sin (2 π x) cos (2 π y) , h(x, y, t) = 1 -b(x, y), u = v = 0, on the computational domain Ω := [0, 1] × [0, 1] with periodic boundary conditions and final time T f := 0.1. In this case, we test the scheme with and without the WB modification. The purpose of this test is twofold. First, with the WB modification, the lake at rest should be exactly preserved (up to machine precision). This will confirm that the WB property is satisfied in this case. Second, without the WB modification, the method should converge to fifth order accuracy. This will verify the correct implementation of the source term. Both convergence trends are presented in Figure 3 and the expected results are obtained. The resolutions of the Cartesian meshes used for this test are 25 2 , 50 2 , 100 2 and 200 2 .

Wet-dry lake at rest

We now present a numerical experiment to show the preservation of a lake at rest steady state in the presence of dry areas. That is to say, the water height will vanish in some parts of the domain. In particular, by virtue of the mPDeC approach, the proposed method is able to deal with such dry states while having a much relaxed CFL constraint compared to traditional high-order techniques. Indeed, we can set CFL ≃ 1 rather than CFL ≃ 1/12. We consider, on the domain Ω := [-5, 5] × [-5, 5], the following bathymetry

b(x, y) :=      exp 1 - 1 1 -r 2 , if r 2 < 1, 0, otherwise,
where mesh size ∥ϵ h (u)∥ h q x q y (b) WB scheme Fig. 5 Wet-dry lake at rest from Section 5.2: convergence test. Left panel: non-WB scheme; we observe second order accuracy. Right panel: WB scheme; we observe an accuracy up to machine precision for each mesh size.

r 2 = x 2 + y 2 . ( 5
This bathymetry represents an island located in the center of the domain. The water height is defined as h(x, y, t) := max(0.7 -b(x, y), 0). In Figure 4, we display the water height (left panel) and then bathymetry (right panel). We observe that dry areas occur in the center of the domain, where the island is located. The test is performed with periodic boundary conditions and final time T f := 1. Just like before, we present the results of a convergence analysis obtained with and without the WB blending. It should be noted that, due to the discontinuity in the derivative of the water height, the non-WB scheme can achieve at most second order convergence, while machine precision is expected by the WB version. The results are reported in Figure 5; they agree with the expected behavior.

Perturbation analysis

Let us consider the computational domain Ω := [-5, 5] × [-2, 2], the bathymetry b defined in Equation (5.1), and the lake at rest steady state characterized by a total water height η 0 := 1.5. Then, we consider the following perturbation of the steady condition

η := η 0 +      0.05 exp 1 - 1 (1 -ρ 2 ) 2 , if ρ 2 < 1, 0, otherwise,
, where we have set

ρ 2 = 9((x + 2) 2 + (x -0.5) 2 ).
We adopt a Cartesian mesh of 100 × 30 elements, with periodic boundary conditions, CFL = 0.8 and a final time T f := 0.375. We still test the non-WB and WB versions of the method, in order to highlight the advantages of the latter setting. However, in this case, we suppress the blending by adopting all the coefficients θ of Section 3.2 equal to 0 in such a way to always use the WB discretization.

The results are displayed in Figure 6. It can be noticed that the evolution of the perturbation is sharply captured by the WB version of the scheme. Instead, in the non-WB case, numerical oscillations, due to the discretization error, propagate from the bathymetry and prevent the proper capturing of the perturbation evolution.

Moving equilibria

In this section, we test the WB properties of the scheme to capture moving equilibria satisfying (2.5). As already specified, they are pseudo-1D states. Therefore, in the context of this section, we focus on the variable s = x, and we drop the dependency on y, being clear that all quantities are constant along the y-direction. At the numerical level, the variable y does not play any role either. Hence, the adopted mesh configurations will be characterized by a uniform distribution of cells along the x-direction, with various values of N x ranging from 25 to 200, and a constant number N y := 5 of cells along the y-direction, with periodic boundary conditions assumed in such direction.

In this 1D frictionless case, moving equilibria are characterized by constant equilibrium variables (2.6). Therefore, although there is no closed-form expression of such steady solutions, they can be computed pointwise, for a given bathymetry, by solving a cubic equation derived from (2.5), see for example [START_REF] Delestre | SWASHES: a compilation of shallow water analytic solutions for hydraulic and environmental studies[END_REF][START_REF] Ciallella | Arbitrary High Order WENO Finite Volume Scheme with Flux Globalization for Moving Equilibria Preservation[END_REF][START_REF] Micalizzi | Novel well-balanced continuous interior penalty stabilizations[END_REF]. The steady flow regime then depends on the prescribed boundary conditions, and is obtained after a transient phase. We focus here on subcritical and supercritical flows, numerically obtained with the initial and boundary conditions described in Table 2, where the final time T f is chosen such that the simulation reaches the steady state (i.e., to make the time residual vanish). We take the following smooth bathymetry b(x) := 0.05 sin(x -12.5) exp(1 -(x -12.5) 2 ), on the computational domain Ω := [0, 25] × [0, 1]. The gravity constant is set here to g := 9.812 as in [START_REF] Ciallella | Arbitrary High Order WENO Finite Volume Scheme with Flux Globalization for Moving Equilibria Preservation[END_REF]. Flow regime T f h(x, 0) q(x, 0) h(0, t) h(25, t) q(0, t) q(25, t)

Subcritical 200 2 -b(x) 0 - 2 4.42 - Supercritical 50 2 -b(x) 0 2 - 24 -
Table 2 Initial and boundary conditions for the subcritical and supercritical flows from Section 5.3. The simulation is always initialized to a lake at rest. Empty cells correspond to Neumann boundary conditions.

Again, we test the WB and the non-WB versions of the scheme. We emphasize that we do not try to merely exactly preserve the steady solution, but to capture it: the WB numerical scheme is expected to converge towards the steady solution with machine accuracy, even after the transient, unsteady phase.

We start by presenting the numerical results obtained for the subcritical steady flow. The solution computed with the WB method, with N x := 200, is presented in Figure 7. We display three quantities: the water height h, the x-discharge q x , and the second component of the equilibrium variables E x h 2 + g(h + b). Recall that both q x and E

(2) x should be constant in this case; we can indeed appreciate the ability of the WB blending to capture constant q x and E

(2)

x . The exact capture of q x and E

(2) x is also visible from the results of the convergence test reported in Figure 8. In particular, we observe that the WB version of the scheme is able to obtain machine precision errors with respect to such variables. Notice that the errors on the water height can be computed following two approaches: the first one, which takes as a reference the exact bathymetry function b(x); the second one, which considers the discrete bathymetry in cell average b i,j . Usually the first approach is employed for classical convergence analysis, however the second one is very common in the field of well-balanced schemes to check whether the scheme is able to preserve the discrete version of the considered equilibrium. While the first method computes the error by using the exact bathymetry evaluated at quadrature points, the second one considers the reconstructed bathymetry to measure the reference equilibrium. For the first convergence test, in line with Remark 1, we expect the error h to scale with first order. However, a second order superconvergence is obtained, due to the exact preservation of q x and E

(2)

x . For the second convergence analysis, the error h (discrete b) provides the proof of the exact preservation of discrete steady states of the WB scheme with machine precision obtained for all meshes. On the other hand, the non-WB scheme produces, as expected, bigger errors which scale with the expected fifth order. Let us notice that a very high level of mesh refinement would be needed in order to obtain, with the non-WB scheme, errors comparable to the ones obtained with the WB version, especially on q x and E

(2)

x . This ensures that, for a given error, the WB method has a much smaller computational cost than the non-WB one.

Similar considerations apply to the supercritical case. The numerical solution computed with the WB method, with N x := 200, is reported in Figure 9. In addition, the convergence plots of both WB and non-WB schemes can be found in Figure 10. Also in this case, we observe the same features and trends as before: the ability of the WB version to capture h with discrete bathymetry and the constant equilibrium variables q x and E

(2)

x up to machine precision, and to obtain much smaller errors with respect to the non-WB version.

Flooding simulations

We finally present the numerical results of flooding simulations performed with the proposed high-order WB positivity-preserving method. While we so far have focused on the validation of the proposed method on standard academic test cases, we now deal with more challenging applications. These applications correspond to waves over

x .

dry areas, and prove the suitability of the proposed approach in the context of realworld situations. We start by presenting a wave over a dry island in Section 5.4.1, and then we move to the simulation of a tsunami over three obstacles in Section 5.4.2.

Wave over a dry island

In this test, we simulate a wave over a dry island. The computational domain is the rectangular region Ω := [-5, 5] × [-2, 2], partitioned into a mesh with 400 × 120 elements. We refer to [START_REF] Ciallella | An arbitrary high order and positivity preserving method for the shallow water equations[END_REF], Section 6.8, for the bathymetry function b(x, y) and the specific initial and boundary conditions. The simulation was run until a final time T f := 5, with a CFL number set to 0.9. The results are presented at various times in Figure 11. The variable η, along with the bathymetry b, have been displayed. Indeed, it allows for a clearer understanding of the underlying physics.

x and h when considering the discrete bathymetry, and second-order accuracy for h when considering the exact bathymetry function.

The simulation starts with a background state moving from left to right at speed u = 1, propelling the wave towards the island. This causes the island to get wet from the left side and to dry from the right side. Thus, the top of the island, initially dry, undergoes multiple wet and dry cycles throughout the whole simulation, without encountering any issue related to negative water height. This is not guaranteed for classical time integration schemes, among which SSPRK schemes, for such high CFL numbers. Various structures are observable in this simulation like vortices and shocks, and the recurring wetting/drying processes are optimally tackled by the proposed scheme.

Tsunami on three obstacles

Finally, the simulation of a tsunami over several obstacles is presented. Simulations of this kind are often performed [START_REF] Guermond | Well-balanced second-order convex limiting technique for solving the Serre-Green-Naghdi equations[END_REF] since they represent a good starting point to move towards the simulation of real coastal engineering problems. In this simulation, we consider a shock impacting three conical obstacles. More specifically, we consider the with b i (x, y) := c(x, y, x i , y i , R i , A i ), where c is a cone function defined as c(x, y, x c , y c , R, A) :=

A R R -(x -x c ) 2 + (y -y c ) 2 , if (x -x c ) 2 + (y -y c ) 2 < R, 0, otherwise.
In particular, we have R i := 0.5 and A i := 3, for all i, and (x 1 , y 1 ) := (1, -1) T , (x 2 , y 2 ) := (1, 1) and (x 3 , y 3 ) := (2, 0). The initial condition is given by

  h u v   (x, y, 0) := [1.5 -b(x, y), 4, 0] T , if x < -3.5, [0, 0, 0] T , otherwise. ( 5
x and h when considering the discrete bathymetry, and secondorder accuracy for h when considering the exact bathymetry function.

The prescribed boundary conditions are

• inflow at the left of the domain, obtained by imposing q(x = -5, y, t) := 3(1 + cos(2πt))e -2t ; • transmissive at the right of the domain;

• solid walls at the top and bottom of the domain.

We remark that, to simulate a realistic configuration, a time-dependent inlet condition has been chosen to represent a series of waves impacting the obstacles after the tsunami. The final time is T f := 3, and we take a CFL condition of 0.8 for added stability.

The results are reported in Figure 12. We start from an initial configuration where the majority of the domain is dry and where the initial tsunami is represented by a discontinuity in the water height, defined in (5.2). Already from the first snapshots in Figures 12(c) and 12(d), we can appreciate the wetting process happening with several structures forming on the right of the three bodies. Thanks to the time-dependent inlet condition, the dynamic of the simulation keeps evolving with shock interactions occurring due to the crushing between new and old wave fronts, while wetting and drying processes keep happening in many parts of the domain.

This simulation best represents the potential of this framework, which is able to retain high-order accuracy, preserve important structures of the model, and deal with complex fluid phenomena. Simulations of this kind are not only challenging but also computationally expensive, due to the accuracy required to capture all flow features. The choice of the considered time-stepping scheme, able to relax the typical severe CFL constraints imposed by positivity preservation, has a huge impact on the computational resources needed to perform these simulations. In fact, keeping the same fifth-order accuracy, we are able to consistently reduce the computational time with respect to classical time integration techniques, provably guaranteeing positivity of the discrete water height. Moreover, the well-balanced procedure is non-intrusive and computationally cheap, and it is able to preserve the equilibrium variables of the model, which are crucial in the context of flooding simulations.

Summary and outlook

In this paper, we presented a high-order, fully well-balanced, unconditionally positivity-preserving framework for flood simulations. The discretization based on the notion of production-destruction terms, presented in [START_REF] Ciallella | An arbitrary high order and positivity preserving method for the shallow water equations[END_REF], has been extended to treat general moving equilibria appearing in shallow water systems. The advantage of this framework lies in the possibility of preserving the positivity of the water height with no constraint on the CFL. This is a real strength with respect to classical time integration schemes, which experience strong CFL reduction as the order of the method increases, and allows for more realistic applications thanks to the huge computational gain. In order to achieve the general WB property, while keeping the productiondestruction formulation, the high-order reconstruction is blended with a WB one, as proposed in [START_REF] Berthon | A very easy high-order well-balanced reconstruction for hyperbolic systems with source terms[END_REF][START_REF] Michel-Dansac | A well-balanced scheme for the shallow-water equations with topography[END_REF]. This allows to achieve structure preservation for moving equilibria reached after a transient simulation, as shown in Section 5.3. On the contrary, when wet-dry simulations are considered, far from existing equilibria, the approach is able to properly perform high-order space and time integration without causing simulation crashes.

There are several perspectives to this work. They and range from deep questions on the numerical analysis and stability of modified Patankar schemes, which is an open research topic [START_REF] Torlo | Issues with positivity-preserving Patankartype schemes[END_REF][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] Izgin | A study of the local dynamics of modified Patankar DeC and higher order modified Patankar-RK methods[END_REF], especially when coupled to space discretizations in the context of PDEs, to the possible development of this approach on unstructured meshes to exploit advanced mesh adaptation algorithms to capture the flow features with even better resolution, and to save even more computational resources.

  PDSs are systems of ODEs characterized by the following structure Nc β=1 p α,β (c) -Nc β=1 d α,β (c), α = 1, . . . , N c , c(0) = c 0 , where c = (c α ) α∈{1,...,Nc} , and where p α,β and d α,β are real non-negative Lipschitzcontinuous functions from R Nc to R * + . More specifically, we are interested in a subfamily of PDSs fulfilling two extra constraints: conservation and positivity. A PDS is said to be conservative if, ∀α, β ∈ {1, . . . , N c } and ∀c ∈ R Nc , we have p α,β (c) = d β,α (c), thus implying

  Nc α=1 c α (t) = Nc α=1 c α (0), ∀t ∈ [0, T f ]. (4.3)

  and c n > 0 =⇒ c n+1 > 0, which are nothing but natural translations, at the discrete level, of the continuous constraints (4.3) and (4.4).

Fig. 2

 2 Fig. 2 Unsteady vortex from Section 5.1: convergence test. Left panel: depiction of the initial condition. Right panel: error lines, showing that the scheme is indeed of fifth order accuracy.

. 1 )Fig. 4

 14 Fig. 4 Wet-dry lake at rest from Section 5.2: depiction of the water height (left panel) and of the bathymetry (right panel).

Fig. 6

 6 Fig. 6 Perturbation analysis of the lake at rest solution: η = h + b isocontours at different times. Top panel, subfigure (a): initial condition. Left panels, subfigures (b), (d) and (f): non-WB scheme; right panels, subfigures (c), (e) and (g): WB scheme.

Fig. 7

 7 Fig. 7 Subcritical flow, test case from Section 5.3. Top left panel: free surface water level η, and bathymetry b rescaled by a factor of 0.15 and shifted by 1.97. Top right panel: discharge qx. Bottom panel: second equilibrium variable E

Fig. 8

 8 Fig. 8 Subcritical flow, test case from Section 5.3: convergence test. Left panel: non-WB scheme; we observe fifth-order accuracy. Right panel: WB scheme; we observe machine precision accuracy for the equilibrium variables qx and E

Fig. 9

 9 Fig. 9 Supercritical flow, test case from Section 5.3. Top left panel: free surface water level η, and bathymetry b rescaled by a factor of 0.45 and shifted by 1.92. Top right panel: discharge qx. Bottom panel: second equilibrium variable E (2)x .

Fig. 10

 10 Fig. 10 Supercritical flow, test case from Section 5.3: convergence test. Left panel: non-WB scheme; we observe fifth-order accuracy. Right panel: WB scheme; we observe machine precision accuracy for the equilibrium variables qx and E

5 Fig. 11

 511 Fig. 11 Wave over a dry island test case from Section 5.4.1: η := h + b and b at different times.

Fig. 12

 12 Fig. 12 Tsunami on three obstacles, test case from Section 5.4.2: η : h + b and b at different times.

+ 0.2x, for x < 0, 1, for 0 ≤ x ≤ 3, 1 + 0.4(x -3), for x > 3,

Acknowledgements M. C. was funded by a postdoctoral fellowship at ENSAM. L. M. was funded by the Schweizerischer Nationalfonds zur Förderung der wissenchaftlichen Forschung (SNF) grant 200020 204917 "Structure preserving and fast methods for hyperbolic systems of conservation laws" and by a Postdoc Fellowship at NCSU. V. M.-D. acknowledges the support of ANR OptiTrust (ANR-22-CE25-0017). P. Ö. was supported by the German Research Foundation (DFG) within SPP 2410, project OE 661/5-1 (525866748) and under the personal grant 520756621 (OE 661/4-1). D. T. was funded by a SISSA Mathematical Fellowship. This work has been developed also in the context of the SHARK-FV conference.

All authors would like to thank Jonatan Núñez for sharing his high-order FV-WENO code on his repository [50]. We have started our work by adapting his code.