T Song

A Main

G Scovazzi

M Ricchiuto

The shifted boundary method for hyperbolic systems: Embedded domain computations of linear waves and shallow water flows

Keywords: Embedded boundary conditions, wave equation, shallow water flows, finite elements, approximate boundary methods conditions aux bords embedded, équations des ondes, équations shallow water, éléments finis, méthodes immergés

We propose a new computational approach for embedded boundary simulations of hyperbolic systems. Applications are shown for the linear wave equations and for the nonlinear shallow water system. The proposed approach belongs to the class of surrogate/approximate boundary algorithms and is based on the idea of shifting the location where boundary conditions are applied from the true to a surrogate boundary. Accordingly, boundary conditions, enforced weakly, are appropriately modified to preserve optimal error convergence rates. This framework is applied here in the setting of a stabilized finite element method, even though other spatial discretization techniques could have been employed. Accuracy, stability and robustness of the proposed method are tested by means of an extensive set of computational experiments for the acoustic wave propagation equations and shallow water equations. Comparisons with standard weak boundary conditions imposed on grids that conform to the geometry of the computational domain boundaries are also presented.

Introduction

In [START_REF] Main | The shifted boundary method for embedded domain computations. Part I: Poisson and Stokes problems[END_REF], the authors introduced an embedded method for the Poisson and Stokes problems, using an approximate (surrogate) domain representation and shifting the location where boundary conditions are imposed from the true to the surrogate boundary. In the present work, we extend those ideas to the case of hyperbolic systems and in particular the case of wave propagation problems in acoustics and shallow water flows.

Part of the motivation of this work is that time-domain acoustics in complex geometry still presents some computational challenges, and similarly for shallow water flows in the presence of very complex boundary/coastline features.

Immersed and embedded boundary methods offer advantages over geometrically conformal computational methods (i.e., methods in which numerical grids conform to the geometrical shape of the physical domains of interest) in that the task of mesh generation is greatly simplified. On the other hand, the enforcement of boundary conditions is more challenging in embedded methods, both from the mathematical and data structure perspectives.

In immersed boundary methods, the flow equations are discretized continuously both inside and outside the physical domain, and a smooth approximation to the Dirac delta function is introduced on the physical boundary, with the purpose of imposing the boundary conditions. Immersed boundary methods are older than embedded methods, dating to the seminal work of Peskin [START_REF] Charles S Peskin | Flow patterns around heart valves: a numerical method[END_REF], and being applied in the finite element context by Boffi and Gastaldi [START_REF] Boffi | A finite element approach for the immersed boundary method[END_REF]. It was observed in [START_REF] Hansbo | An unfitted finite element method, based on nitsche's method, for elliptic interface problems[END_REF] that immersed methods can be interpreted in the context of variational formulations as penalty methods, and for this reason share with them some of the shortcomings.

In embedded methods, the equations are discretized and solved only on the physical domain, and external regions are excluded from the computation. These methods generally do not have the problems mentioned for immersed boundary discretizations, as they employ a "sharp interface approximation." Within the finite element context, this is typically done by means of weak enforcement of boundary conditions through Nitsche's method [START_REF] Nitsche | Uber ein Variationsprinzip zur Losung Dirichlet-Problemen bei Verwendung von Teilraumen, die keinen Randbedingungen unteworfen sind[END_REF] in combination with XFEM strategies, as a way to appropriately construct the solution's approximation spaces [START_REF] Hansbo | An unfitted finite element method, based on nitsche's method, for elliptic interface problems[END_REF].

In contrast to a simple penalty approach, Nitsche's method provides a mechanism to weakly and consistently enforce boundary conditions and does not adversely affect the conditioning of the discretized problem. Unfortunately, the standard Nitsche-XFEM method suffers from instabilities on elements that are cut by the interface in such a way that only a small fraction of them remains inside the physical domain, and consequently their effective mesh size becomes extremely small. As a consequence, these areas of small support require an arbitrarily large penalty parameter, thereby destroying the good conditioning properties of Nitsche's method. A similar challenge exists for embedded boundary methods of finite volume type, which can suffer from a similar though technically different version of the small cut cell problem.

A variety of clever techniques have been attempted to circumvent this difficulty [START_REF] Hashem M Mourad | A bubble-stabilized finite element method for dirichlet constraints on embedded interfaces[END_REF][START_REF] Annavarapu | A robust Nitsche's formulation for interface problems[END_REF][START_REF] Annavarapu | An Efficient Finite Element Method for Interface Problems[END_REF], but in most cases only within the context of interior interface problems. For single material problems, Burman [START_REF] Burman | Ghost penalty[END_REF] introduced the ghost penalty method, in which the variational form is stabilized by introducing a penalization of the solution gradients at the interface separating cut and uncut elements. This method was applied to the Stokes problem in [START_REF] Burman | Fictitious domain methods using cut elements: Iii. a stabilized nitsche method for stokes' problem[END_REF], to the Navier-Stokes equations in [START_REF] Schott | A new face-oriented stabilized xfem approach for 2d and 3d incompressible navier-stokes equations[END_REF], and to two-phase flow in [START_REF] Schott | A face-oriented stabilized nitsche-type extended variational multiscale method for incompressible two-phase flow[END_REF]. The ghost penalty method has some drawbacks, in that the introduction of a fourth order operator in the ghost penalty term may have a delicate implementation in the nonlinear case, may negatively affect the overall condition number, and increases the stencil size when using low order finite elements, with additional complications in the case of parallel computing. An alternative approach introduced for B-spline variational formulations, known as Extended B-splines [START_REF] Höllig | Finite element methods with B-splines[END_REF][START_REF] Höllig | Weighted extended b-spline approximation of dirichlet problems[END_REF], involves eliminating via an extrapolation procedure those cut B-splines with small support. This technique was applied to the Navier-Stokes equations for moving boundary problems in [START_REF] Rüberg | Subdivision-stabilised immersed b-spline finite elements for moving boundary flows[END_REF] and [START_REF] Rüberg | A fixed-grid b-spline finite element technique for fluidstructure interaction[END_REF].

One additional challenge for the classes of methods just mentioned is that they require the geometric construction of the partial elements cut by the embedded boundary, Inria typically a complex and computationally intensive process. Since some sort of adaptive quadrature is used [START_REF] Parvizian | Finite cell method[END_REF][START_REF] Düster | The finite cell method for three-dimensional problems of solid mechanics[END_REF], it is often the case that a non-negligible portion of the overall wall-clock time for a simulation is spent handling the embedded boundary.

The small cut cell problem can be circumvented through the introduction of an approximate domain method, in which the true domain is replaced by a surrogate domain. Boundary conditions are imposed along the frontier of the surrogate domain, whose geometry is chosen to avoid cut cells. The challenge then reduces to designing effective (i.e., accurate and robust) boundary conditions on the surrogate boundary. One of the earliest approximate domain methods was proposed in [START_REF] Fedkiw | A non-oscillatory eulerian approach to interfaces in multimaterial flows (the ghost fluid method)[END_REF] for inviscid multiphase compressible flow. Referred to as the Ghost Fluid Method (GFM), it was later applied to multiphase compressible flow in [START_REF] Farhat | A higher-order generalized ghost fluid method for the poor for the three-dimensional two-phase flow computation of underwater implosions[END_REF][START_REF] Abgrall | Computations of compressible multifluids[END_REF], to compressible fluid-structure interaction in [START_REF] Wang | A computational framework based on an embedded boundary method for nonlinear multi-phase fluid-structure interactions[END_REF], and to multiphase fluid/structure interaction in [START_REF] Wang | Predictive simulation of underwater implosion: Coupling multi-material compressible fluids with cracking structures[END_REF]. An approximate domain approach has also been utilized for viscous incompressible flow in [START_REF] Kang | Prediction of wall-pressure fluctuation in turbulent flows with an immersed boundary method[END_REF], where, however, the approximated domain concept is only harnessed for the continuity equation. Although lacking strong theoretical foundations, the Ghost Fluid Method proved very effective for the simulation of practical problems, due to its ease of implementation and avoidance of the small cut cell instability. Approximate domains methods were also explored in [START_REF] Derek M Causon | Calculation of shallow water flows using a Cartesian cut cell approach[END_REF][START_REF] Derek M Causon | A Cartesian cut cell method for shallow water flows with moving boundaries[END_REF] in the context of the shallow water equations.

In this work, we propose a new approach that falls in the category of embedded finite element methods and leverages a surrogate/approximate boundary strategy. The key feature of the proposed approach is the idea of shifting the location where the boundary conditions are applied from the true to the surrogate boundary, and, most importantly, of appropriately modifying the shifted boundary conditions in order to avoid a reduction in the convergence rates of the overall formulation. In fact, we show that if the boundary conditions originally applied to the true domain are not appropriately modified, only firstorder convergence can be expected. The appropriate (modified) boundary conditions are then applied weakly, using, for example, a Nitsche strategy (other boundary enforcement strategies are also possible, e.g., Lagrange multipliers, etc.). This process yields a method which is simple, robust, accurate and efficient.

We apply our method to the shallow water equations, as a prototypical hyperbolic system. We also consider the linearized limit of the shallow water equations, constituted by the equations of acoustics [START_REF] Scovazzi | Weak boundary conditions for wave propagation problems in confined domains: Formulation and implementation using a variational multiscale method[END_REF][START_REF] Song | A nitsche method for wave propagation problems in time domain[END_REF][START_REF] Scovazzi | A velocity/stress mixed stabilized nodal finite element for elastodynamics: Analysis and computations with strongly and weakly enforced boundary conditions[END_REF][START_REF] Kagawa | Finite-element simulation of nonlinear sound-wave propagation[END_REF].

