THICK INTERFACE COUPLING TECHNIQUE FOR WEAKLY DISPERSIVE MODELS OF WAVES

. The primary focus of this work is the coupling of dispersive free-surface flow models through the utilization of a thick interface coupling technique. The initial step involves introducing a comprehensive framework applicable to various dispersive models, demonstrating that classical weakly dispersive models are encompassed within this framework. Next, a thick interface coupling technique, well-established in hyperbolic framework, is applied. This technique enables the formulation of unified models across different subdomains, each corresponding to a specific dispersive model. The unified model preserves the conservation of mechanical energy, provided it holds for each initial dispersive model. We propose a numerical scheme that preserve the projection structure at the discrete level and as a consequence is entropy-satisfying when the continuous model conserve the mechanical energy. We perform a deep numerical analysis of the waves reflected by the interface. Finally, we illustrate the usefulness of the method with two applications known to pose problems for dispersive models, namely the imposition of a time signal as a boundary condition or the imposition of a transparent boundary condition, and wave propagation over a discontinuous bathymetry.


Introduction
The present research is dedicated to the coupling of dispersive models of free surface flows.While dispersive models, such as the Peregrine equations [43] or Green-Naghdi equations [24], are essential for wave propagation at the free surface, they exhibit less robustness compared to shallow water equations in various scenarios.Specifically, in dry fronts, most numerical strategies for solving dispersive models encounter stability issues.Furthermore, weakly non-linear models (such as Boussinesq-like equations) are not well-posed in these regions.Instabilities also occurs in areas with steep variations in bathymetry (non-Lipschitz bathymetry), where dispersive models lack well-posed solutions [31], or with vertical roofs [12,22].Additionally, understanding the boundary conditions of dispersive models remains incomplete, although some research is paving the way [26,33,34,41].In practical applications, the imposition of flow at the boundary requires source terms on the layers around it [5,47].In all these configurations, the shallow water model proves to be more robust, and various strategies exist.Also for modeling breaking waves, one strategy involves local use of the shallow water model, capable of dissipating energy through shocks [28].The primary goal of this work is to couple dispersive models with shallow water equations to address the mentioned challenges.To achieve this, two tools need to be developed: a robust coupling method and an indicator where the models should be used.The indicator depends on the application (break criterion or a priori estimate of the error) and will be taken into account in future work.The current work focuses on the development of a robust coupling method between dispersive and hyperbolic models.Similarly, we will explore the coupling between weakly dispersive and fully non-linear dispersive models.The main advantage of this coupling is its computational efficiency compared to simulations using the fully nonlinear model throughout the entire domain.Specifically, in regions where water depth is sufficiently large and remains relatively constant, weakly non-linear models yield similar results to fully non-linear models [27] and they are more cost-effective since the matrix used for resolving dispersive terms remains constant throughout the simulation.
This research is grounded in the reformulation of dispersive models as a hyperbolic model projected onto a linear subspace.Subsequently, we employ a coupling strategy commonly used for hyperbolic models [21,23].In Section 2.1, we describe a general framework of dispersive models, called projected hyperbolic models.To our knowledge, this framework has never been described in general but it has been used in [42] to propose a numerical scheme satisfying entropy and in [41] to analyse boundary conditions in a semi-discrete setting.This previous work can be followed to propose the same results for all projected hyperbolic models.It seems that many dispersive models [4,6,13,18,24,38,43,45] obtained by reduction of the water wave model can be written as a projected hyperbolic model, and in Section 2.2 we show some of them.In Section 2.3, we introduce a new model able to couple two projected hyperbolic models within the continuous framework.This model, parametrized by a spatial function , reverts to one of the previously mentioned models when  takes the value of 0 or 1.If the unified model is derived from two mechanical energy-preserving models and uses a single hyperbolic operator, it also preserves mechanical energy.For the models considered in this work, this is notably the case for the coupling between the Green-Naghdi and shallow water equations, and for the coupling between the Yamazaki and shallow water equations.Moving on to Section 3, we consider the discrete counterpart of the previous models.The strategy does not rely on domain decomposition with transmission conditions but involves discretizing the continuous unified model.This approach allows the definition of intermediate areas where  ∈ ]0, 1[, making the coupling smoother.This strategy is generally referred to as thick interface coupling.Applying the numerical strategy proposed in [42], Section 3.1 describes a numerical scheme that preserves the projection structure at the discrete level and as a consequence ensures the dissipation of discrete mechanical energy when verified at the continuous level.Finally, in Section 4, several numerical experiments have been conducted to study the behavior of the coupling.

General setting of the projected hyperbolic models
We focus on models, called later projected hyperbolic models, under the form where the state variables (, ) ∈ R   and  (, ) ∈ R   will be respectively called the potential and the velocity variable, where the variable  ∈ R  refers to space and the variable  ∈ R refers to time.We assume that the matrix (,  ) ∈    +  (R) is diagonalisable so that omitting the right hand side, the model is hyperbolic.Assume the hyperbolic operator exhibits an entropy structure, meaning there exists a pair of entropy-flux functions (, ) dependent on the state variables (,  ) such that for smooth enough solution, the relation (∇ , )   = (∇ , )  (2) holds true.Note that this condition is a common feature in numerous hyperbolic models, as discussed in references such as [7,20].Moreover, assume that the entropy (,  ) = P() + K (,  ) can be decomposed into a potential part dependent solely on the potential variables P() and a kinetic part defined through a weight function  : R   ↦ → R + , i.e.K (,  ) = () 2 | | 2 .Consequently, assuming that (, V) is solution of (1) without the right-hand side Ψ, each component of V are on  2 (() The right-hand side of (1), which will be referred to as the dispersive source term Ψ(, ), is not explicitly defined but is constructed to ensure that the velocity variable satisfies a relation L  ( ) = 0 for a given linear application L  : ( 2 (()))   ↦ → ( 2 (()))  with 0 ≤   ≤   .This linear application will be called the constraint and can be parametrized by the potential variable .In other words, the dispersive source term Ψ act so that the velocity variable  are on the linear subspace A  = Ker(L  ) called the set of admissible functions.At this stage, the model is not yet fully defined since the way to ensure that  ∈ A  is not unique.The dispersive source term Ψ is assumed to belong the dual space A ⊥  with respect the inner product Finally, we incorporate into the model the initial conditions (0, ) =  0 () and  (0, ) =  0 () ∈ A  0 .
Proposition 1.If there exists a pair of entropy-flux functions (, ) satisfying the condition (2), and the dispersive source term Ψ satisfies the orthogonality condition (3) with the inner product associated to the kinetic energy Proof.Multiplying (1) by ∇ , , we have Given that the entropy  depends on the velocity variable only through the kinetic energy, we have ∇   = () .Therefore, by integrating over the spatial variable, the last term vanishes as  ∈ A  and Ψ ∈ A ⊥  .
Note that Proposition 1 gives a global energy conservation law but the flux G (, , Ψ) of a local energy conservation law of the form is still unclear.It's crucial to note that this flux G isn't solely the hyperbolic flux  since the dispersive source term, or more precisely the associated Lagrange multiplier, can also contribute.Refer to Section 2.2 for illustrative instances.The orthogonality property (3) bears a striking resemblance to the duality observed between the pressure gradient and divergence-free velocity fields in incompressible fluid mechanics.More precisely, a significant number of models adhering to the projection structure outlined in Section 2.1 are reductions of the water wave equations [15].

Some examples of projected hyperbolic models
To illustrate the concept, this section provides examples of classical models that satisfy the projection structure defined in Section 2.1 after a reformulation without modification.More precisely, smooth enough solutions the projected hyperbolic model are solution of the classical dispersive model given as an example, and conversely.To explicitly define a model within the family of projected hyperbolic models as outlined in Section 2.1, it is necessary to specify the potential variable , the velocity variable  , the hyperbolic operator , the linear constraint L  (V), and the inner product ⟨V 1 , V 2 ⟩  .Many models can be formulated as projected hyperbolic models, starting with the incompressible Euler model.Both the KdV-BBM equations [4,29] and the Camassa-Holm equations [13] can be demonstrated to fit the projected hyperbolic model framework, as well as high-order dispersive models [37] and the layerwise Green-naghdi model [18].For greater physical relevance, while maintaining a reasonable level of complexity, this work will focus on models with the hyperbolic structure of the shallow water equations extended with some tracers, i.e.
where ℎ(, ) ∈ R + represents the water depth, () ∈ R denotes the bottom depth, and (, ) ∈ R  corresponds to the horizontal velocity (refer to Fig. 1).The bottom depth  can be treated as an unknown, although it satisfies a trivial equation.The tracers  1 ∈ R and  2 ∈ R can be interpreted as degrees of freedom of the vertical velocity.Additional details can be found in [18].Two hyperbolic operators will be considered, where the tracers are either advected with the flow  ad or are steady  st , namely The matrices are not strictly hyperbolic due to multiple eigenvalues involving  and 0. However, these multiple eigenvalues are associated with linearly degenerate fields.Additionally, it is worth noting that for both matrices  ad and  st , the model (1) can be expressed in a conservative form with a source term.Specifically, defining the projected hyperbolic model (1) with the matrices  xx (xx ∈ {ad, st}) can be written as where  st ( ) =  ad ( ) = (0, 0,  0 ∇ 1 , 0, 0)  and Ψ = ((0) 0≤<  , Ψ)  .It is well known that the couple entropy-flux defined by (,  ) = P() + K (,  ) with and the fluxes with and without advection of the tracers respectively reads )︂ ℎ satisfies the condition (2).Obviously, the change of variables (,  ) ↦ →  need to be reversible and in the current cases it reads

Hydrostatic model: the Shallow Water equations
Let us first consider for any V = ( 0 ,  1 ,  2 )  the constraint The constraint (8) has no impact on the first velocity variable  0 and nullifies the other components  1 and  2 .
It is evident that the projected hyperbolic model (1) with the hyperbolic operator  ad or  st defined by ( 5) and the constraint (8) recovers the classical shallow water equations for any inner product ⟨V 1 , V 2 ⟩  .To clarify, we determine the orthogonal set Φ = (  ) 0≤≤2 ∈ (A sw  ) ⊥ from U = (  ) 0≤≤2 ∈ A sw  by computing the inner product ⟨U, Φ⟩  , which vanishes for any  0 if and only if  0 = 0. Consequently, the first two equations of (1) with the hyperbolic operators  ad or  st and the constraint (8) correspond to the shallow water equations, i.e.
(  sw ) The constraint (8) is not the sole constraint leading to the shallow water equations (  sw ).Others may be interested in the constraint L ∅  (V) = ∅ instead, implying that the projection step has no effect also on the tracers  1 and  2 .

Weakly dispersive fully nonlinear models
Now we consider the constraint L fn  (V) = L [ℎ](V) and the inner product As demonstrated in [18], the projected hyperbolic model (1) with the hyperbolic operator  ad defined by ( 5) and the constraint (9) reproduces the classical Green-Naghdi equations [24].To clarify, we determine the orthogonal set Φ = ( This relation should be true for any  0 ∈  2 (()), hence we conclude that Applying this operator to the equation of the velocity variable  in (1), we get where    ∈    ×  (R) and    ∈    (R) are, respectively, the left lower block and the right lower block of the hyperbolic operator .Using the constraint (9) and the hyperbolic operator  ad defined in (5), the Green-Naghdi equations (in the form [30], (5.11)) are recovered, i.e.
where the dispersive operators are defined by with The operators  [h],  fn [h], and  ad [h] are defined to emphasize the contribution of tracer advection in the hyperbolic operator  ad .From (11) and using  st , we can write the projected hyperbolic model (1) with the hyperbolic operator  st defined by (5) and the constraint (9) as follows This model was first introduced in [48], and we refer to it as the Yamazaki model thereafter.
Corollary 1.The models (  fn ad ) and (  fn st ) conserve the energy  defined in (7).More precisely, the local energy balance (4) holds with G = G fn st for (  fn st ) and G = G fn ad for (  fn ad ), which reads )︂ and G fn ad (,  ) =  ad (,  ) + It is possible to interpret the last term of the flow as the vertical integral of the hydrodynamic pressure, see [18].
Proof.Since the inner product used for the projection uses the measure of the kinetic energy, i.e.
, we conclude using Proposition 1 the global conservation of energy.To find the energy flow, simply calculate the product ℎU • Ψ.We have The first term vanishes by definition of R  (Ψ) and the second term is the contribution of the dispersive term to the flux of energy.

Weakly dispersive weakly nonlinear models
While linear for the velocity variable, the subspace (︀ A fn

𝐻
)︀ ⊥ depends on the time dependent water depth ℎ.From a numerical point of view, this implies that the matrix of the linear system has to be recomputed at each iteration, significantly slowing down calculations.A classical simplification of fully nonlinear models involves replacing the water depth ℎ(, ) in the dispersive terms with the bottom depth () under the assumption that the perturbation of the free surface is small enough, i.e. ℎ =  + (), where  is the nonlinearity parameter, as seen in [30].In our framework, this simplification must be applied to the constraint (9).Through similar computations to those in Section 2.2.2, the projected hyperbolic model (1) with the hyperbolic operator  st and the constraint L wn  recovers the Peregrine equations [43].More precisely, the orthogonal set reads (A wn  ) (10).Applying this operator to the equation of the velocity variable  in (1), using the constraint L wn  and the hyperbolic operator  st , the Peregrine equations (in the form [30], (5.23)) are recovered, i.e.
with the dispersive operator defined in (12).Applying the operator R wn  to the equation of the velocity variable  in (1), using the constraint L wn  and the hyperbolic operator  ad , we can write the projected hyperbolic model (1) with hyperbolic operator  ad and the constraint L wn  as follow (  wn ad ) A similar model was first introduced in [32].One consequence of the simplification of the weakly nonlinear model is that the inner product hence the Proposition 1 does not hold.Indeed, it is well-known that the weakly dispersive models (  wn st ) and (  wn ad ) do not conserve the energy  as defined in (7).

Overview of the general strategy
Let us now explore the space-adaptive modeling of projected hyperbolic models, driven by the considerations outlined in Section 1.For simplicity, we focus on the coupling of two models: one represented by   x0 y0 in a subdomain Ω 0 ⊂ R  and the other by   x1 y1 in a different subdomain Ω 1 ⊂ R  with Ω 0 ∩ Ω 1 = 0. Our proposed strategy involves formulating a unified projected hyperbolic model and employing a standard resolution method, as detailed in Section 3.1.Let us introduce the color function (, ) ∈ [0, 1], which acts as a parameter for a model denoted as   x1  x0 y1  y0 such that   x1  x0 y1  y0 recover the model   x0 y0 when  = 0 and   x1 y1 when  = 1.In the point of view of hyperbolic equations, the color function serves as a parameter akin to bathymetry, defining extended variables as  = (ℎ, , ) and  = (,  1 ,  2 )  .As it was done in Section 2.2 for classical model of the literature, the unified projected hyperbolic model is defined by its hyperbolic operator and its constraint.Note that the coupling of  projected hyperbolic models can be achieved by considering the color function as a vector of [0, 1]  −1 , and proceed similarly.
The coupling of nonlinear hyperbolic models is extensively addressed in the literature.While specific cases may pose open questions, the coupling of many models, particularly conservative ones, has well-established solutions.For intricate cases, the series of articles [8][9][10][11] provides in-depth discussions.Regarding the relatively simple hyperbolic operators considered in this work, relevant references include [21,23].
A straightforward coupling strategy involves constructing a unified hyperbolic operator as a linear combination of the hyperbolic operators of the two models being coupled.The advantage of this approach is that if the two models satisfy an entropy relation with the same energy and different fluxes, the unified hyperbolic operators will satisfy an entropy relation by linearity.However, this strategy does not necessarily yield a conservative form (6), and the non-conservative products are typically unclear.Moreover, numerical schemes for approximating solutions of non-conservative hyperbolic models tend to be more complex.Several dispersive models of waves does not satisfy a conservation of energy, such as the Peregrine equations (  wn st ) or the Yamazaki equations (  fn st ), so unified models derived from them cannot satisfy it either.A second coupling strategy, which we will adopt later, was introduced in [23] with the aim of preserving the conservation of variables when both hyperbolic operators exhibit this property.This strategy is based on the continuity of fluxes.The coupled projected hyperbolic model is formulated in the conservative form (6) with the unified conserved variable given by  y1  y0 =  y1 +(1 − ) y0 and the unified flux  y1  y0 =  y1 +(1 − ) y0 , introducing an additional equation for the evolution for .The evolution of  is not crucial, as its value in practice is fixed by a source-term relaxation operator.
It is evident that  y1  y0 ,  y1  y0 tend to  y0 and  y0 where  = 0, and to  y1 and  y1 where  = 1.Then, we can compute the non-conservative form of the unified model (1) with the extended variables  and  .The hyperbolic operator reads .
We can easily verify that the matrix  ad  st is hyperbolic.More precisely, its eigenvalues are  ± √ ℎ (for the shallow water variables ℎ and ), 0 (twice for  and ) and ℎ ĥ (twice for  1 and  2 ).The change of variables  ↦ → (,  ) reads Note that the hyperbolic operator  ad  st does not satisfy the entropy relation (2).
In a similar manner, we introduce the constraint in the entire domain Note that where  = 0, and where  = 1.Finally the coupling of two projected hyperbolic models   x0 y0 and the model   x1 y1 is realized by the unified projected hyperbolic models   x1  x0 y1  y0 defined using the hyperbolic operator  y1  y0 , the constraint L x1  x0  and the inner product

Application to fully non-linear and weakly non-linear unified model
The coupling (  fn  wn ad  st ) between a fully non-linear model (  fn ad ) and a weakly non-linear model (  wn st ) is realized with the help of the constraint defined by (14), which in the present case becomes L fn  wn as defined in (9).We remember that ĥ = as defined by (10).
Applying this operator to the equation of the velocity variable  in (1), using the constraint L fn  wn  and the hyperbolic operator  ad  st , we obtain the unified model with the operators where  0 [h],  0 [h] and  2 [h] are defined in (13).
The coupling (  fn  wn st ) between the Yamazaki model (  fn st ) and the Peregrine model (  wn st ) can be realized in a similar way with the hyperbolic operator  st .We obtain the coupled model The coupling (  fn  wn ad ) between the Green-Naghdi model (  fn ad ) and the Boussinesq-type model (  wn ad ) is realized in a similar way using the hyperbolic operator  ad .We obtain the coupled model As with the classical weakly non-linear models, the coupled models (  fn  wn ad  st ) (  fn  wn st ) and (  fn  wn ad ) do not conserve the energy  in the form defined in (7), since the inner product of the projection is not defined thanks to the kinetic energy of the hyperbolic operator and thus the Proposition 1 cannot be applied.

Application to fully non-linear and hydrostatic models
The coupling (  fn  sw ad ) between the Green-Naghdi model (  fn ad ) and the hydrostatic shallow water model (  sw ) is realized using the constraint defined by (14), which in the present case becomes and the inner product (9).By performing the same calculations as in Section 2.2.2, we obtain that the orthogonal set (A fn  sw Finally, if we apply this operator to the equation of the velocity variable  in (1), use the constraint L fn  sw  and the hyperbolic operator  ad , we obtain the coupled model with the operators where  0 [h],  0 [h] and  2 [h] are defined in (13).
The coupling (  fn  sw st ) between a Yamazaki model (  fn st ) and the hydrostatic shallow water model (  sw ) is realized in a similar way using the hyperbolic operator  st .We obtain the coupled model ), which is Proof.The proof is similar to the proof of Corollary 1.

Application to weakly non-linear and hydrostatic models
The coupling (  wn  sw st ) between the Peregrine model (  wn st ) and the hydrostatic shallow water model (  sw ) is realized using the constraint defined by (14), which in the present case becomes (15) and the inner product . The orthogonal set (A wn  sw , defined by (16).Applying this operator to the equation of the velocity variable  in (1), using the constraint L wn  sw  and the hyperbolic operator  st , we obtain the coupled model with the operator   [] defined in (17).
The coupling (  wn  sw ad ) between a Boussinesq-type model (  wn ad ) and the hydrostatic shallow water model (  sw ) is realized in a similar way with the hyperbolic operator  ad .We obtain the coupled model with the operators defined in (17).
As with the classical weakly nonlinear models, the coupled models (  wn  sw st ) and (  wn  sw ad ) do not conserve the energy  as defined in (7), since the inner product of the projection is not defined thanks to the kinetic energy of the hyperbolic operator and therefore the Proposition 1 cannot be applied.

Discrete framework
The current section is devoted to the numerical resolution of a general projected hyperbolic model.This resolution is based on a splitting between the hyperbolic operator and the dispersive source term.This splitting is often used for the numerical resolution of dispersive models, see [1,14,16,17,32,36,40,44].In the current work, we follow the numerical strategy proposed in [42] because it is based on the projection relation (3) that we used for the coupling, also because it ensures entropy stability at the discrete level, i.e. a discrete counterpart of (4) and stability for a family of boundary conditions see [41].Note that the projection relation is not essential for the numerical resolution of the (1) model, not even for the unified models   x1  x0 y1  y0 .The strategies mentioned above can probably be used.

Overview of the numerical scheme
We consider a tessellation T of the horizontal domain Ω ⊂ R  consisting of Card(T) star-shaped control volumes, see Figure 1.We denote by  ∈ Ta control volume on the tessellation, by F  the set of its faces and by m  its surface area.In addition, for a face  , its length is denoted by m  and the neighbor of  by  is denoted by   , e.g. ∪   =  .The unit normal to  outwards to  is denoted by n    .We consider a finite volume scheme, so that the numerical unknowns    for  ∈ {ℎ, ,  1 ,  2 , Ψ} and the vectorial notations ,  and  are the averaged values in the control volume  at time   of the physical unknowns .
The first parts of the section are devoted to the description of a numerical scheme for a general projected hyperbolic model.Section 3.1.1describes a Godunov-type scheme for the hyperbolic operator, while Section 3.1.2describes the correction step that uses the projection relation (3).Note that the hyperbolic operator can be discretized with any kind of method, in particular explicit schemes are not required.This strategy was already used in [42] for the Green-Naghdi equations.We describe the projection scheme for the coupling of dispersive models in Section 3.2 and for the coupling with the shallow water equations in Section 3.3.

Hyperbolic step
The first step of the numerical scheme is the resolution of the hyperbolic operator.In order to use the most classical hyperbolic schemes, we assume in this paper that the hyperbolic operator can be written under the conservative form (6). We propose to use an explicit conservative Godunov-type scheme, i.e. we set where   =  −1 +    , F   is a numerical approximation of the flux F yy ( ) at time   and through the face  and S   is a numerical approximation of the source term  yy ( ) at time   and on the control volume .In practice, the real unknowns of this step are the unknowns of the shallow water plus the two tracers ℎ, ,  1 ,  2 , since the fluxes of the depth  and the color  vanish.
Several strategies for computing the numerical flux F   and the source term S   of the shallow water equations are described in the literature, see [7,20,35,46] m  if  = 2. (  ;   ) is an upper limit of the wave speed, depending on the numerical flux F used, see [7].In practice, we use an HLL scheme in this work to compute the water depth ℎ   and the horizontal velocity normal to the face    • n    and an upwind scheme (using the mass flux) for the horizontal velocity parallel to the face    • (n    ) ⊥ and the tracers  1 and  2 .The bathymetry source term is treated with the hydrostatic reconstruction [3].At the end of this step, we compute the potential variables  *  = H yy ( *  ) and the velocity variable  *  = U yy ( *  ) from the conservative variable  *  .For a variable  that is discretized on the tessellation, we find that  ⋆ = (  ) ∈T and the product term to term means  ⋆  ⋆ = (    ) ∈T .Assuming that the numerical scheme of the hyperbolic step is entropy satisfying, each component of

Projection step
Since the dispersive source term for the potential variable disappears, we first set  +1 ⋆ =  * ⋆ .As a continuous description, the numerical scheme then only requires the discretization L ⋆ (V ⋆ ) of the contraint L  ( ) and the discretization respectively.Then for any U ⋆ ∈ A ⋆ = Ker(L ⋆ ) we determine the dual discret space A ⊥ ⋆ so that for any Φ ⋆ ∈ A ⊥ ⋆ we have As in the continuous framework, the dual space can be defined thanks to a linear discrete function R ⋆ : − such that A ⊥ ⋆ = Ker(R ⋆ ).We will see in Section 3.2 how to compute in practice the discrete linear function R ⋆ .Finally, the projection step consists of decomposing the discrete functions U * ⋆ into the subspace A  +1 ⋆ and A ⊥  +1

⋆
. More precisely, we obtain the following linear system whose unknowns are Thank to the orthogonality relation (19), we have ( 2 (( ⋆ ))) and the discrete inner product is defined by where () is associated to the kinetic energy.Proof.Using the Pythagorean theorem, we have Since the potential variable is not affected by the second step of the scheme we have and we conclude since the hyperbolic scheme is assumed entropy-satisfying.
The linear system (20) can be interpreted as a mixed formulation, as called for the numerical scheme of incompressible flows, see also [25].In practice, the linear system (20) can be solved by a smaller linear system using the constraint L ⋆ and R ⋆ .For example, if we apply R  +1 ⋆ to the first equation of (20), we get ⋆ ).This system can be further reduced to a system of   −   unknowns using L  +1 ⋆ since the components of  +1 ⋆ are linked.Then the dispersive source term can be explicitly derived by . This strategy is called velocity correction when applied in the context of incompressible flows.Conversely, if we apply L  +1 ⋆ to the first equation of (20), we get    L ⋆ (Ψ +1 ⋆ ) = L ⋆ ( * ⋆ ).This system can be further reduced to a system of   unknowns using R  +1 .As it is the scheme is globally first order.However, it was shown in [42] that following the strategy proposed in [25], a high order scheme can be proposed with only one implicite projection step.

Application to fully non-linear and weakly non-linear unified model
In this section, the projection step for the coupling of weakly nonlinear and fully nonlinear dispersive models is described in detail.We propose to use the simple discretization of the constraints where the centered approximation is used to discretize the divergence and the gradient, i.e.
We will also use the discrete inner products Then we determine the dual subspace (A xx ⋆ ) ⊥ of A xx ⋆ = Ker(L xx ⋆, ).In the present work, we assume a periodic boundary for simplification.In the periodic domain, the following relation applies

Hence for any U
] with the discrete counterpart of the operator (10) ,⋆ to the first equations of the projection scheme (20) we obtain )︀ = 0 the linear system (20) finally resume to the elliptic equation of the horizontal velocity with The projection step can be summarized by solving the linear system (22) to compute

