<figure><img src="image_1.png" /></figure>

# An operational discontinuous Galerkin shallow water model for coastal flood assessment

A.G. Filippini<sup>a,*</sup>, L. Arpaia<sup>b</sup>, V. Perrier<sup>c</sup>, R. Pedreros<sup>a</sup>, P. Bonneton<sup>d</sup>, D. Lannes<sup>e</sup>, F. Marche<sup>f</sup>, S. De Brye<sup>g</sup>, S. Delmas<sup>h</sup>, S. Lecacheux<sup>a</sup>, F. Boulahya<sup>a</sup>, M. Ricchiuto<sup>i</sup>

## ARTICLE INFO

### Keywords:

Discontinuous Galerkin

Coastal hazards

Shallow water equations

Well-balanced scheme

Wet/dry treatment

Entropy viscosity stabilization

## ABSTRACT

Hydrodynamic modeling for coastal flooding risk assessment is a highly relevant topic. Many operational tools available for this purpose use numerical techniques and implementation paradigms that reach their limits when confronted with modern requirements in terms of resolution and performances. In this work, we present a novel operational tool for coastal hazards predictions, currently employed by the BRGM agency (the French Geological Survey) to carry out its flooding hazard exposure studies and coastal risk prevention plans on International and French territories. The model, called UHAINA (wave in the Basque language), is based on an arbitrary high-order discontinuous Galerkin discretization of the nonlinear shallow water equations with SSP Runge-Kutta time stepping on unstructured triangular grids. It is built upon the finite element library AeroSol, which provides a modern C++ software architecture and high scalability, making it suitable for HPC applications. The paper provides a detailed development of the mathematical and numerical framework of the model, focusing on two key-ingredients : (i) a pragmatic $P^0$ treatment of the solution in partially dry cells which guarantees efficiently well-balancedness, positivity and mass conservation at any polynomial order; (ii) an artificial viscosity method based on the physical dissipation of the system of equations providing nonlinear stability for non-smooth solutions. A set of numerical validations on academic benchmarks is performed to highlight the efficiency of these approaches. Finally, UHAINA is applied on a real operational case of study, demonstrating very satisfactory results.

## 1. Introduction

Nowadays, coastal areas are inhabited by around 10% of the world population and are exploited by a huge number of economic activities (ports, industrial and energy installations, tourist activities, etc.). Extreme sea-state conditions can generate highly destructive floods resulting in human casualties and economic losses. Numerical modeling is a fundamental tool for assessing and forecasting coastal hazard and risk. Thanks to their simple structure the shallow water equations can

be a useful modeling tool and provide accurate flooding simulations, in particular when the mean water level exceeds the height of the crest of natural barriers or coastal defenses (overflow). Shallow water codes are well-established for quite some time now in the engineering community and the attention of the scientists is actually mainly focused on using these models to produce huge dataset of flooding simulation to feed either a probabilistic approach to risk evaluation, or warning and prediction systems. In this context, most of the tools fail in reaching

<footnote><sup>a</sup> BRGM, French Geological Survey, Coastal Risk and Climate Change, 3 Av. C. Guillermin 45060 Orleans, France</footnote>

<footnote><sup>b</sup> ISMAR, Istituto di Scienza Marine CNR, Arsenale Castello 2737/F 30122 Venezia, Italy</footnote>

<footnote><sup>c</sup> INRIA, Univ. Pau et Pays de l'Adour, Avenue de l'Universite 64000 Pau, France</footnote>

<footnote><sup>d</sup> EPOC, Univ. Bordeaux CNRS, Allee Geoffroy Saint-Hilaire 33615 Pessac, France</footnote>

<footnote><sup>e</sup> IMB, Univ. Bordeaux CNRS, 351 cours de la Liberation 33400 Talence, France</footnote>

<footnote><sup>f</sup> IMAG, Univ. Montpellier CNRS, Place Eugene Bataillon 34090 Montpellier, France</footnote>

<footnote><sup>g</sup> FLYING WHALES, Paris, France</footnote>

<footnote><sup>h</sup> SUEZ Smart Solutions Limited, Bordeaux, France</footnote>

<footnote><sup>i</sup> INRIA, Univ. Bordeaux, 200 Avenue de la Vieille Tour 33405 Talence, France</footnote>

<footnote>* Corresponding author.</footnote>

<footnote>E-mail addresses: a.filippini@brgm.fr (A.G. Filippini), luca.arpaia@ve.ismar.cnr.it (L. Arpaia), vincent.perrier@inria.fr (V. Perrier), r.pedreros@brgm.fr (R. Pedreros), philippe.bonneton@u-bordeaux.fr (P. Bonneton), david.annes@math.u-bordeaux.fr (D. Lannes), fmarche@math.univ-montp2.fr (F. Marche), sebastien.de.brye@gmail.com (S. De Brye), delmas.simon@gmail.com (S. Delmas), s.lecacheux@brgm.fr (S. Lecacheux), f.boulahya@brgm.fr (F. Boulahya), mario.ricchiuto@inria.fr (M. Ricchiuto).</footnote>
the modern requirements for high resolution, robustness and high performance computing.

When applied to hyperbolic conservation or balance laws, Discontinuous Galerkin (DG) schemes combine many favorable properties of Finite Volume (FV) and finite element methods (Cockburn and Shu, 1998a). Indeed, on one side one can exploit the very high order of accuracy (both in space and time) in the propagation of smooth waves (Giraldo et al., 2002; Wang et al., 2013; Mitchell, 2015; Xu et al., 2018). On the other hand, DG inherits many favorable properties advocated for finite volumes: local cell discrete conservation statements; geometrical flexibility and the possibility of using adaptive unstructured meshes; discontinuous representation of data, which can be exploited to capture bathymetric discontinuities. In addition, the DG scheme provides great parallel scaling properties due to its compact stencil when deployed even on hundreds or thousands of cores on a modern cluster (Salehipour et al., 2013; Cockburn and Shu, 2001; Kelly and Giraldo, 2012; Müller et al., 2019). This list of characteristics makes the DG schemes an advantageous choice for a modern numerical tool for coastal flooding assessment.

DG schemes have become very popular since the early 2000s, starting from the pioneering works of Cockburn and Shu (1998b, 2001). Their application to discretize the shallow water equations has been performed by many authors in the following years, as well on the cartesian (Li and Liu, 2001; Schwanenberg and Harms, 2004; Tassi et al., 2007) and on the spherical domain (Salehipour et al., 2013; Giraldo et al., 2002; Marras et al., 2015), involving arbitrary orders methods on unstructured triangulations (Eskilsson and Sherwin, 2006; Kubatko et al., 2006; Duran and Marche, 2014). The list of references is non-exhaustive.

More recently, several authors have focused on two issues, particularly relevant in flooding applications: the preservation of the motionless steady states over complex varying bathymetries (well-balanced property), and the preservation of the water height positivity in presence of wet/dry fronts, to properly handle flooding and drying events.

To satisfy the well-balanced property, DG schemes take advantage from a vast literature on well balanced high-order schemes (Bermudez and Vazquez, 1994; Greenberg and LeRoux, 1996; LeVeque, 1998; Jin, 2001; Audusse et al., 2004; Xing and Shu, 2005; Castro et al., 2008; Ricchiuto, 2011, 2015a; Berthon and Chalons, 2016; Cheng et al., 2019; Castro and Parés, 2020; Ciallella et al., 2023). Nowadays, very general methodologies exist to preserve arbitrary steady states while benefitting of the very high accuracy of DG schemes (Guerrero Fernández et al., 2022; Mantri et al., 2024). For motionless steady states, the well known hydrostatic reconstruction approach of Audusse et al. (2004) has been generalized in the DG context by Xing and Shu (2006) and Duran and Marche (2014).

Furthermore, wet/dry transitions introduce three distinct numerical difficulties: maintaining positivity of the water height, introduction of artificial free surface gradients generating unphysical waves at the shores, and numerical instabilities due to solution discontinuities. A popular method introduced by Xing et al. (2010) allows to maintain the positivity of cell-average water height by means of a limiter, which rescales the polynomial around the positive average. While this method is claimed to be well-balanced by the authors, it is not unconditionally well-balanced as wet/dry cells are neglected. Some ad-hoc numerical treatment have been developed for linear polynomials (Kesserwani et al., 2010; Bunya et al., 2009; Vater et al., 2019). The issue has been recently tackled by Bonev et al. (2018) on quadrilateral meshes, but their solution, based on finite difference gradient reconstruction on the partially dry cell, does not appear to be directly applicable on triangular meshes. More in general, it can be stated that no satisfactory solution still exist to overcome this problem in a way that does not cause the loss either of the high-accuracy of the scheme or of its mass conservation property. In this work we design, for partially dry cells, a first order scheme (finite volume) that ensures the preservation of steady state (well-balancing), while avoiding unphysical flow over

adverse slope issues (in the sense of Brufau et al. (2002)) and respecting the exact elevation of coastal defense structures, that primarily controls the extent of flooded areas.

As for most CFD problems, discontinuities are likely to develop in numerical hydrodynamics. Indeed, bore like wave profiles generally occur with the use of nonlinear shallow water (SW) equations and are associated to a significant dissipation of energy. Breaking fronts are usually numerically approximated by dissipative shocks. High order finite element methods are known to produce spurious oscillations around such discontinuities, raising a serious stability issue. A specific and adequate stabilization treatment is thus needed for wave breaking to be correctly represented in our discontinuous Galerkin scheme, with no numerical instabilities. In order to prevent the appearance of these non physical oscillations we propose here an entropy viscosity method, inspired by the work of Guermond et al. (2011). This technique is based on shock theory physics (Bonneton, 2007) and it is preferred rather than the use of a purely numerical shock limiter for stability. This method will also help in stabilize wet/dry moving interfaces, where local velocities are prone to grow dramatically when the water height gets very small.

This work combines and extends many of the above mentioned techniques in a single model for operational use in coastal hazard and coastal flooding risk assessment. For completeness, we mention that UHAINA is also formulated in spherical coordinates, embedding the Coriolis terms, as documented in Arpaia et al. (2022b). The same work also shows UHAINA applied at the large scale for tsunami and storm surge hazard predictions.

The paper is organized as follows. In Section 2 we will discuss the numerical scheme used to discretize the shallow water equations, including the main external forcings driving the mean sea water level in a context of coastal flooding (the bottom friction, the atmospheric forcing and the wave setup contributions). In Section 3.2 we will provide detailed information on the implemented wet/dry treatment. Section 4 will be devoted to the entropy based scheme stabilization in shocks. In Section 6 we will validate our model on standard benchmarks for the targeted applications and we will apply the model on a real challenging scenario by simulating the flooding generated by four historical storms inside the Arcachon lagoon (on the France Atlantic coast).

## 2. The shallow water equations for coastal hazard applications

In this work, we use the so-called pre-balanced formulation of the shallow water equation, introduced in Rogers et al. (2001). This set of equation, formulated in terms of free-surface rather than the water depth, is widely use within the FV framework (Kurganov and Levy, 2002; Kurganov and Petrova, 2007; Liang and Borthwick, 2009) and it has been also successfully applied in the DG context (Kesserwani et al., 2010; Duran and Marche, 2014) with the purpose of naturally balancing flux gradient and topography source term, using the approximate Riemann solver of Roe. These equations read as follows :

$$ \frac{\partial U}{\partial t} + \nabla \cdot F(U, b) = S(U, b) \qquad (1) $$

with

$$ U = \begin{pmatrix} \zeta \\ \frac{hu}{hv} \end{pmatrix}; \quad F(U, b) = \begin{pmatrix} hu \\ \frac{1}{2}g(\zeta^2 - 2\zeta b) \\ huv \end{pmatrix}; \quad \frac{hv}{hv^2 + \frac{1}{2}g(\zeta^2 - 2\zeta b)} $$

$$ S(U, b) = \begin{pmatrix} 0 \\ -g\zeta\partial_x b + s_x \\ -g\zeta\partial_y b + s_y \end{pmatrix} $$