We implement the proposed approach for boundary conditions in the context of a stabilized finite element formulation, which takes inspiration specifically from [START_REF] Scovazzi | Weak boundary conditions for wave propagation problems in confined domains: Formulation and implementation using a variational multiscale method[END_REF][START_REF] Song | A nitsche method for wave propagation problems in time domain[END_REF][START_REF] Scovazzi | A velocity/stress mixed stabilized nodal finite element for elastodynamics: Analysis and computations with strongly and weakly enforced boundary conditions[END_REF] and more broadly from [START_REF] Hauke | A symmetric formulation for computing transient shallow water flows[END_REF][START_REF] Hauke | A stabilized finite element method for the saint-venant equations with application to irrigation[END_REF][START_REF] Takase | Space-time SUPG formulation of the shallow-water equations[END_REF][START_REF] Kashiyama | Parallel finite element methods for large-scale computation of storm surges and tidal flows[END_REF][START_REF] Kashiyama | 3-step explicit finite-element computations of shallow-water flows on a massively-parallel computer[END_REF][START_REF] Kashiyama | Adaptive boundary-type finite-element method for wave diffraction in harbors[END_REF][START_REF] Kawahara | Boundary-type finite-element method for wave-propagation analysis[END_REF][START_REF] Kashiyama | A boundary-type finite-element model for watersurface wave problems[END_REF][START_REF] Kashiyama | Boundary-type finite-element method for surfacewave motion based on trigonometric function interpolation[END_REF][START_REF] Kashiyama | Selective lumping finite-element method for nearshore current[END_REF][START_REF] Kawahara | Two-step explicit finite-element method for tsunami wave-propagation analysis[END_REF][START_REF] Smith | Application of finite-element analysis to body wave-propagation problems[END_REF][START_REF] Ricchiuto | Stabilized residual distribution for shallow water simulations[END_REF]. The proposed embedded method is not limited in applicability, however, to the specific stabilized method considered here. For example, the extension of the proposed approach to discontinuous Galerkin methods or residual redistribution schemes is straightforward.

The benefit of an embedded method in the context of shallow water flow is evident when considering the complex morphology of the ocean coastlines. In real scenarios, this geometric complexity may induce considerable meshing costs, which we attempt to significantly reduce with the proposed method. Applications of interest that we consider here are the treatment of complex coastlines as reflective walls in the framework of large scale simulations and the simulation of fine scale urban floods. In both cases, one may actually treat the coast and the constructions as part of the topography, but this will introduce hard requirements on the numerical method and the mesh generation. In fact, the alternative option of adapting the mesh to these geometries often requires very large numbers of cells, even when capturing the fine details of the flow may not be necessary. In this context our method has the benefit of relaxing the requirements on both the numerical scheme and the mesh, while allowing a full second order accuracy. We assess the accuracy, stability and robustness of the proposed method in a battery of tests, which also involve comparisons with finite element simulations on grids that conform to geometric boundaries.

The shallow water equations

The shallow water equations, also known as the de Saint Venant equations, are a simplified version of the Navier-Stokes equations for free-surface flows, and allow for the propagation of nonlinear waves. They are obtained by averaging the Navier-Stokes equations along the depth direction, and can be written as

∂h ∂t + ∂hv 1 ∂x 1 + ∂hv 2 ∂x 2 = 0 , (1a)
∂ ∂t (hv 1 ) + ∂ ∂x 1 hv 2 1 + 1 2 gh 2 + ∂ ∂x 2 (hv 1 v 2 ) = S 1 , (1b)
∂ ∂t (hv 2 ) + ∂ ∂x 1 (hv 1 v 2 ) + ∂ ∂x 2 hv 2 2 + 1 2 gh 2 = S 2 . ( 1c
)
As illustrated in Figure 1, h is the height of the water column, z is the bathymetry of the water bed, free surface level η = h + z and the two-dimensional position and velocity vectors are expressed in Cartesian coordinates as x = {x 1 , x 2 } and v = {v 1 , v 2 }. ∂/∂t indicates derivation with respect to time, and

S = S 1 S 2 = gh S o1 -S f 1 S o2 -S f 2 , (2)
that is, S is the source vector containing the slope of the basin (e.g., river bed, ocean floor) in the ith direction, S oi = -∂z ∂xi and the friction

S f i = f 2 vi √ v 2 1 +v 2 2 h 4/3
is expressed through the Manning's roughness coefficient f . Additional forces can be added to the source term Z if needed, to account, for example, for wind stress or other phenomena. The system (1a)-(1c) expresses the conservation of mass and momentum for a column of water of height h, is hyperbolic in structure and amenable to an analysis by means of the theory of characteristics and Riemann invariants. In the one-dimensional case, there are three characteristic speeds of propagation of information in the system, v -c s , v, and v + c s , where c s = √ gh is the speed of gravity waves. Hence a signal can move with the velocity of the flow (in analogy to the entropy waves for the Euler equations of gas Inria dynamics), or can move with the speed associated with either one of two additional waves, one receding and the other advancing with respect to the flow velocity with the celerity of gravity waves, ±c s , respectively. Equations (1a)-(1c) can also be written in system vector form as

∂ t U + ∇ • G = Z , (3a)
where ∂ t (•) = ∂(•)/∂t, for simplicity, and

U = h hv , G = G v +G h , G v = U ⊗v , G h = (0 2×1 ) T 1/2gh 2 I 2×2
, and Z = 0 S .

(3b) Here U ⊗ v = U v T and ∇ denotes the spatial gradient operator. The divergence operator ∇•, when expressed in Cartesian coordinates, is understood to apply to the second index of the

3 × 2-matrix G, that is ∇ • G = ∂ xj G ij .
The system of equations (3a) is also often expressed as

∂ t U + ∂ xi F i = Z , (4)
where the arrays F i are the columns of the matrix G, and can be decomposed as

F i = F v i + F h i , F v i = v i U , F h i = 1 2 gh 2 0 δ i , (5a)
with

δ i = δ 1i δ 2i . (5b)
Using any set of variables X, it is possible to rewrite the vector form (4) in quasi-linear form as

A 0 ∂ t X + A i ∂ xi X = Z . (6)
Here A 0 = U ,X and A i = F i,X is the ith Euler Jacobian matrix. When the vector of variables X are the conservation variables U , the Euler Jacobians become

A 0 = I 3×3 , (7a)
A 1 =   0 1 0 -v 2 1 + gh 2v 1 0 -v 1 v 2 v 2 v 1   , (7b)
A 2 =   0 0 1 -v 1 v 2 v 2 v 1 -v 2 2 + gh 0 2v 2   . ( 7c
)
Another natural choice for the vector of variables X are primitive variables Y = {h, v T } T , then the corresponding Euler Jacobians are

Â0 =   1 0 0 v 1 h 0 v 2 0 h   , (8a)
Â1 =   v 1 h 0 v 2 1 + gh 2hv 1 0 v 1 v 2 hv 2 hv 1   , (8b)
Â2 =   v 2 0 h v 1 v 2 hv 2 hv 1 v 2 2 + gh 0 2hv 2   . ( 8c
)
Complete specification of the problem requires initial conditions on Ω, namely,

v(t = 0) = v 0 , (9a)
h(t = 0) = h 0 , (9b)
v • n = 0 , on Γ N , (10)
which are typically applied when analyzing flooding around buildings in urban areas, and also when considering a coarse representation of the flow around large-scale coastal areas that are not the main focus of the simulation. The latter case is typical when there is interest in simulating a specific geographic area and it is necessary to also simulate some background, larger areas in order to appropriately characterize the dynamics of the overall flow, but these secondary areas are not of interest in the simulation, and only a coarse approximation is sufficient. Note that this type of boundary conditions conforms to the framework of the general Riemann problem, since the velocity is tangent to the boundary, and therefore only one of the characteristic lines goes across the boundary. In this case then, only one boundary condition is necessary, out of the possible three conditions that can be enforced (i.e., the height of the water column, the normal and tangential components of the velocity).

Open-sea (fixed height) boundary conditions (similar to Dirichlet boundary conditions):

h = η D -z , on Γ D , (11)
where the water height h is given by the free surface level η D and the bathymetry z. These conditions are somewhat artificial, in the sense that a more appropriate boundary condition on the far field of an open domain should be given by non-reflective boundary conditions. These are in general considerably more complicated to implement, and can be constructed starting from the other boundary conditions mentioned here.

Subcritical river inflow boundary conditions (|v

• n| ≤ c s ): h v • n = m I;sub , v • τ = 0 , on Γ I;sub , (12)
which, as the name suggests, are typical of river inflows when the speed of propagation of gravity waves is faster than the normal velocity of the fluid at the inflow, that is |v • n| ≤ c s . Recalling that c s = √ g h, it is easy to see that this situation occurs when the speed of the river inflow is relatively low with respect to the height of the water column. In this case there is always one characteristic line that exists the domain in the upstream Inria direction, and only two Riemann invariants (i.e. conditions) need to be imposed. The boundary condition on the water flux q = h v • n in [START_REF] Düster | The finite cell method for three-dimensional problems of solid mechanics[END_REF] can be replaced by the alternative conditions

v • n = v I;sub , or h = η I;sub -z , on Γ I;sub . ( 13
)
Supercritical river inflow boundary conditions (|v

• n| > c s ): v = v I;sup , h = η I;sup -z , on Γ I;sup , (14)
which are imposed at river inflows when the speed of propagation of gravity waves is slower than the normal velocity of the fluid at the inflow, that is |v • n| > c s . Opposite to the case of subcritical river inflow, this situation occurs for relatively high speed inflows, or relatively low water column heights. In this case no characteristic line can leave the domain upstream, and all three Riemann invariants of the system must be imposed, which means that we need three scalar conditions on the normal and tangential velocity components and the water height.

Subcritical river outflow boundary conditions (|v

• n| ≤ c s ): h v • n = m O;sub , on Γ O;sub . (15)
This condition imposes only one constraint (the mass flow exiting the domain) on the three unknowns at the outflow. Alternatively, other less common boundary conditions can be imposed:

v • n = v O;sub , or h = η O;sub -z , on Γ O;sub . (16)
In this situation there are two characteristics that exit the outflow boundary and only one condition needs to be imposed, on the normal velocity, or the tangential velocity, or the water height.

Supercritical river outflow boundary conditions:

In this case all characteristics exit the outflow, and no condition need to be imposed. The boundary Γ = ∂Ω is partitioned as

Γ D ∪ Γ N ∪ Γ I;sub ∪ Γ O;sub ∪ Γ I;sup ∪ Γ O;

Linearized shallow waters equations: Non-dissipative acoustic wave propagation

In the limit of small disturbances ṽ and h to the hydrostatic state of equilibrium of a fluid, denoted by v = 0 and h, we have

v 1 = ṽ1 1 , (17a)
v 2 = ṽ2 1 , (17b)
h = h + h , h = const. , h 1. (17c)
Substituting [START_REF] Höllig | Finite element methods with B-splines[END_REF] in the shallow water equations (1a)-(1c) and neglecting higher-order terms, we obtain

∂ t h + ∇ • ( h ṽ) = 0 , (18a)
∂ t ( h ṽ) + g h ∇ h = S . (18b)
Removing, for the sake of simplicity, the tilde above v, defining p = g h, dividing the first and second equation by h, denoting by c s = hg the speed of waves, χ = c -2

s
and by b = h-1 S, we have that [START_REF] Wang | A computational framework based on an embedded boundary method for nonlinear multi-phase fluid-structure interactions[END_REF] yield:

χ∂ t p + ∇ • v = 0 , (19a)
∂ t v + ∇p = b , (19b)
which is the mixed system of equations governing the propagation of acoustic waves in a homogenous medium. Typical initial conditions are

v(t = 0) = v 0 , (20a)
p(t = 0) = p 0 , (20b)
and these complement appropriate pressure boundary conditions on the Dirichlet boundary Γ D and normal velocity boundary conditions on the Neumann boundary Γ N . Namely:

v • n| Γ N = v N (x, t) , (21a)
p| Γ D = p D (x, t) . (21b)
Also ( 19) can be written in vector form. Denoting by Y = {p, v T } T the solution state vector, we have

∂ t U + ∂ xi F i = Z , in Ω , (22)
where

U = χp v , Z = 0 b , (23)
and the flux vector F i can be decomposed into two parts, that is, F v i associated with the velocity component of the solution and F p i associated with the pressure component of the solution:

F i = F v i + F p i , (24a)
F v i = v i 1 0 n d ×1 , (24b)
F p i = p                0 δ 1i . . . δ n d i                . ( 24c
)

Inria

Note the similarities in structure between ( 24) and [START_REF] Boffi | A finite element approach for the immersed boundary method[END_REF]. In order to recover the formalism of the quasi-linear shallow water equations found in equation ( 6), we can write

A 0 ∂ t Y + A i ∂ xi Y = Z , (25)
where

A 0 = χ 0 1×n d 0 n d ×1 I n d ×n d , (26a)
A i = F i,Y =         0 δ 1i • • • δ n d i δ 1i . . . 0 n d ×n d δ n d i         . ( 26b
)
Note that the ith Euler Jacobian matrix

A i ∈ R (n d +1)×(n d +1
) are constant and symmetric for i = 0 , . . . , n d .

3 Preliminaries: General notation, the true domain, the surrogate domain and maps e=1 Ω e (where n el is the total number of elements). We denote by h e = h e (Ω e ) the diameter of element e and h = max Ωe∈T h h e . Denoting by ω ⊂ Ω a portion of Ω (e.g., an element domain Ω e ), and by γ a portion of Γ, we define with

(v, w) ω = ω v w and (v , w) ω = ω v • w (27a)
the L 2 (ω)-and (L 2 (ω)) n d -inner products on the interior of ω and with

v, w γ = γ v w and v, w γ = γ v • w (27b)
boundary functionals on γ. In particular, the space L 2 (E ) refers to the set of functions whose traces are square integrable on the interior or exterior boundary set E . Let v 2 = (v, v) ω , and let W k j (Ω) be the Sobolev space with norm Consider now an embedded discretization, in which the computational grid does not conform to the boundary, but is actually overlapped with it (see Fig. 2). The computational grid intersects the true boundary Γ of the domain Ω. We can then introduce a surrogate boundary Γ composed of the edges/faces of the mesh that are the closest to the true boundary Γ. Γ can be constructed, for example, by computing the intersections of the grid and the true boundary Γ and using closest-point projection algorithms to detect the closet face/edge of Γ to Γ. Other choices are of course possible, as long as the overall topology of Γ and Γ are close to each other, that is if Γ has a certain number of holes, the same number of holes needs to also be present also in Γ. To a certain extent, the construction of Γ from Γ is a computational geometry problem, and several techniques can be borrowed from this field. The surrogate boundary Γ encloses the surrogate domain Ω. In particular, ñ indicates the unit outward-pointing normal to the surrogate boundary Γ, to be distinguished from the outward-pointing normal n of Γ. We now define the map

v W k j (Ω) =   α≤k D α x v j L j (Ω)   1/j , 1 ≤ j < ∞ . ( 28
) For j = 2, let H k (Ω) = W k 2 (Ω) and set v W k j = v k .
M : Γ → Γ , (29a) x → x , (29b)
which maps a point x ∈ Γ on the surrogate boundary to a point x ∈ Γ on the true boundary. For example, the map M can also be built by means of the closest-point projection of points in Γ onto Γ, as shown in Figure 2(c).

Remark 1. Note that the closest-point projection, in spite of the segmented/faceted nature of the surrogate boundary Γ is actually a smooth map from points in Γ to points in Γ, and this is the primary reason why we propose it as the back-bone of the map M .

In particular, it will become very important to characterize the map M though a distance vector function

d M (x) = x -x = [ M -I ](x) . (30)
In what follows, for the sake of simplicity, the subscript in the definition of d M will be omitted and we will simply write "d." If the closest-point projection is used, the vector d is aligned with n. This choice is made throughout the rest of this article, and is stated as Inria Assumption 1. The distance vector is defined as d = d n, where the normal n to the true boundary and the normal ñ to the surrogate boundary satisfy

n • ñ ≥ 0 , (31)
The condition n • ñ ≥ 0 means that we require that ñ lies on the half-plane identified by the normal n, a situation that is always verified in practice. Through the map M , it is possible to define the extension ψ on Γ of a function ψ initially defined only on Γ, as

ψ(x) ≡ ψ(M (x)) . ( 32
)
For example, the unit normal vector n and unit tangential vectors τ i (1 < i < n d -1) to the boundary Γ can be extended to the boundary Γ as follows:

n(x) ≡ n(M (x)) , (33a)
τ i (x) ≡ τ i (M (x)) . (33b)
In what follows, with the purpose of simplifying the notation, we will omit the bar from the expressions of n and τ i , whenever this does not cause confusion. Therefore if, in the following, we write n(x) we actually mean n(x), and similarly for τ i (x) and τ i (x). We can also introduce the derivative of a function ψ along the directions n and τ i at a point x ∈ Γ as

ψ , n(x) = ∇ψ(x) • n(x) = ∇ψ(x) • n(M (x)) , (34a)
ψ , τi (x) = ∇ψ(x) • τ i (x) = ∇ψ(x) • τ i (M (x)) . ( 34b
)
Observe also that the following Taylor expansion formula centered at x ∈ Γ holds for a generic field u at x = M (x) ∈ Γ:

u(x) = u(x) + ∇u(x) • (x -x) + O x -x 2 = u(x) + ∇u(x) • (M (x) -x) + O M (x) -x 2 = u(x) + ∇u(x) • d(x) + O d(x) 2 . ( 35
)
The last expression in the chain of equalities can be used to develop a new strategy for the imposition of boundary conditions in the context of embedded methods. This approach is intrinsically only second-order accurate, unless additional terms in the Taylor expansion are included.

Remark 2. Note that the closest-point projection, in spite of the segmented/faceted nature of the surrogate boundary Γ is actually a smooth map from point in Γ to points in Γ, and this is the primary reason why we propose it as the back-bone of the map M .

4 A review of the weak enforcement of boundary conditions

Time-domain linear acoustics

Let us first consider the case of weak enforcement of boundary conditions for conformal (i.e., not embedded) mesh computations of linear acoustic waves. This problem retains all the elements of the hyperbolic system of the shallow water equations, without the inherent complexity of nonlinearities. The discussion that follows is a summary of the approach proposed in [START_REF] Song | A nitsche method for wave propagation problems in time domain[END_REF] for time domain acoustics, and its extension to the context of shallow water equations.

For the sake of clarity and without loss of generality, here we consider only the case of a semi-discrete formulation, in which the equations are discretized in space but not in time. We also consider solution and test spaces that are slightly more regular than the roughest possible case, that is, the following function spaces are chosen for the velocity and pressure

S p = V p =H 1 (Ω) , (36a)
S v = V v =H div (Ω) . ( 36b
)
This last assumption is not a limitation, and the discussion can be extended to the most general case with relative ease. For example, the entire discussion can be adapted with minor changes to the case of discontinuous Galerkin approximation spaces. We consider a conforming finite element discretization, in which the discrete pressure and velocity spaces are subsets of the infinite dimensional ones above. For example, equal-order piecewise polynomial interpolation spaces can be used. Hence we write

S h p ⊂ S p , V h p ⊂ V p , S v ⊂ S v and V v ⊂ V v .
Testing the strong form of the equation agaist the appropriate test functions and integrating by parts, we obtain:

Find p ∈ S h p and v ∈ S h v such that, for all φ ∈ V h p and ψ ∈ V h v , 0 =(φ h , χ ∂ t p h ) Ω -(∇φ , v h ) Ω + φ h , v•n Γ , (37a) 0
=(ψ h , ∂ t v h -b) Ω -(∇ • ψ h , p h ) Ω + ψ h •n , p Γ . ( 37b
)
The terms v•n and p are numerical boundary traces of the solution and, if appropriately defined, can be used to impose boundary conditions weakly. The formulation we propose involves the following choices:

v•n = v • n on Γ D , v N on Γ N , (38a)
p = p D on Γ D , p on Γ N . (38b)
Note that the normal velocity component is enforced directly on the boundary Γ N , while the pressure is free, and the opposite happens on the Dirichlet boundary Γ D . In what follows, we will omit, with the goal of a more readable notation, the superscript h on the components of the solution and the corresponding test functions, in spite of the fact that these fields are assumed fully discretized in space. The previous choice yields the final abstract formulation

0 = (φ , χ ∂ t p) Ω -(∇φ , v) Ω + φ , v N Γ N + φ , v • n Γ D , (39a)
0 = (ψ , ∂ t v -b) Ω -(∇ • ψ , p) Ω + ψ•n , p D Γ D + ψ•n , p Γ N , (39b)
which can be rewritten in system vector form as

0 = (W , A 0 ∂ t Y -Z) Ω -(∂ xi W , F i ) Ω + W , F i n i Γ , (40a)
Inria with W = {φ, ψ T } T . We can then introduce the generalized vector H of boundary conditions as follows:

H(Y ) = F i n i =              v N 1 0 n d ×1 + p 0 n , on Γ N , v • n 1 0 n d ×1 + p D 0 n , on Γ D , (40b)
which, recalling the identities

F p i n i = p 0 n , (41)
F v i n i = v • n 1 0 n d ×1 , (42)
can also be rewritten as

H = H v + H p , where H v (Y ) =        v N 1 0 n d ×1
, on Γ N ,

F v i n i , on Γ D , H p (Y ) =          F p i n i , on Γ N , p D 0 n , on Γ D , (43a)
a form that highlights more clearly the type of boundary condition enforced, and the corresponding data. Then the boundary conditions (21a)-(21b) can be recast as Neumann condition on Γ N :

F i n i = H ⇔ v • n = v N , (44a)
Dirichlet condition on Γ D :

F i n i = H ⇔ pn = p D n . (44b)
Hence, it is possible to reformulate (39a)-(39b) using a stabilized formulation, analogous to the one pursued in [START_REF] Scovazzi | Weak boundary conditions for wave propagation problems in confined domains: Formulation and implementation using a variational multiscale method[END_REF][START_REF] Song | A nitsche method for wave propagation problems in time domain[END_REF][START_REF] Scovazzi | A velocity/stress mixed stabilized nodal finite element for elastodynamics: Analysis and computations with strongly and weakly enforced boundary conditions[END_REF] in the context of hyperbolic wave systems:

0 = (W , A 0 ∂ t Y ) Ω -(∂ xi W , A i Y ) Ω -(W , Z) Ω + W , H Γ + (L W , τ A -1 0 (L t Y -Z)) Ω , (45)
where τ = c τ ∆t 2 , and

L t = A 0 ∂ t + L , L = A i ∂ xi . ( 46
)

Nonlinear shallow water equations

In this case, boundary conditions can be enforced weakly with the following stabilized variational form, inspired again by the work in [START_REF] Scovazzi | Weak boundary conditions for wave propagation problems in confined domains: Formulation and implementation using a variational multiscale method[END_REF][START_REF] Song | A nitsche method for wave propagation problems in time domain[END_REF][START_REF] Scovazzi | A velocity/stress mixed stabilized nodal finite element for elastodynamics: Analysis and computations with strongly and weakly enforced boundary conditions[END_REF], and with close similarities to the variational forms developed in [START_REF] Hauke | A symmetric formulation for computing transient shallow water flows[END_REF][START_REF] Hauke | A stabilized finite element method for the saint-venant equations with application to irrigation[END_REF]:

0 = (W , ∂ t U -Z) Ω -(∂ xi W , F i ) Ω + W , H Γ + (L W , τ Â0 (L t Y -Z)) Ω + (∂ xi W , ν g ij Â0 ∂ xj Y ) Ω , (47)
where again τ = c τ ∆t 2 is the stabilization characteristic time scale. The term 47) is a discontinuity capturing operator, function of the artificial viscosity ν, which is defined as