Application to weakly dispersive and hydrostatic unified model
In this section, the projection step for the coupling of shallow water equations with weakly dispersive models is described in detail.We propose the discretization of the constraints (21).The constraint of the shallow water equations is simply L sw ⋆, (V ⋆ ) = ( 1, ,  2, )  .The constraint of the coupling is therefore Since the shallow water equations do not depend on the inner product used for the projection, we use the inner product of the dispersive model
Then we determine the orthogonal set (A We conclude that to the first equations of the projection scheme (20), we obtain )︀ = 0, the linear system (20) finally goes into the elliptic equation (22) with the parameters The projection step can be summarized by solving the linear system (22) to compute Note that the numerical scheme coupling the Green-Naghdi equations with the shallow water equations   fn  sw ad is entropy-satisfying since it fulfills the assumptions of Proposition 2.

Some practical remarks on the numerical scheme
We would like to emphasize the following points for the numerical computation: (i) The way it is currently presented, the values for the unknowns  and  are computed in the hyperbolic step but remain unchanged.Obviously in practice, the computations related to these variables are not carried out.(ii) Where the shallow water equations are under consideration, i.e. when    = 0, the linear system (22) simplifies to  +1  =  *  .Consequently, the solution of the linear system ( 22) can be confined to the region where    ̸ = 0, resulting in a reduction of computation time.(iii) An advantage of employing the weakly nonlinear models (  wn st ) or (  wn ad ) is that the matrix of the linear system (22) remains independent of time.Inverting this matrix (or factoring it) is advantageous once at the simulation's initialization and subsequently apply this inversion to the right-hand side throughout the simulation.This is still the case for the coupling between weakly nonlinear models and hydrostatic model, i.e. (  wn  sw st ) or (  wn  sw ad ).This advantage is less apparent in cases involving coupling with fully nonlinear models, as seen in (  fn  wn ad  st ), (  fn  wn ad ), (  fn  wn st ).When the color function  remains constant in the part of the domain where the weakly nonlinear model is applicable, i.e.   = 0, there is no need to estimate the coefficients of the matrix of the linear system (22) at each time step.However, the entire system still needs to be solved at each iteration, akin to the fully nonlinear models.Despite this, it is still feasible to independently solve the weakly nonlinear model portion using an iterative block method, such as the block-Jacobi or block-Gauss-Seidel methods.This strategy proves to be advantageous, especially in the context of domain decomposition for large simulations, where the decomposition aligns with the boundaries of the subdomain where   = 0.

