Cipriano Escalante
email: escalante@uma.es

Tomás Morales De Luna

Francisco N Cantero-Chinchilla

Oscar Castro-Orgaz

Vertically averaged and moment equations: modeling of non-hydrostatic free surface flows with an efficient path-conservative scheme

Keywords: dispersive waves, non-hydrostatic models, weighted average, depth-integrated modeling, finite volume, finite difference

come

Introduction

The computation of free-surface water flows in environmental hydraulics entails an accurate and efficient modeling of important physical phenomena, including flows over obstacles, tsunami propagation and run-up and dam break waves. In the shallow water framework, the horizontal length scales are generally larger than the vertical ones. The dispersionless nonlinear shallow water equations (SWE), which are derived by assuming vertical hydrostatic pressure distribution and depth-independent horizontal velocity, are often a reasonable choice in such situations.

However, in many applications, such as for example the propagation of undular bores in an estuary or in the dispersive waves propagating over an obstacle, the hydrostatic pressure hypothesis is no longer valid. In these flows the vertical acceleration is significant, and the pressure distribution is influenced by a significant dynamic component. As a consequence, for a correct simulation of these type of flows, it becomes necessary to explore other models that consider non-hydrostatic effects.

The full three-dimensional (3D) flow modeling is still costly from the computational point of view, specially for large-scale domains and long-time simulations. Therefore, alternative modeling approaches keeping a balance between computational cost and physical accuracy are considered by resorting to the vertically-averaging of the 3D equations. The vertically-averaged modeling of non-hydrostatic flows has been accomplished traditionally by resorting to two main family of models: (i) Boussinesq-type models, and (ii) depth-averaged non-hydrostatic models.

Boussinesq-type models account for the non-hydrostaticity of flow by retaining high-order terms in the Taylor expansion of the velocity potential in the momentum equations. For instance, this is the approach used in [START_REF] Abbott | Accuracy of short-wave numerical models[END_REF][START_REF] Boussinesq | Théorie des ondes et des remous qui se propagent le long d'un canal rectangulaire horizontal, en communiquant au liquide contenu dans ce canal des vitesses sensiblement pareilles de la surface au fond[END_REF][START_REF] Peregrine | Long waves on a beach[END_REF][START_REF] Green | A derivation of equations for wave propagation in water of variable depth[END_REF][START_REF] Madsen | A new form of the Boussinesq equations with improved linear dispersion characteristics. Part 2. A slowly-varying bathymetry[END_REF][START_REF] Nwogu | Alternative form of boussinesq equations for nearshore wave propagation[END_REF][START_REF] Wei | A fully nonlinear boussinesq model for surface waves. part 1. highly nonlinear unsteady waves[END_REF][START_REF] Langtangen | Computational models for weakly dispersive nonlinear water waves[END_REF][START_REF] Kazhyken | Discontinuous galerkin methods for a dispersive wave hydro-morphodynamic model with bed-load transport[END_REF][START_REF] Samii | An explicit hybridized discontinuous galerkin method for serre-green-naghdi wave model[END_REF] among others. This approach allows for an accurate characterization of the dispersive effects from shallow to intermediate water depths. However, the velocity dependent extra-terms retained in the governing equations contains high-order derivatives which results in complex systems that often requires a careful numerical discretization.

Contrary to Boussinesq type models, depth-averaged non-hydrostatic models use first-order derivatives, incorporating dispersive effects by retaining some vertical effects in the depth-averaged process for SWE. More explicitly, the vertical velocity or non-hydrostatic pressure is not assumed negligible, and the pressure is split into a hydrostatic and a non-hydrostatic part. The simplest choice of nonhydrostatic model considers a vertical 1-layer where the horizontal velocity is assumed constant (see for instance [START_REF] Yamazaki | Depth-integrated, non-hydrostatic model for wave breaking and run-up[END_REF][START_REF] Bristeau | An energy-consistent depth-averaged Euler system: Derivation and properties[END_REF][START_REF] Escalante | Non-hydrostatic pressure shallow flows: GPU implementation using finite volume and finite difference scheme[END_REF]). The dispersive properties of non-hydrostatic models can be improved by considering multilayer models, where the vertical domain is decomposed in a prescribed number of layers. Within these layers, the velocity and pressure variables are averaged in the vertical direction. In [START_REF] Stelling | An accurate and efficient finite-difference algorithm for non-hydrostatic free-surface flow with application to wave propagation[END_REF] 1-layer and 2-layer models are presented using a finite-difference algorithm. 2-layer versions of the model are presented in [START_REF] Bai | Depth-integrated free-surface flow with a two-layer non-hydrostatic formulation[END_REF] or [START_REF] Escalante | An Efficient Two-Layer Non-hydrostatic Approach for Dispersive Water Waves[END_REF]. Moreover, multilayer non-hydrostatic models are found in [START_REF] Fernández-Nieto | Layer-averaged approximation of euler equations for free surface flows with a non-hydrostatic pressure[END_REF][START_REF] Fernández-Nieto | A hierarchy of dispersive layeraveraged approximations of Euler equations for free surface flows[END_REF].

As it was shown in [START_REF] Escalante | A general non-hydrostatic hyperbolic formulation for boussinesq dispersive shallow flows and its numerical approximation[END_REF], there is a straight relation between Boussinesq and non-hydrostatic models. Most well-known Boussinesq models may be written in an equivalent non-hydrostatic formulation. Moreover, writing them as non-hydrostatic pressure systems avoids using high-order derivatives that are not easy to treat numerically.

There is a third possible method to modeling non-hydrostatic free surface flows, namely the method of weighted residuals. This technique was originally proposed by Green and Naghdi [START_REF] Green | A derivation of equations for wave propagation in water of variable depth[END_REF][START_REF] Green | Directed fluid sheets[END_REF][START_REF] Green | A direct theory of viscous fluid flow in channels[END_REF][START_REF] Green | A nonlinear theory of water waves for finite and infinite depths[END_REF][START_REF] Green | Further developments in a nonlinear theory of water waves for finite and infinite depths[END_REF], although it has not been properly acknowledged in the literature. They developed the theory of a Cosserat surface or the theory of "directed fluid sheet", where in the mass and momentum conservation laws the velocity vector is expanded in terms of some basis functions depending on the vertical coordinate. The expansions contain proportionality coefficients multiplying the basis functions, to be determined as part of the solution. Given that the 3D flow equations are unable to satisfy these expansions in an arbitrary point of the flow, Green and Naghdi proposed to seek solutions satisfied "at large", e.g., forming integral statements where the residuals of each conservation law are multiplied by a test function. Green and Naghdi theory is therefore a variational method where the dependence of a system of 3D equations on the vertical coordinate is removed. In the usual case the test functions are taken equal to the basis functions, resulting a Galerkin procedure. The lowest order version of their theory is the so-called nowadays Serre-Green-Naghdi theory [START_REF] Serre | Contribution à l'étude des écoulements permanents et variables dans les canaux[END_REF], not requiring to use variational statements using test functions different from unity [START_REF] Castro-Orgaz | Non-linear shallow water flow modelling over topography with depth-averaged potential equations[END_REF]. Interestingly, the method is in reality not original from Green and Naghdi, but it rather originates from the charming work of Kantorovich and Krylov [START_REF] Kantorovich | Approximate Methods of Higher Analysis[END_REF].

They presented in their book exactly the procedure described above with application to the solution of Poisson equation. The method of Kantorovich and Krylov [START_REF] Kantorovich | Approximate Methods of Higher Analysis[END_REF] or the Green-Naghdi theory [START_REF] Green | A derivation of equations for wave propagation in water of variable depth[END_REF][START_REF] Green | Directed fluid sheets[END_REF] is thus a Galerkin procedure suitable to produce vertically-averaged systems of conservation laws for non-hydrostatic flows.

The method was used in open-channel flow applications by Steffler and Jin [START_REF] Steffler | Depth averaged and moment equations for moderately shallow free surface flow[END_REF], who proposed the Vertically Averaged and Moment (VAM) equations model. However, it was not stated that the method is an application of Green-Naghdi theory. In fact, it can be demonstrated that the basis functions used by Steffler and Jin [START_REF] Steffler | Depth averaged and moment equations for moderately shallow free surface flow[END_REF] are the first two shifted Legendre polynomials, while for the test functions they took moments of the residuals with respect to the centroid of the flow section, thus the name "Moment" equations they used for the integrals of the weighted residuals. However, taking moments produce identical statements to using the first shifted Legendre polynomial as test function. Therefore, the VAM model is a particular application of Green-Naghdi theory obtained with a Galerkin procedure using Cartesian coordinates and an orthogonal base formed with shifted Legendre polynomials. The model was applied in [START_REF] Ghamry | Two dimensional vertically averaged and moment equations for rapidly varied flows[END_REF][START_REF] Ghamry | Effect of applying different distribution shapes for velocities and pressure on simulation of curved open channels[END_REF][START_REF] Ghamry | Two-dimensional depth-averaged modeling of flow in curved open channels[END_REF] to steady river flows and in [START_REF] Cantero-Chinchilla | Depth-integrated nonhydrostatic free-surface flow modeling using weighted-averaged equations[END_REF][START_REF] Cantero-Chinchilla | Nearshore coastal flow processes using weighted-averaged equations[END_REF] to unsteady flows in the river and maritime environment. To the date, the only phase resolving model available to solve this formidable set of non-linear PDEs is the finite-volume finite-difference scheme devised in [START_REF] Cantero-Chinchilla | Depth-integrated nonhydrostatic free-surface flow modeling using weighted-averaged equations[END_REF]. Moreover, in [START_REF] Gamero | A MATLAB software platform for modelling vertically-integrated non-hydrostatic flows with moment equations[END_REF], a MATLAB software platform was presented for modeling the VAM model.

However, to the best of the authors' knowledge, the computational efficiency and mathematical stability have not been targeted nor studied for the VAM model so far. These are crucial aspects for an accurate reproduction of the physical phenomena by computation. Furthermore, it would be interesting to study the accuracy and efficiency of the weighted residual models versus other models. This comparison is not only interesting from a physical standpoint, but also form a practical point of view. The main objective of this work is to produce an accurate, robust and efficient numerical approach for VAM model and compare its accuracy predicting dispersive effects against other non-hydrostatic models.

The paper is organized as follows: in Section 2 the VAM governing equations are stated while their derivation process is described. Although the model is not new, the derivation and the formulation presented in this work differs from the original one. The alternative VAM formulation permits us to later develop and propose in Section 3 an efficient numerical scheme based on a two-step approach: (i) the hyperbolic part is firstly solved using a second order path-conservative PVM scheme, and (ii) the hydrostatic terms are then corrected by solving a linear Poisson-like system using an iterative scheme. Later, in Section 4, the numerical method for the VAM model is applied to several test cases as compared with experimental data and the performance of the previous existing numerical scheme for the VAM model [START_REF] Cantero-Chinchilla | Depth-integrated nonhydrostatic free-surface flow modeling using weighted-averaged equations[END_REF], which relies on a hybrid FV-FD scheme along with a time evolution Newton-Raphson (NR) iterative scheme. The results are also compared with other models, including the SGN equations and two-layer non-hydrostatic model. Finally, conclusions of this study are drawn.

Description of the model

For the sake of simplicity, we consider here a two-dimensional flow over an uneven fixed bed, although everything might be generalized to the 3D case. Let (x, z) be the horizontal and vertical directions, respectively, and (u, w) the velocity vector. The Euler system for a fluid with constant density is given by

       ∂ x u + ∂ z w = 0, ∂ t u + u ∂ x u + w ∂ z u + ∂ x p T = 0, ∂ t w + u ∂ x w + w ∂ z w + ∂ z p T = -g. (1)
Remark that here we are neglecting the friction term for the sake of simplicity. We shall consider a decomposition of the pressure into hydrostatic and non-hydrostatic components:

p T = p H + p, with ∂ z p H = -g.
We denote by h(t, x) the fluid depth and we define b(x) the bottom topography. We shall denote by η(t, x) = h(t, x) + b(x) the free surface. As boundary conditions, we take the usual kinematic conditions at the bottom (z = b) and free surface (z = η):