(∂ xi W , ν g ij Â0 ∂ xj Y ) Ω in (
ν = max   0, (L t Y -Z) • Â-1 0 (L t Y -Z) | ∇ξ Y | Â0 1/2 - (L t Y -Z) • τ Â0 (L t Y -Z) | ∇ξ Y | Â0   , (48)
where

g ij = [ξ k ,i ξ k ,j ] -1 , ξ k , k = 1, 2, are the local element coordinates, with | ∇ξ Y | Â0 = Y , ξ0 • Â0 Y , ξ0 + g ij Y , i • Â0 Y , j (49)
and the approximate time derivative in the element parent domain given as

Y , ξ0 = Y (t n+1 ) -Y (t n ) 2 . (50)
See [START_REF] Hauke | A symmetric formulation for computing transient shallow water flows[END_REF][START_REF] Hauke | A stabilized finite element method for the saint-venant equations with application to irrigation[END_REF] for more details. Note also that H = H v + H h , and specifically,

H v (Y ) =                                                              v N h h v , on Γ N , v • n h h v , on Γ D , m I;sub 1 (v • n) n , on Γ I;sub , v I;sup • n h I;sup h I;sup v I;sup , on Γ I;sup , m O;sub 1 v , on Γ O;sub , v • n h h v , on Γ O;sup , H h (Y ) =                                                              1 2 gh 2 0 n , on Γ N , 1 2 gh D 2 0 n , on Γ D , 1
Inria H v (Y ) =                                                              v N h h v , on Γ N , v • n h h v , on Γ D , m I;sub 1 (v • n) n , on Γ I;sub , v I;sup • n h I;sup h I;sup v I;sup , on Γ I;sup , m O;sub 1 v , on Γ O;sub , v • n h h v , on Γ O;sup , H h (Y ) =                                                              1 2 gh 2 0 n , on Γ N , 1 2 g (η D -z) 2 0 n , on Γ D , 1 2 gh 2 0 n , on Γ I;sub , 1 2 g (η I;sup -z) 2 0 n , on Γ I;sup , 1 2 gh 2 0 n , on Γ O;sub , 1 2 gh 2 0 n , on Γ O;sup . (51)
Observe the following equivalences of boundary condition expressions:

H(Y ) v + H(Y ) h H(Y ) -F i n i = 0 ⇔                                                    (v N -v • n) h h v = 0 ⇔ v • n = v N , on Γ N , 1 2 g (η D -z) 2 -h 2 0 n = 0 ⇔ h = η D -z , on Γ D , h v • n = m I;sub and v • τ = 0 , on Γ I;sub , v = v I;sup and h = η I;sup -z , on Γ I;sup , (m O;sub -h v • n) 1 v = 0 ⇔ h v • n = m O;sub , on Γ O;sub , 0 = 0 ⇔ No conditions imposed , on Γ O;sup . ( 53
)

Discrete conservation properties

The variational formulation (47) naturally incorporates a conservation statement. Taking as test function W = 1 j , where 1 j is an array in which the j th entry is unity and all the others are zero, we obtain:

0 = (1 j , ∂ t U -Z) Ω + 1 j , H Γ , (54)
which is a conservation statement of the mass and the various component of the momentum equation. In fact mass or momentum would not change if the integral of the boundary fluxes H vanishes as well as the integral of the source vector Z. A similar statement can be obtained for the acoustic system (54) (recall that the pressure rate equation is effectively a statement of mass conservation).

C-property: Preservation of a lake at rest

When solving the shallow water equations, a main concern is the so-called well-balanced character of the discretization, which consists in the ability of the numerical method of exactly preserving the solution

v = 0 , (55a)
h + z = η 0 (constant) . (55b)
This is the state of a flat free surface in stationary equilibrium and is often encountered as a background or initial state in many applications. It is important to exactly preserve the solution (55), since a small spurious perturbation can be amplified by the bathymetry (or analogous mechanisms). This property is often referred as the C-property (conservation property) and was originally introduced in [START_REF] Bermudez | Upwind methods for hyperbolic conservation laws with source terms[END_REF].

The analysis of the C-property for the stabilized scheme proposed here is very similar to that of the residual based methods studied e.g. in [START_REF] Ricchiuto | An explicit residual based approach for shallow water flows[END_REF]. We will proceed first by omitting the stabilization and discontinuity capturing terms from the discussion. These can be treated in a second stage. Substituting (55) in [START_REF] Scovazzi | A velocity/stress mixed stabilized nodal finite element for elastodynamics: Analysis and computations with strongly and weakly enforced boundary conditions[END_REF] as the solution leads to specific conditions on the discretization to maintain the C-property. Specifically, no inflow or outflow are possible, the integrals involving the term F v i = v i U vanish identically, as well as the volume integrals of the friction source term S f i , as they involve variations of the flow speed. The only terms left are related to the hydrostatic pressure, and the bathymetry source of components gh S oi = -gh ∂z/∂x i . In what follows, we will indicate with Z o the simplified version of Z due to the previous simplifications. Now, if exact quadrature formulas are used for the evaluation of the terms (∂ xi W , F h i ) Ω and (W , Z o ) Ω , with respect to the piecewise linear interpolation spaces used in the approximation of h and z, then the continuity of the approximation and integration by parts over each element leads to

(W , ∂ t U ) Ω = (W , Z o -∂ xi F i ) Ω + W , F h i n i -H Γ . (56)
Let us focus our attention now on the term

(W , Z o -∂ xi F i ) Ω . Note first that the first component of Z o -∂ xi F i ,
corresponding to the mass conservation equation, vanish, and that the if ψ is a vector test function corresponding to the momentum equation, we have

(ψ j , ∂ xj (gh 2 /2) + gh ∂ xj z) Ω = (ψ j , gh∂ xj (h + z)) Ω = (ψ j , gh∂ xj η 0 ) Ω = 0 ,
where the first equality holds because of the assumed exactness of the integration with respect to the linear variation of h within each mesh element, and the second is true as long as z and h are in the same space and (55) holds. We are thus left with the boundary term W , F h i n i -H Γ , which will cancel exactly. Note now that in the case of a stabilized method with a discontinuity capturing operator, the additional terms are functions of the equation residuals, which vanish exactly in the case of the solution (55). Hence, we conclude that the stabilized variational form detailed in [START_REF] Scovazzi | A velocity/stress mixed stabilized nodal finite element for elastodynamics: Analysis and computations with strongly and weakly enforced boundary conditions[END_REF] satisfies the C-property. Remark 3. The above analysis is true as long as no dry areas are present in the domain. If there are elements in which h = 0 at some of the nodes, then the condition ∇(h + Inria z) = ∇η 0 = 0 may be violated. This is due to the fact that the value of the polynomial interpolating the nodal values of h + z may be different than η 0 at the dry node. To cure this issue, one must somehow modify the numerical approximation of the bathymetry to compensate for this unphysical effect. This is particularly necessary in the case in which the dry node is above the wet level in the element. A simple technique for this purpose is suggested in [START_REF] Ricchiuto | Stabilized residual distribution for shallow water simulations[END_REF][START_REF] Ricchiuto | An explicit residual based approach for shallow water flows[END_REF] and consists in modifying the nodal values of the bathymetry as follows:

ẑj = H max if h j = 0 and z(x j ) > H max z(x j ) otherwise
where on each element

H max = max k∈element h k >0 (h k + z(x k ))
Replacing the nodal values of z by ẑ in the computation of the integrals of bathymetric source terms allows to restore the C-property in dry cells. Some authors [START_REF] Brufau | Zero mass error using unsteady wetting-drying conditions in shallow flows over dry irregular topography[END_REF][START_REF] Castro | Numerical treatment of wet/dry fronts in shallow flows with a roe scheme[END_REF] also suggest to couple this correction with a limiter on the mass flux in vicinity of dry areas.

The Shifted boundary method

While weak boundary conditions can be quite effective on conformal grids, their applicability in the case of embedded boundaries is more challenging, for the following reasons:

a. The presence of cut elements of small size can render the overall approach numerically unstable and/or produce poor condition numbers in the ensuing algebraic system.

b. Numerical integration on cut elements can be computationally expensive and/or difficult to implement.

Because both of these issues emanate from the mere existence of cut elements near the boundary, one idea could be to exclude them altogether from the simulation. Excluding these elements has the effect of moving the boundary Γ of the computational domain to the surrogate boundary Γ. Of course, if boundary conditions are naïvely applied on the surrogate boundary, then an O(h) error is introduced. We present next an embedded finite element formulation, which falls under the umbrella of the shifted boundary (SB) method [START_REF] Main | The shifted boundary method for embedded domain computations. Part I: Poisson and Stokes problems[END_REF][START_REF] Main | The shifted boundary method for embedded domain computations. Part II: Linear advection-diffusion and incompressible Navier-Stokes equations[END_REF]. The key idea behind the SB method is not to apply boundary conditions on the true boundary Γ, but, rather, to shift their location to a surrogate boundary Γ. The map M defined in [START_REF] Kim | Mesh type tradeoffs in 2d hydrodynamic modeling of flooding with a godunov-based flow solver[END_REF] and the Taylor formula [START_REF] Parvizian | Finite cell method[END_REF] are instrumental in imposing on Γ a shifted boundary condition that is a second-order accurate approximation to the exact boundary condition on Γ.

Acoustic waves

We start from the hyperbolic system of non-dissipative acoustics described in Section 2.1 since it is a simpler prototype for the nonlinear shallow water equations.

Dirichlet boundary conditions

The goal of this section is to develop a suitable boundary condition on the surrogate Dirichlet boundary ΓD that, up to second-order accuracy, is equivalent to imposing the original boundary condition on the true boundary Γ D . To achieve this, we perform the following Taylor expansion of the pressure field along the direction

d = M (x) -x, from x ∈ ΓD to x = M (x) ∈ Γ D , 0 ≈ p (x) + ∇p (x) • (x -x) -p D (x) + O x -x 2 = p (x) + ∇p (x) • (M (x) -x) -p D (M (x)) + O M (x) -x 2 = p (x) + ∇p (x) • d (x) -p D (M (x)) + O d (x) 2 . ( 57
)
The last expression in the chain of equalities can be used as the modified boundary condition on the surrogate boundary of Dirichlet type ΓD , which preserves the accuracy of the true bounadry condition up to second-order. Note that according to the notation that we introduced in Section 3 we could have written pD (x) = p D (M (x)). Then the shifted boundary condition on surrogate Dirichlet-type boundary ΓD reads

p| ΓD = pD (x) -∇p (x) • d (x) . (58)

Neumann boundary conditions

A similar strategy can be elaborated to derive the boundary condition on the Neumann surrogate boundary ΓN . The main difference with respect to the previous case of Dirichlet boundary conditions is that the normal ñ to the surrogate boundary and the normal n to the true boundary do not coincide. This situation can be resolved by decomposing the unit normal vector ñ at x as ñ = ( ñ•n) n + ( ñ•τ j ) τ j , where n is the normal to Γ N and τ j (j = 1, • • • , n d -1) are the vectors tangent to Γ N , respectively, and we also recall that n(x) = n(M (x)) and τ j (x) = τ j (M (x)) by [START_REF] Hashem M Mourad | A bubble-stabilized finite element method for dirichlet constraints on embedded interfaces[END_REF]. Then we can apply the Taylor expansion to the velocity appearing in the term v(x) • n(x), so that, again for the point x on the surrogate boundary ΓN ,

v(x) • ñ(x) = (v(x) • n) n + v(x) • τ j τ j • ñ ≈ ((v(x) -∇v(x) d) • n) n + (v(x) • τ j ) τ j • ñ = v N -n T ∇v d n • ñ + (v • τ j ) τ j • ñ , (59)
where v N = v N (M (x)). Equation (59) represents the surrogate Neumann boundary condition on ΓN , and requires for consistency the tangential term (v•τ j ) τ j • ñ, a byproduct of the decomposition of ñ in terms of n and τ j .

Variational formulation

Using the approximate (shifted) Dirichlet and Neumann conditions we can modify formulation (37a)-( 37b) as follows:

0 = (φ , χ ∂ t p) Ω -(∇φ , v) Ω + φ , v • ñ Γ , (60a)
0 = (ψ , ∂ t v -b) Ω -(∇ • ψ , p) Ω + ψ• ñ , p Γ , (60b)
where now

v • ñ = v • ñ , on ΓD , v N -n T (∇v) d n • ñ + (v • τ j ) τ j • ñ , on ΓN , (61a)
p = p D -∇p • d , on ΓD , p , on ΓN . ( 61b
)

Inria

Hence, upon substitution, we obtain the variational formulation

0 =(φ , χ ∂ t p) Ω -(∇φ , v) Ω + φ , ((v N -n T (∇v) d)n + (v • τ j )τ j ) • ñ ΓN + φ , v • ñ ΓD , (62a)
0 =(ψ , ∂ t v -b) Ω -(∇ • ψ , p) Ω + ψ• ñ , p D -∇p • d ΓD + ψ• ñ , p ΓN . (62b)
The previous equations can be compared with (39a)-(39b) to highlight differences with respect to the case of conformal grids. Integrating by parts in space leads to the Euler-Lagrange equations,

0 =(φ , χ∂ t p + ∇ • v) Ω -( ñ • n) φ , v + (∇v) d) • n -v N ΓN , (63a)
0 =(ψ , ∂ t v + ∇p -b) Ω -ψ • ñ , p -p D + ∇p • d ΓD (63b)
which show that the SB method enforces the partial differential equations on the interior of surrogate domain Ω and the shifted boundary conditions on the surrogate boundary Γ.

As in the case of conformal grid computations, also the variational SB formulation can be cast in vector form:

0 = (W , A 0 ∂ t Y -Z) Ω -(∂ xi W , A i Y ) Ω + W , F i ñi Γ , (64)
where

F i ñi =                  v N -n T (∇v) d n • ñ + (v • τ j ) τ j • ñ 1 0 n d ×1 + 0 p ñ , on ΓN , v • ñ 1 0 n d ×1 + (p D -∇p • d) 0 ñ , on ΓD . ( 65
)
Denoting H = F i ñi , with a little algebra, somewhat tedious but otherwise straightforward, we obtain

H = F i ñi =                n• ñ   v N    1 0 n d ×1    -n i A v i ∇Y d   + τ j • ñ F v i τ j i + F p i ñi , on ΓN , p D    0 ñ    -ñi A p i ∇Y d + F v i ñi , on ΓD , ( 66
) in which A v i = F v i ,Y and A p i = F p i ,Y , so that ∇F v i d = A v i ∇Y d and ∇F p i d = A p i ∇Y d, respectively.
Then the vector form of (62a)-(62b) is expressed as

0 = (W , A 0 ∂ t Y ) Ω -(∂ xi W , A i Y ) Ω -(W , Z) Ω + W , H Γ + (L W , τ A -1 0 (L t Y -Z)) Ω . ( 67
)

Shallow water equations

In the case of the shallow water equations, it is important to take into consideration the complications that nonlinearities involve. We seek to construct a generalized vector of boundary conditions H, which enforces the modified boundary conditions on the surrogate boundary Γ.

We will decompose this problem into a series of subproblems which will be instrumental in forming the appropriate H for each of the six types of boundary conditions considered here. Let us begin by boundary conditions which enforce the value h = η * -z of the water column height on some portion Γ * of the boundary Γ, where Γ * = Γ D or Γ * = Γ I;sup . Then, on Γ * we will have

1 2 g h 2 (x) ≈ 1 2 g h 2 (x) -gh(x)∇h(x) • d ≈ 1 2 g (η * (x) -z(x)) 2 -gh(x)∇h(x) • d , (68)
so that Hh can take the form

Hh * (x) = H h * (M (x)) -gh(x)∇h(x) • d 0 ñ(x) = 1 2 g h 2 * -gh(x)∇h(x) • d 0 ñ(x) . ( 69
)
Turning now our attention to the term Hv N on Γ N , its design is rather straightforward using (59) and other results already obtained in the acoustics case. A term that requires instead special care is the term Hv * , on the surrogate of Γ * = Γ O;sub , where we impose h v • n = m * . We start from (59) and derive, neglecting quadratic and higher-order terms:

h(x) v(x) • ñ(x) = h(x) ((v(x) • n) n + (v(x) • τ ) τ ) • ñ ≈ h(x) v(x) • n(x) -n T ∇v d n • ñ + h(x) (v • τ ) τ • ñ ≈ h(x) v(x) • n(x) -h(x) n T ∇v d n • ñ + h(x) (v • τ ) τ • ñ ≈ h(x) v(x) • n(x) -(∇h(x) • d) v(x) • n(x) -h(x) n T ∇v d n • ñ + h(x) (v • τ ) τ • ñ = m * -(∇h(x) • d) v(x) • n(x) -h(x) n T ∇v d n • ñ + h(x) (v • τ ) τ • ñ ≈ m * -(∇h(x) • d) (v(x) + (∇v)d) • n(x) -h(x) n T ∇v d n • ñ + h(x) (v • τ ) τ • ñ ≈ m * -(∇h(x) • d) v(x) • n(x) -h(x) n T ∇v d n • ñ + h(x) (v • τ ) τ • ñ , ( 70
)
where we removed the superscript j from τ j since the shallow water equations are inherently two-dimensional. Hence,

Hv * (x) = m * -(∇h • d) v • n -h n T ∇v d n • ñ + h (v • τ ) τ • ñ 1 v , ( 71
)
where as usual, n(x) = n(M (x)), and m * (x) = m * (M (x)).

In the case of a subcritical inlet Γ * = Γ I;sub , one additional boundary condition v • τ = 0 on the tangential component of velocity needs to be imposed on the surrogate boundary. We start from (70) but we consider a revised approximation for v

• τ = 0, in Inria the term h(x) v(x) • ñ(x) ≈ m * -(∇h(x) • d) v(x) • n(x) -h(x) n T ∇v d n • ñ + h(x) (v • τ ) τ • ñ ≈ m * -(∇h(x) • d) v(x) • n(x) -h(x) n T ∇v d n • ñ + h(x) (v(x) • τ (x) -(τ ) T ∇v d) τ • ñ ≈ m * -(∇h(x) • d) v(x) • n(x) -h(x) n T ∇v d n • ñ + h(x) (-(τ ) T ∇v d) τ • ñ = m * , (72)
h(x) v(x) • ñ(x) v(x) ≈ m * (v(x) -∇v d) ≈ m * ((v(x) • n) n + (v(x) • τ ) τ -∇v d) ≈ m * ((v(x) • n) n -∇v d) ≈ m * (((v(x) + ∇v d) • n) n -∇v d) . ( 73
)
Then we obtain the following revised term Hv * :

Hv * (x) = m * -(∇h • d) v • n -h n T ∇v d n • ñ + h (-(τ ) T ∇v d) τ • ñ 1 ((v(x) + ∇v d) • n) n -∇v d , (74)
where quadratic and higher-order terms can be neglected in the products. It remains to consider the case of Hv I;sup on Γ I;sup , which can simply be handled using Taylor expansions. Namely,

Hv I;sup (x) = ((η I;sup -z) -∇h • d) (v I;sup -(∇v)d) 1 v I;sup -(∇v)d . (75)
In conclusion, we obtain the following final variational statement

0 = (W , ∂ t U -Z) Ω -(∂ xi W , F i ) Ω + W , H Γ + (L W , τ Â-1 0 (L t Y -Z)) Ω + (∂ xi W , ν g ij Â0 ∂ xj Y ) Ω , (76a)
with H = Hv + Hh and

Hv (Y ) =                                                                v N -n T ∇v d n • ñ + (v • τ ) τ • ñ h h v , on Γ N , v • ñ h h v , on Γ D , mI;sub    1 ((v(x) + ∇v d) • n) n -∇v d    , on Γ I;sub , (h I;sup -∇h • d) (v I;sup -(∇v)d)    1 v I;sup -(∇v)d    , on Γ I;sup , mO;sub 1 v , on Γ O;sub , v • ñ h h v , on Γ O;sup (76b) RR n°9136
this is not the case, and ∇z = 0 on Γ D or Γ I;sup , then the SB method may not satisfy the C-property. For this reason, in general a corrected form of (68) should be used, namely:

1 2 g h 2 (x) = 1 2 g (η -z) 2 (x) ≈ 1 2 g (η -∇η(x) • d -z) 2 (x) = 1 2 g h 2 * (x) -gh(x)∇η(x) • d + O( d 2 ) , (79)
with the convention that h * (x) = η * (x) -z(x). This modification allows to recover the C-property and is essential when the gradient of the topography/bathymetry is non-zero near boundaries, as in many applications.

Time integration

The integration in time of the equations is performed using the same algorithms presented in [START_REF] Song | A nitsche method for wave propagation problems in time domain[END_REF][START_REF] Scovazzi | A velocity/stress mixed stabilized nodal finite element for elastodynamics: Analysis and computations with strongly and weakly enforced boundary conditions[END_REF], to which the reader can refer for more details, and the analysis of stability and convergence, for body-fitted conformal grids. These algorithms are space-time integrators that are implemented as explicit predictor/multi-corrector with mass lumping. Specifically, we use a second-order Petrov-Galerkin space-time method, which, if only one corrector pass is performed, exactly corresponds to a second-order Runge-Kutta integrator. Full details about its implementation for acoustic wave problems is found in [START_REF] Song | A nitsche method for wave propagation problems in time domain[END_REF].

The extension of this algorithm to the shallow water equations is straightforward.

Numerical results for the wave propagation problem

We present a number of numerical results, to confirm the proposed approach is a robust, stable, and accurate strategy for boundary condition enforcement on embedded boundaries. All computations were run with three corrector passes of the second-order Petrov-Galerkin space-time integrator (similar to RK2) proposed in [START_REF] Song | A nitsche method for wave propagation problems in time domain[END_REF] for a Courant-Lewy-Friedrichs (CFL) condition of 0.5. The density ρ and wave speed c s are both equal to the unit constant. We pick the stability parameter c τ = 0.3 and mass matrices are lumped for all testes. We start with a battery of tests to verify the accuracy of the proposed method. We consider a two-dimensional, radially symmetric, stationary acoustic wave propagation problem, where the computational domain is a disk of radius R = 2.5, as shown in Figure 3, with corresponding conformal and embedded grids. The exact solution is radially symmetric and given by

v r (r , t) = 1 ρ c s J 1 z 1 ,1 R r sin c s z 1 ,1 R t , v θ (r , t) = 0 , p (r , t) = J 0 z 1 ,1 R r cos c s z 1 ,1 R t ,
where v r is the radial component of the velocity, v θ is the tangential component of the velocity, and the scalar z 1 ,1 = 3.83170597020751 is the first root of the first-kind Bessel function J 1 . Zero radial velocity boundary conditions are weakly enforced on the circle of radius R, using the SB method. Specifically, a background rectangular domain [-5, 5] × [-3, 3] is meshed using a fully unstructured triangular mesh, and its size is progressively and hierarchically refined from 0.5 to 0.0312 to check convergence rates. We consider the evolution of the problem until T = 4, that is, the instant of time at which the wave has reflected once against the boundary and is about to return to its original configuration. Figure 4 shows the convergence of the l 1 ([0, T ] : (L 2 (Ω)))-norm of the pressure and velocity for this test. The rates of convergence are very close to the expected second order using lumped mass matrices. The slight degradation of the convergence rates may be a byproduct of a more pronounced accumulation of dispersion error over long time durations, a well-known issue in the case of mass lumping. For comparison, results computed on conformal (body-fitted) grids for the same problem are also presented in Figure 4, and are virtually indistinguishable from the results of the SB method.

Dirichlet boundary condition on an embedded disk

To test the SB method in an acoustic problem with Dirichlet boundary conditions, we consider the same domain geometry of the previous test, and we change the boundary conditions. Consequently, the exact solution has been altered to comply with a zero pressure boundary condition on the embedded disk or radius R,

v r (r , t) = 1 ρ c s J 1 z 0 ,1 R r sin c s z 0 ,1 R t , v θ (r , t) = 0 , p (r , t) = J 0 z 0 ,1 R r cos c s z 0 ,1 R t ,
where z 0 ,1 = 2.40482555769577 is the first root of the first-kind Bessel function J 0 . The solutions are computed until T = 4 with embedded and conformal grids. We give the error of pressure and velocity in the l 1 ([0, T ] : (L 2 (Ω))) norm in Figure 5. Comparing the SB method with the reference Nitsche method on conformal grids, we observe that, although the errors in l 1 -norm associated with the SB are slightly larger, they still converge with order 1.8 and 1.7. We argue that this is actually a small price to pay when weighted against the ease of implementation and simplicity of the method, if compared with strategies that require tedious integration along the true boundary of the computational domain.

Acoustic waves around circle/rectangle

This third test involves the interaction of radial acoustic wave with boundaries of more complex geometry. The problem domain consists of the rectangle [-5, 5] × [-3, 3], in which a circular hole is present, of radius 1.0 and center (3, 0), as well as a square hole of side 2.0, and center (-3, 0). The goal of this test is to validate the robustness of the proposed SB method in the presence of sharp corners, which may induce singularities in the velocity profile. The hump-like initial condition is given by the fields:

v = 0 , p = 1 + A 1 + 2 r 3 R 3 -3 r 2 R 2 , r ≤ 1 , 1 , otherwise ,
where A = 0.01 is the amplitude, R = 1.0 is an intrinsic radius and r = x 2 + y 2 is the magnitude of the radial vector. A zero normal velocity boundary condition is applied on all boundaries. The computation using the SB method was performed on a grid that filled the entire rectangular domain [-5, 5] × [-3, 3] and of element size approximately 0.0312. In this test case, it is easy to observe that the solution is free from spurious oscillations throughout the many reverberations of the wave against the various boundaries, as shown in the pressure plots of Figures 6, 7 and8.