Numerical validation
This section is dedicated to illustrating the coupling strategy proposed in the current work through various numerical experiments.We consider a homogeneous 1D tessellation given by T = [1,  ] ∩ N, where  is the number of control volumes, and the spatial step is   = m  = mΩ  , with m Ω being the domain size.Throughout all simulations, the gravity acceleration  is set to 9.81.
In Section 4.1, we conduct a thorough parameter analysis of the method in a simple scenario, i.e. a traveling wave on flat bottom passing through the coupling interface.We illustrate then how to take advantage of the method with two typical cases, by imposing time-dependent boundary conditions Section 4.2 and for the propagation over a discontinuous bottom Section 4.3.) with some interface thicknesses  = 0,  = 5 and  = 10 with the initial condition (23) with  = 10 −2 ,  = 10 −1 and   = 10 −3 .

Coupling between dispersive and hydrostatic models
We consider the following test case: a solitary wave propagating over a flat bottom in a 1D domain with a length of m Ω = 100.The solitary wave (of the Green-Naghdi equations (  fn ad )) is initially centered at  0 = 40 and travels to the right, i.e.
where  is the shallowness parameter and  the nonlinearity parameter.The bottom depth  = m Ω since the bottom is flat.The color function  is defined by where  is the thickness of the interface and  controls the regularity at  = 60.If  = 0, the interface is sharp such that  = 1 for  < 60 and  = 0 for  > 60 for any .The boundary condition are define as wall at both sides, i.e. (, 0) = (, m Ω ) = 0.The discrete initial conditions are defined as follows: , where   is the spatial step.The additional velocity variables,  0 1, and  0 2, , are determined such that L  0 (︀  0 )︀ = 0. To enforce the discrete wall boundary conditions, a mirror flow approach is employed, i.e. we set , and   , +1 =   , .Let's discuss the results obtained through the coupling of the weakly dispersive model with the shallow water equations.In Figure 2  is slightly smaller than the waves of the uncoupled models, even though the shallow water equations have not yet developed shocks.This reduction in amplitude is attributed to the reflection of a portion of the wave at the interface.
In Figure 3, the difference between the water depth obtained with the coupled model (  fn  sw ad ) and the uncoupled Green-Naghdi equations (  fn ad ) is plotted at time  = 40 √ mΩ with   = 10 −3 ,  = 10 −2 ,  = 10 −1 , and various interface thickness values  = 0,  = 5, and  = 10 with  = 1.As anticipated, the amplitude of the reflected waves diminishes with increasing interface thickness.The reflected waves originate at the boundary of the dispersive model,  = 60, and subsequently propagate to the left.In Figure 4, the same simulations are conducted with  = 10 and varying values of  from 1 to 5. The results demonstrate that the solutions are smoother at  = 60, even if the reflected waves have the same amplitude.Optimal results appear to be achieved for  = 2, where  is regular at  = 60 but without to large variations in the interface.
In this section, we aim to provide a detailed characterization of the amplitude of the reflected waves.region of the domain where  = 1.This quantification is precisely defined by the norm with  out is the number of the output such that the solution at iteration  out is an approximation of the solution at time  = 0.4out √ mΩ .The water depth ℎ yy xx corresponds to the water depth obtained using the numerical scheme for the model   yy xx .For clarity, the reflected waves have been normalized by the amplitude of the initial condition m Ω .
From the left picture of Figure 5, it becomes apparent that the amplitude of the reflected waves converges as   goes to zero, aligning with the expectation that the numerical scheme for the coupled model (  fn  sw ad ) is entropy-satisfying.This observation leads us to conclude that the observed wave is not a numerical artifact but rather an inherent property of the coupled model (  fn  sw ad ), where the color function  undergoes variations.This behavior can be rationalized by the dynamic nature of the wave celerity in (  fn  sw ad ), which changes spatially with the color function .As the interface thickness  increases, the color function  becomes more regular, resulting in a reduction of the amplitude of the reflection waves.However, it's worth noting that the decay of the reflection waves is not particularly rapid, and in practical terms, it remains unclear if increasing the interface thickness is a feasible solution to mitigate this effect.Nevertheless, the amplitude of the reflected wave remains below 1% of the initial wave, even at stiff interfaces.
The center and right pictures of Figure 5 illustrate the dependency of the reflected waves on the physical parameters  and .It is evident that the normalized amplitude of the reflected waves predominantly decreases with both physical parameters, regardless of the interface thickness .Specifically, for sufficiently large interface thickness, it appears that the amplitude of the reflected waves can be expressed as According to [30], the shallow water model (  sw ) is an approximation of order  (︀  2 )︀ of the water waves model, while the Green-Naghdi model (  fn ad ) is an approximation of order  (︀  4 )︀ .From (24), we conclude that the coupling does not introduce an error greater than the modeling error corresponding to the error of the worst-case model.Furthermore, assuming that the shallow water model is used is nearly linear regime area, i.e.  ≈  4 3 , the coupling strategy should leads to result similar to the Green-Naghdi model.This can be use for instead for linear boundary condition forcing Section 4.2.On the other hand, it is possible to use the shallow water model in truly non-linear regions to model breaking waves [28], which leads to a significant difference with the dispersive model since it involves energy dissipation.Similar results are obtained for the models (  fn  sw st ), (  wn  sw st ) and (  wn  sw ad ) as shown respectively on Figure 5 second line, third line and fourth line.