∂ t b + u | z=b ∂ x b = w | z=b , ∂ t η + u |z=η ∂ x η = w |z=η .
Now, we write system (1) in the so-called σ-coordinates. To this end, we consider the change of variable given by

ξ = z -b h , where ξ ∈ [0, 1] for z ∈ [b, η].
Notice that ξ = 0 (respectively ξ = 1) corresponds to the bottom z = b (respectively free surface z = η) level (see Figure 1). Then, denoting by ψ(t, x, ξ) = ψ(t, x, ξh + b), the differential operators read

∂ ξ ψ = h∂ z ψ and h∂ s ψ = ∂ s h ψ -∂ ξ ∂ s (ξh + b) ψ , for s ∈ {t, x}.
Taking into account the new variables, the pressure is written as

p = p H + p, with ∂ ξ p H = -gh,
where it is clear that, assuming that the pressure is null at the surface, we obtain the expression for the hydrostatic pressure p H = -gh (ξ -1) .

Then system (1) in σ-coordinates reads

             ∂ x (h u) + ∂ ξ ( w -u∂ x (ξh + b)) = 0, ∂ t (h u) + ∂ x (h u 2 + h p) + gh∂ x (z b + h) +∂ ξ (ω u -p∂ x (hξ + z b )) = 0, ∂ t (h w) + ∂ x (h u w) + ∂ ξ (ω w + p) = 0,
where ω is given by ω

(t, x, ξ) = w(ξ) -(∂ t (hξ + b) + u(ξ)∂ x (hξ + b)).
Boundary conditions are then given by

ω | ξ =0 = ∂ t b + u | ξ=0 ∂ x b -w | ξ=0 = 0, ω | ξ =1 = ∂ t η + u | ξ=1 ∂ x η -w | ξ=1 = 0.
For the sake of simplicity we drop the symbol (•). We shall assume from now on a rigid bed so that ∂ t b = 0.

Consider now an orthogonal basis of P N [ξ] in the interval [0, 1], consisting of the Legendre polynomials ϕ j : [0, 1] → R, j = 0, 1, . . . , N.

In particular, in this paper, we set by choice N = 2 and the Legendre basis is reduced to

ϕ 0 (ξ) = 1, ϕ 1 (ξ) = 1 -2ξ, ϕ 2 (ξ) = 1 -6ξ + 6ξ 2 .
Notice that these polynomials satisfy ϕ i (0) = 1, and

1 0 ϕ i (ξ)ϕ j (ξ) dξ = µ i δ ij ,
where δ ij is the Kronecker symbol and µ i is the square norm of

ϕ i in P N [ξ].
The VAM model introduced in [START_REF] Cantero-Chinchilla | Depth-integrated nonhydrostatic free-surface flow modeling using weighted-averaged equations[END_REF] considers an approximation u ∈ P 1 [ξ], w, p ∈ P N [ξ], which will be written in the Legendre polynomial basis as

u = u 0 ϕ 0 + u 1 ϕ 1 , w = w 0 ϕ 0 + w 1 ϕ 1 + w 2 ϕ 2 , p = p 0 ϕ 0 + p 1 ϕ 1 + p 2 ϕ 2 .
Since the basis is orthogonal, it results 1 0 f (t, x, z)ϕ j (ξ)dξ = µ j f j (t, x), for f ∈ {u, w, p} and j = 0, 1, 2, where we are setting u 2 = 0.

Multiplying system (2) by ϕ j for j = 0, 1 and integrating in the vertical direction interval [0, 1] will produce

                                 ∂ t h + ∂ x (hu 0 ) = 0, ∂ t (hu 0 ) + ∂ x hu 2 0 + 1 3 hu 3 1 + hp 0 + gh∂ x η = -2p 1 ∂ x b, ∂ t (hw 0 ) + ∂ x hu 0 w 0 + 1 3 hu 1 w 1 = 2p 1 , h∂ x u 0 + 1 3 ∂ x (hu 1 ) + 1 3 u 1 ∂ x h + 2 (w 0 -u 0 ∂ x b) = 0, 1 3 ∂ t (hu 1 ) + 1 3 ∂ x (2hu 0 u 1 + hp 1 ) - 1 3 u 0 ∂ x (hu 1 ) -p 0 - 1 3 p 1 ∂ x h -2(p 0 -p 1 )∂ x b = 0, 1 3 ∂ t (hw 1 ) + 1 3 ∂ x (hu 0 w 1 ) + 1 3 ∂ x u 1 hw 0 + 2 5 hw 2 + 1 3 1 5 w 2 -w 0 ∂ x (hu 1 ) + 2p 0 -p b = 0,
where p b is the non-hydrostatic pressure at the bottom, given by p b = p 0 + p 1 + p 2 . The previous system, consisting in 6 equations for 9 variables is completed with the assumption that the non-hydrostatic pressure vanishes at the surface p 0 -p 1 + p 2 = 0 and the corresponding kinematic conditions

∂ t η + u | ξ=1 ∂ x η -w | ξ=1 = 0, u | ξ=0 ∂ x b -w | ξ=0 = 0.
In view of choice of the Legendre basis, we have

u | ξ=0 = u 0 + u 1 , u | ξ=1 = u 0 -u 1 , w | ξ=0 = w 0 + w 1 + w 2 , w | ξ=1 = w 0 -w 1 + w 2 ,
that used in the kinematic conditions equations gives

∂ t η + (u 0 -u 1 )∂ x η = w 0 -w 1 + w 2 , (u 0 + u 1 )∂ x b = w 0 + w 1 + w 2 .
Subtracting and adding both equations, we write them in the following equivalent way:

h∂ x u 0 + u 1 ∂ x h + 2(u 1 ∂ x b -w 1 ) = 0, h∂ x u 0 + u 1 ∂ x h + 2(w 0 + w 2 -u 0 ∂ x b) = 0.

Compact form and hyperbolicity of the underlying hydrostatic system

The VAM model, as introduced here, can be written as the system

                       ∂ t h + ∂ x (hu 0 ) = 0, ∂ t (hu 0 ) + ∂ x hu 2 0 + 1 3 hu 2 1 + hp 0 + gh∂ x η = -2p 1 ∂ x b, ∂ t (hw 0 ) + ∂ x hu 0 w 0 + 1 3 hu 1 w 1 = 2p 1 , ∂ t (hu 1 ) + ∂ x (2hu 0 u 1 + hp 1 ) -u 0 ∂ x (hu 1 ) -(3p 0 -p 1 ) ∂ x h -6(p 0 -p 1 )∂ x b = 0, ∂ t (hw 1 ) + ∂ x hu 0 w 1 + u 1 hw 0 + 2 5 hw 2 + 1 5 w 2 -w 0 ∂ x (hu 1 ) + 6(p 0 -p 1 ) = 0, (2)
together with the constraints

     h∂ x u 0 + 1 3 ∂ x (hu 1 ) + 1 3 u 1 ∂ x h + 2 (w 0 -u 0 ∂ x b) = 0, h∂ x u 0 + u 1 ∂ x h + 2(u 1 ∂ x b -w 1 ) = 0, h∂ x u 0 + u 1 ∂ x h + 2(w 0 + w 2 -u 0 ∂ x b) = 0 (3)
We write ( 2)-( 3) in a more compact form:

             ∂ t U + ∂ x F (U, w 2 ) + G(U, w 2 )∂ x U + S(U )∂ x b + T (U, ∂ x U, P, ∂ x P, ∂ x b) = 0, I 1 (U, ∂ x U, ∂ x b) = 0, I 2 (U, ∂ x U, ∂ x b) = 0, w 2 = -(w 0 + w 1 ) + (u 0 + u 1 )∂ x b, (4)
with

U =       h hu 0 hu 1 hw 0 hw 1       , P = hp 0 hp 1 , F (U, w 2 ) =              hu 0 hu 2 0 + 1 3 hu 3 1 hu 0 w 0 + 1 3 hu 1 w 1 2hu 0 u 1 hu 0 w 1 + u 1 hw 0 + 2 5 hw 2 0              , G(U, w 2 ) = ghE 2,1 -u 0 E 4,3 - 1 5 w 2 -w 0 E 5,3 S(U ) = ghe 2 ,
where E i,j represents the usual elements of the canonical basis for matrices of size 5 × 5, and e 2 is the second element of the canonical basis in R 5 , that is,

G(U, w 2 )∂ x U =           0 gh∂ x h 0 -u 0 ∂ x (hu 1 ) - 1 5 w 2 -w 0 ∂ x (hu 1 ) 0           , S(U )∂ x b =         0 gh∂ x b 0 0 0 0         .
The non-hydrostatic terms correspond to

T (U, ∂ x U, P, ∂ x P, ∂ x b) =         0 ∂ x (hp 0 ) + 2p 1 ∂ x b -2p 1 ∂ x (hp 1 ) -(3p 0 -p 1 ) ∂ x h -6(p 0 -p 1 )∂ x b 6(p 0 -p 1 ) 0        
, and the constraints

I 1 (U, ∂ x U ) = h∂ x u 0 + 1 3 ∂ x (hu 1 ) + 1 3 u 1 ∂ x h + 2 (w 0 -u 0 ∂ x b) , I 2 (U, ∂ x U ) = h∂ x u 0 + u 1 ∂ x h + 2(u 1 ∂ x b -w 1 ).
Now, let us consider the underlying hydrostatic pressure system given by

∂ t U + ∂ x F (U, w 2 ) + G(U )∂ x U + S(U )∂ x b = 0. ( 5
)
It can be easily check that ( 5) can be written in quasi-linear form as

∂ t U b + A(U, w 2 ) S(U ) 0 0 + ∂ x U b = 0 with A(U, w 2 ) = J F + G(U, w 2 ),
where J F is the Jacobian of the flux F with respect to the conserved variables U. The eigenvalues of the matrix A(U, w 2 ) are

λ 1 = u 0 , λ 2,3 = u 0 ± 1 √ 3 u 1 , λ 4,5 = u 0 ± gh + u 2 1 , λ 6 = 0,
and therefore, the system is hyperbolic and the full eigenstructure can be explicitly computed.

Link with existing dispersive systems in the literature

As presented here, the VAM model differs from the one presented in [START_REF] Cantero-Chinchilla | Depth-integrated nonhydrostatic free-surface flow modeling using weighted-averaged equations[END_REF]. Indeed, in that paper the original variables in [START_REF] Abbott | Accuracy of short-wave numerical models[END_REF] where assumed to have the following polynomial approximations

u(t, x, z) = u 0 (t, x) -u 1 (t, x)(1 -2ξ), w(t, x, z) = w b (t, x)(1 -ξ) + 4 w 2 (t, x)ξ(1 -ξ) + w η (t, x)ξ, p(t, x, z) = p 1 (t, x)(1 -ξ) + 4 p 2 (t, x)ξ(1 -ξ),
where w b and w η are the vertical velocity at the bottom and at the surface respectively and ξ = (z-b)/h. It is then clear that the following relations hold:

w b = w 0 + w 1 + w 2 , w η = w 0 -w 1 + w 2 , u 0 = u 0 , u 1 = -u 1 , p 1 = p | z=b = p 0 + p 1 + p 2 ,
and integrating w and p in the vertical direction

1 3 w 2 = w b + w η 2 -w 0 , 1 3 p 2 = p 1 -p 0 = p 1 + p 2 .
Then, the first three equations in [START_REF] Cantero-Chinchilla | Depth-integrated nonhydrostatic free-surface flow modeling using weighted-averaged equations[END_REF] are obtained by integrating (1) in the vertical direction, that is, considering

f (t, x) = 1 h η b f (t, x, z)dz = 1 0 f (t, x, ξ)dξ = f 0 (t, x),
for any variable f , that is, the first three equations coincide. Then, the other three equations are obtained by considering the weighted residuals, this is, the equations are multiplied by