The same problem is also solved on a conformal grid of almost identical mesh size that is attached to the geometrical shape of the problem domain. As shown in Figures 6,7 and 8, the differences in pressure fields at various times are visually negligible.

An overlay of the results near the rectangular hole from the SB and conformal Nitsche methods is shown in Figure 8, where the solution of the SB method is shown as a solid surface and the solution of the conformal Nitsche method as a black wireframe. The two results show fairly good agreement.

8 Numerical results for shallow water flows Also in the case of shallow water flows, we compare the performance of the proposed SB method against a reference Nitsche method on conformal grids of equivalent resolution. All computations were run with two corrector passes of the second-order Petrov-Galerkin space-time integrator proposed in [START_REF] Song | A nitsche method for wave propagation problems in time domain[END_REF] (very similar to RK2 explicit integrator) with a Courant-Lewy-Friedrichs (CFL) condition of 0.5. Moreover, we chose the stabilization parameter c τ = 0.3, and the primary variables Y = {h, v T } T as solution variables. We also note that the discontinuity capturing operator is active only for the test of Sections 8.3 and 8.4.

C-property check

The C-property of the proposed SB method, discussed in Section 5.4, was validated in three numerical tests involving a body of water initially at rest in the rectangular domain

[-1, 1] × [-1, 1].
In the first test, the bathymetry of the river bed is given by a constant slope in the x-dimension, namely z(x, y) = 1.0 -0.25 (x -x min ) .

(80)

In the second test, the river bed bathymetry has the shape of a hump, located at the center of the computational domain:

z(x, y) = I h 1 0.4 exp -10 x 2 + y 2 , ( 81
)
where I h 1 is the interpolation operator on the space of piecewise linear polynomials defined on the computational grid utilized in simulations. In these two tests, the boundaries are solid walls, and no friction is imposed. A third test is then considered, in which ∇z = 0 and, in particular, z(x, y) = 0.25 + 0.25

I h 1 [cos(π (x -x min ))] , (82)
and fixed height boundary conditions are applied to the left and right hand side boundaries. Specifically, the correction suggested in (79) is used to maintain the C-property for the SB method. The other two sides are considered as solid walls.

For the SB method case, the right wall of the rectangular domain is embedded. The bathymetry of the river bed contributes to nonzero source terms S oi . The water inside the domain is initially set at rest and the bathymetric bed is completely submerged. Ideally, the equilibrium state

h(x, y) + z(x, y) = 1.0 , ( 83
) v = 0 , (84)
should be maintained for all three bathymetries, as well as zero velocities everywhere. The simulations are carried out until T = 100s, using the SB method and the conformal Nitsche method, with a grid of size h = 0.06. The discontinuity capturing operator is switched off since there is no discontinuity forming. The L 2 -norm of height and velocity errors are presented in Table 1 for these three different bottom types. The water height and velocity are unperturbed and the zero velocity is reproduced within machine double precision in the L 2 -norm, for all tests and both numerical methods.

Inria

Table 1: C-property check. L 2 -norm of the error in height and velocity at time T = 100s, for the SB and the conformal Nitsche methods. The C-property is preserved within machine accuracy.

SB method

Conformal Nitsche method height velocity height velocity Constant slope bottom 1.60495e-12 6.2741e-13 6.29469e-14 6.90434e-13 Central hump bottom 4.04608e-12 7.78965e-13 9.81656e-14 9.18004e-13 Cosine shape bottom 6.20764e-14 7.17012e-13 5.56596e-14 7.38136e-13

Convergence tests for shallow water flows inside a rectangular domain

We present convergence tests for the SB method with three different types of boundary conditions: impermeable wall, subcritical river inflow and subcritical river outflow. The computation is confined inside the rectangular domain [0, 8] × [0, 5]. In all simulations, and only for the case of the SB method, the right side of the rectangle -the set {8} × [0, 5] -is embedded, while for all other sides the boundary conditions are enforced weakly with a standard conformal Nitsche approach. Both the SB and the reference conformal Nitsche method were tested on five unstructured triangular grids of size h = 0.625, 0.3125, 0.15625, 0.078125, 0.0390625, respectively. We consider the evolution of the problem until T = 3. The computational domain and the coarsest level of refinement of the conformal and embedded grids are shown in Figure 9.

Impermeable wall boundary conditions

In this case, zero normal velocity boundary conditions are applied to all boundaries of the rectangular domain. We consider the exact solution h = 2 + 0.1 sin (ωt) cos (λx) ,

v 1 = 0.1 sin (λx) , (85a)
v 2 = 0 , (85b)
obtained with the method of manufactured solutions adding to the equations the source terms

S 1 = 0.1 (cos (λx) (2λ + ω cos (ωt)) + 0.1λ cos (2λx) sin (ωt)) , (86a)
S 2 = 0.1 sin (λx) (λ (0.005 -2g + 0.015 cos (2λx)) sin (ωt) + 0.1 cos (λx) (4λ + ω cos (ωt) -gλ sin (ωt) sin (ωt))) ,

S 3 = 0 , (86b)
where we set, in particular, λ = π/L and ω = 1.

Figure 10 shows a comparison of the error convergence rates of the SB method and conformal Nitsche method. Both methods are nearly second order, but with differences in the L 2 -norm of the error that are more noticeable than in the acoustic case. These discrepancies are however within a reasonable range, and the accuracy sacrificed is a small price to pay for the simplicity of the proposed SB method.

River inflow boundary condition

In this case, zero normal velocity boundary conditions are applied to the top and bottom boundaries of the rectangular domain, and an outflow boundary condition on the mass flux h v •n = 0.2 is applied to the left boundary. An inflow boundary condition (h v •n = -0.2) is applied to the right wall, and implemented as embedded in the case of the SB method. The exact solution is given by h = 2 -0.1 cos (ωt) sin (λx) ,

v 1 = -0.1 cos (λx) , (87a)
v 2 = 0 , (87b)
obtained with the source terms:

S 1 = -0.1 (-0.1λ cos(2λx) cos(ωt) -sin(λx)(2λ + ω sin(ωt))) , (88a)
S 2 = -0.1 cos(λx) (λ (-0.005 + 2g + 0.015 cos(2λx)) cos(ωt) -0.1g λ cos(ω) 2 sin(λx) + 0.1 sin(λx) (4λ + ω sin(ωt)) , (88b)
S 3 = 0 , (88c)
where λ = π/(2 L) and ω = 1. Note that the inflow boundary condition is enforced using the two separate conditions

h v • n = 0.2 and v • τ = 0 , (89)
as suggested by the theory of Riemann invariants. As shown in Figure 11, the differences in errors between the conformal Nitsche method and the SB method are small, and both methods yield the expected second-order convergence rate. and is obtained with the source terms

S 1 = -0.1 (-0.1λ cos(2λx) cos(ωt) -sin(λx)(2λ + ω sin(ωt))) , (91a)
S 2 = -0.1 cos(λx) (λ (-0.005 + 2g + 0.015 cos(2λx)) cos(ωt) -0.1g λ cos(ω) 2 sin(λx) + 0.1 sin(λx) (4λ + ω sin(ωt)) ,

S 3 = 0 , (91b)
where λ = π/L and ω = 1. The following three options for the imposition of the outflow boundary conditions are tested: a.) Mass flux: Results with all three previous options are illustrated in Figure 12. There is a discrepancy in the magnitude of errors between the Nitsche and conformal SB method in the case of the fixed height boundary conditions, while in all other cases the matching between the SB and conformal Nitsche method is very tight. In conclusion, we can say that the convergence rates are approximately second order for all implementations of the subcritical outflow boundary conditions.

h v • n = 0.2 , b

Coastline test

This test involves the propagation of a surface gravity-driven wave pulse against a portion of a coastline near Savannah, (GA). The coastline is reconstructed from the digital elevation map (DEM), as presented in Figure 13. The geometry of this problem is consid-erably complicated and may require advanced meshing capability if a discretization using a conformal Nitsche method is pursued. This is problem makes a good candidate for the proposed SB method. The spatial sampling interval of the DEM is about 8.75 meters in x-direction and 10.27 meters in y-direction. A triangulated grid of size 10 meters is used to fill the computation domain of [35200, 39300] × [-58400, -54600], with sufficient resolution. A closeup of the coastline and the embedded grid is presented in Figure 14. The boundary conditions utilized are of the impermeable wall type for the coastline, and represent a simplification, since in general one should admit the possibility of flooding. However, when large areas of coastline that are not the primary focus of a simulation need to be modeled, this is often the boundary condition utilized. Flooding boundary/interface conditions will be considered in future work, due to their specific nature.

The initial solution utilized for this test is the same as the hump given in section 7.2, applied to the water height this time, with parameters Amp = 1, R = 200 and r = (x -37089)

2 + (y -56250) 2 . The magnitude of the velocity and height are presented in Figure 15 and 16 at various times, while elevation plots of the height are presented in Figure 17. No oscillations are visible in this series of plots, in spite of the fact that the complex coastline chosen contains many corners. This test confirms that the proposed approach is a viable and robust strategy in the context of complex shallow water flows.

Toce valley flash flood test: A simplified urban district layout

In [START_REF] Testa | Flash flood flow experiment in a simplified urban district[END_REF], the authors conducted a series of flash flood flow experiments using a scaled model of the Toce river valley (Italy). The main goal of this work was the study of flow patterns in an urban environment, modeled by means of sixteen cubic concrete blocks of side 15cm, arranged in a 4 × 4 grid. Two masonry walls were erected with the purpose of constraining the flow on either side of the urban area. The problem setup is shown in Figure 18. Ten probes measuring water heights are placed at various locations: in front of the first row of houses, in the middle of streets and in the wake of some of the buildings. The initial state of the model is dry, the pump discharge history data of "low flow" is used for the enforcement of upstream boundary condition, and supercritical outflow boundary conditions are imposed at the downstream end of the valley. The frictional coefficient is 0.0162 for the concrete bed. To resolve the possible numerical difficulties associated with the dry bed condition, for which the water level h is zero, a threshold value = 1.0e -4 is applied to h and entropy variables are used as suggested in [START_REF] Hauke | A symmetric formulation for computing transient shallow water flows[END_REF].

In the computations presented here, three modifications and/or assumptions with Inria respect to the experimental setup are adopted to guarantee the well-posedness of numerical simulations. As shown in Figure 19, in the case of SBM, the buildings are embedded and two meshes of different size are used. A conformal mesh with size comparable to the finest SB method grid is used for comparison.

1. A constant inflow boundary condition is utilized at the inlet, since no specific information is available in the experiment. Inria 2. The location of the inlet was shifted to match the position of probe 2, since the bathymetry at the inlet does not seem to be represented accurately in the experimental data provided in [START_REF] Testa | Flash flood flow experiment in a simplified urban district[END_REF]. Consequently, all the results from the numerical simulation reported here are offset by one second, which is the time taken by the water front to cover the distance between probe 1 and probe 2 in the experimental data. The initial flooding time is postponed to t 0 = 4.6s, since until that time, probe 2 shows particularly noisy discharge data.

3. The computational domain is extended to x = 7.5 with a smooth decrease in elevation, to ensure supercritical flow in the outlet area.

The free surface elevations near the urban district at various times are shown in Figure 20 for the SB method simulation. Figure 18 presents the water depths predicted by the SB method from 4.6 to 55 seconds, compared to the laboratory measurements and the Nitsche method results on conformal grids. The sampling time interval is 0.2 seconds. Overall, there is good agreement between the simulations of both the SB method and conformal Nitsche's method with the experiment, with the exception of probe 5, where the computational results match, but are markedly different from the experimental data. However, this fact is well documented in the literature of shallow water flow computations for this test [START_REF] Sanders | Integral formulation of shallow-water equations with anisotropic porosity for urban flood modeling[END_REF][START_REF] Kim | Mesh type tradeoffs in 2d hydrodynamic modeling of flooding with a godunov-based flow solver[END_REF][START_REF] Vincent | Two-dimensional shallow-water model with porosity for urban flood modelling[END_REF], and the discrepancies observed in our computations are very similar to the ones observed in many other numerical simulations of this problem reported in the literature.

Figure 1 :

 1 Figure 1: Definition of water height h and the bathymetry z.

  Consider the domain Ω, an open set in R n d with Lipschitz boundary Γ = ∂Ω, where n d is the number of space dimensions. Let T h be a shape-regular triangulation (in the sense of Ciarlet) constituted by a family of non-overlapping n d -partitions/elements Ω e of Ω (e.g., triangles/quadrilaterals for n d = 2 or tetrahedra/hexahedra for n d = 3) such that Ω = n el

( a )

 a The surrogate domain Ω and true domain Ω. (b) The surrogate boundary Γ and true boundary Γ. (c) The distance vector d, and the unit vectors normal and tangent to the true boundary Γ.

Figure 2 :

 2 Figure 2: The surrogate domain, the boundary, and the vector d.

Figure 3 :

 3 Figure 3: Linear acoustic wave convergence test. Computational domain geometry and grids at the coarsest level of refinement: conformal grid (left) and embedded grid (right). For the embedded grid, light blue indicates Ω, that is the active elements inside the computational domain, dark blue indicates the elements outside, and red and orange indicate the elements intersecting the circular boundary (in white).

Figure 4 :Figure 5 :

 45 Figure 4: Convergence of the l 1 ([0, T ] : (L 2 (Ω)))-norm for the acoustic pulse propagation problem on an embedded disk domain with zero normal velocity boundary conditions. On the left, pressure, and on the right, velocity. A comparison between SB method and the conformal Nitsche method is shown for the errors of each solution field.

  (a) t = 0.0: Conformal Nitsche method (left) and SB method (right). (b) t = 3.5: Conformal Nitsche method (left) and SB method (right).(c) t = 6.0: Conformal Nitsche method (left) and SB method (right).

Figure 6 :

 6 Figure 6: Comparison between pressure contours obtained with the conformal and SB methods, for the propagation of an acoustics wave in a complex two-dimensional domain at various times.

Inria

  (a) t = 0: Conformal Nitsche method (left) and SB method (right). (b) t = 3.5: Conformal Nitsche method (left) and SB method (right).

  (c) t = 6.0: Conformal Nitsche method (left) and SB method (right).

Figure 7 :

 7 Figure 7: Comparison between pressure elevation plots obtained with the conformal and SB methods, for the propagation of an acoustics wave in a complex two-dimensional domain at various times.

  (a) mesh size ∆x = 0.125. (b) mesh size ∆x = 0.0625. (c) mesh size ∆x = 0.0312.

Figure 8 :

 8 Figure 8: Zoom of the elevation plot of the pressure solution near the rectangular boundary for various meshes. The results from SB method are shown as a solid contour, and are overlaid with the results of the conformal Nitsche method, shown as a black wireframe. The results of the conformal Nitsche and SB methods are difficult to tell apart.

Inria

Figure 9 :

 9 Figure 9: Nonlinear shallow water convergence test. Computational domain geometry and grids at the coarsest level of refinement: conformal grid (left) and embedded grid (right). For the embedded grid, light blue indicates Ω, that is the active elements inside the computational domain, and red and orange indicate the (inactive) elements intersecting the circular boundary (in white).

Figure 10 :

 10 Figure 10: Convergence of the l 1 ([0, T ] : (L 2 (Ω)))-norm for the shallow water flows with zero normal velocity boundary condition on the embedded vertical right side of the rectangular domain. On the left, height, and on the right, velocity. Comparison between the SB and conformal Nitsche methods are presented.

Figure 11 :

 11 Figure 11: Convergence of the l 1 ([0, T ] : (L 2 (Ω)))-norm for the shallow water flows with an inflow boundary condition on the embedded vertical right side of the rectangular domain. On the left, height, and on the right, velocity. A comparison between the SB and the conformal Nitsche methods is shown.

Figure 12 :

 12 Figure 12: Convergence of the l 1 ([0, T ] : (L 2 (Ω)))-norm for the shallow water flows with outflow boundary condition on the embedded vertical right side of the rectangular domain. On the left, height, and on the right, velocity. Three different choices, namely fixed height, fixed velocity flux, and fixed mass flux, for outflow boundary conditions are tested and verified. A comparison between SB method and the conformal Nitsche method is shown.

Inria

8. 2 . 3 Figure 13 :

 2313 Figure 13: The coastline near Savannah (GA).

Figure 14 :

 14 Figure 14: A closeup of the coastline and the embedded grids utilized in the computation. White solid lines represent the zero elevation isoline. Light blue indicates the active elements inside the computational domain, dark blue indicates the elements outside, and red and orange indicate the elements intersected by the coastline.

Figure 17 :

 17 Figure 17: Height elevation plots for the propagation of a tsunami along the coastline near Savannah, (GA).

Inria

  (a) Modified valley geometry and building layout with locations of probes. (b) Height of water at the probe locations 3 through 10, as depicted in Figure 18(a).

Figure 18 :

 18 Figure 18: Toce valley [52] flash flood test. Comparison of water depth numerical predictions with experimental data for various probes. The horizontal axis represents time in seconds and the vertical axis represents water depth in centimeters.

Figure 19 :

 19 Figure 19: Toce valley flash flood test. Zoomed views of the embedded and conformal grids near the buildings in the urban area: (left) embedded grid (coarse), (middle) embedded grid (fine), (right) conformal grid. In the left and middle figures, light blue indicates the elements inside the computational domain, dark blue indicates the elements outside, and red and orange indicate the elements that intersect with the embedded boundaries.

(a) T = 6 (

 6 b) T = 7.5 (c) T = 9.8

Figure 20 :

 20 Figure 20: Toce valley [52] flash flood test. Free surface elevation plots for the SB method solution.

  sup , in which all intersections (if any) of the partitions Γ D , Γ N , Γ I;sub , Γ O;sub , Γ I;sup and Γ O;sup are empty.

RR n°9136

Publisher Inria Domaine de Voluceau -Rocquencourt BP 105 -78153 Le Chesnay Cedex inria.fr ISSN 0249-6399

Acknowledgments

Guglielmo Scovazzi and the other members of his research team at Duke University (Ting Song and Alex Main) would like to gratefully acknowledge the support of the U.S. Department of Energy, Office of Science, Advanced Scientific Computing Research under Early Career Research Program Grant SC0012169, and the Army Research Office (ARO) under Grant W911NF-13-1-0452. Mario Ricchiuto has been partially funded by the TANDEM contract, reference ANR-11-RSNR-0023-01 of the French Programme Investissements d'Avenir. The support of the Associated Team HAMSTER between Duke University and the CARDAMOM INRIA team is also acknowledged. The authors would also like to thank Dr. Rodrigo Pedreros (at BRGM Orleans, France) and Dr. Chris Kees and Dr. Matthew Farthing (both at the Coastal and Hydraulics Laboratory, US Army ERDC, in Viscksburg, Mississippi) for valuable discussions.

and

where

Discrete conservation properties

Also in the case of the SB method statement of conservation can be derived for the variational formulations (67) and (76), by using the same test function W j defined in Section 4.3, but this time supported over Ω only. Namely, for (76), we have,

and a similar statement holds for (67).

C-property

The analysis of the C-property for the SB method is very similar to the one for the base conformal method, but with one important difference, in that now the boundary condition values are approximated with Taylor expansions, and for this reason the boundary terms

analogous to (56), may not always simplify exactly. Comparing with (76b)-(76c), we have that perfect cancellations always occur on the boundaries Γ N , Γ I;sub , Γ O;sub , and Γ O;sup , which do not need further care. However, on the boundary ΓD and ΓI;sup , the shifting of boundary conditions produces terms of the type gh∇h • d.

If the domain is defined in a way that ∇z = 0 in correspondence of Γ D and Γ I;sup , then ∇h = ∇(h + z) = ∇η 0 = 0, and we recover, in the end, (W , A 0 ∂ t Y ) Ω = 0. If however Inria