Application to boundary conditions
In this section, we illustrate how to take advantage of the coupling strategy to impose a temporal signal as a boundary condition.To achieve this, we consider a 1D domain with m Ω = 100, initially at rest with ℎ 0 = 1 and  0 = 0, and study the propagation of a plane wave given by with  representing the frequency and  denoting the wave number.Due to the non-linearity, the waves (25) are not solutions of the models.Additionally, owing to the dispersion relation of the models, the frequency and wave number are interrelated.For all the models considered above, the dispersion relation is expressed as except for the shallow water equations (  sw ) where () = .To impose the signal at the boundaries, we employ different methods depending on the models: -For the shallow water equations (  sw ), we enforce the water depth at the left boundary of the domain, ℎ(, 0) = ̃︀ ℎ(, 0).At the discrete level, this is achieved by setting ℎ  0 = ̃︀ ℎ(  , 0), and the velocity is computed to preserve the outward Riemann invariant: )︀ see [20].To maintain the constraint L sw  ( ) at the boundaries, we set   ,0 = 0.For the right boundary, a transparent (free) boundary condition can be imposed by considering homogeneous Neumann boundary conditions for all the unknowns.At the discrete level, this involves setting ℎ   +1 = ℎ   ,    +1 =    , and   , +1 =   , .-For the dispersive models (  fn ad ), (  fn st ), (  wn st ), or (  wn ad ), we enforce the discharge and the hydrodynamic pressure at the left boundary of the domain as proposed in Section 5.2 of [41].This strategy involves fixing the discharge at the boundary, the hydrodynamic pressure, and the vertical velocities.In this document, the hydrodynamic pressure at the boundary is neglected, and the vertical velocity   1 is computed using the constraint L xx  ( ) of the model.For the right boundary, we use the transparent boundary condition as proposed in Section 5.1.2 of [41].
-One application of the coupled models (  fn  sw ad ), (  fn  sw st ), (  wn  sw st ) or (  wn  sw ad ) is to impose boundary conditions similar to those for the shallow water equations (  sw ) by setting  = 0 close to the boundaries.In practice, we set   = 0 only for the two control volumes closest to each bound, i.e.  ∈ {1, 2,  − 1,  }, and   = 1 in the rest of the domain.The thickness of the shallow water domain corresponds to the stencil of the linear system (22).It is noteworthy that certain boundary conditions of the dispersive models cannot be specified using this technique, especially those involving non-vanishing hydrodynamic pressure at the boundary.However, for many practical cases, specifying time series of water depths at the boundaries suffices, and this aligns with the current scenario.
The results of simulations have been plotted in Figure 6 with   = 10 −2 , the wave number  = 10 −1 , the frequency  computed with (26) (even for (  sw )), and  = 10 −2 .On the first line, we display the results at time 30, before the wave reaches the right boundary.It is observed that the results obtained with (  fn ad ) and (  fn  sw ad ) are nearly identical, validating the boundary condition on the right.The main advantage of imposing the boundary condition with the coupled model lies in its simplicity.On the second line, the results at time 100 are plotted, after several periods have left the domain.Due to a small wave reflected at the right boundary, the results of (  fn ad ) are shifted, and the mean water depth is not correct, as it is not fixed by the left boundary condition.A better transparent boundary condition for the Green-Naghdi equations could potentially correct this drawback (see [26]).With the coupled model (  fn  sw ad ), the reflection at the right boundary almost disappears, and the mean water depth is fixed at the left boundary condition regardless.Higher non linear waves can be simulated as well, as shown in Figure 7, where the simulations were plotted with the same parameters except for , which was set to 10 −1 .Similar results were observed with the other models (  fn  sw st ), (  wn  sw st ), or (  wn  sw ad ).