being $\zeta$ the free surface elevation, $h = \zeta - b$ the total water depth, $b$ the bathymetry and $u$ the depth-averaged horizontal velocity vector of components $u$ and $v$ in the $x-$ and $y-$ direction respectively.

$s_x$ and $s_y$ represent the external forces that model bottom friction, atmospheric pressure and wind stress, all of paramount importance for
coastal flooding. In this work we use a Manning model for the bottom friction and a quadratic model for the wind stress. The source term then writes:

$$s_x = -C_F(h, \mathbf{u})\mathbf{hu} - \frac{h}{\rho_0}\partial_x P_{atm} + C_D \frac{\rho_{air}}{\rho_0}|u_{10}|u_{10}$$

$$s_y = -C_F(h, \mathbf{u})\mathbf{hv} - \frac{h}{\rho_0}\partial_y P_{atm} + C_D \frac{\rho_{air}}{\rho_0}|u_{10}|v_{10}$$

where $C_F$ is the Manning friction coefficient, $\rho_{air}$ and $\rho_0$ represent respectively the density of air and water, $P_{atm}$ stays for the atmospheric pressure, $u_{10} = [u_{10}, v_{10}]$ is the wind velocity field ten meters above the ground and $C_D$ the wind stress coefficient proposed by Charnock (1955). The value of the Charnock parameter is set to 0.022 after the calibration of Muller et al. (2014).

Besides the effects of wind and pressure, another key factor which increases the mean water level in coastal zones is short-waves breaking, an effect commonly referred to as wave setup. Previous studies revealed that wave-induced setup can be a significant contribution to the storm surge under energetic wave conditions and can even become the dominant one (Kennedy et al., 2012; Lerma et al., 2017; Pedreros et al., 2018; Lavaud et al., 2020). Few studies have shown that the depth-varying circulation in surf zones can increase the maximum wave setup at the coast (Apotsos et al., 2007; Guérin et al., 2018). However, 2DH models embedding a coupling for short waves can provide sufficiently accurate solutions in operational context compared to 3D models (as stressed in the work Lavaud et al., 2020). In order to account for wave-current interactions in our model, we adopt the extensively used approach of Mei (1989). We modify the source term adding the following terms:

$$+ \frac{1}{\rho_0} \nabla \cdot \mathbf{R} + \nabla \cdot \left( v_t \nabla(h\mathbf{u}) \right), \quad (2)$$

where $\mathbf{R}$ is the waves radiation stress tensor (Phillips, 1977) and the turbulent eddy viscosity $v_t$ associated with lateral mixing term (Reniers et al., 2004) is defined following the approach of (Battjes, 1975):

$$v_t = M h \left( \frac{D_w}{\rho_0} \right)^{1/3},$$

with $M$ a constant coefficient and $D_w$ the rate of energy loss due to wave breaking.

We compute $\mathbf{R}$ and $D_w$ by using the spectral wave model WAVE WATCH III (Tolman, 2014). This model solves the equations for the conservation of the wave action (Komen et al., 1994) to simulate the generation and propagation of wind-waves. Waves characteristics (significant wave height, wave direction, wavelength and waves-to-ocean energy flux) are then used to compute $D_w$ and the radiation stresses components for Eq. (2), following the relations of the linear theory (Mei, 1989; Dingemans, 1997) and without accounting for surface rollers to limit model complexity at this stage.

### 3. Well-balanced discontinuous Galerkin discretization

In this work, we discretize system (1) with an arbitrary high-order discontinuous Galerkin method. In order to do this, we take a computational domain partition $\mathcal{M}_h$ of non-overlapping elements $K$. The approximated vector solution $U_h$ and the bathymetry $b_h$ belong to a finite dimensional space $\mathcal{V}_h$ and are expressed as polynomials of order $k$ on each element $K$:

$$U_h(\mathbf{x}, t) = \sum_{i=1}^{N\mathrm{DOFK}} U_i(t) \varphi_i(\mathbf{x}), \quad b_h(\mathbf{x}, t) = \sum_{i=1}^{N\mathrm{DOFK}} b_i \varphi_i(\mathbf{x}) \quad (3)$$

where $\{\varphi_i\}_{i=1}^{N\mathrm{DOFK}}$ is a Lagrangian polynomial expansion basis, with $N\mathrm{DOFK} = (k+1)(k+2)/2$. A weak formulation of the problem is obtained by multiplying (1) by $\varphi_i$ and integrating the result on a given element $K$ of $\mathcal{M}_h$:

$$\int_K \frac{\partial U_h}{\partial t} \varphi_i(\mathbf{x}) d\mathbf{x} - \int_K F(U_h, b_h) \cdot \nabla \varphi_i(\mathbf{x}) d\mathbf{x} + \quad (4)$$

$$\int_{\partial K} F(U_h, b_h) \cdot \hat{\mathbf{n}} \varphi_i(\mathbf{x}) d\mathbf{x} = \int_K S_h(U_h, b_h) \varphi_i(\mathbf{x}) d\mathbf{x}$$

where $\partial K$ designate the element boundary and $\hat{\mathbf{n}}$ the outward boundary normal of unitary length. For the sake of simplicity in the notation, we will drop the subscript $h$ in the following paragraphs, but quantities should be understood as discrete ones unless explicitly stated. The integrals are computed with Gaussian quadrature of precision degree equal to $3k-1$ for the volume integrals and $3k$ for the boundary integrals. This allows to integrate exactly the hydrostatic terms which is a key ingredient to retain the well-balanced property (Xing et al., 2010).

The term $F(U, b) \cdot \hat{\mathbf{n}}$ in (4) is evaluated through a numerical flux $H(U^-, U^+, b^-, b^+)$, which depends on the reconstructed states $U^\pm, b^\pm$ at each side of the interface, denoting with $^-$ and $^+$ respectively the interior and exterior traces on a face, with respect to element $K$. The reconstruction process is also fundamental for the well-balanced property of the scheme. For this, we apply the hydrostatic reconstruction of Audusse et al. (2004) and Duran and Marche (2014). We first define a single valued term for the topography at the interface:

$$\check{b} = \max(b^-, b^+) \quad (5)$$

Then, we use the following relations:

$$\check{h}^- = \max(0, \zeta^- - \check{b}) \quad ; \quad \check{h}^+ = \max(0, \zeta^+ - \check{b}) \quad (6)$$

$$\check{\zeta}^- = \check{h}^- + \check{b} \quad ; \quad \check{\zeta}^+ = \check{h}^+ + \check{b} \quad (7)$$

$$\check{h}u^- = \check{h}^- u^- \quad ; \quad \check{h}u^+ = \check{h}^+ u^+ \quad (8)$$

The numerical flux $H$ is thus defined as:

$$H = H_{LF} (\check{U}^-, \check{U}^+, \check{b}) + H_{\sigma} (\check{U}^-, \check{U}^+, b^-, b^+, \check{b}) \quad (9)$$

where $H_{LF}$ is the global Lax-Friedrichs flux:

$$H_{LF} = \frac{1}{2} (F(\check{U}^-, \check{b}) + F(\check{U}^+, \check{b})) \cdot \hat{\mathbf{n}} - \frac{\alpha_{LF}}{2} (\check{U}^+ - \check{U}^-) \quad (10)$$

with $\alpha_{LF} = \max_{\partial K \in \Sigma} \max_{q \in \partial K}(\lambda_q)$ the maximum eigenvalue of the flux Jacobian $\partial F / \partial U$ among all the faces quadrature points $q$. $H_\sigma$ is a well-balanced correction term, which is necessary to obtain the numerical preservation of a lake at rest solution (C-property), in analogy with a finite volume scheme. The goal is to ensure that the flux and the source terms are numerically balanced, such as to preserve exactly the steady state solution ($\mathbf{u} = 0$ and $\zeta = \text{const}$) in presence of bathymetry variations. $H_\sigma$ reads:

$$H_{\sigma} = \begin{pmatrix} g \check{\zeta}^- & 0 & 0 \\ 0 & g \check{\zeta}^- & 0 \\ 0 & 0 & g \check{\zeta}^- \end{pmatrix} \cdot \hat{\mathbf{n}} \quad (11)$$

Still, more efforts are needed to preserve the steady state in presence of wet/dry fronts, and especially when using high order polynomial approximations. The numerical treatment of wet/dry regions will be discussed in detail in the paragraph 3.2.

#### 3.1. Positivity-preserving limiter

To guarantee positivity of the cell-averaged water depth, we apply the well-established limiter of Xing et al. (2010). This positivity-preserving limiter rescales the solution around the cell average water height, while preserving mass conservation, via the simple formula:

$$\tilde{\mathbf{W}} = \overline{\mathbf{W}} + \Theta (\mathbf{W} - \overline{\mathbf{W}}) \quad (12)$$

where $\mathbf{W} = [h, hu]$ is the vector of conservative variables and $\Theta$ is defined as:

$$\Theta = \min \left( 1, \frac{\overline{h}}{h_{min} - \overline{h}} \right) \quad (13)$$

with $\overline{h}$ the cell-averaged water depth and $h_{min} = \min_{x \in K} h(x)$ the minimum water depth over the element, computed taking into account
the values at the degrees of freedom and at the special set of Gauss-Lobatto quadrature points (Xing et al., 2010). The free surface elevation is then constructed via its definition: $\tilde{\zeta}(x) = \tilde{h}(x)+b(x)$. The so-computed variables $\tilde{U}^n$ are then used instead of $U^n$ into the DG scheme (4). It can be noted that this limiter does not take any effect in wet regions, being active only in nearly dry areas.

## 3.2. Wet/dry treatment

Handling wet/dry frontiers in a well-balanced manner can be particularly tricky especially for high order DG schemes, where high order continuous polynomials must be used to approximate a solution with discontinuous derivatives at the shoreline. First of all, in order to detect dry regions we introduce two threshold parameters $\epsilon_h^{wd}$ and $\epsilon_u^{wd}$, acting independently on water height and velocity respectively. $\epsilon_h^{wd}$ is used to avoid working with very small water depth quantities. If $h < \epsilon_h^{wd}$, the degree of freedom is considered as dry and $h = 0$ is set. $\epsilon_h^{wd}$ is small and variable with the mesh size. The second parameter $\epsilon_u^{wd}$ is used to avoid division by zero, when the velocity values are recovered from the momentum. If $h < \epsilon_u^{wd}$ in a degree of freedom, the velocity and the momentum are set to zero there. Note that $\epsilon_h^{wd} < \epsilon_u^{wd}$ should hold. In the context of discontinuous Galerkin discretization, the degrees of freedom inside an element can be all wet, all dry or a mix of the two. In the latter case, the element is said partially dry and it is crossed by the wet/dry frontier. When such case holds and high order polynomials are used to approximate the solution, negative water heights and artificial pressure gradients may appear, leading to instabilities. These are due to the discontinuity in the derivative of the state variables present inside the element. In order to recover both positivity and well-balancedness, additional corrections must be performed and will be discussed hereafter.

As pointed out by Bonev et al. (2018), the outcome $\tilde{U}^n$ of the positivity preserving limiter is not well-balanced. Fig. 1(a)-(b) depict in one dimension the lake at rest solution and its numerical approximation using a parabolic polynomial function ($k = 2$), in a partially dry cell K. As $h_{min} < 0$ in the cell, the positivity limiter (12) activates, computing a new $\tilde{h}(x)$ solution over the cell (Fig. 1(c)) which satisfies the positivity property. As consequence, a non-physical sea surface gradient is generated in the cell, together with non-zero fluxes at the two interfaces. This globally leads to spurious waves generation at the shores which pollute the entire solution. An ad-hoc computation of the free-surface gradient may be deployed in partially dry cells in order not to spoil the well-balanced property in presence of dry areas (the so-called “dry lake at rest”): a correction of the bathymetry (Brufau et al., 2002, 2004) in order to recover a flat free-surface gradient, or a computation of the free-surface gradient only using the wet degrees of freedom (Bonev et al., 2018).

### 3.2.1. P⁰ at wet/dry cells