1 h η b 2 z -b h - 1 2 f (t, x, z) dz = 1 0 (2ξ -1) f (t, x, ξ)dξ = 1 0 ϕ 1 (ξ) f (t, x, ξ)dξ,
so that the other three equations also coincide, although they are written by using the unknowns f for each variable f .

Dispersive properties of linear waves

In this Subsection, the dispersive features of the model are presented for a complete description of the VAM model. In particular, we focus on the linear dispersion relation for the wave and group velocities and the linear shoaling gradient.

As usual, the governing equations are linearised around the steady state solution corresponding to a constant water thickness H and zero velocity given by

b = -H, u 0 = u 1 = w 0 = w 1 = 0,
where we are considering the perturbations of free surface η with respect to the zero level. We then consider the following asymptotic expansion h = H + ϵη (1) + O(ϵ 2 ), f = ϵf (1) + O(ϵ 2 ), ∀f ∈ {u 0 , u 1 , w 0 , w 1 }.

In the following, we neglect the notation (1) for the sake of simplicity. Using this linearisation in (2)-(3), we shall neglect O 2 terms and keep the system at first order:

           ∂ t η + H∂ x u 0 = 0, ∂ t u 0 + ∂ x p 0 + g∂ x η = 0, H∂ t w 0 = 2p 1 , ∂ t u 1 + ∂ x p 1 = 0, H∂ t w 1 + 6(p 0 -p 1 ) = 0, (6)
together with the restrictions

     H∂ x u 0 + H 3 ∂ x u 1 + 2w 0 = 0, H∂ x u 0 -2w 1 = 0, H∂ x u 0 + 2(w 0 + w 2 ) = 0. ( 7
)
We then conduct a standard Stokes-type Fourier analysis were the linearised variables are supposed to be of the form

f (x, t) = f • e i(ωt-kx) , ∀ f ∈ {η, u 0 , u 1 , w 0 , w 1 },
where ω is the wave frequency and k the wave number. By substituting the previous expression in ( 6)-( 7), we get a linear system for the unknowns f . Non-trivial solutions for the aforementioned linear system are given when the following linear dispersion relation for the wave celerity, C = ω/k, holds:

C 2 √ gH = (kH) 2 12 + 1 (kH) 4 144 + 5 (kH) 2 12 + 1
.

Then, following a general procedure similar to the one described in [START_REF] Escalante | Nonhydrostatic layer-averaged Euler system with layerwise linear horizontal velocity[END_REF] .

An interesting index that measures the shoaling of impinging waves on the presence of a slope, is the shoaling gradient coefficient, proposed by Madsen et al (see [START_REF] Madsen | A new form of the Boussinesq equations with improved linear dispersion characteristics. Part 2. A slowly-varying bathymetry[END_REF]). The shoaling coefficient γ is defined as the quantity that fulfills the relation

∂ x η η = -γ ∂ x H H .
As for the group velocity, we follow the general procedure for the computation of γ from the linear dispersion relation for C (see [START_REF] Escalante | Nonhydrostatic layer-averaged Euler system with layerwise linear horizontal velocity[END_REF]), and it yields .

In Figure 2 we show comparisons of the described linear dispersion relations with respect to the Airy linear theory, which establishes

C 2 Airy gH 0 = tanh(kH 0 ) kH 0 , C 2
g,Airy

gH 0 = (2kH 0 + sinh(2kH 0 )) 2 2kH 0 (2 sinh(2kH 0 ) + sinh(4kH 0 )) , and
γ Airy = kH 0 tanh(kH 0 ) (1 -kH 0 tanh(kH 0 )) 1 -tanh 2 (kH 0 ) tanh(kH 0 ) + kH 0 1 -tanh 2 (kH 0 ) 2 .
As it is seen in Figure 2, the proposed VAM model ensures a minimal error of less than 5 % in an extended range of kH ∈ [0, 4] for the wave celerity and group velocity. system that fits better. In order to take into account a whole range in kH, following the ideas in [START_REF] Escalante | An Efficient Two-Layer Non-hydrostatic Approach for Dispersive Water Waves[END_REF] let us consider the weighed L 1 -norm relative errors for a fixed interval [0, (kH) max ], that is, we define

ϵ C ((kH) max ) = (kH)max 0 1 kH C 2 (kH) -C 2 Airy (kH) C 2 Airy (kH) dkH, ϵ Cg ((kH) max ) = (kH)max 0 1 kH C 2 g (kH) -C 2 g,Airy (kH) C 2 g,Airy (kH) dkH, ϵ γ ((kH) max ) = (kH)max 0 1 kH γ(kH) -γ Airy (kH) γ Airy (kH) dkH. (8)
These L 1 -norm errors are shown in Figure 3 for all the above mentioned systems. As we see, when considering the whole range kH ∈ [0, 6] for the wave celerity, the VAM system performs better than the other systems. It is for large ranges of wavenumber that the Two-Layer system will give a better result. Nevertheless, as shown in Figure 2, for large values of kH the errors rapidly increase. Similar results are obtained for the group velocity, where now VAM is better up to (kH) max = 5, or for the shoaling gradient, where VAM is again better up to (kH) max = 3.5.

Numerical scheme

Let us describe now the numerical scheme used to discretize the system (4). The numerical method is based on a two-step projection-correction method. First, we shall solve the underlying non-conservative hydrostatic and hyperbolic system given by [START_REF] Bristeau | An energy-consistent depth-averaged Euler system: Derivation and properties[END_REF].

Then, in a second step, non-hydrostatic terms will be taken into account. System (5) is numerically approximated by means of a second order finite volume Polynomial Viscosity Matrix method (from now on PVM); it is positive-preserving, well-balanced, path-conservative method based on segments [START_REF] Castro Díaz | A class of computationally fast first order finite volume solvers: PVM methods[END_REF]. As usual, we consider a set of N finite volume cells the cell average of the function U (x, t) on cell I i at time t. Regarding non-hydrostatic terms, we consider the mid-point x i of each cell I i and denote the point values of the function P at time t by

I i = [x i-1/2 , x i+1/2 ] with constant lengths ∆x and define U i (t) = 1 ∆x I i U (x, t) dx,
P i (t) = P (x i , t).
Non-hydrostatic terms will be approximated by second order compact finite-differences.

Time stepping

Assume given time steps ∆t n , and denote t n = j≤n ∆t j . To obtain second order accuracy in time, a two-stage second-order TVD Runge-Kutta scheme [START_REF] Gottlieb | Total variation diminishing Runge-Kutta schemes[END_REF] is adopted. At the kth stage, k ∈ {1, 2}, the two-step projection-correction method is given by

                                   U ( k) -U (k-1) ∆t + ∂ x F (U (k-1) , w (k-1) 2 ) + G(U (k-1) , w (k-1) 2 )∂ x U (k-1) + S(U (k-1) )∂ x b = 0, U (k) -U ( k) ∆t + T (U (k) , ∂ x U (k) , P (k) , ∂ x P (k) , ∂ x b) = 0, I 1 (U (k) , ∂ x U (k) , ∂ x b) = 0, I 2 (U (k) , ∂ x U (k) , ∂ x b) = 0, w (k) 2 = -(w (k) 0 + w (k) 1 ) + (u (k) 0 + u (k) 1 )∂ x b, (9a) (9b) (9c) (9d)
where U (0) is U at the time level t n , U ( k) is an intermediate value in the two-step projection-correction method that contains the numerical solution of the hyperbolic (system (9)) at the corresponding kth stage of the Runge-Kutta, and U (k) is the kth stage estimate. Finally, the value of the solution at the t n+1 time level is obtained as: 2) .