Application to discontinuous bathymetry
In this section, we illustrate how to take advantage of the coupling strategy to simulate the propagation of a wave over a discontinuous bottom.As indicated in [31], dispersive models are only well-posed with a bottom that is slightly more regular than Lipschitz, i.e.  1+,∞ with  > 0. In our framework, this limitation can be understood by the fact that with a non-Lipschitz bottom, the constraint L [h](V) defined by ( 9) does not define a linear subspace.However, this is not the case for the constraint L  (V) defined in (8).It is noteworthy that the shallow water model is not clearly derived over a discontinuous bottom [19]; it is not wellposed due to non-conservative product [7], and its numerical approximation can lead to different results using different schemes [2].Nevertheless, many numerical schemes for the shallow water model are stable and used for practical applications.The objective of this section is to propose a stable numerical strategy to perform numerical simulations of dispersive models over a discontinuous bottom.
Let us consider in a 1D domain with m Ω = 100 with wall boundary conditions and initially a solitary wave given by (23) with  0 = 20,  = 10 −2 and  = 10  .This ensures that the thickness of the domain where  = 0 is minimized, yet sufficient for the projection step not to encounter the discontinuity of the bottom.
In Figure 8, we show the results obtained with (  sw ) (red lines), (  fn ad ) (green lines) and (  fn  sw ad ) (blue lines) with a small discontinuity  = 10 −1 .After time  = 4, the results of (  fn  sw ad ) become significantly not relevant.A discontinuity of the velocity is observed at the discontinuity of the bottom, positive at the left of the discontinuity and negative at the right.This leads to the formation of a hole in the water that grown in time until time  = 8 where the water depth vanishes.Note that the simulation of (  fn ad ) with the scheme described in Section 3 can simulate the propagation of waves over continuous bottom without any problem as shown in [42].The results obtained with (  sw ) and (  fn  sw ad ) do not seem inconsistent over the discontinuity.For larger discontinuities of the bottom, we observed that the wave above the discontinuity of the bottom become discontinuous for both (  sw ) and (  fn  sw ad ), see Figure 9 for  = 9 • 10 −1 .Although probably unphysical, this discontinuity does not interfere with the simulation and the transmitted and reflected waves seem consistent with what we expect.Obviously, the relevance of the numerical results required more comparisons with finer model like 3D Navier-Stokes or observations.
In Figure 8, we present the results obtained with (  sw ) (red lines), (  fn ad ) (green lines), and (  fn  sw ad ) (blue lines) over a small discontinuity  = 10 −1 .Beyond time  = 4, the results of the uncouple model (  fn ad ) become significantly less relevant.A discontinuity in velocity is observed at the bottom's discontinuity, positive to the left and negative to the right.This leads to the formation of a hole in the water, growing over time until it vanishes at  = 8 where the water depth becomes zero.It is noteworthy that the simulation of (  fn ad ) with the scheme described in Section 3 can successfully simulate wave propagation over continuous bottoms, as shown  in [42].The results obtained with (  sw ) and (  fn  sw ad ) do not appear inconsistent over the discontinuity.For larger bottom discontinuity, we observed that the wave becomes discontinuous for both (  sw ) and (  fn  sw ad ), as seen in Figure 9 for  = 9 • 10 −1 .Although potentially unphysical, this discontinuity does not disrupt the simulation, and the transmitted and reflected waves seem consistent with expectations.It is essential to acknowledge that the relevance of numerical results requires further comparison with more refined models, such as 3D Navier-Stokes, or validation through experimental data.