In this work we use a simplified treatment which consists in dropping the polynomials order to zero order in wet/dry cells, thus reverting the scheme to first-order Finite Volume on such cells. Although the global order of the solution is spoiled and we loose a finer scale representation of the bathymetry, this allows to verify in a simple way both well-balancing and positivity in presence of dry area. This correction is needed only for cells with adverse slope where a non-physical slope in the free surface is introduced by the polynomial representation. These particular configuration of wet/dry cells over an adverse slope are detected as:

$$K = K_{P0} \quad \text{if} \quad \forall j = 1, N \text{ DOF}_K \quad \exists j, h_j < \epsilon_h^{wd} \text{ and } b_j + \epsilon_h^{wd} \ge \zeta_w^{max} \quad (14)$$

with $\zeta_w^{max} = \max_{x_i \in DOF_K, h_i > \epsilon_h^{wd}} \zeta_i$ the maximum free-surface of wet nodes on the cell itself. Namely, for a wet/dry cell with adverse slope we set:

$$\varphi_i = 1, \Theta = 0 \quad \forall i = 1, N \text{ DOF}_K \quad \text{if} \quad K = K_{P0}$$

Thus, we impose:

$$\zeta_j^* = \zeta_w^{max} \quad \forall j = 1, N \text{ DOF}_K \quad \text{if} \quad K = K_{P0} \quad (15)$$

and we compute a constant bathymetry that conserves the mass, see Fig. 1(d):

$$b_j^* = \zeta_w^{max} - \bar{h} \quad \forall j = 1, N \text{ DOF}_K \quad \text{if} \quad K = K_{P0} \quad (16)$$

With such treatment, the wet/dry frontier matches the edges separating wet P⁰ cells from adjacent dry ones. It is worth to mentioned that, every time a new solution has been computed, the original bathymetry representation is reintroduced in all cells and (14) has to be evaluated again.

### 3.2.2. Numerical flux correction at wet/dry edges

In addition to the polynomial degree limiting described in the previous section, on partially dry cells we also modify the right and left states used to compute the numerical fluxes across two cells as follows (Audusse et al., 2004):

$$\tilde{b}^- = b^* - \max(0, b^* - \zeta^{*-}) \quad (17)$$

$$\tilde{h}^- = \max(0, \zeta^{*-} - b^*) \quad ; \quad \tilde{h}^+ = \max(0, \zeta^+ - b^*)$$

$$\tilde{\zeta}^- = \tilde{h}^- + \tilde{b}^- \quad ; \quad \tilde{\zeta}^+ = \tilde{h}^+ + \tilde{b}^- \quad (18)$$

$$\tilde{h}u^- = \tilde{h}^- u^- \quad ; \quad \tilde{h}u^+ = \tilde{h}^+ u^+ \quad (19)$$

It is easy to check that such corrections, when applied at a wet/dry edge, allow to:

* nullify both free-surface and bathymetry jumps, balancing the right-hand-side in the “dry lake at rest” case.

* nullify advective fluxes (if presents), ensuring a zero flux condition.

The latter point is particularly relevant for overflowing applications. In fact, apart from the preservation of the lake at rest state and positive water depth, another important physical constraint is preserved at a discrete level. The wet/dry front can only advance when the water depth gradient in the wet/dry cell is larger than the bathymetry gradient. This numerical artifact is also known in the literature as flow over adverse slope effect (Brufau et al., 2002) and, by enforcing such constraint, we avoid non-physical overtopping. In a DG method this is translated into a zero numerical advective flux if at a wet/dry edge (e.g. with $K^-$ wet and $K^+$ dry) the free-surface on the wet side is lower then the bathymetry on the dry side, $\eta_h^- < b_h^+$.

Further to this point, in flooding applications, a crucial aspect is the preservation of the coastal defenses height (like dykes) in the discrete bathymetry $b_h$: coupled with the aforementioned flow over adverse slope condition, it allows an accurate computation of the flooded area. The DG method is well suited to this end, since a discontinuous representation of the bathymetry allows to maintain, in presence of the corrections (16) and (17), the original height value on the dry side of the edge. This may not work properly in situations in which both the elements adjacent to an edge are partially dry cells. This may occur when a defense structure is represented in the model by a single high point in the bathymetry (see en example in Fig. 2(a)), but the same can also quite often occurre as a transitory situation in a flooding process. In this case, in fact, (16) and (15) will produce the numerical representation described in Fig. 2(b), allowing water to be transferred from one side of the dyke to the other, even if the sea water level on both is below the structure crest. To overcome this problem, dry faces are identified at the beginning of any time step computation. Then, the numerical flux across the edge is computed as in 3.2.2, by considering the internal state obtained by the polynomial degree limiting, while keeping the dry state values for the external state (namely $\zeta^+ = b^+$ and $h^+u^+ = 0$).
<figure><img src="image_2.png" /><figcaption>Fig. 1. (a) One dimensional sketch of an element containing a wet/dry interface in the case of a constant water free surface; (b) quadratic DG polynomial representation inside the element; (c) effect of the application of the positivity limiting on the partially dry cell; (d) result of the polynomial degree limiting approach: first order polynomial correction for well-balancedness.</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>Fig. 2. (a) One dimensional sketch of a situation in which a dyke is surrounded by wet/dry cells; (b) numerical representation resulting from the P⁰ treatment of wet/dry cells.</figcaption></figure>

### 3.3. Time discretization

High-order temporal accuracy is achieved with a Strong-Stability-Preserving Runge-Kutta method (Gottlieb et al., 2001) up to third order, while standard Runge-Kutta scheme is used to go beyond. Time step is computed according to:

$$ \Delta t = CFL \min_K \left( \frac{\rho_K}{2 \alpha_{LFK}} \right) \qquad (20) $$

with $\rho_K$ the element inscribed circle radius and $\alpha_{LFK}$ the local Lax-Friedrich parameter. For linear stability reason, the CFL value is set to $\frac{1}{2k+1}$ (being k the polynomial order of the approximation).

In order to avoid very small time-step, the bottom friction term has been implemented with a fractional step method after the RK update (4). A first order implicit time scheme is used for the friction step:

$$ hu_j^{n+1} = \frac{hu_j^*}{1 + \Delta t C_F (h_j^{n+1}, u_j^*)} \qquad j = 1, N \text{ DOF}_K \qquad (21) $$

where $hu_j^*$ is the momentum computed at the jth degree of freedom after evolving (4) with RK.

### 4. Entropy based stabilization

As recall by Pasquetti (2017), the SW system admits a convex entropy E which satisfies the following relation:

$$ \frac{\partial E}{\partial t} + \nabla \cdot \mathbf{G} \le 0 \qquad (22) $$

with E = $\frac{h}{2}u^2 + \frac{g}{2}h^2 + ghb$ actually the shallow water energy and G = (E + $\frac{g}{2}h^2$) u the energy flux. The energy satisfies the exact conservation

equation (22) in smooth regions, while satisfying the strong inequality in shocks. As a direct consequence, the energy dissipation D can then be easily calculated from the ensuing entropy residual equation:

$$ -D = \frac{\partial E}{\partial t} + \nabla \cdot \mathbf{G} \qquad (23) $$

We should expect D ≈ 0 in smooth regions and, locally, D ≥ 0 in discontinuities, which makes it an efficient criterion to detect shocks by the way and a useful indicator for breaking fronts tracking.

The idea is to stabilize the numerical scheme by defining an artificial viscosity $\mu$ based on the residual of Eq. (23): $v = v(D)$. Following Guermond et al. (2011), a viscosity term is then added to the
hyperbolic system (1):

$$ \frac{\partial U}{\partial t} + \nabla \cdot F(U, b) = S(U, b) + \nabla \cdot (\nu \nabla U) \quad (24) $$

For the discretization of this additional diffusion term, only the cell integral, resulting from the integration by parts applied to the variational formulation of Eq. (24), will be kept. Indeed, there is no concern about the numerical consistency of this term as we are essentially interested in adding some dissipation to the system. The final form of the discretized system will then result in:

$$ \int_K \frac{\partial U_h}{\partial t} \varphi_h(x) \, dx - \int_K F_h(U_h, b_h) \cdot \nabla \varphi_h(x) \, dx + \int_{\partial K} F_h(U_h, b_h) \cdot \hat{n} \, \varphi_h(x) \, dx = \int_K S_h(U_h, b_h) \varphi_h(x) \, dx - \nu_K(D_h) \times \int_K \nabla U_h \cdot \nabla \varphi_h(x) \, dx \quad (25) $$

### 4.0.1. Entropy based viscosity computation

Based on Guermond et al. (2011) we compute the artificial viscosity $\nu$ from the energy dissipation $D$ with the following formula for the element $K$:

$$ \nu_K = \kappa \alpha \nu_{\max} \min \left( 1, \frac{\delta x |\bar{D}_K|}{\tau D_{th}} \right) \quad (26) $$

where $\kappa(k)$, $\alpha$ and $\tau$ are calibration parameters of the method. $\delta x$ is the mesh size of the cell $K$, $\bar{D}_K$ is the averaged value of energy dissipation over the cell. $\nu_{\max}$ is given by:

$$ \nu_{\max} = \delta x \max_K \left( |\mathbf{u} \cdot \mathbf{n}| + \sqrt{g\bar{h}} \right) \quad (27) $$

$D_{th}$ corresponds to a scaling factor. The idea behind formula (26) relies on willing to impose some diffusion to stabilize the solution only in shock regions. The minimum is used to control the maximum amount of dissipation that is injected to the solution.

Ideally we would like to introduce at most an amount of dissipation corresponding to the numerical diffusion inherent to a first order scheme, which is known to be always stable in shocks. Then, to calibrate $\alpha$, the scaling dissipation $D_{th}$ was initially set close to 0 and $\alpha$ was increased till removing all spurious oscillations, until the value $\alpha = \frac{1}{2}$ was selected.

In case of breaking waves, the scaling dissipation $D_{th}$ is given by the 1D water wave shock theory (Bonneton, 2007):

$$ D_{th} = \frac{g}{4} \left( \frac{g(h_1 + h_2)}{2h_1h_2} \right)^{\frac{1}{2}} (h_2 - h_1)^3 \quad (28) $$

By defining $\bar{h} = \frac{h_1+h_2}{2}$, $H = h_2 - h_1$ and assuming $H \sim \gamma \bar{h}$, expression (28) can be simplified and reduced to:

$$ D_{th} \sim \beta_{th} g^{\frac{3}{2}} \bar{h}^{\frac{5}{2}} \quad (29) $$

with $\beta_{th} = \frac{\tau^3}{2\sqrt{4-\gamma^2}}$. The resulting implemented formula will thus read:

$$ \nu_K = \kappa \alpha \nu_{\max} \min \left( 1, \frac{\delta x |\bar{D}_K|}{\beta g^{\frac{3}{2}} \bar{h}^{\frac{5}{2}}} \right) \quad (30) $$

with $\beta = \tau \beta_{th}$. The parameter $\beta$ is controlled by the theoretical parameter $\beta_{th}$. Usually, for breaking waves in field conditions $\gamma \in [0.6, 0.8]$, so that $\beta_{th} \in [0.057, 0.14]$. In academic test cases, $\gamma$ is rather close to 0.6, hence the initially retained value is $\beta_{th} = 0.057$. Furthermore, the evaluation of $\tau$ was performed by reducing its value till removing unwanted oscillations around breaking wave fronts. After some tests, the value $\tau = \frac{1}{10}$ was observed to be the best compromise to get stabilized shocks with a minimal imposed viscosity. Finally, as $\nu$ should have a lower value for higher order methods, the value $\kappa = \frac{1}{k}$ was tested and gave satisfactory results in simulated test cases.

### 4.0.2. Solution of the energy equation

The last missing step for employing formula (26) is to compute $\bar{D}_K$ via (23). With this in mind, we first compute the entropy $E_h$ on our $k$th order DG nodal basis on element $K$, via its definition (22). For this, we have to project the energy expression on our $k$th order DG nodal basis on cell $K$ and the operation requires a mass matrix inversion. Since we are interested in a cell-averaged value of the dissipation we use the following cell-centered finite volume scheme for the energy equation (23):