U n+1 = 1 2 U n + 1 2 U (
Note that, equations (9a-9c) requires, at each stage of the calculation respectively, to solve a Poissonlike equation for each one of the variables contained in P (k) . This will be described below. Finally, once the non-hydrostatic pressure is known, the values of U (k) and w (k) 2 are updated by means of equations (9a) and (9d) respectively.

The usual Courant-Friedrichs-Lewy (CFL) restriction should be considered for the computation of the time step ∆t.

Finite volume discretization for the underlying hyperbolic system

A second order path-conservative PVM scheme for the discretization in space of the system (9), at the kth stage of the Runge-Kutta, reads as follows (see [START_REF] Castro Díaz | A class of computationally fast first order finite volume solvers: PVM methods[END_REF][START_REF] Castro | High order finite volume schemes based on reconstruction of states for solving hyperbolic systems with nonconservative products. applications to shallow-water systems[END_REF] for further details):

U ( k) -U (k-1) ∆t = - 1 ∆x D (k-1),- i+1/2 + D (k-1),+ i-1/2 + VI i , (10)
where, neglecting the time dependence k for for the sake of simplicity in the notation,

D - i+1/2 = D - i+1/2 (U - i+1/2 , U - i+1/2 , w - 2,i+1/2 , w - 2,i+1/2 , b + i+1/2 , b - i+1/2 ) = 1 2 F (U + i+1/2 , w + 2,i+1/2 ) -F (U - i+1/2 , w - 2,i+1/2 ) +G(U i+1/2 , w 2,i+1/2 ) U + i+1/2 -U - i+1/2 + S(U i+1/2 )(b + i+1/2 -b - i+1/2 ) ± 1 2 Q i+1/2 U + i+1/2 -U - i+1/2 + A -1 i+1/2 S(U i+1/2 )(b + i+1/2 -b i+1/2 ) -, VI i = F (U - i+1/2 , w - 2,i+1/2 ) -F (U + i-1/2 , w + 2,i-1/2 ) + G(U i , w 2,i ) U - i+1/2 -U + i-1/2 + S(U i ) b - i+1/2 -b + i-1/2
where

U i+1/2 = (h i+1/2 , (hu 0 ) i+1/2 , (hu 1 ) i+1/2 , (hw 0 ) i+1/2 , (hw 1 ) i+1/2 , b i+1/2 ) T , with h i+1/2 = h + i+1/2 + h - i+1/2 2 , b i+1/2 = b + i+1/2 + b - i+1/2 2 , (hυ) i+1/2 = h i+1/2 υ i+1/2 , υ i+1/2 = υ + i+1/2 h + i+1/2 + υ - i+1/2 h - i+1/2 h + i+1/2 + h - i+1/2
, for υ = u 0 , u 1 , w 0 , w 1 , w 2 .

The values U ± i+1/2 and w ± 2,i+1/2 are defined by a reconstruction procedure on the variables to the left ( -) and right ( + ) of the inter-cell x i+1/2 . This reconstruction procedure is done using a MUSCL reconstruction operator (see [START_REF] Leer | Towards the ultimate conservative difference scheme. V. A second-order sequel to Godunov's method[END_REF]), combined with a minmod limiter. Q i+1/2 corresponds to the viscosity matrix associated to the numerical method. In general, Q i+1/2 is obtained by a polynomial evaluation of the Roe Matrix. Here, Q i+1/2 is defined as

Q i+1/2 = β 0 Id + β 1 A i+1/2 , being β 0 = S R |S L | -S L |S R | S R -S L , β 1 = |S R | -|S L | S R -S L ,
where S L , S R are estimates of the smallest and largest wave speeds, respectively, at the interface x i+1/2 , as usual. Here, we use

S L = min u 0,i+1/2 -gh i+1/2 + u 2 1,i+1/2 , u 0,i -gh i + u 2 1,i+1/2 , S R = max u 0,i+1/2 + g ĥi+1/2 + u 2 1,i+1/2 , u 0,i+1 + gh i+1 + u 2 1,i+1/2 .
A i+1/2 denotes the matrix

A i+1/2 = J F (U i+1/2 , w 2,i+1/2 ) + G(U i+1/2 , w 2,i+1/2 ),
being J F the Jacobian matrix of the flux F.

The scheme described previously corresponds to a segment path-conservative extension of the HLL scheme [START_REF] Harten | On upstream differencing and godunov-type schemes for hyperbolic conservation laws[END_REF] for non-conservative systems.

Finite difference discretization for the non-hydrostatic terms

Now we describe the discretization of the non-hydrostatic terms. For this task, we shall substitute U (k) from (9a) in equations (9b)-(9c), resulting the Poisson-like equations

a 1 ∂ xx p (k) 0 + a 2 ∂ x p (k) 0 + a 3 p (k) 0 + a 4 ∂ xx p (k) 1 + a 5 ∂ x p (k) 1 + a 6 p (k) 1 = RHS 1 b 1 ∂ xx p (k) 0 + b 2 ∂ x p (k) 0 + b 3 p (k) 0 + b 4 ∂ xx p (k) 1 + b 5 ∂ x p (k) 1 + b 6 p (k) 1 = RHS 2 (11)
where the coefficients a j , b j as well as the Right-Hand-Sides RHS j depends on U ( k) and w ( k) 2 . Equations [START_REF] Castro-Orgaz | Non-linear shallow water flow modelling over topography with depth-averaged potential equations[END_REF] are discretized using second order finite differences. In order to obtain point value approximations of the non-hydrostatic pressure variables, terms p

(k) 0,i , p (k)
1,i and RHS j,i , j ∈ {1, 2}, i ∈ {1, . . . , N } will be approximated at every mid point x i of each cell I i . Further, the corresponding space derivatives will be approximated using compact second order central finite differences

∂ x P (k) i = P (k) i+1 -P (k) i-1 2∆x , ∂ xx P (k) i = P (k) i+1 -2P (k) i + P (k) i-1 ∆x 2 (12)
To compute the coefficients a j , b j , RHS j , which contain terms depending on U ( k) as well as its first derivatives, we will use the averaged values on the cell I i as second order point value approximations at the center of the cells. Similarly, the space derivatives are computed from the averaged values using compact second order finite differences

∂ x U ( k) i = U ( k) i+1 -U ( k) i-1 2∆x . ( 13
)
After replacing ( 12) and ( 13) in [START_REF] Castro-Orgaz | Non-linear shallow water flow modelling over topography with depth-averaged potential equations[END_REF], ones obtain a linear system, that is solved following an iterative method similar to the one described in [START_REF] Escalante | An Efficient Two-Layer Non-hydrostatic Approach for Dispersive Water Waves[END_REF].

Remark 1. In practice one could use any independent linear combination of equations ( 9b), (9c) and (9d). Nevertheless, although this is equivalent at continuous level, it plays a relevant role at discrete level when solving iteratively the final linear system. In practice, we have empirically found that better results are obtained when we use

I 1 (U (k) , ∂ x U (k) , ∂ x b) + I 2 (U (k) , ∂ x U (k) , ∂ x b) = 0, I 1 (U (k) , ∂ x U (k) , ∂ x b) -I (U (k) , ∂ x U (k) , ∂ x b) = 0,
and then we proceed in a similar way as stated before. Doing so we have found that the iterative algorithm for the linear system converges 1.6 times faster.

Final numerical scheme

For the sake of clarity, a guideline of the final numerical scheme will be highlighted as follows. For every kth stage k ∈ {1, 2} of the Runge-Kutta method, the problem ( 4) is split into two parts. A two-step projection-correction method is used:

• Finite volume step (solving the hydrostatic system): From [START_REF] Castro Díaz | A class of computationally fast first order finite volume solvers: PVM methods[END_REF], solve explicitly the hyperbolic system at the kth stage of the Runge-Kutta:

U ( k) -U (k-1) ∆t = - 1 ∆x D (k-1),- i+1/2 + D (k-1),+ i-1/2 + VI i ,
by means of a PVM path-conservative finite volume scheme combining a MUSCL reconstruction operator to obtain the intermediate value U ( k) in the two-step projection-correction method that contains the numerical solution of the underlying hyperbolic system.

• Finite difference step (non-hydrostatic pressure correction):

-Solve the Poisson-like equations ( 11) to obtain the non-hydrostatic pressures P (k) . To do so, compact centred second order finite differences are used for the discretization of the derivatives that appear in [START_REF] Castro-Orgaz | Non-linear shallow water flow modelling over topography with depth-averaged potential equations[END_REF], and a linear system is solved to obtain P (k) . -With the computed non-hydrostatic pressure terms P (k) , then hu 0 , hu 1 , hw 0 , hw 1 can be updated from (9a) and w 2 from (9d). To do that, a second order point value approximation in the centre of the cell will be used to compute the non-hydrostatic contribution.

Finally, 2) .

U n+1 = 1 2 U n + 1 2 U (
Concerning the boundary conditions, and the wetting and drying treatment, we follow procedure applied in [START_REF] Escalante | An Efficient Two-Layer Non-hydrostatic Approach for Dispersive Water Waves[END_REF].

Numerical tests

In this section, we highlight the ability of the proposed VAM numerical model to simulate a wide variety of complex situations involving dispersive water waves. Comparisons with some standard numerical tests for dispersive water wave systems are shown, including analytical solutions and laboratory data. We shall as well compare the results with the classic and well-known SGN as introduced in [START_REF] Green | A derivation of equations for wave propagation in water of variable depth[END_REF][START_REF] Serre | Contribution à l'étude des écoulements permanents et variables dans les canaux[END_REF][START_REF] Lannes | Derivation of asymptotic two-dimensional time-dependent equations for surface water wave propagation[END_REF][START_REF] Fernández-Nieto | A hierarchy of dispersive layeraveraged approximations of Euler equations for free surface flows[END_REF].

In what follows, the gravity acceleration is set to g = 9.81, and, unless stated otherwise, all the quantities are expressed in units of measure of the International System of Units.

Steady flow over fixed bedforms

We consider a steady non-hydrostatic flow over a fixed bedform in this first test case. As the flow passes over the bottom bump, the pressure drops below the hydrostatic value due to vertical acceleration. We intend to accurately reproduce the net drag force or form drag determined by the non-hydrostatic bottom pressure distribution.

The numerical test reproduced here was experimentally investigated in [START_REF] Sivakumaran | Steady shallow flow over curved beds[END_REF], where symmetric and asymmetric bed-forms were considered, and experimental measurements for the free-surface level and the bed pressure were obtained. A flume of length 9.15 m, width 0.3 m and 0.65 m deep was used. Constant inflow discharge was imposed. Two types of experiments are considered:

• A Gaussian profile with length 1.2 m and maximum elevation equal to 0.2 m and an inflow discharge q = 0.11197 m 2 /s.

• An asymmetrical hump 0.3 m high and 1.5 m long and inflow q = 0.11165 m 2 /s

For the first test case, we consider a computational domain given by [-1.5, 1.5], whereas [-0.5, 1.7] is used for the second test case. For all the test cases, the computational domain is covered with constant cells of length ∆x = 0.01 m, the CFL number is set to 0.4, and the simulation time was t = 50 s.

In the case of the first experiment, to simulate the inflow, water is initially at rest and a closed gate is located at x = 1. The water depth is h = 0.015 m, downstream of the gate and h = 0.34 m upstream. The rest of the flow variables are set initially to zero. A constant discharge boundary condition as used in the experiment was prescribed at the upstream end of the domain, and free-outflow boundary condition is imposed at the downstream part. The gate was removed instantaneously from the flume, and a dam-break-like flow was generated. In the case of the second experiment, the inflow is simulated by setting a water depth level to h = 0.02 m downstream of the gate and h = 0.42 m upstream. Friction term is neglected, given that the flow resistance is not relevant in these tests.

The results are shown in Figure 4. The numerical approach presented here for the VAM model correctly reproduces both the free surface and the bed pressures for the experimental data.

Solitary wave solution

A usual test for non-hydrostatic free-surface flows corresponds to the simulation of solitary waves. We consider here the general expression given in [START_REF] Bristeau | An energy-consistent depth-averaged Euler system: Derivation and properties[END_REF]:

h(t, x) = H * + A sech 2 1 H * γA 2(A + H * ) (x -ct) , u(t, x) = c 1 - H * h(t, x) , w(t, x) = - cA h(t, x) γA 2(A + H * ) sech 3 1 H * γA 2(A + H * ) (x -ct) sinh 1 H * γA 2(A + H * ) (x -ct) (14)
where A and H * are constant fixed values and c = g(A + H * ). We denote by η = h -H * the corresponding free surface.

The particular case γ = 3/2 corresponds to the SGN system [START_REF] Serre | Contribution à l'étude des écoulements permanents et variables dans les canaux[END_REF][START_REF] Green | A derivation of equations for wave propagation in water of variable depth[END_REF], and that value will be used throughout this paper. In that case, we shall consider a large-amplitude solitary wave using H * = 1 m and A = 0.55 m.

The considered computational domain is [-50, 700], and is covered with a set of constant cells of ∆x = 0.1 m length. The CFL is equal to 0.9, and the simulation time was t = 150 s. Free-outflow boundary conditions are imposed on both sides of the domain. A constant bottom is set to b = H * = 1 m; the initial water depth variable is initially given by h as in [START_REF] Escalante | A general non-hydrostatic hyperbolic formulation for boussinesq dispersive shallow flows and its numerical approximation[END_REF]; the mean of the horizontal and vertical velocities for the VAM model, u 0 and w 0 , are given by u, and w as in [START_REF] Escalante | A general non-hydrostatic hyperbolic formulation for boussinesq dispersive shallow flows and its numerical approximation[END_REF], respectively; the rest of flow variables for the VAM model are initially set to zero.

The solitary wave selected is not an exact solitary wave solution for the VAM model. Therefore we expected it to be perturbed and to change its shape. However, a solitary wave solution for the VAM model is reached after some time. That can be seen in Figure 5, where amplitude increases until a stable solitary wave solution of amplitude A = 0.65.

Similarly, we compare in Figure 6 the profiles given by the computed solitary wave obtained with the VAM model; the analytical solution [START_REF] Escalante | A general non-hydrostatic hyperbolic formulation for boussinesq dispersive shallow flows and its numerical approximation[END_REF]; and the solitary wave profile obtained by Tanaka solving the exact boundary value problem [START_REF] Tanaka | The stability of solitary waves[END_REF], taken from [START_REF] Carter | The kinematics and stability of solitary and cnoidal wave solutions of the serre equations[END_REF]. It can be observed a good fitting of the solitary wave given by the VAM model with the Tanaka solution.

Solitary wave run-up on a plane beach

A classic test for non-hydrostatic shallow flows is the one corresponding to the experimental setup carried out by Titov and Synolakis and described in [START_REF] Titov | Modeling of breaking and nonbreaking long-wave evolution and runup using VTCS-2[END_REF]. Incident solitary waves of multiple relative amplitudes were simulated to study propagation, breaking and run-up over a planar beach with a slope of 1 : 19.85. Experimental data are available in [START_REF] Titov | Modeling of breaking and nonbreaking long-wave evolution and runup using VTCS-2[END_REF] for surface elevation at different times.

In the test case used here, the still water level is H * = 1 m. The bathymetry of the problem is given in Figure 7. The initial condition for the VAM model is provided by a solitary wave of amplitude A = 0.3 m centred at point x = 25 m, and that is done using [START_REF] Escalante | A general non-hydrostatic hyperbolic formulation for boussinesq dispersive shallow flows and its numerical approximation[END_REF] as in the previous Subsection.

In this case, a Darcy coefficient of n = 0.01 s • m -1/3 was used to account for the glass surface roughness effects appearing in the experiments. The computational domain [-15, 50] is covered with cells of constant length ∆x = 0.1 m. The CFL number is set to 0.4, and free-outflow boundary conditions were imposed everywhere. The simulation time is t = 15 s. Figure 8 shows snapshots at different times, t g/H * = t 0 . A good agreement between experimental and simulated data is seen. We remark that wave breaking is observed at t g/H = 20 and 25 during the experiment. That wave breaking effect is well reproduced by the VAM model without any special treatment. Moreover, the run-up is in good agreement with the data and the wet-dry transitions are correctly treated by the numerical scheme. Now that we concluded that the VAM system correctly reproduces the experimental data, let us now compare the results given by VAM system with the ones obtained using the well-known SGN system. The comparative results are shown in Figure 9. We remark that initially both systems give similar results. Nevertheless, as time evolves and the wave front travels up the plane beach, they began to differ. First, the shock front position is misplaced and SGN gives a slower speed for this shock. Second, we see in the SGN simulations the characteristic overshoot in the surface. It is well-known that in such situations a breaking-wave mechanism is needed (see [START_REF] Escalante | Non-hydrostatic pressure shallow flows: GPU implementation using finite volume and finite difference scheme[END_REF] and references therein). This is not the case for VAM system, where the breaking-wave is correctly simulated without any further treatment. We remark as well that the run-up for SGN does not coincide with VAM.

If we compare the wall-clock times needed for each of them, they are quite similar: 3.261s for SGN vs 3.464s. Therefore, we see that the gain on the quality and precision on the solution for VAM are remarkable, with very little increase in computational effort.

Propagation of sinusoidal waves over submerged bar

The transformation of periodic waves approaching a submerged bar into deep water waves as described in [START_REF] Beji | Numerical simulation of nonlinear wave propagation over a bar[END_REF][START_REF] Dingemans | Comparison of computations with boussinesq-like models and laboratory measurements[END_REF] is also a popular test case for phase-resolving models.

For this test, the spatial domain [0, 30] with a submerged trapezoidal obstacle is considered as in Figure 10. The domain is covered with cells of the constant length of ∆x = 0.02 m. An incident sinusoidal wave train for the free-surface η is imposed as a boundary condition at the left-hand side of the domain x = 0 m. That is done as in [START_REF] Escalante | An Efficient Two-Layer Non-hydrostatic Approach for Dispersive Water Waves[END_REF] for the free-surface η (Figure 10). The results at the different gauge points are given in Figure 11. First, we see an excellent agreement for VAM model with the experimental data. We remark that we recover similar observations stated in the literature, such as [START_REF] Ma | Shock-capturing non-hydrostatic model for fully dispersive surface wave processes[END_REF] where σ-coordinates are used, or [START_REF] Escalante | An Efficient Two-Layer Non-hydrostatic Approach for Dispersive Water Waves[END_REF] where an enhanced two-layer version of the non-hydrostatic pressure multilayer system LDN H 0 is used. The results in [START_REF] Chazel | Numerical simulation of strongly nonlinear and dispersive waves using a Green-Naghdi model[END_REF] with a three-parameter Green-Naghdi model optimized for uneven bottoms show the same level of agreement. Therefore, we can obtain excellent agreement with a model that does not need the calibration of parameters (such as [START_REF] Escalante | An Efficient Two-Layer Non-hydrostatic Approach for Dispersive Water Waves[END_REF] or [START_REF] Chazel | Numerical simulation of strongly nonlinear and dispersive waves using a Green-Naghdi model[END_REF]) nor the use of complex vertical discretizations [START_REF] Ma | Shock-capturing non-hydrostatic model for fully dispersive surface wave processes[END_REF].

In the same figure we see as well the results obtained with SGN. The comparison with experimental data for gauge points located after the submerged bar is not good. When comparing the wall-clock times for each system, VAM requires 62.694s while SGN needs 57.235s for the whole simulation. This means that the results are largely improved with similar computational effort.

Moreover, in order to highlight the performance of VAM model, we show as well the results obtained with the model LDN H 0 introduced in [START_REF] Escalante | An Efficient Two-Layer Non-hydrostatic Approach for Dispersive Water Waves[END_REF][START_REF] Escalante | Numerical simulations of a dispersive model approximating free-surface euler equations[END_REF] for the case of two layers. The LDN H 0 is a multilayer non-hydrostatic model which tends to the dispersive relations provided by Airy linear theory as the number of layers increases. Of course, the counterpart of increasing the number of layers is that the computational cost increases as well. In this case we see that two layers are sufficient to obtain good results at all gauges. In particular we see that VAM model fit experimental data as well as a two-layer models. Therefore, we obtain comparable results without the extra computational cost of increasing the number of layers.

Dam-break flood waves

With the subsequent two proposed numerical tests, we check the shock-capturing capabilities of the numerical model by considering a dam-break flood wave. Both situations will produce a rarefaction wave with an edge moving faster than the long-wave celerity and an undular shock front composed of a solitary-like leading wave and the second train of cnoidal-like called Favre waves. The experiments are described in [START_REF] Favre | Étude théorique et expérimentale des ondes de translation dans les canaux découverts[END_REF][START_REF] Soares-Frazão | A second-order semi-implicit hybrid scheme for one-dimensional boussinesq-type waves in rectangular channels[END_REF], where a flume 15.24 m long, 0.4 m wide, and 0.4 m high is used.

For the first test case, the domain [-10, 10] is discretized with ∆x = 0.01 m. The bottom is set to The corresponding results for the free surface are given in Figure 12. The results generally agree with the experimental data again: the rarefaction wave profile is well reproduced, and the shock front position is well predicted, especially at t = 0.22 s and t = 0.32 s. Nevertheless, the amplitude of the leading solitary-type wave of the bore front is overestimated for t = 0.76 s. We remark again that no empirical wave breaking mechanism is considered here.

We further study the simulation of Favre waves generated in a laboratory flume using the experiments presented in [START_REF] Frazao | Undular bores and secondary waves -experiments and hybrid finite-volume modelling[END_REF][START_REF] Soares-Frazão | A second-order semi-implicit hybrid scheme for one-dimensional boussinesq-type waves in rectangular channels[END_REF]. In this case the domain [0, 30] is discretized using constant cells of ∆x = 0.01 m length. The bottom is set to b = 0, and the water height is h = 0.251. The rest of the flow variables were initially set to zero. A bore with a Froude number F = 1.104, where Results are shown in Figure 13 where again inviscid computations show excellent agreement with experimental data for all water level gauges.

F 2 = h u (h u + h d ) 2h 2

Conclusions

An efficient and new numerical approach to solve a vertically averaged and moment equations non-hydrostatic system (VAM model) is presented. This numerical approach is based on a two-step algorithm: first, the hyperbolic part of the system is solved with a second-order path-conservative PVM scheme. Second, non-hydrostatic terms are corrected by solving a linear Poisson-like system using an iterative method. The numerical approach has proven to be efficient, and the results are accurate and provide good to experimental data. Moreover, compared to the well-known SGN system, the scheme described here and applied to the VAM system has a similar computational effort, while the results are largely improved. Furthermore, the results given by VAM system are comparable to those obtained with a two-layer non-hydrostatic system, while being less expensive from the computational point of view. Therefore, using the VAM system combined with an efficient numerical technique is a reasonable choice to deal with non-hydrostatic simulations of free surface flows.

Figure 1 :

 1 Figure 1: Sketch of the definition of 2D free surface flow and σ-coordinate transformation

Figure 2 Figure 2 :

 22 Figure 2: (a), (b) Relative errors (in percentage) for the wave celerity C 2 and group velocity C 2 g respectively. (c) Comparison with the shoaling gradient γ.

Figure 3 :

 3 Figure 3: (a), (b), (c) weighed L 1 -norm errors as defined by (8) (in percentage) for the wave celerity C 2 , group velocity C 2g , and shoaling γ respectively for a given maximum wavenumber (kH) max .

Figure 4 :

 4 Figure 4: Steady flow over fixed bedforms. Blue solid lines correspond to numerical results, and star values represent experimental data. A Gaussian profile is given on the left column, while an asymmetric profile is on the right. The upper figures represent the numerical free surface flow over the bottom (solid black line). Lower figures correspond to the bed pressure.

Figure 5 :Figure 6 :

 56 Figure 5: Solitary wave propagation test for times: t = 25 s (upper left), t = 50 s (upper right), t = 100 s (lower left), and t = 150 s (lower right).

Figure 7 :Figure 8 :

 78 Figure 7: Sketch of the topography

Figure 9 :

 9 Figure 9: Comparison VAM (blue) and SGN (magenta) surface evolution at times t g/H = 15, 20.30, 21.87, 30.

Figure 10 :

 10 Figure 10: Periodic waves breaking over a submerged bar. Sketch of the topography and layout of the wave gauges

Figure 11 :

 11 Figure 11:Comparison of data time series (red star points) and numerical values at wave gaugesG 1 , G 2 , G 3 , G 4 , G 5 , G 6 , G 7 , G 8 .

d,

  was generated by imposing h = 0.285989 m on the left-hand side of the domain. Free-outflow boundary conditions for the rest of the flow variables were imposed everywhere.The evolution of the undular bore was obtained for water level gauges positioned at the following distances from the gate placed at x = 0: C 0 (x = 6.15 m), C 1 (x = 13.15 m), C 2 (x = 15.65 m), C 3 (x = 18.15 m), C 4 (x = 20.65 m), and C 5 (x = 23.15 m).

Figure 12 :

 12 Figure 12: Dam-break flow for at times t = 0.22 s (upper left), t = 0.32 s (upper right), t = 0.52 s (lower left), and t = 0.76 s (lower right).

  0

Figure 13 :

 13 Figure 13: Favre wave test for F = 1.104. Temporal evolution of water levels at gauges C 0 (x = 6.15 m), C 1 (x = 13.15 m), C 2 (x = 15.65 m), C 3 (x = 18.15 m), C 4 (x = 20.65 m), and C 5 (x = 23.15 m). Starred values correspond to experimental measurements and blue solid lines to numerical results.

  the group velocity C g = C +k∂ k C can be obtained from the linear dispersion relation for C, and it yields

	C 2 g √ gH	=	kH 14 35831808	+	kH 12 186624	kH 8 1296 + 31kH 10 + kH 6 108 82944 + 233kH 8 + kH 4 12 20736	+ +	kH 2 3 233kH 6 + 1 1728 +	31kH 4 48	+	4kH 2 3	+ 1

Acknowledgments

This work is partially supported by proyects RTI2018-096064-B-C2(1-2) and PID2020-114688RB-I00 funded by MCIN/AEI/10.13039/501100011033 and "ERDF A way of making Europe".


domain: phys spi

text:
Quentin Bennehard
email: quentin.bennehard@onera.fr

Jacques Peter
email: jacques.peter@onera.fr

Marco Carini
email: marco.carini@onera.fr

COMPARISON OF EFFICIENT UNCERTAINTY QUANTIFICATION TECHNIQUES ON A TRANSONIC AIRFOIL, USING RANS COMPUTATIONS AND ADJOINT METHOD

Keywords: transonic airfoil, CFD, geometric uncertainty, polynomial chaos, compressed sensing, adjoint gradient-enhancement

published or not. The documents may come

Comparison of efficient uncertainty quantification

techniques on a transonic airfoil, using rans computations and adjoint method

INTRODUCTION

During the last decades, computing power has largely increased, contributing to make Computational Fluid Dynamics (CFD) a quite mature tool for industrial design, widely exploited for various applications, especially in the aeronautical sector. However, the awareness of discrepancies between the ideal conditions of numerical simulations and the real ones has motivated the increasing attention towards sensitivities of classical optimal shapes to uncertain parameters robust design techniques [START_REF] Ghisu | Affordable Uncertainty Quantification for Industrial Problems: Application to Aero-Engine Fans[END_REF] and robust design techniques [START_REF] Sabater | Robust Design of Transonic Natural Laminar Flow Wings under Environmental and Operational Uncertainties[END_REF].

These uncertainties are usually of two kinds: epistemic and aleatory. The first one originates in our lack of knowledge in the modelling of the physical phenomena. In this context, Uncertainty Quantification (UQ) can help making our physical models more reliable. The second type is related to the variability in shape or flow conditions and although it cannot be controlled, it still needs to be considered in the design process. For example, deviations of real aerodynamic shapes from their reference design are often encountered, not only due to manufacturing tolerances but also to temporary and permanent degradation of aerodynamic surfaces along their lifespan.

Consequently, to be able to provide a robust aerodynamic shape by numerical optimization, uncertainties must be integrated in the design process. The state of the art of shape optimization provides, with the adjoint approach, the capability to tackle high dimensional design space [START_REF] Peter | Discrete adjoint method for shape optimization and mesh adaptation in the elsA code[END_REF]. However, on the other hand, the number of uncertain inputs can also become significantly large and, today, for such complex cases treated with High-Fidelity (HiFi) CFD, designers are struggling with what is known as the "curse of dimensionality". This curse represents the main bottleneck for the widespread application of UQ techniques in the industrial framework.

The main focus of this paper will be on efficient polynomial chaos techniques which can afford a reduced number of expensive CFD computations: compressed sensing (Least Angle Regression [START_REF] Blatman | Adaptive sparse polynomial chaos expansion based on least angle regression[END_REF], Basis Pursuit Denoising [START_REF] Chen | Basis pursuit[END_REF]) and adjoint-gradient enhanced variants [START_REF] Mura | Least squares approximation-based polynomial chaos expansion for uncertainty quantification and robust optimization in aeronautics[END_REF], [START_REF] Ghisu | Affordable Uncertainty Quantification for Industrial Problems: Application to Aero-Engine Fans[END_REF] of the standard least-square approximation. Their implementation will mainly rely on the functionalities already available in the open-source toolboxes OpenTURNS [START_REF] Baudin | OpenTURNS: An Industrial Software for Uncertainty Quantification in Simulation, Handbook of Uncertainty Quantification[END_REF] (http://openturns.github.io/) and eQuadrature [START_REF] Seshadri | Effective-Quadratures (EQ): Polynomials for Computational Engineering Studies[END_REF] (https://equadratures.org/). These techniques will be compared on the a realistic design test case represented by a transonic airfoil whose shape design parameters are considered uncertain with focus on efficiently estimating the statistics associated with its aerodynamic coefficients.

METHODOLOGY

Polynomial Chaos Expansion is a well establish approach [START_REF] Lüthen | Sparse polynomial chaos expansions: Literature survey and benchmark[END_REF] to derive a surrogate model approximating a function of interest 𝑓(𝒙) in the form of multivariate polynomials 𝐻 𝒋 (𝒙):

𝑓(𝒙) = ∑ 𝑐 𝒋 𝐻 𝒋 (𝒙), ∞ 𝒋 (1)
where 𝒙 is a d-dimensional vector of independent random variables, j the multi-index associated with univariate polynomials and 𝑐 𝒋 the polynomial coefficients. The polynomial expansion {H j (x)} j=0 ∞ is chosen as a complete orthogonal basis with respect to the inner product associated with the joint Probability Density Function (PDF) of 𝒙, 𝜌(𝒙):

〈𝐻 𝒋 , 𝐻 𝒌 〉 = ∫ 𝐻 𝒋 (𝒙) 𝜌(𝒙)𝐻 𝒌 (𝒙) = 𝛿 𝒋,𝒌 ,

with (𝒙) = ∏ 𝜌 𝑖 (𝑥 𝑖 ) 𝑑 𝑖=1

, 𝜌 𝑖 (𝑥 𝑖 ) being the marginal probability density function characterizing the i th uncertain variable 𝑥 𝑖 . Different truncation rules, e.g. total order or a hyperbolic cross [START_REF] Lüthen | Sparse polynomial chaos expansions: Literature survey and benchmark[END_REF], can be applied to Eq. ( 1), resulting in a different definition of the multi-index set 𝐽 spanned by j:

𝑓(𝒙) = ∑ 𝑐 𝒋 𝐻 𝒋 (𝒙). 𝐽 𝒋 (3)
In the present paper, the total order rule is adopted by fixing a common maximum degree 𝑝 for all the polynomial terms:

|𝒋| = ∑ 𝑗 𝑖 𝑑 𝑖=1 ≤ 𝑝, (4)
which leads to a total number n of polynomial terms (and hence of coefficients) to be determined given by n = (p + d)! (p! d!) ⁄ . The PCE coefficients can be computed by various methods. In particular, when considering non-intrusive techniques, the different approaches can be essentially divided into two families: i) Spectral projection methods, where the coefficients are computed by the numerical approximation of the projection of 𝑓(𝒙) on the orthogonal polynomial basis using classical Gauss quadrature and their more efficient sparse versions [START_REF] Gerstner | Numerical integration using sparse grids[END_REF]. ii) Collocation methods, where the PCE coefficients are computed by fitting Eq. ( 1) to a given number of collocation points in the uncertain space, to best reproduce the exact function value:

[ 𝐻 𝟎 (𝒙 1 ) … 𝐻 𝒏 (𝒙 1 ) ⋮ ⋱ ⋮ 𝐻 𝟎 (𝒙 𝑁 ) ⋯ 𝐻 𝒏 (𝒙 𝑁 ) ] ( 𝑐 𝟎 ⋮ 𝑐 𝒏 ) = ( 𝑓(𝒙 1 ) ⋮ 𝑓(𝒙 𝑁 )
), [START_REF] Baudin | OpenTURNS: An Industrial Software for Uncertainty Quantification in Simulation, Handbook of Uncertainty Quantification[END_REF] that can be written as 𝐴𝒄 = 𝒃 in compact form. Oversampling up to 2n is often employed to satisfy the full-rank condition for the above matrix, resulting in a classical Least Square Approximation (LSA) of the PCE solution. This is equivalent to seek an approximation for the PCE coefficients in the form of an l2-minimization problem:

𝒄 = 𝑎𝑟𝑔𝑚𝑖𝑛 ‖𝐴𝒄 -𝒃‖ 2 . (6)

Gradient-enhanced LSA

When gradient information is available, it can be used to introduce additional equations for the considered collations points while reducing their number in order to achieve full rank. In this case, the LSA system in Eq. ( 6) assumes the following form:

[ 𝐴 (0) ⋮ 𝐴 (𝑑) ] ( 𝑐 𝟎 ⋮ 𝑐 𝒏 ) = ( 𝒃 (0) ⋮ 𝒃 (𝑑)
), [START_REF] Destarac | Drag extraction from numerical solutions to the equations of fluid dynamics: the far-field philosophy in 43ème Congrès[END_REF] with:

𝐴 (𝑗) = [ 𝐻 𝟎 (𝑗) (𝒙 1 ) ⋯ 𝐻 𝒏 (𝑗) (𝒙 𝑁 ) ⋮ ⋱ ⋮ 𝐻 𝟎 (𝑗) (𝒙 1 ) … 𝐻 𝒏 (𝑗) (𝒙 𝑁⋯ ) ],
𝒃 (𝑗) = ( 𝑓 (𝑗) (𝒙 1 ) ⋮ 𝑓 (𝑗) (𝒙 𝑁 )

), [START_REF] Mura | Least squares approximation-based polynomial chaos expansion for uncertainty quantification and robust optimization in aeronautics[END_REF] Where 𝐴 (0) = 𝐴 and𝐴 (0) = 𝐴 𝐻 𝒌 (𝑗) = 𝜕𝐻 𝒌 𝜕𝑥 j ⁄ for 𝑗 ≥ 1𝑗 ≥ 1. Since the gradient at each collocation point provides d additional scalar information, the total number of samples can be ideally reduced by the same factor with respect to the standard LSA formulation. However, attention has to be paid to accuracy issues in computing the derivative information, as often occurs for adjoint solvers employed in CFD applications. In order to reduce the sensitivity from numerical noise affecting the gradient accuracy, a null-space method has been proposed in [START_REF] Ghisu | Gradient-enhanced least-square polynomial chaos expansions for uncertainty quantification and robust optimization[END_REF] showing how this approach can achieve improved accuracy, especially for standard deviation, on a turbomachinery test case. The same method as implemented in the eQuadratures toolbox is employed here to compute the solution of the Gradient-Enhanced LSA problem (LSA-GE) presented in Eq. ( 7).

Compressed sensing techniques

When dealing with very limited amount of data (e.g. as a direct consequence of expensive function evaluations) resulting in an underdetermined system for Eq. ( 1), an alternative solution can be sought in a sparse form, i.e. by minimizing the number of non-zero entries of 𝒄:

𝒄 = 𝑎𝑟𝑔𝑚𝑖𝑛 ‖𝒄‖ 0 𝑠𝑢𝑐ℎ 𝑡ℎ𝑎𝑡 ‖𝐴𝒄 -𝒃‖ 2 ≤ 𝜀, (9)
where 𝜀 ≥ 0 is a given tolerance on the interpolation condition and ‖𝒄‖ 0 denotes the l0 norm of 𝒄. Indeed, although real PCE models are not truly sparse, they are expected to be compressible, i.e. to feature a rapid decay in the magnitude of the coefficients at increasing order of the expansion, with most of the variance being captured by a few terms. This provides the rational behind the compressed sensing approach. Unfortunately, the cost of solving Eq. ( 9) grows exponentially with d and a convex relaxation based on l1-minimization is often considered in practice. When 𝜀 = 0 this is also referred to as Basis Pursuit Denoising (BPDN) which leads to solve the following problem:

𝒄 = 𝑎𝑟𝑔𝑚𝑖𝑛 ‖𝒄‖ 1 𝑠𝑢𝑐ℎ 𝑡ℎ𝑎𝑡 ‖𝐴𝒄 -𝒃‖ 2 = 0, (10)
using convex optimization algorithm. Besides convex optimization solvers, greedy methods have also been proposed to find a sparse PCE solution [START_REF] Lüthen | Sparse polynomial chaos expansions: Literature survey and benchmark[END_REF]. In these methods regressors are added to the model one by one according to some selection criterion, to find a heuristic solution to the burdensome l0 minimization problem. An example is represented by LARS [START_REF] Blatman | Adaptive sparse polynomial chaos expansion based on least angle regression[END_REF] where the regressors are added according to their correlation with the current residual and the PCE coefficients are then updated using a least angle strategy. It should be noted that LARS with the LASSO (Least Absolute Shrinkage and Selection Operator) modification can be interpreted as a l1-optimization solver. Both techniques, BPDN and LARS will be considered in the present study. In particular, for LARS, the implementation available in the open-source library Open-TURNS [START_REF] Baudin | OpenTURNS: An Industrial Software for Uncertainty Quantification in Simulation, Handbook of Uncertainty Quantification[END_REF] is employed.

THE RAE2822 TEST CASE

The considered test case is represented by the compressible viscous flow around the RAE2822 airfoil at transonic conditions: Mach number of 0.725, Reynolds number of 6.5e6, and angle of attack of 2.55°: these conditions corresponds to the AGARD test case 7 [START_REF] Cook | Aerofoil RAE2822 -Pressure Distributions, and Boundary Layer and Wake Measurements[END_REF]. The RANS solution is computed by using the elsA software [START_REF] Cambier | The Onera elsA CFD Software: Input from Research and Feedback from Industry[END_REF] (ONERA-SAFRAN property) on a structured grid with a CH topology featuring a total of ~8e5 cells and near-wall resolution of 𝑦 + ~ 1.0. The CFD solution is converged down to 7 orders of magnitude in the decrease of the residual. The employed mesh and the computed Mach field are illustrated in Figure 1. A far-field drag analysis [START_REF] Destarac | Drag extraction from numerical solutions to the equations of fluid dynamics: the far-field philosophy in 43ème Congrès[END_REF] is applied to CFD results in order to extract the different phenomenological sources of drag. In particular, the far-field estimation of the total drag coefficient CDff and the wave drag coefficient, CDw, are considered as Quantity of Interest (QoI) for the present UQ analysis, in addition to the lift, CL, and pitching moment, CMy, coefficients. Corresponding values for the reference airfoil geometry and the considered flight conditions are reported in Table 1.

CL

Uncertain shape parameterization

The airfoil shape is parameterized by means of both camber and thickness parameters at five different positions uniformly distributed along the chord, i.e. at 16.7%, 33.3%, 50%, 66.7%, 83.3%, as illustrated in Figure 2. The deviation of each parameter from its corresponding baseline value is assumed to follow a zero-mean β distribution. The convention used to define the probability density function of the β distribution is recalled Eq. ( 11):

𝑓 𝑋 (𝑥) = (𝑥-𝑎) 𝛼-1 (𝑏-𝑥) 𝛽-1 (𝑏-𝑎) 𝛼+𝛽-1 𝐵(𝛼,𝛽)
with 𝑥 ∈ [𝑎, 𝑏] and 𝐵 being the Euler's β function [START_REF] Peter | Discrete adjoint method for shape optimization and mesh adaptation in the elsA code[END_REF] For all camber parameters the same value of variance is employed. The same assumption applies to the thickness parameters. Two different levels of amplitude, namely Ω1 and Ω2 are considered for the uncertain perturbations: the associated sets of β-PDF parameters are reported in Table 2, while a representation of the envelope of the resulting airfoil shape is illustrated in Figure 3 and Figure 4 for Ω1 and Ω2 cases, respectively. Not surprisingly, the envelope is wider for Ω2 and as a consequence, the pressure coefficient distribution shows quite large deviations from the reference results, especially affecting the shock position on the suction side of the airfoil, as illustrated in Figure 4. For each realization of the uncertain design parameters, a mesh deformation process is used to adapt the baseline airfoil mesh and re-run CFD computations for both primal and adjoint analysis.

UQ RESULTS AND DISCUSSION

The comparison of UQ methods presented in this section is limited to second-order PCE approximation, as often adopted in robust optimization study for realistic applications [START_REF] Ghisu | Gradient-enhanced least-square polynomial chaos expansions for uncertainty quantification and robust optimization[END_REF]. For the considered 10 uncertain parameters, this results in a total of 66 coefficients to be computed. The standard LSA solution is employed here as a reference to asses the accuracy of the results obtained by LSA-GE, LARS and BPDN with a limited amount of data as well their potential gain in terms of efficiency. The reduced number of samples is defined as a minimum needed to attain a prescribed precision computed by the Q2 criteria. The efficiency gain is heuristically defined as the ratio between the number of samples used by LSA and their reduced number required by the other techniques, also including adjoint computations. In particular, we assume that the cost of computing a gradient for a single aerodynamic coefficient equals the one of a standard CFD run as done by [START_REF] Ghisu | Gradient-enhanced least-square polynomial chaos expansions for uncertainty quantification and robust optimization[END_REF] and [START_REF] Mura | Least squares approximation-based polynomial chaos expansion for uncertainty quantification and robust optimization in aeronautics[END_REF]. However, for some complex applications, this assumption could lead to underestimate the effective cost of adjoint CFD computations, which are known to be less robust than primal ones.

For the two sets of uncertain amplitudes, Ω1 and Ω2, a DoE is generated using 1024 samples distributed according to a Sobol sequence. For each airfoil shape realization, a CFD evaluation of CDff, CDw, CL, CMy and the associated adjoint computations are carried out, resulting in a total of 5024 CFD runs. All the 1024 function evaluations are employed to feed the LSA problem. The values of mean and standard deviation estimated by using standard LSA for the different aerodynamic coefficients are reported in Table 3.

QoI

Mean and variance estimation

The comparison of mean and variance estimations obtained by the enhanced LSA methods for CDff, CDw, CL and CMy is summarized in Table 4-8, 5-9, 6-10 and 7-11, respectively. All the meta-models show very good accuracy on the mean estimations while the accuracy is reduced for the variance up to a difference of ~7% with respect to the LSA reference value. For the Ω1 case, efficiency gains are greater for CDff and CMy while almost halve for CDw and CL. LARS method achieves the largest gain for each coefficient except for the wave drag. The gain is almost doubled for CDff and CL with respect to LSA-GE. BPDN achieves a performance similar to that of LARS for CDw and CMy. A similar behavior is also observed for the Ω2 case, with LARS and BPDN showing a noticeably better performance than LSA-GE, with an efficiency gain which is improved by a factor ranging from ~2 up ~4 with respect to LSA-GE. The resulting PDFs are also compared in Figure 5 and Figure 6 for case Ω1 and Ω2, respectively. For case Ω1 we can observe that all the PDF functions are substantially symmetric except for the pitching moment coefficient. The differences among the different methods are quite limited. The obtained best using analytical PDF laws are reported in Table 8, using the Open-TURNS notation. Differently from case Ω1, for case Ω2 CDff and CDw PDFs are characterized by a non-negligible skewness while those of CL and CMy are almost symmetric as confirmed by fitting results shown in Table 13. In addition, compared to case Ω1, larger discrepancies are observed among the different methods, especially in CDff and CDw. For CDw, we also observe a non-physical tail of the PDF at values lower than 0.

PCE method

Sensitivity analysis

The sensitivity analysis based on PCE first order and total Sobol indices is presented in Figure 7 and Figure 8, for case Ω1 and Ω2, respectively. These results indicate that the camber variable at 83% of the chord has the greatest influence on the variance of CL and CMy while all thickness variables have a negligible effect. Conversely, for both drag coefficients, the most influent design variables are represented by thickness and camber at 17% followed by those at the 33% of the chord. A lower contribution is also observed for the camber at the 83% of the chord. All thickness variables, except the last one, show a non-negligible contribution to the variance of CDw, which is not surprising from physical viewpoint. When analyzing interaction effects, we can observe that for case Ω1, they are almost negligible: minor deviations of total Sobol index from first order ones are only observed for the wave drag coefficients. The differences between the considered PCE methods are also very limited. On the contrary, for case Ω2, relevant interaction effects are observed also for the far-field drag coefficient.

CONCLUSIONS

This work treats two uncertainty quantification problems in dimension d=10. Four methods of computing a polynomial chaos expansion have been compared. The main purpose is to try to minimize the computational cost required to accurately estimate the mean and variance of four aerodynamic coefficients on a transonic airfoil. It is shown that 'gradient enhanced' and 'compress sensing' methods present significant performance gain when compared to the standard least square approximation. More precisely, the LARS method, which looks for a sparse PCE, achieve the most significant gain across all the test cases. In addition, it was found that the gain provided by LSA-GE method was sometimes lower than d, especially when dealing with larger variance of the uncertain parameter in test case Ω2.

The achieved reduction in the size of the DoE by means of these advanced PCE techniques pave the way to a robust shape optimization using HiFi computations, building a polynomial chaos expansion at each step of the optimization. Further improvements are necessary, e.g. by combining these techniques with the dimension reduction of the uncertain design space, in order to further reduce the cost of the UQ surrogate model.

Figure 1 :

 1 Figure 1: Computed Mach field for the transonic flow around the RAE2822 (right) and employed structured mesh (left).

Figure 2 :

 2 Figure 2: RAE2822 baseline airfoil geometry (light blue line) and position of the control points used to modify the camber and thickness laws (orange makers and line).

Figure 3 :

 3 Figure 3: Airfoil shape envelope (left) and corresponding variations of the pressure coefficient distribution (right) for the DoE associated with the uncertain set Ω1.

Figure 4 :

 4 Figure 4: Airfoil shape envelope (left) and corresponding variations of the pressure coefficient distribution (right) for the DoE associated with the uncertain set Ω2.

Figure 5 :

 5 Figure 5: Comparison of the histogram of the training data and the PDFs resulting from the different PCE methods for the case Ω1.

Figure 6 :

 6 Figure 6: Comparison of the histogram of the training data and the PDFs resulting from the different PCE methods for the case Ω2.

Figure 7 :

 7 Figure 7: Comparison of the Sobol indices obtained from the different PCE methods for the case Ω1.

Figure 8 :

 8 Figure 8: Comparison of the Sobol indices obtained from the different PCE methods for the case Ω2.

Table 1 :

 1 Computed values of the aerodynamic coefficients for the RAE2822 airfoil at considered transonic conditions. Drag coefficient values are reported in drag counts.

		CDff	CDw	CMy
	0.8076	139.37	24.20	-0.09465

Table 2 :

 2

β-PDF parameters for camber and thickness uncertain design variables.

Table 3 :

 3 Reference values of mean and standard deviation obtained by using standard LSA.

Table 4 :

 4 Comparison of the different methods on the estimation of the CDff statistics for uncertain set Ω1.

		Mean	Mean dif-ference	Sigma	Sigma differ-ence	Q2	DoE size	Perfor-mance gain
	LSA-reference	139.92 -	3.99 -	-		1023 -
	LSAGE	139.97	0.03%	4.21	5.51% 0.979	30	17.1
	LARS	139.79	-0.09%	4.02	0.86% 0.977	33	31.0
	BPDN	139.76	-0.11%	3.96	-0.84% 0.975	42	24.4
	PCE method	Mean	Mean dif-ference	Sigma	Sigma differ-ence	Q2	DoE size	Perfor-mance gain
	LSA-reference	25.12 -	2.88 -	-		1023 -
	LSAGE	25.18	0.25%	2.88	0.18% 0.934	66	7.8
	LARS	25.07	-0.21%	2.75	-4.34% 0.926	188	5.4
	BPDN	25.14	0.10%	2.81	-2.37% 0.924	182	5.6

Table 5 :

 5 Comparison of the different methods on the estimation of the CDw statistics for uncertain set Ω1.

	PCE method	Mean	Mean dif-ference	Sigma	Sigma difference	Q2	DoE size	Perfor-mance gain
	LSA-reference	0.809 -	6.02E-03 -	-	1023 -
	LSAGE	0.809	0.01%	5.86E-03	-2.68% 0.958	66	7.8
	LARS	0.808	-0.03%	5.98E-03	-0.66% 0.957	69	14.8
	BPDN	0.809	0.01%	6.11E-03	1.48% 0.949	132	7.8

Table 6 :

 6 Comparison of the different methods on the estimation of the CL statistics for uncertain set Ω1.

	PCE method	Mean	Mean dif-ference	Sigma	Sigma differ-ence	Q2	DoE size	Perfor-mance gain
	LSA-reference	-0.095 -	1.84E-03 -	-		1023 -
	LSAGE	-0.095	0.05% 1.86E-03	0.84% 0.964	35	14.6
	LARS	-0.095	-0.08% 1.71E-03	-7.22% 0.963	41	25.0
	BPDN	-0.095	-0.03% 1.75E-03	-5.33% 0.951	43	23.8

Table 7 :

 7 Comparison of the different methods on the estimation of the CMy statistics for uncertain set Ω1.

Table 8 :

 8 Parameters of the fitted distributions obtained by the LSA reference PCE method for the case Ω1

	QoI	Mean	Sigma		Skewness Fitted PDF Fitted PDF Parameters
	CDff	139.92	3.99	0.203 LogNormal μ =-5.23, σ=7.44e-2, γ=8.64e-3
	CDw	25.12	2.88	0.232 Dirichlet	α=[75.98, 30174.1]
	CL	0.8087 6.02E-03	0.08103 Beta	α=5.51, β=5.79, a=0.79, b=0.83
	CMy	-0.09495 1.84E-03	-0.224 Beta	α=6.78, β=4.51, a=-0.103, b=-0.0897
	PCE method	Mean	Mean dif-ference	Sigma	Sigma differ-ence	Q2	DoE size	Perfor-mance gain
	LSA-reference	146.80 -		24.13 -	-		1024 -
	LSAGE		146.52		-0.19%	24.45	1.29% 0.985	66	7.8
	LARS		146.22		-0.40%	23.74	-1.63% 0.975	92	11.1
	BPDN		146.46		-0.24%	24.24	0.45% 0.976	93	11.0

Table 9 :

 9 Comparison of the different methods on the estimation of the CDff statistics for uncertain set Ω2.

	PCE method	Mean	Mean dif-ference	Sigma	Sigma differ-ence	Q2	DoE size	Perfor-mance gain
	LSA-reference	29.97 -	16.76 -	-		1024 -
	LSAGE	29.50	-1.58%	16.42	-2.04% 0.939	66	7.8
	LARS	29.12	-2.86%	16.48	-1.65% 0.937	54	19.0
	BPDN	29.76	-0.71%	16.48	-1.65% 0.930	77	13.3

Table 10 :

 10 Comparison of the different methods on the estimation of the CDw statistics for uncertain set Ω2.

	PCE method	Mean	Mean differ-ence	Sigma	Sigma differ-ence	Q2	DoE size	Perfor-mance gain
	LSA-reference	0.803 -		3.69E-02 -	-		1024 -
	LSAGE	0.803	-0.01%	3.51E-02	-4.91% 0.965	66	7.8
	LARS	0.805	0.31%	3.70E-02	0.26% 0.953	30	34.1
	BPDN	0.805	0.29%	3.60E-02	-2.43% 0.954	27	37.9

Table 11 :

 11 Comparison of the different methods on the estimation of the CL statistics for uncertain set Ω2.

	PCE method	Mean	Mean dif-ference	Sigma	Sigma differ-ence	Q2	DoE size	Perfor-mance gain
	LSA-reference	-0.095 -	1.12E-02 -	-		1024 -
	LSAGE	-0.095	0.31%	1.13E-02	0.24% 0.952	29	17.7
	LARS	-0.095	0.02%	1.09E-02	-2.95% 0.984	35	29.3
	BPDN	-0.095	-0.05%	1.05E-02	-6.59% 0.958	35	29.3

Table 12 :

 12 Comparison of the different methods on the estimation of the CMy statistics for uncertain set Ω2.

Table 13 :

 13 Parameters of the fitted distributions obtained on the LSA reference PCE method for the case Ω2.

	QoI	Mean	Sigma	Skewness	Fitted Dis-tribution	Distribution Parameters
	CDff	146.79	24.13	0.748 LogNormal μ =-4.86, σ=0.295, γ=6.578e-3
	CDw	29.96	16.76	0.692 Beta	α=2.71, β=6.43, a=-4.65e-4, b=1.12e-2
	CL	0.8029 3.68E-02	-0.0527 Beta	α=5.12, β=4.95, a=0.678, b=0.923
	CMy	-0.094812 1.12E-02 0.01716 Beta	α=5.53, β=5.03, a=-0.135, b=-0.0584

ACKNOLEDGMENTS

This study is supported by the European project NEXTAIR, ONERA's internal project HORUS and previous work at ONERA, studying efficient quadrature on analytical test cases [14]. The NEXTAIR project has received funding from the European Union's Horizon Europe re-search and innovation program under grant agreement No 101056732. Views and opinions expressed are however those of the author(s) only and do not necessarily reflect those of the European Union. Neither the European Union nor the granting authority can be held responsible for them.


Maïte ´vicre ´1

Vincenzo Lionetti
email: vincenzo.lionetti@uniroma1.it

Lal Prem

Kashyap

Plant cell wall in pathogenesis, parasitism and symbiosis, Volume II

Keywords: cell wall remodeling, symbiosis, cell wall integrity, plant immunity, plant parasitism, plant cell wall

A wide range of organisms that interact with plants must interface with the plant cell wall (CW) [START_REF] Lionetti | Plant cell wall in pathogenesis, parasitism and symbiosis[END_REF]. The view of the CW as only a static cellular barrier in these interactions is outdated. Cell wall polysaccharides, phenolic compounds, and proteins, in addition to regulating important growth and development processes, are also sources of elicitors that activate cell signaling pathways [START_REF] Nguema-Ona | Arabinogalactan proteins in root-microbe interactions[END_REF]. Surveillance mechanisms detect CW contacts with other organisms, and specific signaling pathways and responses are activated [START_REF] Swaminathan | Plant cell wall integrity perturbations and priming for defense[END_REF] (Figure 1). During biotic interactions and abiotic stresses, the structure and composition of plant CW can be regulated at the biosynthetic level and through precise, continuous post-synthetic remodeling. As a consequence, the CW must be understood as a strategic space between organisms where intelligent and dynamic molecular strategies are implemented to overwhelm a fight or cooperate for specific physiological processes [START_REF] Bacete | Plant cell wall-mediated immunity: cell wall changes trigger disease resistance responses[END_REF][START_REF] Castilleux | Cell wall extensins in root-microbe interactions and root secretions[END_REF][START_REF] De Lorenzo | Cell wall traits that influence plant development, immunity, and bioconversion[END_REF].

Cell wall enzymes and their inhibitors play key roles in apoplastic metabolism [START_REF] Rui | A wall with integrity: surveillance and maintenance of the plant cell wall under stress[END_REF]. Invertases (INVs) and pectin methylesterases (PMEs) play essential roles in carbohydrate metabolism, stress responses, and sugar signaling [START_REF] Bellincampi | Plant cell wall dynamics and wall-related susceptibility in plant-pathogen interactions[END_REF][START_REF] Tauzin | Sucrose and invertases, a part of the plant defense response to the biotic stresses[END_REF][START_REF] Del Corpo | AtPME17 is a functional arabidopsis thaliana pectin methylesterase regulated by its PRO region that triggers PME activity in the resistance to botrytis cinerea[END_REF]. In this collection, Coculo and Lionetti reviewed the roles of invertase inhibitors (INVI) and pectin methylesterase inhibitors (PMEI) belonging to the "Plant Invertase/Pectin Methylesterase Inhibitor Superfamily" [START_REF] Lionetti | Three pectin methylesterase inhibitors protect cell wall integrity for arabidopsis immunity to botrytis[END_REF]). An updated overview of the specific activity of the characterized isoforms, their specific functions in plant physiology, and their applications in biotechnology is provided. After the pioneering work in 2007 [START_REF] Lionetti | Overexpression of pectin methylesterase inhibitors in arabidopsis restricts fungal infection by botrytis cinerea[END_REF], several piece of evidence supported the role of PMEIs in plant resistance to stresses [START_REF] An | Pepper pectin methylesterase inhibitor protein CaPMEI1 is required for antifungal activity, basal disease resistance and abiotic stress tolerance[END_REF][START_REF] Liu | A pectin methylesterase inhibitor enhances resistance to verticillium Wilt1[OPEN[END_REF]. With a genome-wide analysis and transcriptomics of the PMEI genes in Brassica napus, Wang et al. identified several BnPMEIs as resistance gene candidates in response to Sclerotinia sclerotiorum, suggesting them as possible tools to breed new and improved genotypes more resistant to Sclerotinia stem rot. Specific CW changes can occur in the vicinity of plasmodesmata during viral infection [START_REF] Lionetti | How do pectin methylesterases and their inhibitors affect the spreading of tobamovirus?[END_REF][START_REF] Stavolone | Extracellular matrix in plants and animals: hooks and locks for viruses[END_REF]. Viruses can modify pectin, callose, and structural proteins near plasmodesmata by increasing their size exclusion limit [START_REF] Kozieł | Plant cell wall as a key player during resistant and susceptible plant-virus interactions[END_REF]. Loṕez-Gonzaĺez et al. addressed the little-explored impact of virus infections on secondary CW and plant development. The authors found a correlation between the developmental alterations induced in Arabidopsis by specific strains of turnip mosaic virus and specific changes in xylan and lignin biosynthesis. Although callose deposition is a response to both elicitors and pathogens, the mechanisms involved in its biosynthesis and degradation remain to be unraveled (German et al., 2023). Li et al. summarized the research progress on plant callose and its synthesizing enzymes in plant physiology.

Microfungi are also important etiological agents for plants [START_REF] Doehlemann | Callose metabolism and the regulation of cell walls and plasmodesmata during plant mutualistic and pathogenic interactions[END_REF]. The fungal pathogen Sphaerulina musiva causes stem canker with the consequent mortality of Populus trees. Populus deltoides can induce a lignified periderm to contain the pathogen, but the precise characterization of lignin changes in response to S. musiva infection is unknown. Bryant et al.identified a higher syringyl:guaiacil ratio, a higher Klason lignin content and lower p-hydroxybenzoate content in Septoria-infected P. deltoides trees compared to the healthy plant. This knowledge can favor biotechnological approaches aimed at improving the resilience and increasing the biomass yield of Populus for biofuel production.

Interesting contributions in this collection concern plant-parasite interactions. A fine-tuned re-arrangement of host CW is induced in response to infection by both plant-parasitic cyst nematodes and root-knot nematodes [START_REF] Zhang | The complex cell wall composition of syncytia induced by plant parasitic cyst nematodes reflects both function and host plant[END_REF][START_REF] Bozbuga | Host-specific signatures of the cell wall changes induced by the plant parasitic nematode, meloidogyne incognita[END_REF][START_REF] Meidani | Cell wall modifications in giant cells induced by the plant parasitic nematode meloidogyne incognita in wild-type (Col-0) and the fra2 arabidopsis thaliana katanin mutant[END_REF]. Veronico et al. found that drought stress affects CW metabolism in tomato roots, limiting feeding site development and reproduction of the nematode Meloidogyne incognita. Parasitic plants, such as Cuscuta species, severely damage economically important crops [START_REF] Jhu | Parasitic plants: an overview of mechanisms by which plants perceive and respond to parasites[END_REF]. These green parasites absorb resources through an invasive organ called the haustorium, which differentiates into vascular hyphae that establish a connection with the host plant's vasculature. The degradation and modification of host CWs allow haustorium to effectively invade host tissues. Yokoyama et al. propose that Cuscuta campestris APETALA2/ETHYLENE RESPONSE FACTORS (ERFs) can activate the transcription of the CW enzymatic genes in haustorium to favor its invasion of tobacco and Arabidopsis plants.

Plant CW is a field of molecular dialogues and agreements with symbiotic microbes to establish intimate interfaces for developmental coordination and nutrient exchange [START_REF] Balestrini | Cell wall remodeling in mycorrhizal symbiosis: a way towards biotrophism[END_REF]. Plants can establish mutualistic symbiosis with arbuscular mycorrhizal fungi, and phosphorus transfer across the CW specialized interfacial compartment is an important process in the mycorrhizal pathway [START_REF] Begum | Role of arbuscular mycorrhizal fungi in plant growth regulation: implications in abiotic stress tolerance[END_REF]. Exploiting the Rhizophagus irregularis-Lotus japonicus interaction, Nguyen and Saito showed that polyphosphate in fungal CWs and apoplastic phosphatases play an important role in phosphorus transfer at the symbiotic interface in arbuscules. Mycorrhizal fungi can be involved in mutualistic interactions during orchid seed germination [START_REF] Põlme | Host preference and network properties in biotrophic plant-fungal associations[END_REF]. Chen et al. identified several genes codifying CW structural proteins such as epidermisspecific secreted glycoprotein, proline-rich receptor-like protein, and leucine-rich repeat (LRR) extensin-like protein, which are particularly involved in the symbiosis of Tulasnella and Serendipita fungi with Dendrobium officinale.

Bioinoculants represent an environmentally-friendly agricultural practice to alleviate drought stress in crops [START_REF] Kour | Drought adaptive microbes as bioinoculants for the horticultural crops[END_REF] These important contributions advance our understanding of the relationships between plants and the environment at the CW interface, which will be helpful to engineer biotechnological strategies for agriculture and bioenergy fields. As a closing remark, we are grateful to the authors and reviewers for their invaluable contributions to this Research Topic.

  . The work presented by Wilmowitz et al. indicates that the inoculation of maize seeds with Glomus sp. and Bacillus sp. can help to cope with drought stress, preventing inhibition of photosynthesis and disruption of redox balance. Glomus sp. and Bacillus sp. can modify pectin methylesterification and hemicellulose content of maize leaves, possibly leading to alleviation of the negative effects of drought.

FIGURE 1

 1 FIGURE 1 Schematic representation of the main systems of perception and responses related to the plant cell wall (CW) in the different biotic interactions. CW can play a central role in combined biotic and abiotic stress. Symbiont-induced CW modifications can improve plant development, nutrition, and tolerance to abiotic stresses. Abiotic stresses can alter CW composition, affecting efficient parasitism. M/H/DAMPS=Microbe/Herbivore/Damage Associated Molecular Patterns.

Frontiers in Plant Science frontiersin.org

Funding

The research was supported by Sapienza University of Rome, Grants RM120172B78CFDF2, RM11916B7A142CF1, RG12117A898EABE0, RM122181424F1F42 and -Rome Technopole" -European Union Next-GenerationEU PNRR -RT12218451D2C509 to VL.

Author contributions

VL draft the editorial text. MV revised and approved the final version of the editorial. All authors approved the submitted version.

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher's note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.