Conclusion
In this work, we provide a general overview of the projected hyperbolic models, encompassing several existing dispersive models from the literature.This generalization allows us to extend previous results, including the study of boundary conditions [41] and the development of entropy-satisfying numerical schemes [42], to all models respecting this structure.Additionally, the general framework facilitates the formulation of continuous and discrete couplings of dispersive models while preserving essential properties like energy conservation.We illustrate the usefulness of the method through two applications known to present challenges for dispersive models, namely the imposition of a time signal as a boundary condition, and wave propagation over a discontinuous bottom.
The focus of this work was primarily on classical weakly dispersive models, either coupled together or with a shallow water model.A natural extension of this work would involve high-order dispersive models [37,39].Specifically, coupling weakly dispersive models with high-order dispersive models could be explored to optimize computation time by restricting the use of high-order models to regions where they are most needed.Another intriguing extension could involve investigating time-dependent functions  such as those used for modeling breaking waves [28].A time-dependent function  could also be appropriately introduced to achieve adaptive coupling between the models through an a priori estimator.
then the smooth enough solutions of the projected hyperbolic model (1) satisfy to the following energy conservation law   ∫︁ R  (,  ) d = 0.

Figure 1 .
Figure 1.Illustration of the notations: Left: interpretation of the unknowns in the vertical plan.Right: Finit volume discretization in the horizontal plan.

2 .
The models (  fn  sw ad ) and (  fn  sw st ) conserve the energy  defined in (7).More precisely, the local energy balance (4) holds with G = G fn  sw ad applies to (  fn  sw ad ) and G = G fn  sw st applies to (  fn  sw st

⋆
since the components of Ψ +1 ⋆ are linked.Then  +1 ⋆ =  * ⋆ −    Ψ +1 ⋆.This strategy is called pressure correction when applied in the context of incompressible flows.At the end of this step, we compute the conservative variables  +1