$$ \bar{D}_K = -\frac{1}{|K|} \left( \int_K \frac{\partial E_h}{\partial t} \, dx + \int_{\partial K} \mathbf{G}_h \cdot \hat{\mathbf{n}} \, dx \right) \quad (31) $$

A Lax-Friedrich numerical flux $\mathbf{G}_h$ is again used to evaluate $\mathbf{G}_h \cdot \hat{\mathbf{n}}$ on element faces:

$$ \mathbf{G}_h = \frac{1}{2} (\mathbf{G}_h^{-} + \mathbf{G}_h^{+}) - \alpha_{LF} (E^{+} - E^{-}) \quad (32) $$

Finally, the unstationary term $\frac{\partial E_h}{\partial t}$ is computed by means of a $(k+1)$th order Backward Differentiation Formula (BDF) time scheme, in accord with the $k$th polynomial approximation order of the DG spatial discretization used. Of course, this means that lower order BDF schemes are needed for the start up in the first time steps. $\bar{D}_K$ is computed once per time step and is then used to compute the solution at each intermediate stage of the $(k+1)$th Runge-Kutta time scheme.

## 5. Implementation in the AeroSol library

UHAINA relies on a few libraries developed at the INRIA center, such as the C++ finite element library AeroSol (AeroSol, 2011), the tools for parallel mesh partitioning PaMPA (PaMPA, 2017) and SCOTCH (SCOTCH, 2009). These libraries provide a modern environment to work with finite element schemes and to handle data management for distributed memory parallel computation. AeroSol is a library dedicated to the solution of advection-diffusion problems with arbitrary high-order finite element discretizations on hybrid meshes (triangles and quadrilaterals). It provides an advanced programming environment designed by the purpose of performance and HPC. It employs runtime polymorphism with interface classes to switch between different numerical and physical combinations, e.g. physical models, numerical fluxes or time integrators for a given spatial operator. Among the AeroSol classes that are mandatory, we mention a Mesh class that collects all the cells and the cell boundary faces along with geometric and topological properties of the mesh, e.g. how are the cells connected and where are their vertices. A Finite Element class describes the properties of the finite element space on the unit cell. This includes, for example, how many degrees of freedom are present and where they are located as well as the values and gradients of shape functions. The Quadrature class takes care of the loops over quadrature points, storing the location and the weights of quadrature points on the unit cell. A Model class is dedicated to the implementation of the advective and diffusive fluxes and of the source term; a NumericalFlux class implements the formula for numerical flux. Finally a Scheme class implements the spatial operator: it performs all the loops over the cells and cell boundaries, it computes the integrals associated to the various finite element terms and it inverts the finite element mass-matrix.

UHAINA consists of several files containing derived classes that implement the specific choices described in the previous sections by overriding methods for the Model, the NumericalFlux and of course for the Scheme class. We mention that, in the latter, both the interpolations at quadrature points and the integrals are computed by local matrix-matrix products. The steady matrices that only depends on stationary data (base functions and geometrical transformations for the quadrature formulas), are precomputed at initialization time leaving only the task of performing the matrix-matrix product into the residual computation.
<figure><img src="image_4.png" /><figcaption>Fig. 3. Well-Balancing: 3D view of the lake at rest solution over nonlinear bathymetry in presence of beaches.</figcaption></figure>

Finally we mention the AeroSol class for input/output, the IO class, which allows a standardization of the template for input/output files. UHAINA makes use of the UGRID format, which stores unstructured (or flexible mesh) model data in the Unidata Network Common Data Form (NetCDF) form, adding conventions for specifying the topology for unstructured (e.g. triangular) grids. NetCDF files are binary formats (machine-independent) that are commonly used in climatology, meteoro-logy, and oceanography applications. It is an input/output for many GIS applications and it is used for general scientific data exchange.

# 6. Results

Hereafter, errors with respect to a reference solution are computed with the following L2-norm definition:

$$ \|\epsilon(U)\|_{L_2} = \sqrt{\sum_K \int_K (U_h - U_{ref})^2 dx} $$

## 6.1. Well balanced property with wet/dry cells

A first important property to be satisfied by a hydrodynamic code for coastal flooding applications is the preservation of a motionless steady state, also known as C-property, in the presence of an irregular bottom and wet/dry cells. In this test case we deploy a nonlinear bathymetric profile, involving a bump and a hollow of the same dimensions and defined by the following analytical expression:

$$ b = (h_0 - 0.2) + \begin{cases} 0.5 \exp^{-25(x-1.2)^2-50(y-0.7)^2} & \text{if } (x > 0.68) \\ -0.5 \exp^{-50(x-0.45)^2-100(y-0.4)^2} & \text{else} \end{cases} \quad (33) $$

The computational domain is the rectangle [0, 2] × [0, 1] m, and we use an unstructured mesh of 1868 elements with a typical mesh size of 0.05 m. The initial solution is the lake at rest state defined by $\zeta = h_0$ and $u = 0$ m/s. Fig. 3 shows a 3D view of the solution for $h_0 = 0.2$ m.

The numerical scheme must be able to maintain a stable and unperturbed solution for as long as desired. Thanks to the use of the pre-balanced formulation of the shallow water equations, our numerical scheme is exactly well-balanced if the initial free surface level of the water plan is set to $\zeta_0 = 0$ m (that is the $h_0 = 0$ m case). By this, we mean that the computed error in this configuration is exactly zero because all terms contributing to the equilibrium balance are zero precisely for as long as wish. Our numerical strategy to handle wet/dry frontiers allows to extend this property even in presence of beaches and at any order. Moreover, this special configuration of $\zeta_0$ is not binding for applications, since a negative offset can always be applied to the bathymetry without impacting the simulation results.

In the case a different value for $h_0$ is used instead, not trivial integrals computations must balance to maintain the equilibrium state in (4). Exact integration of each term is thus mandatory, but it may not be sufficient to guarantee that the error will not grow over time. In order to check this possibility, we ran the code for 1400 s, which represents the time that a wave (with phase velocity $\sqrt{gh_0} = 1.4$ m/s)

<figure><img src="image_5.png" /><figcaption>Fig. 4. Well-Balancing: time-history of the $L_2$ errors on the momentum for all k-orders.</figcaption></figure>

would take to travel 1000 times through the computational domain. The total number of iterations to reach the final time was 504 000 for the P0 scheme, 1680 000 for P1, 2548 000 for P2 and 5 082 000 for P3, due to the different CFL restrictions. Table 1 reports the last values of the residuals computed by the code, together with the L2-norm of the solution errors with respect to the initial solution. All quantities are very close to the zero machine, but a certain growth happens, driven by the error on the momentum (see Fig. 4), albeit over very long time and it is not relevant for our applications. Indeed, we have tried to run an operational case on the Arcachon lagoon (more details will be discussed in Section 6.8) in the lake-at-rest configuration, observing only negligible errors for the usual duration of a storm (36 h): $< 1e-12$ on the free surface elevation and $< 1e-11$ on the momentum.

## 6.2. Flow over adverse slope

With the following test case, we show that our scheme does not produce numerical artifacts in wet/dry cells in presence of a non-zero flow and that it does not produce artificial overtopping of protection structures. To this end, we have constructed a very simple test case in a rectangular domain of dimension $[-80, 80] \times [0, 30]$ m with a linearly varying bathymetry along the x-coordinate:

$$ b = \begin{cases} -15 & x < 0 \\ -15 + 17/30 x & 0 \le x < 30 \\ 2 - 1/10 (x - 30) & 30 \le x < 40 \\ 1 & x > 40 \end{cases} \quad (34) $$

The bathymetry is characterized by a steep gradient close to the coastline, culminating in a dyke that protects a flat area inland. A cosinussoidal signal is imposed at the seawall boundary of the domain with amplitude $a = 2.5$ m and period $\tau = 1200$ s using a Dirichlet condition, while a wall condition is imposed to the other boundaries. The top-left picture of Fig. 5 gives a representation of the initial solution for this test case.

Figs. 5 and 6 show the evolution of the numerical solution computed using the polynomial approximations $P^1$ and $P^2$ respectively. In both plots, the magnitude of the velocity has been expanded by a factor 10, for sake of clarity. It should be noted that no unphysical flooding is generated in the two simulations, which corroborates our treatment of wet/dry interfaces, discussed in Sections 3.2 and 3.2.2. Our scheme respects the topographical dimension of the dyke, even in the challenging case where its crest is described by a single mesh node, when some other codes need instead to double the point introducing a very small element, with consequent impact on the computational time-step. This property is of the utmost importance for an operational flood simulation tool, as directly related to the volume of flooding water.
<table><caption>**Table 1**

Well-Balancing: Computed residuals and errors after *t* = 1400 s in the test case configuration case in which *h*₀ = 0.2 m.</caption><thead><tr><th rowspan="2">Approx.</th><th colspan="3">Residuals</th><th colspan="3">Errors</th></tr><tr><th>Mass</th><th>Momentum x</th><th>Momentum y</th><th>ζ</th><th>hu</th><th>hu</th></tr></thead><tbody><tr><td>P0</td><td>1.72e-17</td><td>6.94e-18</td><td>6.94e-18</td><td>3.16e-18</td><td>2.12e-16</td><td>4.69e-17</td></tr><tr><td>P1</td><td>2.02e-17</td><td>1.39e-17</td><td>1.04e-17</td><td>1.09e-15</td><td>1.13e-13</td><td>1.03e-13</td></tr><tr><td>P2</td><td>4.51e-17</td><td>1.73e-17</td><td>2.08e-17</td><td>7.22e-15</td><td>1.57e-12</td><td>1.41e-12</td></tr><tr><td>P3</td><td>8.57e-17</td><td>5.12e-17</td><td>3.51e-17</td><td>1.26e-16</td><td>3.39e-12</td><td>3.46e-12</td></tr></tbody></table>

<figure><img src="image_6.png" /><figcaption>**Fig. 5.** Flow over adverse slope: evolution of the solution computed using a $P^1$ polynomial approximation. On the left: 3D zoom around the dyke. On the right: 1D extraction along the center line of the computational domain. The y-axis of the plot has to be intended as m (on the left) for the free surface and bathymetry and as m/s (on the right) for the velocity.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>**Fig. 6.** Flow over adverse slope: evolution of the solution computed using a $P^2$ polynomial approximation. On the left: 3D zoom around the dyke. On the right: 1D extraction along the center line of the computational domain. The y-axis of the plot has to be intended as m (on the left) for the free surface and bathymetry and as m/s (on the right) for the velocity.</figcaption></figure>

Indeed, it can be appreciated that flooding only occurs after the water level has reached the height of the defense structure and, similarly, a thin film of water remains trapped on land after the sea water retired. This ability of the code is exploited in practical situations, as shown hereafter in the study on the Arcachon's lagoon 6.8.

### 6.3. Carrier and Greenspan's nonbreaking sloping planar beach problem

This test case, introduced in Carrier and Greenspan (1958), describes the physical process in which a monochromatic wave is let
<figure><img src="image_8.png" /><figcaption>Fig. 7. Carrier and Greenspan’s nonbreaking sloping planar beach problem: computed free surface elevation at $3T+dt$, $dt \in \{T/5, 2T/5, 3T/5, 4T/5, T\}$ (left) and shoreline evolution (right) for the fourth order scheme on the mesh with $\Delta x = 0.125$ m.</figcaption></figure>

run-up and run-down on a plane beach. The periodic wave of dimensionless amplitude $\tilde{A}$ and frequency $\tilde{\omega}$ travels shoreward on a plane beach without breaking and is than reflected out to sea, generating a standing wave in the nearshore region. The case reproduces interesting conditions to assess the robustness of the scheme’s wet/dry treatment. In Carrier and Greenspan (1958), the authors derived an analytical solution in the context of the shallow water equations. The solution makes use of two dimensionless variables (denoted with tilde notation in the following) $\tilde{\sigma}$ and $\tilde{\lambda}$ defined as:

$$ \tilde{\lambda} = 2(\tilde{t} + \tilde{u}) \quad \text{and} \quad \tilde{\sigma} = 2\sqrt{\tilde{\xi} - \tilde{x}} \qquad (35) $$

Being $L$ the typical length scale of the problem and $\alpha$ the beach slope, nondimensional variables are computed by:

$$ \tilde{x} = x/L, \quad \tilde{\xi} = \zeta/(\alpha L), \quad \tilde{u} = u/\sqrt{g\alpha L}, \quad \tilde{t} = t/\sqrt{L/g\alpha} \qquad (36) $$

The analytical solution is then formulated as follows :

$$ \left\{ \begin{array}{l} \tilde{u} = -\frac{\tilde{A}J_1(\tilde{\sigma}) \sin(\tilde{\lambda})}{\tilde{\sigma}} \\ \tilde{\xi} = \frac{\tilde{A}}{4} J_0(\tilde{\sigma}) \cos(\tilde{\lambda}) - \frac{\tilde{u}^2}{4} \end{array} \right. \qquad (37) $$

where $J_0$ and $J_1$ stand for the Bessel functions of zero and first order. We consider here the solution obtained by setting $\tilde{A} = 0.6$, with length scale $L = 20$ m and beach slope $\alpha = 1/30$. The system solution for $t = 0$ s is used as initial solution, while the boundary signal is given also by solving (35) at the seaward boundary point. We let the solution oscillate for three full periods and then look at the solution at times $3T + dt$, for $dt \in \{T/5, 2T/5, 3T/5, 4T/5, T\}$. On Fig. 7 we show the results provided by a simulation that uses the fourth order scheme ($k = 3$) and a mesh with resolution $\Delta x = 0.125$ m. On Fig. 7 (left) the computed free surface elevation is compared against the analytical solution provided by system (35). On the right, analogous comparison is made in therm of the shoreline time evolution. Overall, the obtained results are satisfying and attest to the effectiveness of the wet/dry technique employed in this work with arbitrary order of polynomial approximation.

To push further our analysis, the analytical solution provided by system (35) is used to perform a mesh convergence study and to assess hp-convergence of the scheme in the runup. Fig. 8 shows the $L_2$ norms of the error obtained by halved four times the mesh size of the computation. As expected, due to the P0-approximation in the wet/dry elements, all schemes exhibits first order slopes, but the computed error is progressively lower as $k$ increases, confirming the benefit of using high order polynomials.

## 6.4. Solitary wave run-up on a planar beach

Next, we check the ability of the model to reproduce breaking and run-up with a benchmark by the NOAA Center for Tsunami Research:

the solitary wave-run up on a plane beach. Synolakis (1987) carried out laboratory experiments for incident solitary waves of multiple relative amplitudes over a planar beach. Detailed description of the test case, along with the initial conditions, can be found for example in Filippini et al. (2016), Li and Liu (2001), Delis et al. (2008), Gottlieb et al. (2001), Castro et al. (2008) among many others. Here we consider a wave characterized by the ratio amplitude over depth of $a/h = 0.28$, with still water depth $h = 1$ m. This wave breaks strongly both in the run-up and the run-down phases of the motion. For this reason it is often used to validate wave-breaking strategies inside Boussinesq codes (Filippini et al., 2016; Bonneton et al., 2010; Chazel et al., 2011; Kazolea and Ricchiuto, 2018) and many others. If dispersive properties would be required to well represent wave propagation and shoaling, the shallow water equations still provide an accurate representation of the subsequent wave breaking, run-up and run-down processes. The computational domain is $x \in [-20, 100]$ m with $\Delta x = 0.2$ m and a wall boundary condition imposed at both grid extremes. We perform the test by using four different polynomial degrees in the DG finite element representation, from P0 (constant) to P3 (cubic). The artificial viscosity parameter has been tuned by using the value $\beta_{th} = 0.057$, according to the wave breaking conditions of the experiment. Finally, according to the literature, a Manning coefficient of $0.01$ s/m$^{1/3}$ is used to define the glass surface roughness.

Fig. 9 compares the free surface results obtained with different degrees of polynomial approximation and against the data provided by the laboratory experiment of Synolakis (1987) at different times. The image clearly shows that all the schemes perform in good agreement with the measurements and are comparable with the results found in the literature (for example in Delis et al. (2008) and Li and Raichlen (2002)). In general, the first order P0 scheme is much more dissipative than the other three, as expected, while P2 and P3 produce results very close to each other, except in the run up process. It can also be noted that the wet/dry mechanism, described in Section 3.2, works very efficiently regardless of the order of the schemes. In addition, the entropy viscosity based shock limiter 4 proves to be very effective in avoiding spurious oscillations in the solution without dissipating too much of both the shock that originates in the run-up (around $t\sqrt{g/h} = 20$) and in the run-down (around $t\sqrt{g/h} = 60$).

## 6.5. Runup on a conical island

Another common case of study used in the literature to assess the quality of shallow water and Boussinesq models is the solitary wave runup over a conical island (see Nikolos and Delis, 2009; Kazolea et al., 2012; Ricchiuto, 2015b among many others). It is based on the laboratory experiments performed by Liu et al. (1995). It consists in a solitary wave traveling over a 7.2 m diameter conical island. The island is 0.625 m tall with 1/4 side slopes and it is positioned in the center
<figure><img src="image_9.png" /><figcaption>Fig. 8. Carrier and Greenspan's nonbreaking sloping planar beach problem: h-convergence (left) and p-convergence (right).</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Fig. 9. Solitary wave run-up on a planar beach: computed free surface elevation at different increasing times of the simulation (from top-left to bottom-right) for different DG approximations.</figcaption></figure>

of a 30 m wide by 25 m long flat bottom basin. Fig. 10 (left) shows the computational grid used for this test, which is refined close to the cone and characterized by a mesh size spacing from 0.5 m to 0.075 m. A solitary wave with $a/h_0 = 0.2$ is initially placed at $x_s = 0$ m and propagates rightward in a water depth of $h_0 = 0.32$ m. The solitary wave shape imposed at the left boundary of the domain is defined by the free surface perturbation :

$$d\zeta = a h_0 \operatorname{sech}((x - x_s)\sqrt{3a/(4h_0^3)}) \quad (38)$$

with a corresponding initial velocity derived by the linearized shallow water theory $u = (\sqrt{g/h_0}, 0)d\zeta$. Fig. 10 (right) shows snapshots of the computed free surface solution at four different times. The incoming

solitary wave shoals and breaks on the front side of the island ($t = 5.6$ s), in correspondence of gauge 9. Then, secondary waves run around the cone and meet behind it giving a rear side run up ($t = 9.68$ s), while a circular reflected wave moves away from the island. Finally, the rear peak splits again into smaller waves propagating back around the island, while the major wave front moves away from it.

The simulation has been performed with different degrees of the DG polynomials, from $P^0$ (constant) to $P^3$ (cubic). Due to the presence of a breaking wave front, the artificial viscosity term was activated and tuned with $\beta_{th} = 0.0015$, based on the $\gamma$ (cf Section 4.0.1) computed from the experimental records. The results are compared in terms of free surface elevation at four wave gauges: $g_6 = (-3.6, 0)$ m, $g_9 =$
<figure><img src="image_11.png" /><figcaption>Fig. 10. Conical island: view of the mesh grid (left) and snapshots of the solution at different simulation time (right).</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>Fig. 11. Conical island: gauge time series of the simulated free surface using different polynomial approximation orders.</figcaption></figure>

($-2.6, 0$) m, $g_{16} = (0, -2.58)$ m and $g_{22} = (2.6, 0)$ m, located as reported in Fig. 10 (left). In general, the physics of the test is well represented and the results are comparable with those of other authors in the literature using shallow water codes (Nikolos and Delis, 2009; Ricchiuto, 2015b). The oscillations recorded in the backwash phase, which are not reproduced by the model, are driven by wave dispersion and require the use of a Boussinesq model in order to be caught. The four pictures of Fig. 11 clearly show that better accuracy in the representation of peaks and troughs in all stations is achieved by increasing the polynomial degree of the simulation, in particular at gauge 22, situated on the back of the island.

## 6.6. Okushiri tsunami experiment

To further confirm the soundness of the wetting/drying procedure and of our shock limiting technique, we consider the second benchmark of the third international workshop on long wave runup models: Tsunami runup onto a complex three-dimensional beach. The case is thoroughly described on the web pages (Benchmark, 2004; Tsunami) and in Liu et al. (2008) to which we refer for details. The test is a scaled down laboratory reproduction of the tsunami wave that hit the Okushiri island in Japan in 1993. The bathymetry data consist in a 400 times scaled down geometry represented in Fig. 12 (left). The site is
<figure><img src="image_13.png" /><figcaption>**Fig. 12.** The Okushiri tsunami experiment : mesh view with gauges positions (left); 3D sketch of the bathymetry with contour lines of the 0 m and 32 m levels. Dashed white line highlights the area of the Monai valley, when maximum runup is expected.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>**Fig. 13.** The Okushiri tsunami experiment : 3d views of the solution at time $t = 14.5$ s, $t = 16.5$ s and $t = 17.5$ s using the second order scheme.</figcaption></figure>

characterized by a small 10 m height island (Muen island), positioned in front of the coast, and by a steep and narrow valley (the Monai valley), where the maximum runup of 32 m has been observed in the historical event (highlighted by the white dashed circle in the picture). The computational domain is a rectangular bassin $[0, 5.448] \times [0, 3.402]$ m and it has been discretized with the mesh shown in Fig. 12 (right). The mesh is adapted to the bathymetry variations, it contains 42 892 elements with maximum and minimum mesh sizes given roughly by 0.1 m and 0.025 m. Please note that the mesh size recommended in Benchmark (2004), Tsunami and Liu et al. (2008) for this test is of 0.014 m which would give a mesh with approximately 189 000 triangles in absence of mesh adaptation, roughly four times more than in the mesh used here. The wave signal displayed in Fig. 14(a) is introduced from the left boundary of the domain, while wall boundary conditions are imposed to the three other boundaries. During the experiment, probes have been set to measure the water height history in three locations shown in Fig. 12.

Fig. 13 (left) shows the initial water retirement preceding the arrival of the tsunami. After hitting the beach, the tsunami wave reflects, and a large wave travels toward the right to hit the steepest slopes in the region of the Monai valley Fig. 13 (right). The reflected wave eventually reaches and inundates the Muen island. After the first reflection at the coast, the computational domain is crossed by multiple breaking fronts of significative height. To run the simulations, the entropy based shock limiter has been tuned using the value $\beta_{th} = 0.08$, based on the theoretical $\gamma$ computed from the greatest recorded shock by the gauges.

The comparison of the computed sea surface height with the experimental data provided is shown on Fig. 14(b), (c) and (d). Simulations have been run with $k$ from 0 to 3, but $k = 2$ is not reported in the figures because little differences where observed with respect to the P3 computation. The agreement between measured and computed heights is satisfactory, with shocks fronts being well captured by the scheme. The P3 simulation gives the closer results to the data, with the highest peaks and steepest shock profiles.
<figure><img src="image_15.png" /><figcaption>Fig. 14. The Okushiri tsunami experiment: inlet boundary wave (a) and comparisons with experiments in gauges 5 (b), 7 (c) and 9 (d).</figcaption></figure>

As a last verification, Fig. 15 shows the maximum runup plot obtained by the three different simulations. As a reference, we have reported in the picture the contour lines corresponding to the bathymetry elevation at 0 m and 32 m, in dashed green and purple respectively. The picture clearly show that only the P3 simulation obtain an accurate results in the maximum runup region. In particular, the P0 and P1 runs underestimate the maximum runup by roughly 8 m, while the P3 run overshoots the line of the 32 m giving a more conservative prediction of the maximum runup of about 38 m. As already remarked, these results were obtained by employing a grid with four times less elements than a uniform mesh following the prescriptions of Benchmark (2004), Tsunami and Liu et al. (2008) would contain, stressing the advantage on the use of an adaptive unstructured mesh.

## 6.7. The SeaSide experiment

The following benchmark simulates a tsunami wave striking an urban area. The test case reproduces a laboratory experiment carried out at the Oregon State University and presented in Park et al. (2013). The experiment takes place in a rectangular basin 48.8 m long, 26.5 m wide and 2.1 m deep, involving a complex topography that includes a seawall and several blocks, inspired by the actual buildings of the Seaside city, Oregon, at a scale 1 : 50 (Fig. 16 left). The incoming wave enters from the left boundary and wall boundary conditions are imposed in the remaining boundaries of the computational domain. In our case, the interest in performing this test is double: to test the model behavior in a quasi-operational context with a large mesh; to evaluate the effect of increasing the formal approximation order of the scheme in flooding applications on irregular topography. An unstructured numerical grid of 172 854 elements is used, refined in the

<figure><img src="image_16.png" /><figcaption>Fig. 15. The Okushiri tsunami experiment: maximum runup plot.</figcaption></figure>

area where the buildings are placed and with a minimum reference size of $h_{ref} = 0.03$ m. The individual structures and buildings are approximated as steep-sided topography and not by imposing some wall boundary conditions to the scheme. For the breaking front, the value of $\beta_{th} = 0.0785$ has been used, computed from the data recorded by the gauges. Finally, a Manning coefficient of $0.01 \text{ s/m}^{1/3}$ is also used.
<figure><img src="image_17.png" /><figcaption>Fig. 16. Seaside: Contour lines of the topography (left) and gauge measurement locations (right).</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>Fig. 17. Seaside: 3D view snapshots of the flooding evolution for increasing times of the simulation.</figcaption></figure>

Fig. 17 shows 3D snapshots of the computed water field. The incoming tsunami wave advances and breaks on the city's seafront (Fig. 17 (top-left)). Flooding occurs, while a part of the wave is reflected backward towards the sea. The water is slowed by buildings and accelerates in the streets and empty spaces perpendicular to the coastline (Fig. 17 (top-right)). The flow, then, reaches the end of the domain, and is reflected by the basin walls (Fig. 17 (bottom-left)). A secondary front is formed, which floods the city from inland (Fig. 17 (bottom-right)).

The results of the simulations can be compared with those of the experiment in terms of free surface elevation and velocity at a large number of locations (16 right). Fig. 18 shows the free surface evolution for $P^1$ and $P^2$ simulations at some of the gauges of the onshore region.

As the tsunami wave approaches, it undergoes wave shoaling and breaking before rushing into the city. The shallow water equations are inadequate to reproduce such phenomena driven by dispersion, for which Boussinesq-type models, with an appropriate breaking strategy, should be used instead. For this reason the first peak of the incoming wave is underestimated at gauge A1 and B1, on the sea front. Despite this, it can be noted that the model performs well overall and that the wet/dry treatment is robust, whatever the order of the polynomial approximation used. Our results are good compared to those performed

by other authors in the literature (Park et al., 2013; Kazolea et al., 2023). The results provided by the $P^2$ approximation globally fit better the experimental data, especially for the most onshore gauges, albeit with a large increase in the computational cost of the simulation: 3 h for the $P^1$ case and 11 h for $P^2$, both on 128 cores.

## 6.8. Example of operational application: sea surface level simulation on the Arcachon lagoon

Our model has been employed by the BRGM agency under operational conditions to study flooding hazard on real sites. To show the full potential of UHAINA in such conditions, we present here the simulation of the sea state and storm surge generated by the most impacting recent storms on the Arcachon lagoon (south-west of France) (Lecacheux et al., 2023b,a). The modeling strategy is based on chaining the WW3 spectral wave model and the UHAINA hydrodynamic model, as described in Section 2. This chain simulates the propagation of waves and water levels from the open sea to the interior of the lagoon, the generation of wave setup, the wind-induced tilting of the free surface, and finally the overflow and propagation of water on land. Both models operate on the same unstructured computational grid, covering the entire Arcachon lagoon and part of ocean in front of the lagoon entry, describing a
<figure><img src="image_19.png" /><figcaption>Fig. 18. Seaside : Free surface elevation measured in some gauges locations.</figcaption></figure>

<figure><img src="image_20.png" /><figcaption>Fig. 19. View of the computational domain and of the mesh.</figcaption></figure>

semi-circle with a 17 km radius (see Fig. 19). The mesh has 426 000 elements with a variable size ranging from 600 m resolution offshore to 10 m resolution on land. In the lagoon the grid resolution is driven by the bathymetry gradient and hessian, with a particular focus on the two channels that connect the basin to the ocean. In Section 6.2, we have assessed the code's capacity to model water flow over coastal defenses represented by a single high point. Here, we exploit this property by making extensive use of constraint lines in order to force the edge of the mesh triangles to follow the main topo-bathymetric features, that are likely to influence the water flow. These include coastal structures, main roads and rivers or channels (see mesh zooms on the right of Fig. 19). In this way, the model respects the presence of obstacles on the ground and the value of their elevation (extracted from the data), but with a sufficiently coarse mesh to limit the computational cost.

Four historical storms have been simulated: Klaus (2009), Xynthia (2010), Emma (2018) and Justine (2021). These four events were chosen for the variety of situations they represent : high variability in water level components (tidal regime, atmospheric surge and wave setup). Moreover, observations were available to determine the model forcing and validate the results of the simulations. The Arcachon lagoon is characterized by a great temporal variability in the shape of the inlet channels and in the bathymetry of the shallow areas, due to erosion-sedimentation processes. The shape of the entrance paths to the lagoon, in particular, is the main driver for tidal propagation and distorsion, therefore it must be described as accurately as possible. As

a consequence, two different bathymetric datasets have been used to simulate Klaus and Xynthia on one side and Emma and Justine on the other (Lecacheux et al., 2023b,a). Concerning the bottom friction, a heterogeneous grid for the Manning coefficient was constructed, based upon the land coverage. Its values range from 0.02 s/m¹/³ offshore up to 0.33 s/m¹/³ in high-density urbanized areas.

Given the large computational cost of employing high order polynomials and the limited impact on the results of the realistic flood test, the simulations were performed with the $P^1$ scheme. Seaward boundary conditions consist of water level time series for the UHAINA model and wave conditions for WW3. Wind forcing is then applied homogeneously within the computational domain of both models, while WW3 also receive the information of a variable water level within the domain from a preliminary UHAINA computation without waves. The time series of the water level imposed at the sea boundary are shown in Fig. 20 (red lines) for the four storms played. These signals are obtained by superposing the tidal level and the storm surge. The tidal level has been computed using FES2014, to which an offset has been applied to account for the difference between the mean water level and the zero level of the topography. The storm surge component of the offshore water level is the result of large-scale regional simulations (Pedreros and Paris, 2012 for the Xynthia case) or reanalysis (Arpege model (Météo-France) for the other storms). Informations on waves are issued from the wave buoy located off the cape of “Cape-Ferret” (https://candhis.cerema.fr). Finally, wind observations from the “Cape-Ferret”
<figure><img src="image_21.png" /><figcaption>Fig. 20. Comparison of simulated water level with respect to tide gauge signal at the Arcachon station for 4 historical events.</figcaption></figure>

meteorological station (provided by Météo-France) are used to define
the wind time-series for each historical event.

For the four storms, the simulation results are compared with the
measurements from the Eyrac tide gauge at the entrance to the lagoon
(cf. Fig. 19). Fig. 20 shows that, for the four storms, the observed levels
are very close to the simulated ones, with peak errors of less than 5 cm.
A comparison of the blue and red lines in the images shows that the
distortion and amplification of the water signal, as it enters the lagoon,
are also well reproduced.

A detailed mapping of the flooded areas is only available for
the Xynthia storm, enabling a comparison with the simulations (see
Fig. 21). In general, the results are very satisfactory, with the simulated
limits close to those observed on the field. As the maximum sea level
remains relatively moderate (with a return period of around ten years),
the areas affected are mainly the lowest-lying areas around the oyster-
farming and port districts. There are, however, a few discrepancies,
highlighted in the picture with letters (a) and (b). In (a), the seafront
neighborhood of Andernos-les-Bains village is not flooded in the simu-
lations, although it was particularly affected in 2010 by overtopping
waves. The present paper only focuses on the overflowing type of
coastal flooding and the representation of flooding by overtopping is
left for further studies to come, since it will require wave dispersion
effects. It is also worth to mention that the digital terrain model and the
mesh take into account the new seafront developments, built after 2010
to provide better protection against extreme events. As a consequence,
even overflowing flooding in this sector would be described differently
respect to historical observations before 2010. In (b), the observed
flooded area in the ornithological reserve sector is noticeably larger
than that computed by the model. This difference can be explained by
the fact that river discharge in the Leyre river has not been taken into
account in the simulation, while having certainly been one of the main
drivers for flooding in the area.

The 4 historical simulations just shown were used in Lecacheux
et al. (2023b,a) to validate the model on the site. Subsequently, the
model was used to produce an atlas of 300 simulations of variable water

level, waves characteristics and wind intensity for forcasting purposes,
showing the robustness of the solutions proposed in this paper.

6.8.1. *Code scaling*

The previous case on the Arcachon lagoon, is used here to evaluate
the performance of our model when used in parallel. We perform a
scaling test by keeping constant the size of the problem (namely the
total number of degrees of freedom times the number of the problem
unknowns) and linearly increasing the number of computation cores
used. The test has been performed by simulating 5000 time-steps with
the same spatial grid described in the previous section. Different degree
of polynomial approximation have been used, corresponding to 426 000
DOFs in the P0 case, 1.28 millions DOFs for P1, 2.56 millions DOFs for
P2 and 3.84 millions DOFs for P3.

Fig. 22 shows, on a log–log scale, the results of the scaling analysis
in terms (a) of the measured speedup and (b) of the computed simulated
time per core-hours as function of the Number of Cores (NC). For
NC = 1, we have assumed the CPU time to be eight times slower
than the elapsed time for NC = 8. The authors benefitted from the
use of the cluster at the Centre de Calcul Scientifique en région Centre-
Val de Loire (Cascimodot). The nodes used for the computations are
bi-processors AMD Epyc 7702 2 GHz, 128 Gb RAM, that is to say
a NUMA (Non-Uniform Memory Access) architecture. Each processor
consists of 64 cores. The measured scalability is very satisfying. A loss
of performance must be highlighted when passing from NC = 64 to
NC = 128, where only an 80% speedup is measured (possibly due to
the computational architecture based on processors with 64 cores each;
multicore scaling on NUMA hardware Mercier and Jeannot, 2011; Chai
et al., 2006). Moreover we observe that the code scalability increases
with k, probably due to a better ratio between the number of operations
and the time for communications. The main factor driving these good
performances is the very compact communication stencil required by
a DG discretization, as outlined by other authors (Salehipour et al.,
2013; Cockburn and Shu, 2001; Kelly and Giraldo, 2012; Müller et al.,
2019). In a modern context in which high resolution and high number
<figure><img src="image_22.png" /><figcaption>**Fig. 21.** Comparison of the submerged areas in the simulations (blue shading) with respect to the ground surveys (black shading). In the picture, $H_{max}$ is the maximum water height (in m) computed during the simulation.</figcaption></figure>

<figure><img src="image_23.png" /><figcaption>**Fig. 22.** Scaling performance of the model.</figcaption></figure>

of simulations are the key ingredients to determine coastal hazards and risks, an unstructured model with good scaling properties appears as mandatory choice among others in the literature.

## 7. Conclusions

In this work we have presented a new arbitrary high-order tool for coastal hazards predictions on practical scenarios. Our tool is based on a DG discretization of the shallow water equations in a pre-balanced formulation. Particular attention in this work has been paid to the numerical treatment of partially dry elements, developing a pragmatic treatment of the wet/dry interface. The algorithm, based upon a polynomial degree limiting of the solution, has been proved to be capable to ensure well-balancedness, robustness together with the respect of

topographic features (coastal defenses, roads, channels, rivers, etc.) at an arbitrary order of global DG representation. These are features of fundamental importance for an operational tool. At the same time, the $P^0$ description on partially dry cells limits the beneficial impact on the solution of using high-order polynomials. It can be improved, for example, by making use of finite volume subcells as proposed in Meister and Ortleb (2016) and more recently in Haidar et al. (2022).

In practical flooding applications, the mesh size is constrained by the need to represent the main topographic features driving the flooding. Therefore, the use of linear polynomials with finer meshes seems the best compromise between accuracy and an affordable computational cost. With such limitation high order polynomials could become effective when coupled with subcell resolution (Haidar et al., 2022; Gaburro and Dumbser, 2021). Some p-adaptation strategy, may also allow to better leverage the use of high order polynomial. Certainly, this
will be most effective when combined with some adaptation strategy based on mesh movement (Arpaia and Ricchiuto, 2018, 2020; Arpaia et al., 2022a), remeshing (Wallwork et al., 2020; Barral et al., 2017; Arpaia et al., 2022a).

A novel approach to ensure the stability of the scheme in nearshore regions and on shocks by an artificial viscosity term has been proposed. The value of the viscosity coefficient is directly linked to the physical dissipation of the equations, scaled on the amount of dissipation provided by the 1D shallow water shock theory. Some calibration parameters have been introduced to give the user more control over the amount of dissipation injected into the scheme, but numerical tests on academic benchmarks have shown that the theoretical value performs well in different context and at all *k*-orders.

The overall numerical scheme has been validated on several academic test cases from the literature, showing its capability in dealing with several numerical issues at an arbitrary high-order of DG polynomial representation and providing globally good results. UHAINA also performed very satisfactorily in flooding simulations on a real scenario and, proving the robustness of the solutions proposed in this article, it has already been used to produce an atlas of 300 simulations of variable water level, waves characteristics and wind intensity for forecasting purposes (Lecacheux et al., 2023b,a).

The performances of the model have also been investigated, showing good scaling at all *k*-orders. However, there is room for improvement, for example by at the same time avoiding the computation on dry elements and dynamically re-allocating the CPU cores on the wet part of the domain, following the shoreline movement. In addition, optimized time schemes for DG, as the ones proposed in Kubatko et al. (2014), can be used instead of the classic explicit SSP Runge-Kutta schemes to relax the CFL condition and reduce the CPU time. Finally, the extension of this work in the context of the Boussinesq-type models, for an exhaustive modeling of coastal flooding, which also includes overtopping waves, is left for further studies and developments.

### CRediT authorship contribution statement

A.G. Filippini: Writing – original draft, Software, Methodology, Conceptualization. L. Arpaia: Writing – original draft, Software, Methodology, Conceptualization. V. Perrier: Software. R. Pedreros: Validation, Supervision, Funding acquisition. P. Bonneton: Writing – review & editing, Funding acquisition, Conceptualization. D. Lannes: Writing – review & editing, Conceptualization. F. Marche: Conceptualization. S. De Brye: Software. S. Delmas: Software. S. Lecacheux: Validation, Funding acquisition. F. Boulahya: Software. M. Ricchiuto: Writing – original draft, Supervision, Funding acquisition, Conceptualization.

### Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

### Data availability

Data will be made available on request.

### Acknowledgments

This work has been financially supported by Région Nouvelle Aquitaine, France, IDEX Université de Bordeaux, France, CNRS, France, INRIA, and Inter-CARNOT, France. Other sources of funding have been the Carib-Coast INTERREG, France project (https://www.interreg-car.aibes.fr/caribcoast), the RDI Gironde mission, France and R/D BRGM funds, France. All the developments presented have been implemented in the platform UHAINA co-developed by EPOC, IMAG, IMB, INRIA and BRGM. We are very grateful to Météo-France, SHOM and SIBA for providing us with the data used for simulating historical storm flooding in the Arcachon lagoon (i.e. topo-bathymetry, meteorological and oceanic measurements, flooded land extension).

### References

AeroSol: https://team.inria.fr/cardamom/aerosol/.

Aptsos, A., Raubenheimer, B., Elgar, S., Guza, R.T., Smith, J.A., 2007. Effects of wave rollers and bottom stress on wave setup. J. Geophys. Res.: Oceans 112 (C2).

Arpaia, L., Beaugendre, H., Cirrottola, L., Froehly, A., Lorini, M., Nouveau, L., Ricchiuto, M., 2022a. H- and r-adaptation on simplicial meshes using MMG tools. In: Sevilla, R., Perotto, S., Morgan, K. (Eds.), Mesh Generation and Adaptation: Cutting-Edge Techniques. Springer International Publishing, Cham, pp. 183–208.

Arpaia, L., Ricchiuto, M., 2018. R-adaptation for shallow water flows: conservation, well-balancedness, efficiency. Comput. & Fluids 160, 175–203, URL https://www.sciencedirect.com/science/article/pii/S0045793017303900.

Arpaia, L., Ricchiuto, M., 2020. Well balanced residual distribution for the ALE spherical shallow water equations on moving adaptive meshes. J. Comput. Phys. 405, 109173, URL https://www.sciencedirect.com/science/article/pii/S0021999119308782.

Arpaia, L., Ricchiuto, M., Filippini, A.G., Pedreros, R., 2022b. An efficient covariant frame for the spherical shallow water equations: Well balanced DG approximation and application to tsunami and storm surge. Ocean Model. 169, 101915.

Audusse, E., Bouchut, F., Bristeau, M.O., Klein, R., Perthame, B., 2004. A fast and stable well-balanced scheme with hydrostatic reconstruction for shallow water flows. SIAM J. Sci. Comput. 25, 2050–2065.

Barral, N., Olivier, G., Alauzet, F., 2017. Time-accurate anisotropic mesh adaptation for three-dimensional time-dependent problems with body-fitted moving geometries. J. Comput. Phys. 331, 157–187, URL https://www.sciencedirect.com/science/article/pii/S0021999116306155.

Battjes, J.A., 1975. Modelling of turbulence in the surfzone. In: Proceedings of 2nd Annual Symposium on WaterWays, Harbours and Coastal Eng. ASCE on Modeling Techniques. pp. 357–375.

Benchmark problem 2, Tsunami runup onto a complex three-dimensional beach, The Third International Workshop on Long-Wave Runup Models, URL https://nctr.pmel.noaa.gov/benchmark/Field/Field_Okushiri/.

Bermudez, A., Vazquez, M.E., 1994. Upwind methods for hyperbolic conservation laws with source terms. Comput. & Fluids 23 (8), 1049–1071.

Berthon, Christophe, Chalons, Christophe, 2016. A fully well-balanced, positive and entropy-satisfying Godunov-type method for the shallow-water equations. Math. Comp. 85 (299), 1281–1307.

Bonev, B., Hesthaven, J.S., Giraldo, F.X., Kopera, M.A., 2018. Discontinuous Galerkin scheme for the spherical shallow water equations with applications to tsunami modeling and prediction. J. Comput. Phys. 362, 425–448.

Bonneton, P., 2007. Modelling of periodic wave transformation in the inner surf zone. Ocean Eng. 34, 1459–1471.

Bonneton, P., Barthelemy, E., Carter, J.D., Chazel, F., Cienfuegos, R., Lannes, D., Marche, F., Tissier, M., 2010. Fully nonlinear weakly dispersive modelling of wave transformation, breaking and runup. arXiv preprint arXiv:1004.3480.

Brufau, P., García-Navarro, P., Vázquez-Cendón, M.E., 2004. Zero mass error using unsteady wetting-drying conditions in shallow flows over dry irregular topography. Internat. J. Numer. Methods Fluids 45 (10), 1047–1082.

Brufau, P., Vázquez-Cendón, M.E., García-Navarro, P., 2002. A numerical model for the flooding and drying of irregular domains. Internat. J. Numer. Methods Fluids 39 (3), 247–275.

Bunya, S., Kubatko, E.J., Westerink, J.J., Dawson, C., 2009. A wetting and drying treatment for the Runge-Kutta discontinuous Galerkin solution to the shallow water equations. Comput. Methods Appl. Mech. Engrg. 198 (17–20), 1548–1562.

Carrier, G.F., Greenspan, H.P., 1958. Water waves of finite amplitude on a sloping beach. J. Fluid Mech. 4 (1), 97–109.

Castro, Ma., Gallardo, J.M., López-García, J.A., Parés, C., 2008. Well-balanced high order extensions of Godunov's method for semilinear balance laws. SIAM J. Numer. Anal. 46 (2), 1012–1039.

Castro, Manuel J., Parés, Carlos, 2020. Well-balanced high-order finite volume methods for systems of balance laws. J. Sci. Comput. 82 (2), 48.

Chai, L., Hartono, A., Panda, D.K., 2006. Designing high performance and scalable MPI intra-node communication support for clusters. pp. 1–10.

Charnock, H., 1955. Wind stress on a water surface. Q. J. R. Meteorol. Soc. (81), 639–640.

Chazel, F., Lannes, D., Marche, F., 2011. Numerical simulation of strongly nonlinear and dispersive waves using a Green–Naghdi model. J. Sci. Comput. 48, 105–116.

Cheng, Y., Chertock, A., Herty, M., Kurganov, A., Wu, T., 2019. A new approach for designing moving-water equilibria preserving schemes for the shallow water equations. J. Sci. Comput. 80, 538–554.

Ciallella, Mirco, Torlo, Davide, Ricchiuto, Mario, 2023. Arbitrary high order WENO finite volume scheme with flux globalization for moving equilibria preservation. J. Sci. Comput. 96 (2), 53.

Cockburn, B., Shu, C.W., 1998a. The local discontinuous Galerkin method for time-dependent convection-diffusion systems. SIAM J. Numer. Anal. 35 (6), 2440–2463.

Cockburn, B., Shu, C.W., 1998b. The Runge-Kutta discontinuous Galerkin method for conservation laws V: multidimensional systems. J. Comput. Phys. 141 (2), 199–224.

Cockburn, B., Shu, C.W., 2001. Runge-Kutta discontinuous Galerkin methods for convection-dominated problems. J. Sci. Comput. 16, 173–261.
Delis, A.I., Kazolea, M., Kampanis, N.A., 2008. A robust high-resolution finite volume scheme for the simulation of long waves over complex domains. Internat. J. Numer. Methods Fluids 56 (4), 419-452.

Dingemans, M., 1997. Water wave propagation over uneven bottoms.

Duran, A., Marche, F., 2014. Recent advances on the discontinuous Galerkin method for shallow water equations with topography source terms. Comput. & Fluids 101, 88–104.

Eskilsson, C., Sherwin, S.J., 2006. Spectral/hp discontinuous Galerkin methods for modelling 2D Boussinesq equations. J. Comput. Phys. 212 (2), 566–589.

Filippini, A.G., Kazolea, M., Ricchiuto, M., 2016. A flexible genuinely nonlinear approach for nonlinear wave propagation, breaking and run-up. J. Comput. Phys. 310, 381–417.

Gaburro, E., Dumbser, M., 2021. A posteriori subcell finite volume limiter for general PNPM schemes: Applications from gasdynamics to relativistic magnetohydrodynamics. J. Sci. Comput. 86 (3), 1–41.

Giraldo, F.X., Hesthaven, J.S., Warburton, T., 2002. Nodal high-order discontinuous Galerkin methods for the spherical shallow water equations. J. Comput. Phys. 181 (2), 499–525.

Gottlieb, S., Shu, C.W., Tadmor, E., 2001. Strong stability-preserving high-order time discretization methods. SIAM Rev. 43, 89–112.

Greenberg, J.M., LeRoux, A.Y., 1996. A well-balanced scheme for the numerical processing of source terms in hyperbolic equations. SIAM J. Numer. Anal. 33 (1), 1–16.

Guérin, T., Bertin, X., Coulombier, T., De Bakker, A., 2018. Impacts of wave-induced circulation in the surf zone on wave setup. Ocean Model. 123, 86–97.

Guermond, J.L., Pasquetti, R., Popov, B., 2011. Entropy viscosity method for nonlinear conservation laws. J. Comput. Phys. 230, 4248–4267.

Guerrero Fernández, Ernesto, Escalante, Cipriano, Castro Díaz, Manuel J., 2022. Well-balanced high-order discontinuous Galerkin methods for systems of balance laws. Mathematics 10 (1).

Haidar, A., Marche, F., Vilar, F., 2022. A posteriori finite-volume local subcell correction of high-order discontinuous Galerkin schemes for the nonlinear shallow-water equations. J. Comput. Phys. 452, 110902.

Jin, S., 2001. A steady-state capturing method for hyperbolic systems with geometrical source terms. ESAIM Math. Model. Numer. Anal. 35 (4), 631–645.

Kazolea, M., Delis, A.I., Nikolos, I.K., Synolakis, C.E., 2012. An unstructured finite volume numerical scheme for extended 2D Boussinesq-type equations. Coast. Eng. 69, 42–66.

Kazolea, M., Filippini, A.G., Ricchiuto, M., 2023. Low dispersion finite volume/element discretization of the enhanced Green–Naghdi equations for wave propagation, breaking and runup on unstructured meshes. Ocean Mod. 182.

Kazolea, M., Ricchiuto, M., 2018. On wave breaking for Boussinesq-type models. Ocean Model. 123, 16–39.

Kelly, J.F., Giraldo, F.X., 2012. Continuous and discontinuous Galerkin methods for a scalable three-dimensional nonhydrostatic atmospheric model: Limited-area mode. J. Comput. Phys. 231 (24), 7988–8008.

Kennedy, A.B., Westerink, J.J., Smith, J.M., Hope, M.E., Hartman, M., et al., 2012. Tropical cyclone inundation potential on the Hawaiian islands of Oahu and Kauai. Ocean Mod. 52.

Kesserwani, G., Liang, Q., Vazquez, J., Mosé, R., 2010. Well-balancing issues related to the RKDG2 scheme for the shallow water equations. Int. J. Numer. Methods Fluids 62 (4), 428–448.

Komen, G.J., Cavaleri, L., Donelan, M., Hasselmann, K., Hasselmann, S., Jansen, P.A.E.M., 1994. Dynamics and Modelling of Ocean Waves. Cambridge University Press, Cambridge, U.K.

Kubatko, E.J., Westerink, J.J., Dawson, C., 2006. Hp discontinuous Galerkin methods for advection dominated problems in shallow water flow. Comput. Methods Appl. Mech. Engrg. 196 (1), 437–451.

Kubatko, E.J., Yeager, B.A., Ketcheson, D.I., 2014. Optimal strong-stability-preserving Runge–Kutta time discretizations for discontinuous Galerkin methods. J. Sci. Comput. 60, 313–344.

Kurganov, A., Levy, D., 2002. Central-upwind schemes for the Saint-Venant system. ESAIM Math. Model. Numer. Anal. 36 (3), 397–425.

Kurganov, A., Petrova, G., 2007. A second-order well-balanced positivity preserving central-upwind scheme for the Saint-Venant system. Commun. Math. Sci. 5 (1), 133–160.

Lavaud, L., Bertin, X., Martins, K., Arnaud, G., Bouin, M., 2020. The contribution of short-wave breaking to storm surges: The case Klaus in the Southern Bay of Biscay. Ocean Mod. 156.

Lecacheux, S., Filippini, A.G., Pedreros, R., Rohmer, J., Arpaia, L., Louisior, J., Nicolae Lerma, A., Idier, D., Quique, R., Paradis, D., Rose, F., 2023a. Maps of potential coastal flood zones on the arcachon lagoon for crisis management support. Law Hum. Behav..

Lecacheux, S., Pedreros, R., Filippini, A.G., Arpaia, L., Rohmer, J., Louisior, J., Nicolae Lerma, A., Quique, R., Pragout, J., 2023b. Appui à la Mission RDI de Gironde : Elaboration de Cartographies ZIP Submersion Marine Sur le Bassin d'Arcachon. Rapport BRGM/RP-73272-FR.

Lerma, A.N., Pedreros, R., Robinet, A., Senechal, N., 2017. Simulating wave setup and runup during storm conditions on a complex barred beach. Coast. Eng. 123.

LeVeque, R.J., 1998. Balancing source terms and flux gradients in high-resolution godunov methods: the quasi-steady wave-propagation algorithm. J. Comput. Phys. 146 (1), 346–365.

Li, H., Liu, R., 2001. The discontinuous Galerkin finite element method for the 2D shallow water equations. Math. Comput. Simul. 56 (3), 223–233.

Li, Y., Raichlen, F., 2002. Non-breaking and breaking solitary wave run-up. J. Fluid Mech. 456, 295–318.

Liang, Q., Borthwick, A.G.L., 2009. Adaptive quadtree simulation of shallow flows with wet-dry fronts over complex topography. Comput. & Fluids 38, 221–234.

Liu, P.L.F., Cho, Y.S., Briggs, M.J., Kanoglu, U., Synolakis, C.E., 1995. Runup of solitary waves on a circular island. J. Fluid Mech. 302, 259–285.

Liu, P.L.F., Yeh, H.H., Synolakis, C., 2008. Advanced Numerical Models for Simulating Tsunami Waves and Runup. Vol. 10.

Mantri, Yogiraj, Öffner, Philipp, Ricchiuto, Mario, 2024. Fully well-balanced entropy controlled discontinuous Galerkin spectral element method for shallow water flows: Global flux quadrature and cell entropy correction. J. Comput. Phys. 498, 112673, URL https://www.sciencedirect.com/science/article/pii/S0021999123007684.

Marras, S., Kopera, M.A., Giraldo, F.X., 2015. Simulation of shallow-water jets with a unified element-based continuous/discontinuous Galerkin model with grid flexibility on the sphere. Q. J. R. Meteorol. Soc. 141 (690), 1727–1739.

Mei, C.C., 1989. Applied Dynamics of Ocean Surface Waves. World Scientific.

Meister, A., Ortleb, S., 2016. A positivity preserving and well-balanced DG scheme using finite volume subcells in almost dry regions. Appl. Math. Comput. 272, 259–273.

Mercier, G., Jeannot, E., 2011. Improving MPI applications performance on multicore clusters with rank reordering. In: Recent Advances in the Message Passing Interface: 18th European MPI Users' Group Meeting, EuroMPI 2011, Santorini, Greece, September 18-21, 2011. Proceedings 18. pp. 39–49.

Mitchell, William F., 2015. How high a degree is high enough for high order finite elements? Procedia Comput. Sci. 51, 246–255, URL https://www.sciencedirect.com/science/article/pii/S1877050915010431. International Conference On Computational Science, ICCS 2015.

Müller, A., Kopera, M.A., Marras, S., Wilcox, L.C., Isaac, T., Giraldo, F.X., 2019. Strong scaling for numerical weather prediction at petascale with the atmospheric model NUMA. Int. J. High Perform. Comput. Appl. 33 (2), 411–426.

Muller, H., Pineau-Guillou, L., Idier, D., Ardhuin, F., 2014. Atmospheric storm surge modeling methodology along the French (Atlantic and English channel) coast. Ocean Dyn. (64(11)), 1671–1692.

Nikolos, I.K., Delis, A.I., 2009. An unstructured node-centered finite volume scheme for shallow water flows with wet/dry fronts over complex topography. Comput. Methods Appl. Mech. Eng. 198 (47–48), 3723–3750.

PaMPA : https://project.inria.fr/pampa/ (Release v2.0.0).

Park, H., Cox, D.T., Lynett, P.J., Wiebe, D.M., Shin, S., 2013. Tsunami inundation modeling in constructed environments: A physical and numerical comparison of free-surface elevation, velocity, and momentum flux. Coast. Eng. 79, 9–21.

Pasquetti, R., 2017. Viscous stabilizations for high order approximations of Saint-Venant and Boussinesq flows.

Pedreros, R., Idier, D., Muller, H., Lecacheux, S., Paris, F., et al., 2018. Relative contribution of wave setup to the storm surge: observations and modeling based analysis in open and protected environments (Truc Vert beach and Tubuai island). J. Coast. Res. 85.

Pedreros, R., Paris, F., 2012. Expertise Xynthia. Rapport BRGM/RP-61078-FR, p. 70.

Phillips, O.M., 1977. The dynamics of the upper ocean, 2ndedition. J. Fluid Mech. 88 (4), 793–796.

Reniers, A.J.H.M., Roelvink, J.A., Thornton, E.B., 2004. Morphodynamic modeling of an embayed beach under wave group forcing. J. Geophys. Res. 109, C01030.

Ricchiuto, M., 2011. On the C-property and generalized C-property of residual distribution for the shallow water equations. J. Sci. Comput. 48, 304–318.

Ricchiuto, M., 2015a. An explicit residual based approach for shallow water flows. J. Comput. Phys. 80, 306–344.

Ricchiuto, M., 2015b. An explicit residual based approach for shallow water flows. J. Comput. Phys. 280, 306–344.

Rogers, B., Fujihara, M., Borthwick, A.G.L., 2001. Adaptive Q-tree Godunov-type scheme for shallow water equations. Internat. J. Numer. Methods Fluids 35 (3), 247–280.

Salehipour, H., Stuhne, G.R., Peltier, W.R., 2013. A higher order discontinuous Galerkin, global shallow water model: Global ocean tides and aquaplanet benchmarks. Ocean Mod. 69, 93–107.

Schwanenberg, D., Harms, M., 2004. Discontinuous Galerkin finite-element method for transcrirical two-dimensional shallow water flows. J. Hydraul. Eng. 130 (5), 412–421.

SCOTCH : https://gitlab.inria.fr/scotch/scotch.

Synolakis, C.E., 1987. The runup of solitary waves. J. Fluid Mech. 185, 523–545.

Tassi, P.A., Bokhove, O., Vionnet, C.A., 2007. Space discontinuous Galerkin method for shallow water flows—kinetic and HLLC flux, and potential vorticity generation. Adv. Water Resour. 30 (4), 998–1015.

Tolman, H.L., 2014. User Manual and System Documentation of WAVEWATCH III Version 4.18. Technical Report NOAA / NWS / NCEP / MMAB Technical Note, 316.

Tsunami runup onto a complex three-dimensional beach; Monai Valley, Benchmarks of the NOAA Center for Tsunami Research, URL https://nctr.pmel.noaa.gov/benchmark/Laboratory/Laboratory_MonaiValley/index.html.
Vater, S., Beisiegel, N., Behrens, J., 2019. A limiter-based well-balanced discontinuous Galerkin method for shallow-water flows with wetting and drying: Triangular grids. Internat. J. Numer. Methods Fluids 91 (8), 395–418.

Wallwork, J.G., Barral, N., Kramer, S.C., Ham, D.A., Piggott, M.D., 2020. Goal-oriented error estimation and mesh adaptation for shallow water modelling. SN Appl. Sci. 2.

Wang, Z.J., Fidkowski, K., Abgral, R., Bassi, F., Caraeni, D., Cary, A., Decon-
inck, H., Hartmann, R., Hillewaert, K., Huynh, H.T., Kroll, N., May, G.,
Persson, P.O., van Leer, B., Visbal, M., 2013. High-order CFD methods: current
status and perspective. Internat. J. Numer. Methods Fluids 72 (8), 811–
845, arXiv:https://onlinelibrary.wiley.com/doi/pdf/10.1002/fld.3767. URL https:
//onlinelibrary.wiley.com/doi/abs/10.1002/fld.3767.

Xing, Y., Shu, C.W., 2005. High order finite difference WENO schemes with the exact conservation property for the shallow water equations. J. Comput. Phys. 208 (1), 206–227.

Xing, Yulong, Shu, Chi-Wang, 2006. High order well-balanced finite volume WENO schemes and discontinuous Galerkin methods for a class of hyperbolic systems with source terms. J. Comput. Phys. 214 (2), 567–598.

Xing, Y., Zhang, X., Shu, C.W., 2010. Positivity-preserving high order well-balanced discontinuous Galerkin methods for the shallow water equations. Adv. Water Resour. 33, 1476–1493.

Xu, H., Cantwell, C.D., Monteserin, C., Eskilsson, C., Engsig-Karup, A.P., Sherwin, S.J., 2018. Spectral/hp element methods: Recent developments, applications, and perspectives. J. Hydrodyn. 30, 1–22.