F Chazel

M Benoit

A Ern

S Piperno

A double-layer Boussinesq-type model for highly nonlinear and dispersive waves

Keywords: Water waves, Boussinesq-type models, multi-layer models, velocity potential formulation, nonlinear dispersive waves, deep water, Padé approximants, velocity profiles, linear shoaling

We derive and analyze in the framework of the mild-slope approximation a new double-layer Boussinesq-type model which is linearly and nonlinearly accurate up to deep water. Assuming the flow to be irrotational, we formulate the problem in terms of the velocity potential thereby lowering the number of unknowns. The model derivation combines two approaches, namely the method proposed by Agnon et al. (

Introduction

During the past two decades, the original Boussinesq (1872) model for flat bottom and its uneven bottom version derived by [START_REF] Peregrine | Long waves on a beach[END_REF] have been widely studied and extended to tackle more and more realistic physical problems. Consequently, Boussinesq-type models have emerged as an attractive and commonly used tool for coastal applications and engineering purposes. The derivation of such models is based on a polynomial approximation of the vertical profile of the velocity field, which allows to reduce the problem size from three to two space dimensions. These models are usually formulated as conservation equations for mass and momentum, including spatial and temporal derivatives of the free surface elevation and the velocity. In practice, their range of applicability is measured in terms of kh, where k is the wavenumber and h the water depth.

The conventional Boussinesq model for uneven bottoms [START_REF] Peregrine | Long waves on a beach[END_REF], which employs a quadratic polynomial approximation for the vertical flow distribution, is a depth-averaged model based on two fundamental assumptions, namely weak nonlinearity and weak dispersion. Its range of applicability is limited to kh < 0.75, as stated in [START_REF] Madsen | A new Boussinesq method for fully nonlinear waves from shallow to deep water[END_REF][START_REF] Madsen | Boussinesq-type formulations for fully nonlinear and extremely dispersive water waves: derivation and analysis[END_REF], so that this model has poor dispersion properties in intermediate depths. Moreover, the weakly nonlinear assumption limits the largest wave height that can be modelled accurately. As a result, substantial efforts have been devoted to extend the linear and nonlinear range of applicability of Boussinesq-type models. The first historical improvement, proposed by [START_REF] Nwogu | An alternative form of the Boussinesq equations for nearshore wave propagation[END_REF], consists in using a reference velocity at a specified depth, allowing the resulting model to be linearly applicable at intermediate depths. Similar models for short-amplitude and long waves have been more recently proposed and rigorously justified [START_REF] Bona | Boussinesq equations and other systems for smallamplitude long waves in nonlinear dispersive media, I: Derivation and linear theory[END_REF][START_REF] Bona | Long waves approximations for water waves[END_REF][START_REF] Chazel | Influence of bottom topography on long water waves[END_REF]). An effort similar to the one of [START_REF] Nwogu | An alternative form of the Boussinesq equations for nearshore wave propagation[END_REF] was pursued by [START_REF] Madsen | A new form of the Boussinesq equations with improved linear dispersion characteristics. II. A slowly varying bathymetry[END_REF], which was followed by the works of [START_REF] Liu | Model equations for wave propagation from deep to shallow water[END_REF] and [START_REF] Wei | A fully nonlinear Boussinesq model for surface waves: Part I. Highly nonlinear unsteady waves[END_REF] in which the authors efficiently removed the weak nonlinearity assumption, allowing the model to simulate wave propagation with strong nonlinear interaction. According to the reviews proposed by [START_REF] Madsen | Higher order Boussinesq-type equations for surface gravity waves-derivation and analysis[END_REF], 1999), these new Boussinesq-type models allow to extend the linear range of applicability to kh = 6, but it turned out that a similar improvement on the nonlinear characteristics was very difficult to reach. Then, so-called high-order Boussinesq-type models were derived to further enhance the deep-water linear and nonlinear accuracy using a higher-order (at least fourth-order) polynomial approximation for the vertical flow distribution. One such example is the formulation of [START_REF] Gobbi | A fully nonlinear Boussinesq model for surface waves. Part 2. Extension to O((kh) 4 )[END_REF] which uses a fourth-order polynomial approximation: this model exhibits excellent linear properties up to kh = 6 for the dispersion relation and up to kh = 4 for the vertical profiles of orbital velocities, whereas nonlinear behaviour is fairly well captured up to kh = 3. The price to pay for such an improvement is a significant increase in computational cost, since this model includes up to fifth-order derivatives and therefore requires a complex numerical scheme.

Over the past decade, two parallel approaches have emerged, one aiming at extending even more the range of applicability of the model of [START_REF] Gobbi | A fully nonlinear Boussinesq model for surface waves. Part 2. Extension to O((kh) 4 )[END_REF] into very deep water without increasing too much the numerical complexity, and another aiming at lowering the numerical cost of the latter while at least preserving the linear and nonlinear properties. The first approach has been extensively investigated by Madsen and co-workers with a first breakthrough in 1999 [START_REF] Agnon | A new approach to high-order Boussinesq models[END_REF]. In this work, the authors presented a new procedure by which it was possible to achieve the same accuracy on both linear and nonlinear properties. The main idea is to obtain approximate solutions to the Laplace equation (combined with the exact nonlinear free surface and bottom conditions) through truncated series expansions. While the formulation of [START_REF] Agnon | A new approach to high-order Boussinesq models[END_REF] involves velocity variables evaluated at the still-water level and is limited to kh = 6, the extended method proposed Article submitted to Royal Society by [START_REF] Madsen | A new Boussinesq method for fully nonlinear waves from shallow to deep water[END_REF][START_REF] Madsen | Boussinesq-type formulations for fully nonlinear and extremely dispersive water waves: derivation and analysis[END_REF] completely removes the conventional shallow-water limitation, allowing for modelling fully nonlinear and highly dispersive waves up to kh ≈ 40, i.e. in very deep water. This extended approach is based on velocity variables taken at optimized levels and optimal expansions through the use of Padé approximants. An extension to rapidly varying bathymetry has been proposed recently [START_REF] Madsen | A Boussinesq-type method for fully nonlinear waves interacting with a rapidly varying bathymetry[END_REF]). Although a few numerical approaches based on this model have been proposed [START_REF] Fuhrman | Numerical solutions of fully non-linear and highly dispersive Boussinesq equations in two horizontal dimensions[END_REF][START_REF] Engsig-Karup | Nodal DG-FEM solution of high-order Boussinesq-type equations[END_REF][START_REF] Engsig-Karup | DG-FEM solution for nonlinear wave-structure interaction using Boussinesq-type equations[END_REF], the counterpart for this wide range of applicability is the numerical complexity of the underlying model which includes derivatives up to fifth-order and consists of more equations (and more unknowns) than the alternative model of [START_REF] Gobbi | A fully nonlinear Boussinesq model for surface waves. Part 2. Extension to O((kh) 4 )[END_REF]. An interesting alternative approach has been chosen by [START_REF] Jamois | A numerical study of nonlinear wave run-up on a vertical plate[END_REF], where the authors propose to use a velocity potential formulation and to truncate the infinite series expansions of the solutions to the Laplace equation at a lower order: the resulting model is linearly and nonlinearly accurate only up to kh = 10, but entails a much lighter numerical complexity with less equations and with derivatives up to fourth-order only.

The second approach that has been studied is the double-layer approach, as proposed among others by Lynett & Liu (2004a, b) and [START_REF] Audusse | A multilayer Saint-Venant system: derivation and numerical validation[END_REF]. This approach is based on the idea of trading fewer unknowns and higher spatial derivatives for more unknowns and lower spatial derivatives. The multi-layering concept developed in the above references can be seen as an efficient mathematical tool to reduce the order of derivatives in any model, while increasing its linear and nonlinear range of applicability. The double-layer modelling proposed by Lynett & Liu (2004a, b) is purely conceptual since the two layers have the same density. However, the resulting model, which is based on classical depth-integrated Boussinesq-type equations, allows to model accurately wave propagation up to kh = 6, both linearly and nonlinearly. This offers a very interesting alternative to high-order models, since this double-layer model is less complex (including derivatives up to third-order only) and more accurate in deep water.

The present work is mainly inspired by these two approaches, namely the one of [START_REF] Madsen | A new Boussinesq method for fully nonlinear waves from shallow to deep water[END_REF][START_REF] Madsen | Boussinesq-type formulations for fully nonlinear and extremely dispersive water waves: derivation and analysis[END_REF] and the one of Lynett & Liu (2004a, b). The primary goal of this paper is to offer an efficient alternative to the models of [START_REF] Madsen | A new Boussinesq method for fully nonlinear waves from shallow to deep water[END_REF][START_REF] Madsen | Boussinesq-type formulations for fully nonlinear and extremely dispersive water waves: derivation and analysis[END_REF] by mixing their procedure, the simplifications of [START_REF] Jamois | A numerical study of nonlinear wave run-up on a vertical plate[END_REF], and the double-layer concept of Lynett & Liu (2004a, b). We aim at deriving a model which is 1) applicable to complex domains such as coastal areas, islands or estuaries and 2) linearly and nonlinearly accurate up to deep water, but with lower complexity than the previous models (i.e. lower order of derivatives and lower number of equations). The model proposed herein satisfies all these conditions as it exhibits excellent linear and nonlinear dispersive properties up to kh = 10, consists of four equations in both one and two horizontal dimensions (denoted by 1DH and 2DH, respectively), and includes up to second-order spatial derivatives only.

The present model hinges on a static Dirichlet-Neumann operator and its approximation using a double-layer technique and Padé approximants. The advantage of using the static operator (that is, defined on a fixed domain) as opposed to the usual Dirichlet-Neumann operator defined at the free surface is that the static operator (or its approximation) can be computed once and for all. The Dirichlet-Neumann operator has been extensively investigated over the past two decades. Exact expressions of the static operator, thereby leading to exact dispersion rela-tions, can be found in the work of [START_REF] Craig | Numerical simulation of gravity waves[END_REF], [START_REF] Dommermuth | A high-order spectral method for the study of nonlinear gravity waves[END_REF], [START_REF] Smith | An operator expansion formalism for nonlinear surface waves over variable depth[END_REF], and in the work of [START_REF] Matsuno | Nonlinear evolution of surface gravity waves over an uneven bottom[END_REF] and its extension to very general bathymetries by Artiles & Nachbin (2004a, b). However, the application of the above approaches to complex 2DH domains has not been reported yet. Thus, with an eye toward coastal engineering applications, we prefer to base our approach on approximating the static Dirichlet-Neumann operator. Furthermore, we mention the new promising approach of [START_REF] Ablowitz | On a new non-local formulation of water waves[END_REF] based on an exact integral representation of the usual Dirichlet-Neumann operator where no approximation is needed. Finally, we observe that the double-layer technique used to derive the approximate static Dirichlet-Neumann operator helps reducing the order of the derivatives while improving the accuracy of the model.

The paper is organized as follows. In §2, our model is formulated in terms of a static Dirichlet-Neumann operator, and we derive in §3 an approximation to this operator. A linear analysis of the model is presented in §4, including linear dispersion, vertical profiles of velocities, and linear shoaling. These properties are optimized based on Stokes linear wave theory, and it is shown that the model is accurate even for deep water conditions. Finally, in §5, numerical simulations are developed in 1DH to assess the nonlinear properties of the model for flat bottom conditions.

Derivation of the double-layer formulation (a) Governing equations and boundary conditions

We aim at formulating a double-layer Boussinesq-type model for the threedimensional irrotational flow of an inviscid and incompressible fluid with a free surface. We focus here on so-called gravity waves or water waves, i.e. the evolution of a fluid under the only influence of gravity. The capillary effects due to the presence of surface tension are neglected. Moreover, we assume constant atmospheric pressure at the free surface of the fluid. We adopt a Cartesian coordinate system, where we denote by X = (x, y) the horizontal coordinates and by z the vertical one, the z-axis pointing upwards. The time-dependent fluid domain is bounded from below by the static sea bottom and from above by the time-dependent free surface. We restrict this study to the case where the bathymetry and the free surface elevation are single-valued continuous functions, i.e. they can be described by the graph of two functions X → z(X) = -h(X) and (t, X) → η(t, X) respectively. The level z = 0 corresponds to the still-water level. As shown in figure 1, the fluid is divided into two layers by an interface z = ẑ(X) = -σh(X) where σ is an arbitrary parameter in ]0, 1[. Thus, the thickness of the two layers are constant fractions of the stillwater depth and do not depend on the free surface elevation. Unless the bottom is flat, the interface level ẑ is therefore spatially (but not temporally) variable. The upper layer of the fluid is denoted by Ω 1 and the lower layer by Ω 2 , namely

Ω 1 (t) = { (X, z) ; ẑ(X) ≤ z ≤ η(t, X)} and Ω 2 = { (X, z) ; z(X) ≤ z ≤ ẑ(X)}.
We point out that this fluid division into two layers is purely conceptual since both layers have the same density.

As far as the bathymetry is concerned, we assume in this work that the stillwater depth h verifies |∇h| 1, which corresponds to the classical mild-slope approximation. This approximation consists in neglecting all the quadratic (and

0 Ω1 Ω2 z=η(t,X) z=ẑ(X) z=z(X)
Figure 1. Representation of the fluid domain higher order) terms in ∇h along with the derivatives of h of order greater than one. Physically, this approximation means that we assume the wavelength of the free surface waves to be shorter than the distance over which the bathymetry (and thus the still-water depth) varies appreciably. We point out that, in the mild-slope framework, the overall amplitude of bottom topography levels can still be large. The mild-slope approximation plays an important role in the derivation and linear optimization of the present model, and it seems quite arduous to incorporate higher-order bathymetric terms without significantly increasing the complexity of the model. For very general bathymetries in 1DH, we refer to Artiles & Nachbin (2004a, b). For clarity, all the equations derived using this mild-slope approximation are indicated in this work by the symbol m ≈ instead of the equality symbol. Since the flow is assumed irrotational, there exists a velocity potential φ such that u ≡ ∇φ, w ≡ ∂ z φ, where u denotes the horizontal velocity of the fluid, w the vertical velocity, and ∇ the horizontal gradient operator (∂ x , ∂ y ) T . We define the velocity potentials φ i and the vertical velocities w i in each layer by

φ i = φ | Ω i , w i = ∂ z φ i
where the subscript i ∈ {1, 2} denotes the layer index.

The fluid motion in each layer is governed by the following equations written in terms of the velocity potential φ i and the vertical velocity w i ,    (2.1b) where P i is the reduced pressure field in each layer, g the gravitational acceleration, and R the Bernoulli constant. Equation (2.1a) corresponds to the Laplace (or continuity) equation and (2.1b) corresponds to the Bernoulli (or momentum) equation. The Bernoulli constant R only depends on time. Therefore, up to a time-dependent shift of the velocity potential, we can take this constant equal to P atm where P atm is the constant atmospheric pressure at the free surface. At the free surface z = η(t, X), the following boundary conditions are enforced:

∆φ i + ∂ 2 z φ i = 0 , (X, z) ∈ Ω i , (2.1a) ∂ t φ i + 1 2 |∇φ i | 2 + 1 2 w 2 i + gz + P i = R , (X, z) ∈ Ω i ,
∂ t η + ∇η • ∇φ 1 -w 1 = 0 , at z = η , (2.2a) P 1 = P atm , at z = η , (2.2b)
where (2.2a) is the usual kinematic free surface condition expressing that the free surface is a bounding surface, i.e. no fluid particle can cross it. At the interface z =

Article submitted to Royal Society ẑ(X) between the layers, the following natural continuity conditions are enforced:

φ 1 = φ 2 , w 1 = w 2 , P 1 = P 2 , at z = ẑ . (2.3)
Observe that φ 1 = φ 2 and w 1 = w 2 at z = ẑ imply ∇φ 1 = ∇φ 2 and thus u 1 = u 2 at z = ẑ, hence recovering the continuity conditions on the horizontal and vertical velocities enforced by Lynett & Liu (2004a). Finally, the system of equations is closed by the usual kinematic condition at the sea bottom z = z(X):

∇h • ∇φ 2 + w 2 = 0 , at z = z , (2.4)
which expresses that the sea bottom is a bounding surface. We now introduce

     φ 1 (t, X) = φ 1 (t, X, z = η(t, X)) , w 1 = w 1 (t, X, z = η(t, X)) , φ i (t, X) = φ i (t, X, z = ẑ(X)) , w i = w i (t, X, z = ẑ(X)) , φ 2 (t, X) = φ 2 (t, X, z = z(X)) , w 2 = w 2 (t, X, z = z(X)) .
Following [START_REF] Zakharov | Stability of periodic waves of finite amplitude on the surface of a deep fluid[END_REF], we can reformulate the equations as

∂ t φ 1 + 1 2 |∇ φ 1 | 2 - 1 2 w 1 2 (1 + |∇η| 2 ) + gη = 0 , (2.5a) ∂ t η + ∇η • ∇ φ 1 -w 1 (1 + |∇η| 2 ) = 0 , (2.5b)
where (2.5a) is the Euler equation expressed at the free surface and (2.5b) the kinematic condition at the free surface,

∆φ 1 + ∂ z w 1 = 0 , ẑ ≤ z ≤ η , 1 2 (2.5c) ∆φ 2 + ∂ z w 2 = 0 , z ≤ z ≤ ẑ , 1 (2.5d)
where (2.5c) and (2.5d) are the Laplace equations in each layer, and

φ 1 = φ 2 , 1 2 (2.5e
)

w 1 = w 2 , 1 |2 (2.5f) w 2 + ∇h • ∇φ 2 m ≈ 0 , 1 (2.5g)
where (2.5e) and (2.5f) correspond to the continuity conditions at the interface, and (2.5g) to the kinematic condition at the bottom, where we used the mildslope hypothesis to neglect the |∇h| 2 term. We point out that we work with a velocity potential formulation, unlike [START_REF] Madsen | A new Boussinesq method for fully nonlinear waves from shallow to deep water[END_REF][START_REF] Madsen | Boussinesq-type formulations for fully nonlinear and extremely dispersive water waves: derivation and analysis[END_REF] who formulated the governing equations in terms of the velocity variables u and w. This choice stems from our will to minimize the total number of equations in the model: in 2DH and in the present double-layer framework under the irrotational flow assumption, our velocity potential formulation allows to consider two less equations than with a velocity formulation.

In the sequel, equations (2.5a) and (2.5b) are left unchanged as they define the fully nonlinear time-stepping problem. We focus on the Laplace equations and the remaining boundary conditions to close the time-stepping problem by expressing the Article submitted to Royal Society vertical velocity at the free surface w 1 in terms of the velocity potential φ 1 at that surface, the free surface η, and the bathymetry h. This relation corresponds to the well-known Dirichlet-Neumann operator. The next subsections and §3 are devoted to the crucial construction of an accurate, yet computationally cheap approximation to this operator.

w 1 = G[η, h] φ 1 .
(2.6)

This Dirichlet-Neumann operator is at the heart of the derivation of Boussinesqtype models since the structure and accuracy of these models essentially depend on the method used to construct an approximation to this operator. Once this approximation is derived, there are two options. The first one is to eliminate the vertical velocity variable w 1 from the equations by plugging (2.6) into (2.5a), (2.5b). This method has been classically used in Boussinesq-type models and has the advantage of lowering the number of equations to solve at each time step, but significantly increases their complexity. The other option has been used for instance by [START_REF] Madsen | A new Boussinesq method for fully nonlinear waves from shallow to deep water[END_REF][START_REF] Madsen | Boussinesq-type formulations for fully nonlinear and extremely dispersive water waves: derivation and analysis[END_REF] and consists in keeping w 1 in the equations, which entails to solve (2.5a), (2.5b) and then to compute w 1 through the use of the Dirichlet-Neumann operator G[η, h] at each time step. This is the method we have chosen to use to keep equations complexity to a minimum. The main difficulty in finding an approximation to the Dirichlet-Neumann operator is that it involves solving the Laplace equations (2.5c) and (2.5d) along with the boundary conditions (2.5e)-(2.5g) and φ 1 = φ 1 at z = η, on a timedependent domain bounded from above by the free surface z = η. Keeping w 1 in the equations involves constructing an approximation to G[η, h] at each time step, which can be a serious computational problem as it increases the numerical cost. An interesting work-around to this issue consists in constructing an alternative Dirichlet-Neumann operator expressed at the still-water level, and then finding a closure between the unknown functions at the free surface z = η and the ones at the still-water level z = 0. To this end, we introduce φ 0 (t, X) = φ 1 (t, X, z = 0) and w 0 (t, X) = w 1 (t, X, z = 0) and denote by G 0

[h] the Dirichlet-Neumann operator corresponding to η = 0, i.e. G 0 [h] = G[0, h] . (2.7) The translated operator G 0 [h] is such that G 0 [h]ψ = ∂ z φ 1 | z=0 = w 0 where (φ 1 , φ 2 ) solves the boundary value problem            ∆φ 1 + ∂ 2 z φ 1 = 0 , ẑ ≤ z ≤ 0 , ∆φ 2 + ∂ 2 z φ 2 = 0 , z ≤ z ≤ ẑ , φ 1 = ψ , at z = 0 , φ 1 = φ 2 , ∂ z φ 1 = ∂ z φ 2 , at z = ẑ , ∂ z φ 2 + ∇h • ∇φ 2 = 0 , at z = z .
(2.8)

Article submitted to Royal Society

The main advantage of this translated Dirichlet-Neumann operator is that the boundary value problem (2.8) is posed on a static domain bounded from above by the still-water level z = 0 and from below by the sea bottom z = z.

Remark 2.1. We point out that φ 0 is not defined everywhere since the bounding free surface η can take negative values. This issue can be solved by extending the solutions to the Laplace equation above z = η when η < 0, using the fact that both φ 1 and its normal derivative w 1 are continuous at the free surface. This mathematical trick allows to artificially define φ 0 and w 0 when η < 0. This tool has been implicitly used by many authors such as [START_REF] Nwogu | An alternative form of the Boussinesq equations for nearshore wave propagation[END_REF], [START_REF] Wei | A fully nonlinear Boussinesq model for surface waves: Part I. Highly nonlinear unsteady waves[END_REF], [START_REF] Gobbi | A fully nonlinear Boussinesq model for surface waves. Part 2. Extension to O((kh) 4 )[END_REF] or [START_REF] Madsen | A new Boussinesq method for fully nonlinear waves from shallow to deep water[END_REF][START_REF] Madsen | Boussinesq-type formulations for fully nonlinear and extremely dispersive water waves: derivation and analysis[END_REF] who derived models based on an horizontal velocity (or potential) variable taken at free-surface independent levels, thus allowing these levels to exceed the bounding value z = η for η negative enough.

(c) Closure relation and model formulation

Our goal is twofold. Firstly to construct an approximation to the translated Dirichlet-Neumann operator G 0 [h], and secondly to look for closure relations between the unknowns φ 1 , w 1 , φ 0 , and w 0 . The second objective can be readily achieved via a Taylor expansion of φ 1 and w 1 at the still-water level z = 0. Indeed, combining the MacLaurin expansions of φ 1 and w 1 at respectively the fourth and third order (see Remark 2.2) and the Laplace equation ∆φ 1 = -∂ z w 1 at z = 0 yields the desired closure relations, namely

φ 1 = (1 -α 1 ∆) φ 0 + ( β 1 -γ 1 ∆) w 0 , w 1 = -β 1 ∆ φ 0 + (1 -α 1 ∆) w 0 ,
(2.9)

where

α 1 = η 2 2 , β 1 = η, and γ 1 = η 3 6 .
We can now state our double-layer model as follows:

                     ∂ t φ 1 + 1 2 |∇ φ 1 | 2 - 1 2 w 1 2 (1 + |∇η| 2 ) + gη = 0 , ∂ t η + ∇η • ∇ φ 1 -w 1 (1 + |∇η| 2 ) = 0 , 1 - η 2 2 ∆ + (η - η 3 6 ∆) G app 0 [h] φ 0 = φ 1 , w 1 = -η∆ + (1 - η 2 2 ∆) G app 0 [h] φ 0 .
(2.10)

where G app 0

[h] is an approximation to the static Dirichlet-Neumann operator G 0 [h] which is detailed in the following section.

The main advantages of this model are that 1) it consists of only four equations both in 1DH and 2DH, 2) it can be used in complex domains, and 3) the approximate Dirichlet-Neumann operator G app 0

[h] can be computed once and for all at t = 0 since this operator is static. Furthermore, we will see in §3 that it only includes at most second-order horizontal derivatives. This is a major improvement in comparison with high-order Boussinesq-type models such as those of Jamois et Article submitted to Royal Society al. (2006) and [START_REF] Madsen | A new Boussinesq method for fully nonlinear waves from shallow to deep water[END_REF][START_REF] Madsen | Boussinesq-type formulations for fully nonlinear and extremely dispersive water waves: derivation and analysis[END_REF] which contain respectively fourth-and fifth-order derivatives, and consist of respectively five equations in 1DH and 2DH, and five equations in 1DH and seven in 2DH.

Remark 2.2. The truncation orders used in (2.9) can be motivated by a dimensional analysis. We scale the vertical coordinate z and the free surface η with the typical amplitude a of the surface waves, the horizontal coordinate X by the typical wavelength λ and introduce the mean depth h 0 . The truncation errors in the two equations of (2.9) are respectively of order O( ε 4 µ 2 , ε 5 µ 2 ) andO(ε 4 µ 2 , ε 3 µ 2 ), where ε = a/h 0 and µ = h 2 0 /λ 2 correspond respectively to the nonlinearity and dispersion parameters. Thus, the truncated terms are almost third and fourth powers of the steepness parameter s defined by s = ε √ µ, whose maximum value is admittedly s max ≈ 0.142 for a stable wave [START_REF] Williams | Limiting gravity waves in water of finite depth[END_REF]. Combining this value and the fact that we consider fully nonlinear waves, for which ε is of order O(1), motivates the truncation order in (2.9).

An approximate static Dirichlet-Neumann operator (a) Theoretical solutions to the Laplace equations

The first step in the derivation of the approximate static Dirichlet-Neumann operator G app 0 [h] is to look for solutions to the Laplace equations

∆φ i + ∂ 2 z φ i = 0 , (X, z) ∈ Ω i , (3.1)
where we have redefined the upper-layer domain as Ω

1 = { (X, z) ; ẑ(X) ≤ z ≤ 0}.
To this end, we follow the generalized Boussinesq procedure introduced by [START_REF] Madsen | A new Boussinesq method for fully nonlinear waves from shallow to deep water[END_REF][START_REF] Madsen | Boussinesq-type formulations for fully nonlinear and extremely dispersive water waves: derivation and analysis[END_REF] which consists in looking for a solution under the form of an infinite Taylor series in the vertical coordinate. The main difference between this method and the classical Boussinesq procedure is that in the latter, one looks for a finite series solution in the vertical coordinate, i.e. a low-order polynomial in the variable z. The generalized method of [START_REF] Madsen | A new Boussinesq method for fully nonlinear waves from shallow to deep water[END_REF][START_REF] Madsen | Boussinesq-type formulations for fully nonlinear and extremely dispersive water waves: derivation and analysis[END_REF] allows to find exact infinite series solutions instead of approximate solutions. We first introduce two arbitrary expansion levels z 1 and z 2 in each layer, namely z 1 (X) = -σ 1 h(X) with 0 < σ 1 < σ and z 2 (X) = -σ 2 h(X) with σ < σ 2 < 1, and the associated unknowns φi and wi such that φi = φ i (X, z = z i (X)) and wi = w i (X, z = z i (X)) for i ∈ {1, 2}. We now look for solutions to the Laplace equations in the form of infinite Taylor series in the vertical variable (zz i ),

φ i (X, z) = n≥0 (z -z i ) n φ (n) i (X) , (3.2)
where the choice of the vertical variable (zz i ) instead of z actually allows to save one step compared to the procedure of [START_REF] Madsen | A new Boussinesq method for fully nonlinear waves from shallow to deep water[END_REF][START_REF] Madsen | Boussinesq-type formulations for fully nonlinear and extremely dispersive water waves: derivation and analysis[END_REF]. Injecting (3.2) into the Laplace equations (3.1) and using the mild-slope assumption leads to the recurrence relation

∆φ (n) i -2(n + 1)∇z i • ∇φ (n+1) i + (n + 2)(n + 1)φ (n+2) i m ≈ 0 . (3.3) Article submitted to Royal Society Observing that ∆(∇z i • ∇φ (k) i ) m ≈ ∇z i • ∆∇φ (k)
i for all k yields the expression of φ (2n) i and φ

(2n+1) i in terms of φi and wi ,        φ (2n) i m ≈ (-1) n (2n)! ∆ n φi + (-1) n-1 (2n -1)! ∇z i • ∆ n-1 ∇ wi , φ (2n+1) i m ≈ (-1) n (2n + 1)! ∆ n wi + (-1) n 2n (2n + 1)! ∇z i • ∆ n ∇ φi .
Plugging these expressions into the ansatz (3.2) provides the desired expressions for the velocity potentials φ i and vertical velocity w i in terms of φi and wi for all (X, z) belonging to Ω i ,

φ i (X, z) = C(z -z i ) φi + S(z -z i ) wi + ∇z i • Γ φi , w i (X, z) = -S(z -z i )∆ φi + C(z -z i ) wi + ∇z i • Γ wi , (3.4)
where C and S are infinite-series pseudo-differential operators defined by

C(λ) = n≥0 (-1) n λ 2n (2n)! ∆ n , S(λ) = n≥0 (-1) n λ 2n+1 (2n + 1)! ∆ n , (3.5)
and where the slope terms Γ φi and Γ wi are given by

   Γ φi = (z -z i ) C(z -z i )∇ φi + S(z -z i )∇ wi -S(z -z i )∇ φi , Γ wi = (z -z i ) -S(z -z i )∇∆ φi + C(z -z i )∇ wi + S(z -z i )∇ wi .
(3.6)

The expression (3.4) of φ i provides a theoretical formulation of an exact solution to the Laplace equations (3.1). Strictly speaking, we can verify that they are in fact solutions to (3.1) with residuals of order O(|∇h| 2 , ∆h) which are negligible within our mild-slope approximation framework.

(b) Truncation of the Taylor series

The previous solutions to the Laplace equations (3.1) are purely theoretical since they involve infinite-series pseudo-differential operators. To deal with this problem, we obviously need to truncate the series at a finite order, and this raises the question of choosing the order of truncation. Through this choice, we have to reach a compromise between the accuracy of the truncated expression and the numerical complexity of the final model. In fact, the truncation order is a key factor for the domain of validity of the model: the higher the truncation order is, the better dispersive effects are reproduced, so that the model is applicable in deeper water. We recover here the common paradigm encountered in works based on asymptotic expansions of the solutions to (3.1) where retaining higher order terms increases the domain of validity in intermediate or deep water.

Within our double-layer framework, the increase of the number of unknowns allows us to lower the truncation order in comparison with the works of [START_REF] Madsen | A new Boussinesq method for fully nonlinear waves from shallow to deep water[END_REF][START_REF] Madsen | Boussinesq-type formulations for fully nonlinear and extremely dispersive water waves: derivation and analysis[END_REF], and even of [START_REF] Jamois | A numerical study of nonlinear wave run-up on a vertical plate[END_REF]. We take advantage of this possibility, and truncate the operators C and S by retaining only the first two terms of the series, which leads to the following approximations

C(λ) = 1 - λ 2 2 ∆ + O(λ 4 ∆ 2 ) , S(λ) = λ - λ 3 6 ∆ + O(λ 5 ∆ 2 ) . (3.7)
Article submitted to Royal Society

We first plug these approximations into (3.4), which leads to the following truncated expressions of φ i and

w i                φ i (X, z) = (1 -α i ∆) φi + (β i -γ i ∆) wi + ∇z i • β i (∇ -α i ∇∆) φi +(β i ∇ -γ i ∇∆) wi -(β i ∇ -γ i ∇∆) φi , w i (X, z) = (-β i ∆ + γ i ∆ 2 ) φi + (1 -α i ∆) wi + ∇z i • β i (-β i ∇∆ +γ i ∇∆ 2 ) φi + (∇ -α i ∇∆) wi + (β i ∇ -γ i ∇∆) wi , (3.8)
where

α i (z) = (z -z i ) 2 2 , β i (z) = z -z i , γ i (z) = (z -z i ) 3 6 . (3.9)
We now reformulate these expressions by applying the operators P i , defined for all smooth enough scalar-valued function u by

P i u = u -β i ∇z i • ∇u, which yields the following expressions for all (X, z) in Ω i ,      P i φ i m ≈ (1 -α i ∆) φi + (β i -γ i ∆) wi -∇z i • (β i ∇ -γ i ∇∆) φi , P i w i m ≈ (-β i ∆ + γ i ∆ 2 ) φi + (1 -α i ∆) wi + ∇z i • (β i ∇ -γ i ∇∆) wi .
(3.10)

The advantage of this new formulation will become clear in the sequel.

Remark 3.1. The chosen order of truncation can be motivated as in remark 2.2 by a dimensional analysis in the case of a flat bottom. Scaling here z and the expansion levels z i with h 1 = σh 0 in the upper-layer and h 2 = (1σ)h 0 in the lower-layer, the dimensionless form of (3.4)-(3.5) is obtained by replacing λ 2 by µ i in (3.5), where µ i = h 2 i /λ 2 corresponds to the dispersion parameter in each layer. Supposing that all the derivatives are of order O(1), we can analyze the order of the third terms of each series in deep water, for instance kh 0 = 10, where µ i is considerably higher than in shallow water. If we restrict σ to the range [0.25, 0.75], this leads to the following estimates for n = 2 in deep water: µ n i /((2n)!) ≈ 0.08, µ n i /((2n + 1)!) ≈ 0.01 and µ n+1 i /((2n+1)!) ≈ 0.02. Thus, the third term of each infinite series in (3.4)-(3.5) is very small in deep water: these terms and all the subsequent ones can be neglected. The same kind of analysis performed for an uneven bottom leads to the same result. This motivates the truncation order chosen in (3.8).

(c) Padé approximants

Using the truncated expression (3.10), it is possible to construct an approximation to the static Dirichlet-Neumann operator, but involving up to fourth-order differential operators at first order in h and fifth-order differential operators in the slope terms. Consequently, we now present a method to lower the maximum order of the derivatives in (3.10) based on Padé approximants. We follow the strategy introduced by [START_REF] Madsen | A new Boussinesq method for fully nonlinear waves from shallow to deep water[END_REF][START_REF] Madsen | Boussinesq-type formulations for fully nonlinear and extremely dispersive water waves: derivation and analysis[END_REF] and expand the variables φi and wi in terms of auxiliary variables φ * i and w * i through the relations

φi = M i (z i ∇) φ * i , wi = M i (z i ∇) w * i , M i (z i ∇) = 1 + p i z 2 i ∆ + q i z i ∇z i • ∇ , (3.11)
Article submitted to Royal Society where p i and q i are arbitrary coefficients to be determined. We now plug this ansatz into (3.10) and conduct the same dimensional analysis as previously, which yields the following expressions:

                           P i φ i m ≈ 1 -(α i -p i z 2 i )∆ φ * i + β i -(γ i -β i p i z 2 i )∆ w * i +∇z i • (q i z i -β i )∇ + (γ i -β i p i z 2 i -α i (q i + 4p i )z i )∇∆ φ * i + (β i q i z i ∇ -γ i (q i + 4p i )z i ∇∆) w * i , P i w i m ≈ -β i ∆ + (γ i -β i p i z 2 i )∆ 2 φ * i + 1 -(α i -p i z 2 i )∆ w * i +∇z i • (-β i (q i + 4p i )z i ∇∆) φ * i + (β i + q i z i )∇z 2 i -(γ i -β i p i z 2 i + α i (q i + 4p i )z i )∇∆ w * i ,
(3.12)

where we have again kept the first two terms in each modified truncated series, except the fifth-derivative of φ * i appearing in the slope terms of P i w i . This choice is motivated by the mild-slope approximation, but it clearly unbalances the global structure of the slope terms and does have a negative impact on the linear shoaling. However, we present in §4 a remedy to this problem.

We now aim at lowering the maximum order of the derivatives in (3.12) while preserving the overall accuracy of the truncated expressions (3.10). This goal can be achieved by choosing the coefficients p i and q i in order to introduce Padé approximants in the equations. In [START_REF] Madsen | A new Boussinesq method for fully nonlinear waves from shallow to deep water[END_REF][START_REF] Madsen | Boussinesq-type formulations for fully nonlinear and extremely dispersive water waves: derivation and analysis[END_REF], the authors use Padé approximants as a way to improve truncation accuracy without increasing the order of the derivatives. In the present work, we rather view the Padé approximants as a means to cancel high-order derivatives in (3.12) while preserving the accuracy of the truncated expressions (3.10).

Practically, lowering the maximum derivative order in (3.12) means requiring that the factors γ iβ i p i z 2 i and q i + 4p i (respectively in front of the fourth-order and third-order derivatives) vanish in each layer, via an appropriate choice of the constants p i , q i , and σ i (that define the expansion levels z i ). Since the quantity γ i -β i p i z 2 i depends on the vertical variable z through γ i and β i , it is impossible that this factor vanishes over the whole still-water column. Nevertheless, the truncated expressions (3.10) of P i φ i and P i w i need to be evaluated only at some levels, namely z = 0 for the upper boundary of Ω 1 , z = ẑ for the interface and z = z for the sea bottom. Consequently, requiring that the quantity γ 1β 1 p 1 z 2 1 vanishes at the stillwater level and at the interface, and that the quantity γ 2β 2 p 2 z 2 2 vanishes at the interface and at the sea bottom eliminates all the fourth-order derivatives from the model. Using (3.9), this yields

       p 1 = 1 6 , p 2 = 1 6 1 -σ 1 + σ 2 , σ 1 = σ 2 , σ 2 = σ + 1 2 .
Then, taking

q i = -4p i yields q 1 = - 2 3 and q 2 = - 2 3 1 -σ 1 + σ 2
. We observe that the expansion levels z 1 and z 2 , respectively defined by z 1 = -σ 1 h and z 2 = -σ 2 h, Article submitted to Royal Society are thus taken respectively at the middle of the upper layer at rest and at the middle of the lower one.

We can now plug the previous values for p i , q i , and z i into the expressions (3.12) of P i φ i and P i w i , and evaluate them at the three boundary levels z = 0, z = ẑ, and z = z. We first obtain the following expressions at the still-water level:

     1 + σ 2 β * 1 ∇h •∇ φ 0 = (1 -α * 1 ∆) φ * 1 + β * 1 w * 1 + ∇h • γ * 1 ∇φ * 1 -δ * 1 ∇w * 1 , 1 + σ 2 β * 1 ∇h •∇ w 0 = -β * 1 ∆φ * 1 + (1 -α * 1 ∆) w * 1 -ε * 1 ∇h •∇w * 1 , (3.13)
where

α * 1 = σ 2 12 h 2 , β * 1 = σ 2 h , γ * 1 = σ 2 12 h , δ * 1 = σ 3 12 h 2 , ε * 1 = 5σ 2 12 h . (3.14)
At the interface z = ẑ, we use the mild-slope approximation to obtain

                 1 -σ 2 β * 1 ∇h •∇ φ 1 m ≈ (1 -α * 1 ∆) φ * 1 -β * 1 w * 1 + ∇h •[-ε * 1 ∇φ * 1 + δ * 1 ∇w * 1 ], 1 -σ 2 β * 1 ∇h •∇ w 1 m ≈ β * 1 ∆φ * 1 + (1 -α * 1 ∆) w * 1 + ∇h •[γ * 1 ∇w * 1 ], 1 + σ+1 2 β * 2 ∇h•∇ φ 2 m ≈ (1 -α * 2 ∆) φ * 2 + β * 2 w * 2 + ∇h•[γ * 2 ∇φ * 2 -δ * 2 ∇w * 2 ], 1 + σ+1 2 β * 2 ∇h •∇ w 2 m ≈ -β * 2 ∆φ * 2 + (1 -α * 2 ∆) w * 2 -∇h •[ε * 2 ∇w * 2 ], (3.15)
where

α * 2 = (1 -σ) 2 12 h 2 , β * 2 = 1 -σ 2 h , γ * 2 = 5σ + 1 12 (1 -σ)h , δ * 2 = (1 -σ) 3 12 h 2 , ε * 2 = σ + 5 12 (1 -σ)h .
(3.16)

Finally, the expressions at the sea bottom z = z are

     1 -σ+1 2 β * 2 ∇h •∇ φ 2 m ≈ (1 -α * 2 ∆) φ * 2 -β * 2 w * 2 + ∇h •[δ * 2 ∇w * 2 -ε * 2 ∇φ * 2 ], 1 -σ+1 2 β * 2 ∇h •∇ w 2 m ≈ β * 2 ∆φ * 2 + (1 -α * 2 ∆) w * 2 + ∇h •[γ * 2 ∇w * 2 ].
(3.17)

(d) Formulation of the approximate static Dirichlet-Neumann operator

The final step in deriving our approximate static Dirichlet-Neumann operator G app 0 [h] is to reformulate the boundary conditions at z = ẑ and z = z in terms of P i φ i and P i w i instead of φ i and w i . At the interface, the two operators acting on φ 1 and w 1 on one side, and φ 2 and w 2 on the other side, differ from each other. A simple reformulation of the continuity condition

φ 1 = φ 2 is 1 - σ 2 β * 1 ∇h • ∇ φ 1 m ≈ 1 - h 4 ∇h • ∇ 1 + σ + 1 2 β * 2 ∇h • ∇ φ 2 .
Consequently, the continuity conditions (2.5e), (2.5f) take the form

     1 - σ 2 β * 1 ∇h • ∇ φ 1 = 1 - h 4 ∇h • ∇ 1 + σ + 1 2 β * 2 ∇h • ∇ φ 2 , 1 - σ 2 β * 1 ∇h • ∇ w 1 = 1 - h 4 ∇h • ∇ 1 + σ + 1 2 β * 2 ∇h • ∇ w 2 .
(3.18)

Article submitted to Royal Society

Finally, applying the operator 1-((σ+1)/2)β * 2 ∇h•∇ to (2.5g) leads to the following reformulation of the kinematic condition at the bottom

1 - σ + 1 2 β * 2 ∇h • ∇ w 2 + ∇h • ∇ 1 - σ + 1 2 β * 2 ∇h • ∇ φ 2 m ≈ 0 . (3.19)
Gathering all the previous results, we are able to construct the system of five equations on the six unknowns φ 0 , w 0 , φ * 1 , w * 1 , φ * 2 , and w * 2 that defines the approximate operator G app 0 [h] linking w 0 to φ 0 . The first one corresponds to the reformulated Dirichlet condition (P 1 φ 1 )| z=0 = (P 1 | z=0 )φ 0 , i.e. the first equation of (3.13). The second and third equations correspond to the continuity conditions at the interface z = ẑ, and are obtained by plugging (3.15) into (3.18). The fourth one is the condition at the sea bottom derived by plugging (3.17) into (3.19), and the last one corresponds to the Neumann condition (P 1 w 1 )| z=0 = (P 1 | z=0 )w 0 as expressed by the second equation of (3.13). These five equations can be recast as

                       M 11 M 12 0 0 M 21 M 22 M 23 M 24 M 31 M 32 M 33 M 34 0 0 M 43 M 44             φ * 1 w * 1 φ * 2 w * 2       =       P 0 φ 0 0 0 0       , P 0 w 0 = -β * 1 ∆φ * 1 + (1 -α * 1 ∆ -ε * 1 ∇h • ∇) w * 1 , (3.20)
with the differential operators

                                       P 0 = 1 + σ 2 β * 1 ∇h • ∇, M 11 = 1 -α * 1 ∆ + γ * 1 ∇h • ∇ , M 12 = β * 1 -δ * 1 ∇h • ∇ , M 21 = 1 -α * 1 ∆ -ε * 1 ∇h • ∇ , M 22 = -β * 1 + δ * 1 ∇h • ∇ , M 23 = -(1 -α * 2 ∆) -∇h • (γ * 2 - h 4 )∇ + h 4 α * 2 ∇∆ , M 24 = -β * 2 + (δ * 2 + h 4 β * 2 )∇h • ∇ , M 31 = β * 1 ∆ , M 32 = 1 -α * 1 ∆ + γ * 1 ∇h • ∇ , M 33 = β * 2 ∆ - h 4 β * 2 ∇h • ∇∆ , M 34 = -(1 -α * 2 ∆) + ∇h • (ε * 2 + h 4 )∇ - h 4 α * 2 ∇∆ , M 43 = β * 2 ∆ + ∇h • (∇ -α * 2 ∇∆) , M 44 = 1 -α * 2 ∆ + (γ * 2 -β * 2 )∇h • ∇ . (3.21)
Incidentally, we observe that it is possible to eliminate all the third-order derivatives from the operators M 23 , M 33 , M 34 , and M 43 in these equations. Indeed, applying the operator ∇h • ∇ to respectively the fourth, first, and second equation of (3.20), using the mild-slope approximation and combining the results yields

α * 2 ∇h • ∇∆w * 2 m ≈ ∇h • ∇w * 2 + β * 2 ∇h • ∇∆φ * 2 , α * 1 ∇h • ∇∆φ * 1 m ≈ ∇h • ∇φ * 1 + β * 1 ∇h • ∇w * 1 -∇h • ∇φ 0 , α * 2 ∇h • ∇∆φ * 2 m ≈ 2β * 1 ∇h • ∇w * 1 + ∇h • ∇φ * 2 + β * 2 ∇h • ∇w * 2 -∇h • ∇φ 0 .
(3.22)

Article submitted to Royal Society

Plugging these results into (3.20), we formulate our approximate static Dirichlet-Neumann operator G app 0 [h] as follows:

                       M 11 M 12 0 0 M 21 M 22 M 23 M 24 M 31 M 32 M 33 M 34 0 M 42 M 43 M 44             φ * 1 w * 1 φ * 2 w * 2       =       P 0 φ 0 Q 1 φ 0 Q 2 φ 0 Q 3 φ 0       , P 0 w 0 = -β * 1 ∆φ * 1 + (1 -α * 1 ∆ -ε * 1 ∇h • ∇) w * 1 , (3.23)
where the differential operators P 0 , M 11 , M 12 , M 21 , and M 31 are given by (3.21) and the new operators are defined by

                       M 22 = -β * 1 + (δ * 1 - h 2 β * 1 )∇h • ∇, M 23 = -(1 -α * 2 ∆) -γ * 2 ∇h • ∇, M 24 = -β * 2 + δ * 2 ∇h • ∇, M 32 = 1 -α * 1 ∆ + (γ * 1 - 3 σ 1 -σ h)∇h • ∇, M 33 = β * 2 ∆ - 3 1 -σ ∇h • ∇, M 34 = -(1 -α * 2 ∆) + (ε * 2 - 3 2 h)∇h • ∇, M 42 = -2β * 1 ∇h • ∇, M 43 = β * 2 ∆, M 44 = 1 -α * 2 ∆ + (γ * 2 -2β * 2 )∇h • ∇, Q 1 = - h 4 ∇h • ∇ , Q 2 = 3 σ -1 ∇h • ∇ , Q 3 = -∇h • ∇ .
(3.24) We denote by M = (M ij ) 1≤i,j≤4 the matrix differential operator linking φ * 1 , w * 1 , φ * 2 , and w * 2 to φ 0 in (3.23). We can write each of the operators M ij as a sum of a first-order (in h) operator P ij and a mild-slope operator ∇h • Q ij , yielding

M = P + ∇h • Q , ( 3.25)
where

P =      1 -α * 1 ∆ β * 1 0 0 1 -α * 1 ∆ -β * 1 -1 + α * 2 ∆ -β * 2 β * 1 ∆ 1 -α * 1 ∆ β * 2 ∆ -1 + α * 2 ∆ 0 0 β * 2 ∆ 1 -α * 2 ∆      , Q =      γ * 1 ∇ -δ * 1 ∇ 0 0 -ε * 1 ∇ (δ * 1 -h 2 β * 1 )∇ -γ * 2 ∇ δ * 2 ∇ 0 (γ * 1 -3 σ 1-σ h)∇ -3 1-σ ∇ (ε * 2 -3 2 h)∇ 0 -2β * 1 ∇ 0 (γ * 2 -2β * 2 )∇      .
We denote by U the vector (φ * 1 , w * 1 , φ * 2 , w * 2 ) T and by F the right-hand side (P 0 φ 0 , Q 1 φ 0 , Q 2 φ 0 , Q 3 φ 0 ) T . The differential system in (3.23) then takes the form

(P + ∇h • Q) U = F .
(3.26)

A Fourier analysis of the differential operator P shows that it can be inverted in the case of a flat bottom. Considering an uneven bottom, we can write P = P 0 + P H where P 0 = P(h 0 ), h 0 being the mean depth, and where P H is of order O(∇h). This ensures that the differential operator P is invertible for ∇h small enough. Finally, the differential operator R defined by

R = P -1 (Id -∇h • QP -1 ) , ( 3
     N 1 = -β * 1 ∆ , N 2 = 1 -α * 1 ∆ -ε * 1 ∇h • ∇ , Q 0 = 1 - σ 2 β * 1 ∇h • ∇ m ≈ P -1 0 , (3.29)
and plugging the expressions of φ * 1 and w * 1 into the last equation of (3.23), so as to obtain the explicit relation between w 0 and φ 0 , and thus the explicit expression of our approximate static Dirichlet-Neumann operator

G app 0 [h] G app 0 [h] = Q 0      N 1 0 0 0 0 N 2 0 0 0 0 0 0 0 0 0 0      R      P 0 Q 1 Q 2 Q 3      . (3.30)
This expression completes the formulation (2.10) of our double-layer Boussinesqtype model and will be further improved in §4 f to tighten the model shoaling properties. Once again, we stress that the major advantage of the operator G app 0

[h] is that it is static. Hence, we can construct it at t = 0 once and for all.

Once we have computed φ 0 , we can compute the values for φ * 1 , φ * 2 , and w * 2 using (3.28). Therefore, we can recover the vertical profiles of the velocity potentials φ 1 , φ 2 and the vertical velocities w 1 , w 2 over the whole water column using a generalization of (2.9) for any z ∈ (0, η), and plugging the computed values of p i and q i into (3.12). We point out that we have neglected the third-and fourth-order derivatives in the following expressions, to obtain only second-order derivatives. Of course, there is a price to pay for this choice, as discussed in the linear analysis of the vertical profiles in §4 e. We use the expressions (3.31) in the vertical region z ∈ (0, η) and the expressions

     φ 1 (t, X, z) = (1 - z 2 2 ∆) φ 0 + (z - z 3 6 ∆) w 0 , w 1 (t, X, z) = -z∆ φ 0 + (1 - z 2 2 ∆) w 0 ,
               φ i (t, X, z) m ≈ (1 + β i ∇z i • ∇) 1 -α ‡ i (z)∆ φ * i + β ‡ i (z) -γ ‡ i (z)∆ w * i -β ‡ i (z)∇z i • ∇φ * i , w i (t, X, z) m ≈ (1 + β i ∇z i • ∇) -β ‡ i (z)∆φ * i + 1 -α ‡ i (z)∆ w * i +β ‡ i (z)∇z i • ∇w * i , (3.32) Article submitted to Royal Society for z ∈ [ẑ, min(0, η)] if i = 1 and for z ∈ [z, ẑ] if i = 2, where            α ‡ 1 (z) = z 2 (z + σh) + σ 2 12 h 2 , α ‡ 2 (z) = 1 2 (z + h)(z + σh) + (1 -σ) 2 12 h 2 , β ‡ 1 (z) = z + σ 2 h , β ‡ 2 (z) = z + σ + 1 2 h , γ ‡ 1 (z) = z 6 (z + σ 2 h)(z + σh) , γ ‡ 2 (z) = 1 6 (z + h)(z + σ + 1 2 h)(z + σh) .
(3.33) Using (3.31) in the region between the still-water level and the free surface and (3.32) elsewhere instead of applying (3.32) everywhere seems to provide a more accurate description of the nonlinear profiles, as specified by [START_REF] Madsen | A new Boussinesq method for fully nonlinear waves from shallow to deep water[END_REF][START_REF] Madsen | Boussinesq-type formulations for fully nonlinear and extremely dispersive water waves: derivation and analysis[END_REF]. This property has also been observed during the nonlinear simulations performed on the present model in §5.

Linear analysis of the double-layer model

The goal of this section is to analyze the linear properties of the model (namely the phase and group velocities, the vertical profiles of velocity potential and vertical velocity, and the linear shoaling) and to optimize their accuracy in relation to the results of Stokes linear theory.

(a) Linearization of the governing equations

In order to investigate these linear properties, we restrict the analysis to the onedimensional case. We linearize the governing equations (2.10) around steady-state, which yields φ 1 = φ 0 and w 1 = w 0 , and leads to the linearized model

     ∂ t φ 0 + gη = 0 , (4.1a) ∂ t η -w 0 = 0 , (4.1b) w 0 = G app 0 [h]φ 0 , (4.1c) or equivalently, up to O(h 2 x ) terms,                          ∂ t φ 0 + gη = 0 , (4.2a) ∂ t η -w 0 = 0 , (4.2b)      M 11 M 12 0 0 M 21 M 22 M 23 M 24 M 31 M 32 M 33 M 34 0 M 42 M 43 M 44           φ * 1 w * 1 φ * 2 w * 2      =      P 0 φ 0 Q 1 φ 0 Q 2 φ 0 Q 3 φ 0      , (4.2c) P 0 w 0 = -β * 1 ∂ 2 x φ * 1 + 1 -α * 1 ∂ 2 x -ε * 1 h x ∂ x w * 1 , (4.2d)
where h x is the bottom slope and the differential operators M ij are left unchanged, except that they are here 1-D operators. We point out that in this linearized model, the slope terms are kept in order to investigate the linear shoaling properties. For convenience, we apply the operator P 0 to equations (4.2a), (4.2b) to obtain

                         ∂ t φ 0 + gN = 0 , (4.3a) ∂ t N -W 0 = 0 , (4.3b)      M 11 M 12 0 0 M 21 M 22 M 23 M 24 M 31 M 32 M 33 M 34 0 M 42 M 43 M 44           φ * 1 w * 1 φ * 2 w * 2      =      φ 0 Q 1 φ 0 Q 2 φ 0 Q 3 φ 0      , (4.3c) W 0 = -β * 1 ∂ 2 x φ * 1 + 1 -α * 1 ∂ 2 x -ε * 1 h x ∂ x w * 1 , ( 4.3d)
where φ 0 = P 0 φ 0 , N = P 0 η, W 0 = P 0 w 0 , and where we have plugged the relation

P -1 0 m ≈ 1 -(σ/2)β * 1 h x ∂ x into φ 0 = P -1 0 φ 0 to obtain Q i φ 0 = Q i P -1 0 φ 0 m ≈ Q i φ 0 for i ∈ {1, 3}, since Q i are differential operators of order O(h x ).
Plugging the expression of φ 0 and W 0 in terms of φ * 1 and w * 1 (given respectively by the first line of the differential system (4.3c) and by (4.3d)) into equations (4.3a) and (4.3b) leads to the final reformulation of the linearized model

             M 11 ∂ t φ * 1 + M 12 ∂ t w * 1 + gN = 0 , ∂ t N + β * 1 ∂ 2 x φ * 1 -1 -α * 1 ∂ 2 x -ε * 1 h x ∂ x w * 1 = 0 , [M 21 -Q 1 M 11 ] φ * 1 + [M 22 -Q 1 M 12 ] w * 1 + M 23 φ * 2 + M 24 w * 2 = 0 , [M 31 -Q 2 M 11 ] φ * 1 + [M 32 -Q 2 M 12 ] w * 1 + M 33 φ * 2 + M 34 w * 2 = 0 , -Q 3 M 11 φ * 1 -Q 3 M 12 w * 1 + M 43 φ * 2 + M 44 w * 2 = 0 . (4.4)
We now look for solutions of the classical form

           η(x, t) = Ae iθ , θ = ωt -kx , φ * 1 = -i(B 1 + ih x B 2 )e iθ , w * 1 = i(C 1 + ih x C 2 )e iθ , φ * 2 = -i(D 1 + ih x D 2 )e iθ , w * 2 = i(E 1 + ih x E 2 )e iθ , (4.5)
where

A, B 1 , B 2 , C 1 , C 2 , D 1 , D 2 , E 1
, and E 2 are slowly spatially-varying functions (i.e. of the general form F (νx) with ν 1), k is the wavenumber, and ω the wave frequency. The complex conjugate parts of these expressions have been left out for brevity. As stated in [START_REF] Madsen | A new Boussinesq method for fully nonlinear waves from shallow to deep water[END_REF][START_REF] Madsen | Boussinesq-type formulations for fully nonlinear and extremely dispersive water waves: derivation and analysis[END_REF], the B 2 , C 2 , D 2 , and E 2 contributions are necessary because of the bottom slope, and since the velocity potential variables are not in phase with the free surface at the lowest order in h x , but are so at the next order.

(b) Linear dispersion relation

To determine the linear properties of the two-layer model, we substitute the desired form of solutions (4.5) into the linear formulation (4.4) and collect terms at the lowest order in h x . Thus, we obtain a linear system of five homogeneous equations in A, B 1 , C 1 , D 1 , and E 1 . This system has non-trivial solutions if its determinant vanishes, which yields the following dispersion relation

c 2 gh = ω 2 ghk 2 = 1 + a 2 (kh) 2 + a 4 (kh) 4 + a 6 (kh) 6 1 + b 2 (kh) 2 + b 4 (kh) 4 + b 6 (kh) 6 + b 8 (kh) 8 , ( 4.6)
Article submitted to Royal Society where c is the wave celerity and where the (a i ) and (b i ) coefficients are given by

     a 2 = 2S + 1 12 , a 4 = S(2S + 1 12 ) , a 6 = S 3 , b 2 = 2S + 5 12 , b 4 = 3S 2 + 2 3 S + 1 144 , b 6 = S 2 (2S + 5 12 ) , b 8 = S 4 ,
where S = σ(1σ)/12. This dispersion relation is compared in §4 e to the exact linear dispersion relation given by Stokes linear theory, namely

c 2 s gh = ω 2 ghk 2 = tanh kh kh , ( 4.7)
and to its [6, 8] Padé approximation which has the same rational form as (4.6).

(c) Linear vertical profiles

Coming back to the previous linear system in A, B 1 , C 1 , D 1 , and E 1 , we can now express each of the unknowns B 1 , C 1 , D 1 , and E 1 in terms of A, which leads to tedious expressions that are not given here for brevity. For a flat bottom, the expressions of φ i and w i on the whole water column are given by (3.32), (3.33) since

   φ i (t, X, z) = 1 -α ‡ i (z)∆ φ * i + β ‡ i (z) -γ ‡ i (z)∆ w * i , w i (t, X, z) = -β ‡ i (z)∆φ * i + 1 -α ‡ i (z)∆ w * i .
(4.8)

Plugging the ansatz (4.5) for a flat bottom (i.e. without the mild-slope contributions) into the previous expressions and using the computed values of B 1 , C 1 , D 1 , and E 1 in terms of A leads to the expressions of φ 1 (z), w 1 (z) in the upper layer and φ 2 (z), w 2 (z) in the lower layer, in terms of k, h, ω, A, and σ. Finally, we recover the linear vertical profiles over the whole water column using

φ(z) = φ 1 (z) for ẑ ≤ z ≤ 0 , φ 2 (z) for z ≤ z ≤ ẑ , and w(z) = w 1 (z) for ẑ ≤ z ≤ 0 , w 2 (z) for z ≤ z ≤ ẑ . (4.9)
The resulting vertical profiles will be compared in §4 e to the theoretical linear profiles φ s (z) and w s (z) coming from Stokes linear theory, namely

       φ s (z) = Ag ω cosh k(z + h) cosh kh sin(ωt -kx) , w s (z) = Agk ω sinh k(z + h) cosh kh sin(ωt -kx) .
(4.10)

(d) Linear shoaling

We now aim at determining the linear shoaling gradient γ 0 of the double-layer model defined by

A x A = γ 0 h x h . (4.11)
ansatz (4.5), we then collect terms at the next order, i.e. terms proportional to the first derivatives of all the variables. Doing this leads to a new inhomogeneous system of linear equations on the unknowns A x , B 2 , C 2 , D 2 , and E 2 involving the first derivatives of k and h (namely k x and h x ) and the first derivatives of B 1 , C 1 , D 1 , and E 1 . Differentiating the previously computed expressions of B 1 , C 1 , D 1 , and E 1 in terms of A, the derivatives of B 1 , C 1 , D 1 , and E 1 can be expressed only in terms of A x , k x , and h x . Then, differentiating the linear dispersion relation (4.6) allows to express k x in terms of k, h, h x , and σ. Plugging all these relations into the inhomogeneous system on A x , B 2 , C 2 , D 2 , and E 2 , we are able to eliminate all the unknowns but A x and express it in terms of A, h x , kh, and σ, thereby yielding the linear shoaling gradient γ 0 of our double-layer model. The detailed analytic expression for γ 0 is not reported here for brevity.

The computed shoaling gradient will be compared in §4 f to the exact shoaling gradient γ s , which was derived by [START_REF] Madsen | A new form of the Boussinesq equations with improved linear dispersion characteristics. II. A slowly varying bathymetry[END_REF] using energy flux conservation combined with Stokes linear theory, namely

γ s = 2kh sinh 2kh + 2k 2 h 2 (1 -cosh 2kh) (2kh + sinh 2kh) 2 . (4.12) (e) Optimization of

linear properties

The goal is now to optimize the linear properties of our double-layer model by minimizing the errors between these properties and the exact linear properties coming from Stokes linear theory. To this end, we tune the free parameter σ which defines the interface level ẑ = -σh. The different errors between the model linear properties and the theoretical ones are computed as follows. We respectively measure the errors on the phase celerity, the vertical profiles of the velocity potential and the vertical velocity, and the linear shoaling gradient as

E α (K, σ) = K 0 1 kh E 2 α (σ, kh) d(kh) , (4.13)
with α ∈ {c, φ, w, γ}, K being a reference relative water depth, and (4.14) where c, c s , (φ, w), (φ s , w s ), γ 0 , and γ s come respectively from (4.6), (4.7), (4.9), (4.10), (4.11), and (4.12). We point out that in all these errors, the weighting by 1 kh helps keeping the errors to a minimum for low wave numbers, like in [START_REF] Madsen | A new Boussinesq method for fully nonlinear waves from shallow to deep water[END_REF][START_REF] Madsen | Boussinesq-type formulations for fully nonlinear and extremely dispersive water waves: derivation and analysis[END_REF]. Doing this, we sacrifice some accuracy at very high wavenumbers (i.e. in deep water), but we reinforce the model accuracy in shallow water. This weighting by 1 kh is especially well-suited for the shoaling gradient errors, which are far more critical in shallow water than in deep water.

       E 2 c (σ, kh) = c -c s c s 2 , E 2 φ (σ, kh) = 1 h 0 -h φ(z) -φ s (z) φ s (0) 2 dz , E 2 w (σ, kh) = 1 h 0 -h w(z) -w s (z) w s (0) 2 dz , E 2 γ (σ, kh) = (γ 0 -γ s ) 2 ,
Article submitted to Royal Society At this point, we could have minimized each of these errors individually, but doing this leads to quite different optimal values for σ, ranging from 0.2 to 0.5. Furthermore, the minimization of the shoaling gradient error is quite problematic: we will see later that whatever value we choose for σ, the range of validity in kh is limited. We thus choose to minimize E c , E φ , and E w simultaneously to infer the optimal value for σ, and then optimize the shoaling gradient error E γ differently. We start with the minimization of the errors E c , E φ , and E w through the average error

E tot (K, σ) = 1 3 (E c (K, σ) + E φ (K, σ) + E w (K, σ)
), and we compute the optimal value of σ for several typical values of K: the shallow-water value K = π/2, the intermediate depth value K = π, and the deep water values K = 2π and K = 10. In this work, we do not optimize σ for larger values as the vertical profiles have systematically shown an error peak of at least 2% within the range kh ∈ [0, 10] for larger K, for instance K = 15 or K = 20. Table 1 summarizes the optimal values σ opt for each value of K.

Figure 2 plots c/c s to assess the dispersion error on the phase celerity. The upper panel compares the errors obtained for each value of σ opt . The price to pay for extending the linear range of validity towards deep water values is the growth of a small error peak around kh = 3. Indeed, we can see that a 2% error is reached at the very deep water value kh = 24 for σ opt = 0.473 with a very small 0.01% error peak at kh ≈ 3, whereas the same error is reached at kh = 28 for σ opt = 0.314, but with a 0.04% error peak at kh ≈ 3. However, such an error is not significant, and the overall accuracy of the double-layer model for σ opt = 0.314 appears to be excellent up to very deep water. In the same way, the lower panel of figure 2 compares the error on the phase celerity of our double-layer model with σ opt = 0.314 with the errors obtained for the Padé [6, 8] approximation, the model of [START_REF] Jamois | A numerical study of nonlinear wave run-up on a vertical plate[END_REF], and the one of [START_REF] Madsen | A new Boussinesq method for fully nonlinear waves from shallow to deep water[END_REF][START_REF] Madsen | Boussinesq-type formulations for fully nonlinear and extremely dispersive water waves: derivation and analysis[END_REF]. We remark that our double-layer model accuracy is far better in deep water than what is achieved with the Padé [6,8] approximation and the model of [START_REF] Jamois | A numerical study of nonlinear wave run-up on a vertical plate[END_REF]: a 2% error is reached at the very deep water value kh ≈ 28 for the double-layer model, whereas the same error is already reached at kh ≈ 18 for the Padé [6, 8] approximation and kh ≈ 12 for the model of [START_REF] Jamois | A numerical study of nonlinear wave run-up on a vertical plate[END_REF]. In comparison with these results, Lynett & Liu (2004a, b) showed that their double-layer model reaches the same 2% error (not plotted here) at kh = 8. As for the model derived by [START_REF] Madsen | A new Boussinesq method for fully nonlinear waves from shallow to deep water[END_REF][START_REF] Madsen | Boussinesq-type formulations for fully nonlinear and extremely dispersive water waves: derivation and analysis[END_REF], a 2% error is reached at kh = 30, i.e. at a slightly greater value than the one reached by our model with σ = 0.314.

Figure 3 plots the depth-averaged errors E φ (upper panel) and E w (lower panel) on the vertical profiles of the velocity potential and the vertical velocity. We remark that the difference between the errors obtained with each value of σ opt remains very small in shallow water. On the contrary, the benefit for taking σ = 0.314 clearly appears for both the velocity potential and the vertical velocity in deep water: for the vertical velocity profile, a 1% error is reached at kh ≈ 4 and a 2% error at et al. (2006), and the dash-dotted line the model of [START_REF] Madsen | A new Boussinesq method for fully nonlinear waves from shallow to deep water[END_REF][START_REF] Madsen | Boussinesq-type formulations for fully nonlinear and extremely dispersive water waves: derivation and analysis[END_REF].

kh ≈ 8. As far as the velocity potential is concerned, a 1% error is reached at k ≈ 6.5 and a 2% error is reached at kh ≈ 10. By comparison, the model derived by [START_REF] Madsen | A new Boussinesq method for fully nonlinear waves from shallow to deep water[END_REF][START_REF] Madsen | Boussinesq-type formulations for fully nonlinear and extremely dispersive water waves: derivation and analysis[END_REF] yields a 2% error at kh = 12 for both horizontal and vertical velocity profiles. The difference between the errors on the velocity potential and the vertical velocity component can be attributed to the fact that we have neglected the fourth-order derivative term in the expression (4.8) of w i (z). As mentioned earlier, this choice entails to sacrifice some accuracy on the profile of the vertical velocity. Nevertheless, the global accuracy for both vertical profiles is still very good, up to the deep water value kh = 10. This analysis of the phase celerity and velocity profile errors does not exhibit any major advantage for the choice σ = 0.473 instead of σ = 0.314 in shallow water.

The additional errors made with σ = 0.314 are at most of order 0.2% in shallow and intermediate water. On the other hand, the advantage of this value clearly appears in deep water as it appreciably extends the linear range of validity of the model, especially for the vertical profiles. Therefore, we decide to adopt the value σ = 0.314 in the sequel.

Article submitted to Royal Society (f ) Improved model with tightened shoaling properties

We now consider the linear optimization of the shoaling gradient properties. As shown on figure 4, the model linear shoaling gradient γ 0 (dashed line) only matches the exact linear shoaling gradient γ s (solid line) up to kh = 3 and swiftly departs from it beyond that value. We have attempted to optimize this shoaling gradient individually, but no value of σ makes the two curves fit beyond kh = 3. Therefore, a more subtle optimization is needed. This can be achieved using the following method. Going back to the full formulation of the approximate static Dirichlet-Neumann operator G app 0 [h] (3.23), we introduce a new constant parameter r and apply the differential operator 1 + rh∇h • ∇ to the last equation of (3.23), which yields (4.15) where

P * 0 w 0 m ≈ -β * 1 ∆ -rhβ * 1 ∇h • ∇∆ φ * 1 + 1 -α * 1 ∆ -(ε * 1 -rh)∇h • ∇ -rhα * 1 ∇h • ∇∆ w * 1 ,
P * 0 = 1 + σ 2 β * 1 + rh ∇h • ∇.

Article submitted to Royal Society

This new formulation does not allow further optimization yet, since the additional terms cancel out in the shoaling analysis. However, an interesting option is to neglect the third-order derivative on w * 1 in the slope terms. By doing this, we greatly improve the shoaling gradient properties, as we will see in the sequel. For the moment, we rewrite (4.15) as follows:

P * 0 w 0 = -β * 1 ∆ -rhβ * 1 ∇h • ∇∆ φ * 1 + 1 -α * 1 ∆ -(ε * 1 -rh)∇h • ∇ w * 1 . (4.16)
We then take advantage of (3.22) to eliminate the third-order derivative from the previous equation, and obtain

P * 0 w 0 = -β * 1 ∆ - 6r σ ∇h • ∇ φ * 1 + 1 -α * 1 ∆ -(ε * 1 + 2rh)∇h • ∇ w * 1 + 6r σ ∇h • ∇ φ 0 .
(4.17)

This implies to redefine the approximate operator G app 0 [h] as follows:

G app 0 [h] = Q * 0     N 1 0 0 0 0 N 2 0 0 0 0 N 3 0 0 0 0 0     R     P 0 Q 1 Q 2 Q 3     . (4.18) where    N 1 = -β * 1 ∆ - 6r σ ∇h • ∇ , N 2 = 1 -α * 1 ∆ -(ε * 1 + 2rh)∇h • ∇ , N 3 = 6r σ ∇h • ∇ , Q * 0 = 1 - σ 2 β * 1 + rh ∇h • ∇ m ≈ (P * 0 ) -1 . (4.19)
Our modified model is thus identical to (2.10) but with G app 0

[h] as redefined above. Starting from it, the new linearized model remains essentially the same, except that the second equation of (4.4) now is

∂ t N * + β * 1 ∂ 2 x + 6r σ h x ∂ x φ * 1 -1 -α * 1 ∂ 2 x -(ε * 1 + 2rh)h x ∂ x w * 1 - 6r σ h x ∂ x φ 0 m ≈ 0 ,
with N * = P * 0 η. This new formulation does not modify the phase celerity and the vertical profiles, but allows to further minimize the shoaling gradient error. We can now compute the new shoaling gradient γ 0 . Optimizing the parameter r so that the error E γ (K, σ) is minimized for K = 10 and σ = 0.314 yields r = 0.0076. Remark 4.1. A dimensional analysis shows that this small value of r is coherent with our earlier choice to neglect the third-order derivatives of w * 1 in (4.16).

Figure 4 displays the optimized shoaling gradient for the previous value for r. The improvement is quite impressive since the new shoaling gradient exhibits a very good agreement up to kh ≈ 12. This accuracy is the same as that reached by the model of [START_REF] Madsen | A new Boussinesq method for fully nonlinear waves from shallow to deep water[END_REF][START_REF] Madsen | Boussinesq-type formulations for fully nonlinear and extremely dispersive water waves: derivation and analysis[END_REF].

To sum up, our final double-layer Boussinesq-type model consists of (2.10) and (4. properties: the phase celerity is accurate up to kh = 28, the vertical velocity profiles are accurate up to kh = 10 for the velocity potential and up to kh = 8 for the vertical velocity component, and the shoaling gradient is accurate up to kh = 12. We emphasize that these properties are not affected by a slight variation of σ, which makes the model robust towards the parameter σ. These results are quite similar to those obtained by [START_REF] Madsen | A new Boussinesq method for fully nonlinear waves from shallow to deep water[END_REF][START_REF] Madsen | Boussinesq-type formulations for fully nonlinear and extremely dispersive water waves: derivation and analysis[END_REF], but the main advantage of the present model is that it contains lower-order derivatives and fewer equations, especially in 2DH. The present double-layer approach is hence a very good alternative to the most advanced high-order Boussinesq models as it offers almost the same linear properties with a lower complexity.

Numerical simulations: nonlinear behaviour

On the basis of the model (2.10) derived in §2 (we do not use here the version derived in §4 f since we consider a flat bottom), a classical finite-difference scheme is developed to study numerically some nonlinear properties of the model in 1DH.

We consider the propagation of two-dimensional periodic and regular nonlinear waves, without change of form, over a flat bottom and without any ambient flow field. For this situation, numerical reference solutions can be obtained by the so-called stream function method, or more precisely the Fourier approximation of the stream function [START_REF] Dean | Stream function representation of nonlinear ocean waves[END_REF][START_REF] Rienecker | A Fourier approximation for steady water waves[END_REF]. Unlike analytical wave theories (such as Stokes or cnoidal wave theories), this numerical approach is applicable whatever the relative water depth and steepness are, and very accurate solutions can be obtained by increasing the number of terms in the Fourier series (e.g. 10, 20, or 50 if necessary for very steep waves). This method was previously implemented in a software called Stream_HT by one of the authors [START_REF] Benoit | Near-bottom kinematics of shoaling and breaking waves: experimental investigation and numerical prediction[END_REF]. For the selected application, the domain of interest covers one wave-length (L = 64 m ; k = 2π/L = 0.098 rad m -1 ) and periodic conditions are imposed at the two lateral boundaries. We consider a still water depth of h = 96 m, so that the relative water depth is h/L = 3/2 yielding kh = 3π ≈ 9.425, which corresponds to deep water conditions. The wave height is chosen as H = 6.4 m, so that the Article submitted to Royal Society steepness is H/L = 0.1 or kH/2 = π/10 ≈ 0.314, i.e. about 70% of the theoretical maximum value of the steepness for a stable wave [START_REF] Williams | Limiting gravity waves in water of finite depth[END_REF]). These conditions correspond to highly dispersive and very nonlinear waves. For this case, the period computed with the stream function approach (at order 20) is T = 6.094 s, yielding a wave celerity of C = L/T ≈ 10.502 m s -1 . The solution obtained with the stream function approach for the free surface elevation and the free surface potential is imposed as the initial condition in our simulations.

The time integration scheme is a classical fourth-order four-stage explicit Runge-Kutta scheme, which is known to possess a wide stability region. However, owing to the nonlinear nature of the considered test case, this scheme can develop some high-frequency instabilities. To avoid such instabilities, an 8th-order Savitsky-Golay smoothing filter is applied twice after each time step to η, φ 1 , and w 1 . The price to pay for this filter is a negligible loss of accuracy of the model. As far as spatial discretization is concerned, all derivative operators in (2.10) are replaced by centered fourth-order finite difference approximations combined with periodic boundary conditions. The covered domain (of one wave-length) is discretized with 32 cells of constant size (equal to 2 m) and a time-step of 0.122 s (corresponding to T /50) is used during the simulations. We have verified that using a refined mesh of 128 cells and 200 time steps per period does not yield any significant improvement.

Numerical integration of the double-layer model (2.10) is performed over a duration of 25 T . Simulations have been performed with the deep-water (kh = 10) optimal value σ = 0.314, even if any value in [0.28, 0.36] would lead to very similar results. Outside this range, the results quickly deteriorate. Results obtained after durations of 10 T and 25 T are plotted on figure 5 for the free surface elevation and the free surface velocity potential. The results are compared to the reference solution obtained with the stream function approach (which propagates at constant Article submitted to Royal Society speed and without change of form). The results appear to be very good since the two curves fit very well until t ≈ 20 T , after which small discrepancies become observable. Since grid convergence has been verified, this difference can be attributed to the approximation in (2.9), where we have neglected fourth-order (and higher) nonlinear terms. However, we remark that the global aspect of the model curve still corresponds to that of the reference solution at t = 25 T . Furthermore, the free surface elevation computed with our model only shows a phase shift error with the reference solution: the forms are the same and the amplitude of the waves are equal. An interesting remark is that we can use these results to compute the nonlinear phase celerity error approximatively. Indeed, taking for instance the free surface elevation results and measuring the distance between the crests of the two curves yields an approximate value of the difference of celerity between these curves. This value provides us a measure of the nonlinear celerity error of the model. We found that our model with σ opt = 0.314 exhibits a nonlinear phase celerity error of about 0.08%, which is an impressive result. To conclude, the model shows an excellent nonlinear behaviour, and we can expect its nonlinear range of validity to reach up to kh = 10, at least for flat bottom conditions.

  (b) A translated Dirichlet-Neumann operator The Dirichlet-Neumann operator associated to the problem (2.5c)-(2.5g) is denoted by G[η, h] and is defined by G[η, h]ψ = ∂ z φ 1 | z=η for any smooth enough function ψ, where (φ 1 , φ 2 ) solves the boundary value problem composed of equations (2.5c)-(2.5g) along with the Dirichlet condition φ 1 = ψ at z = η. One can thus simply write the closure between the unknowns w 1 , φ 1 , and η as

Figure 2 .

 2 Figure 2. Comparison of linear phase celerity with the exact Stokes result. The top figures compare the errors obtained for our model with σ = 0.314 (solid line), σ = 0.365 (dashed line), σ = 0.428 (dash-dotted line), and σ = 0.473 (dotted line). The top right figure is a zoom on the region kh ∈ [0, 12]. The bottom figure compares the errors for different models: the solid line represents our double-layer model with σ = 0.314, the dotted line the Padé [6, 8] approximation, the dashed line the model of[START_REF] Jamois | A numerical study of nonlinear wave run-up on a vertical plate[END_REF], and the dash-dotted line the model of[START_REF] Madsen | A new Boussinesq method for fully nonlinear waves from shallow to deep water[END_REF][START_REF] Madsen | Boussinesq-type formulations for fully nonlinear and extremely dispersive water waves: derivation and analysis[END_REF].

Figure 3 .

 3 Figure 3. Depth-averaged errors on the vertical profile of the velocity potential (top) and vertical velocity (bottom) for σ = 0.314 (solid line), σ = 0.365 (dashed line), σ = 0.428 (dash-dotted line), and σ = 0.473 (dotted line).

Figure 4 .

 4 Figure 4. Shoaling gradient: model with (r = 0.0076; dotted line) and without (r = 0; dashed line) optimization with σ = 0.314, and exact shoaling gradient (solid line).

  Figure5. Snapshots of the free surface elevation (top) and the free surface velocity potential (bottom) at t = 10 T (left) and t = 25 T (right). The solid line is our model for σ = 0.314 and the dotted line is the reference solution computed using the stream function approach.

  .27)Article submitted to Royal Societywhere Id is the identity operator, is thus an approximate inverse of the operator P + ∇h • Q up to O(|∇h| 2 ) terms. Hence,

	U	m ≈ RF ,	(3.28)
	which yields the explicit expressions of φ * 1 , w * 1 , φ * 2 , and w * 2 in terms of φ 0 .
	The very last step consists in introducing the operators

Table 1 .

 1 Optimal values for σ

	K	π/2	π	2π	10
	σopt	0.473	0.428	0.365	0.314

Article submitted to Royal Society


Gwenael Chaillou
email: gwenael.chaillou@lirmm.fr

Philippe Maurine
email: philippe.maurine@lirmm.fr

Nadine Azemard
email: nadine.azemard@lirmm.fr

Jean-Marc Galliere
email: jean-marc.galliere@lirmm.fr

Iterative Method for Performance Prediction Improvement of Integrated Circuits

Keywords: Design flow, IR Drop, libraries, performances

Currently, the design flow of digital integrated circuits relies on the use of pre-characterized libraries. These libraries, pre-calculated outside of any context, are used at each step of the design flow to predict the performance of the circuit. They are necessary for the calculation of a wide variety of characteristics, including the IR Drop. However, the behavior of the logic gates described in these libraries is not representative of the real power supply of the standard-cells and this difference is even more important as the number of logic gates in the circuits increases. In this paper, we will show that with the current design flow, the use of these libraries leads to an underestimation of the circuit performances. To remedy this, we propose an iterative approach based on a new use of voltageadaptive libraries that allows to compute the IR Drop more precisely and thus to better predict the performances of the circuits.

I. INTRODUCTION

Nowadays, for designing a digital integrated circuit, the expert must verify the correct operation of the circuit at each step of the design flow and with regard to the different parameters. The main external parameters are called PVT and refer to the variations of manufacturing Process (P), supply Voltage (V) and Temperature (T). Ever-increasing number of components in circuits, manufacturing transistor fluctuations, power supply and temperature, which are stochastic phenomenons, move towards an increasing impact on integrated circuit performances especially in advanced semiconductor technology nodes [1].

These parameters are included in the design flow through the use of pre-characterized libraries. This is called the corner method. This conventional design method consists in using simple best-and worst-case combinations of the PVT parameters to predict the circuit performances. Indeed, with the help of CAD tools, the designer implements the application by ensuring, at each step of the flow, the achievement of the specifications but also the correct circuit operation regarding both transistor manufacturing fluctuations and circuit environment (power supply, temperature). However, this worst-case modelling is considered unrealistic and tends to be rejected [START_REF] Rout | Process corner variation aware design of low power current starved VCO power[END_REF]. The necessity to adopt a better approach is demonstrated in [1].

Moreover, the logic gates in a circuit are not supplied directly by the external power supply (VDD), know generally at ±2%, but by the same source delivered through a power grid. This power grid is resistive, and the resistivity increases significantly as the size of the technologies is reduced [3]. As the process reduced from 28nm to 7nm, conductor resistance through the power grids has increased by factor 10. This resistivity induces voltage Drop and, as a result, each logic gate has a proper local supply voltage (Vi) lower than the one available at the regulator's output. It is well known that power grids can cause timing errors. For this purpose, IR Drop analysis [START_REF] Nithin | Dynamic voltage (IR) drop analysis and design closure: Issues and challenges[END_REF] helps the designer to counteract the negative effects of the latter. IR Drop is included in common design closure by a fixed limit of 10% of external supply. In our case, IR Drop estimations are made with Voltus IC Power Integrity Solution [5]. As a rule of the thumb, if the IR Drop of a circuit is over 10%, it is admitted that its impact on the delays will be too important and therefore the circuit, as it stands, is not validated. This margin is fixed in this way, and its application has shown that it is intended to be pessimistic. The value of IR Drop is calculated from libraries, but we will see in this article that there is a misuse of those libraries which leads to an overestimation of the IR Drop.

The libraries describe the behavior of each standard-cell of the used technology. These files describe mainly the consumption and the delays of the standard-cells and are classified according to the PVT parameters. Then, using these libraries impose to work with a fixed supply voltage, temperature and manufacturing process for all the gates in the circuit. This procedure is not suitable for working with unknown and variable Vi. Therefore, we propose in this paper a new approach which assign to each gate a corresponding part of library. For this study, we focused on variations of voltage supply, manufacturing process and IR Drop. Temperature is considered constant.

The remainder of this paper is organized as follows. In sect. II, we will give more details about libraries and IR Drop in actual design flow. Then, we propose a new iterative approach in sect. III. The results will be described in sect. IV. Conclusion and perspectives of this article are provided in sect. V.

II. LIBRARIES AND IR DROP IN DESIGN FLOW

In this section we will explain the use of libraries and IR Drop in the current design flow. We will also demonstrate why we consider libraries to be inaccurate and how their use leads to an over-estimation of the IR Drop.

A. Libraries use in the current design flow

The logic gates of the libraries are pre-characterized a priori, i.e., before their integration within a logic path. The libraries are composed of look-up tables. These look-up tables allow to describe the behavior of the logic gates according to the parameters of the library (PVT) and also the input net transition and output net capacitance. At each step of the design flow, the designer uses these look-up tables to ensure the achievement of the specifications and the correct circuit operation (see Fig. 1). The tools used by designers are becoming more and more precise and complex in their calculations. This is in order to have results that are closer to reality. But whatever the level of calculation, the data used are those found in the look-up tables. However, for calculations of more complex phenomena, we notice that these data are not adapted anymore.

As seen on Fig. 1, the IR Drop is calculated at almost the end of the design flow. The designer uses the library corresponding to the corner with the worst IR drop, i.e., high voltage, best process, and high temperature to estimate the IR Drop in the circuit. The result is then compared to the fixed margin of 10% to determine if the circuit is validated or not. At this stage, the supply grid has been considered. Each logic gate actually sees a proper supply voltage Vi provided by the power grid and depending on the IR Drop generated in its area. However, the library used to deal with this the IR Drop, at gate-level, is still the same and characterized for only one value of VDD (high). Moreover, only one library is used for all the gates in the circuit. This does not correspond to the fact that Vi varies over the entire surface of the circuit, and that each standard-cell has its own power supply Vi. In this paper, it will be shown that it is possible to use several libraries for the same simulation with a single circuit.

B. Libraries misuse

The use of libraries in the usual design flow is not adapted to the variations of IR Drop. Indeed, the designer calculates the IR Drop using a library that describes the behavior of the logic gates in the circuit for a given power supply. However, as an activity necessarily induces IR Drop, it is impossible for a logic gate to receive all the expected power supply. Its real performances will be lower than those described in the library. In other words, it is impossible to precharacterize the performance of a logic gate as a function of VDD. This is because it depends on the IR Drop, generated by the gate itself and also by the standard cells surrounding it. Then, as the IR Drop is overestimated in the actual design flow, a more realistic estimation of the IR Drop in the circuit will give a lower value. This concept is illustrated on Fig. 2. As it is shown, with precharaterized standard-cell libraries, the IR Drop is overestimated compared to reality. With this overestimation, the 10% rule is troublesome. Indeed, it shows that a circuit could appear to be below the 10% margin when in reality it is not the case. This prediction error is also linked to the process. In fast process, the logic gates commute faster and consume more power. So, the IR Drop generated will be higher. And the error between the theoretical supply voltage used in the libraries and the actual supply voltage will also be greater.

Then, the overestimation of IR Drop calculated in the usual way will also be higher. However, it is precisely in this case (fast process) that the IR Drop is calculated with the actual design flow to be compared to the 10% limit. Finally, the risk is again to have to rework unnecessarily the design.

By using libraries adapted to the context of each standardcell, we could therefore gain in precision on the estimation of the IR Drop and obtain a more realistic IR Drop value. This decrease of the IR Drop leads to a gain on the temporal performances of the circuit and could allow to improve the efficiency, the power consumption or the surface of the circuit. The objective of the second part is to describe the method that will allow to estimate more precisely the IR Drop.

III. ITERATIVE APPROACH

In this part, a new approach to the use of libraries is proposed as a solution to the problems outlined above. This method is iterative. Indeed, it consists in repeating the classical design flow several times. We will describe the modifications to be made to the flow and how to link these cycles together.

A. Working environment

The first step is to design the circuit with CAD tools. The most important point is the library file. It must contain the libraries corresponding to the usual corners. For example, the library corresponding to the worst-case corner used for setup timing calculations is the low voltage, worst process and high temperature library: wc_LV_HT.lib. The difference with the current design flow (fig. 1) is that the libraries must be available on several power supply values. The purpose is to have libraries that can cover all possible values of Vi, the value of the supply voltage delivered to the standard cell from the power grid, i.e., the value after deducting the IR Drop. The power supply values of the libraries must go from VDD-Dropmax to VDD-Dropmin with a step that will depend on the needs. The work in this paper shows that a step of about 10 mV seems to be consistent. In order to replace permute active libraries during the design flow, the logic gates in the different libraries must not have the same name; the names of the logic gates of a library must be specific to this library. In this way, two different libraries can be assigned to two physically identical gates. For instance, NAND2X4 can be declined to NAND2X4_0.9 and NAND2X4_0.8 if they belong to libraries characterized for 0.9V and 0.8V respectively .

The second step is to correctly set up these logic gates and libraries in the initialization files. It must contain all the names of the standard-cells for all Vi values, but their physical LEF description remains the same. The corners also need to be assigned to the corresponding libraries. Once the environment is set up, the circuit can be placed and routed the usual way with design tools. In this paper, the method is applied with Genus [6] and Innovus [7].

B. Method

Once all the libraries are ready for use and the circuit is placed and routed, the IR drop is calculated using Voltus. So, Fig. 3 shows the result of an IR Drop analysis. Each square represents a standard cell, so this figure shows the IR Drop at the supply point of each cell. This first analysis is done with the usual way. For this example, the IR drop varies from 88mV to 99mV. These results must then be retrieved in a log file. If the method is applied on Voltus, the report_power_rail_result command can be used to obtain a file giving the IR Drop value for each instance in the circuit. The next step is to process this file. In this paper, the process is done with Python. The objective is to create an Innovus tcl command file that will replace the standard cells of each instance by those corresponding to their new supply voltage value after deduction of the IR Drop. This step, as well as all the steps of the method, can be seen in the Fig. 4. Processing the IR Drop report in Python allows to both list all the standard-cells and assign them their new power supply value, which is closer to reality.

For instance, the IR Drop analysis indicates that a NAND2X4 logic gate instantiated as g1516 has an IR Drop value of 98.1 mV. The supply value of the new library associated with g1516 will therefore be 801.9mV (VDD -98.1 = 801.9). If there is a step of 5 mV between the libraries, the library used will be the one with the closest value. Then, we will find the following command in the Innovus tcl command file:

-cell NAND2X4_0.800 -inst g1516

Once the file is fully written, it can be sourced in Innovus to reassign all standard-cells. Fig. 5 shows an exemple of libraries switching on the whole surface of a circuit. The libraries are now reassigned and the design flow can be followed as if it were a new circuit. The IR drop of the circuit is then re-calculated from the new libraries. For this second cycle, the values of IR Drop decrease significantly. Indeed, the libraries are characterized at the real supply voltages of the circuit, whereas for the first cycle the supply voltages were ideal. During the second cycle the standards will have less activity and will generate less IR Drop. With this second IR Drop result, the library assignment operation can be repeated to start a new cycle. It is then necessary to perform as many cycles as necessary to stabilize the IR Drop value. The work presented in this paper shows that about 3 cycles are needed to stabilize the IR Drop value. However, the exact number of cycles always depends on the context and the desired accuracy.

C. Delays estimations in the method

In the usual design flow, delays are calculated after placement and routing and before taking the IR Drop into account. These delays are inaccurate because they are calculated for an ideal power supply. We can say that these delays are not IR-aware. The link between the delays and the IR Drop is made only through the 10% margin, despite the inaccuracies induced. However, with the presented method, the delays are re-calculated after the last assignment of libraries. As the libraries are re-assigned based on the IR Drop, the obtained delays are now IR-aware. This method allows us to obtain a more precise idea of the values of the delays. Moreover, if these delays are IR-aware, so the 10% IR Drop margin becomes obsolete. The removal of this pessimistic margin could lead to a gain in performance.

In summary, this iterative method (Fig. 4) requires the availability of libraries covering a wider range of supply values. Indeed, it consists in performing several IR Drop analysis by reassigning the library corresponding to the correct power value to each instance of the circuit between each analysis. A priori, 3 IR Drop analysis are necessary. Performing these 3 iterations allows to obtain a stable and more accurate IR Drop value. Knowing the IR Drop of the circuit and being able to integrate it with the use of the libraries allows to calculate the IR-aware delays.

IV. APPLICATION AND RESULTS

In this section, we apply the iterative method to a simple circuit and compare the results with actual design flow. A 32bit-multiplier is taken as an example (see Fig. 6). This multiplier is made from 65nm technology and composed by 7k logic gates. Its area is 1 mm² and it has a clock period of 3600ps. This multiplier has been designed so that its IR Drop is relatively large and around the 10% margin, in order to be able to correctly illustrate the remarks. The results are specific to this circuit, but it can be extended to any type of circuit.

A. Evolution of IR Drop

The IR Drop is calculated for each standard-cell at each iteration as described in the method. It is therefore possible to plot the Vi seen by a standard-cell as a function of the iterations (see Fig. 7). Iteration 0 is equivalent to the current design flow. The value of Vi at iteration 0 is 992mV. The value of VDD used for IR Drop analysis is 1.1V. It's corresponding to the high voltage for this technology. So, the IR Drop at iteration 0 is 108mV or 9.81% of the external supply voltage. With the classical design flow, the circuit operation seems to be bordeline, but acceptable. It doesn't exceed the limit of 10% IR Drop. On the second cycle, the IR Drop significantly decreased, as expected. With the iterative approach, Vi is stabilized at 997mV, i.e. an IR Drop of 103mV. IR Drop has thus decreased by 4.6% of its value, and it now reaches 9.3% of the external supply value. The value of IR Drop has decreased slightly. It is still quite close to the 10% limit.

B. Evolution of delay estimations

One of the advantages of the iterative method is to be able to calculate IR-aware delays. Thus, the uncertainty about the 10% rule will be eliminated. Indeed, it is very interesting to remove this uncertainty in the case of the multiplier because the value of IR Drop is very close to the limit. As for the IR Drop, the delays at iteration 0 are the same as the delays calculated with the classical design flow. As a reminder, these delays are not IR-aware. With the iterative method, it is possible to recalculate them at each iteration, i.e., at each time we change the libraries. These delays are then calculated according to the IR Drop value of the previous iteration. Finally, as the IR Drop is constant on the last iteration, the delays of the circuit are also stabilized. Fig. 8 compares the final delays (last iteration) to the initial delays of the logical paths for the multiplier. The setup critical path delay has been slowed down by 16%. The setup critical path delay is used as operating condition of the circuit as it represents the longest path. If this delay is correct, all the delays are validated. This increase in delays was expected as this delay is IR-aware, unlike to the initial delay. This increase in delay goes hand in hand with a 5% decrease in overall circuit power. However, the clock period being 3.6ns and the critical path delay being 3.9ns the timing slack is negative. The timing slack being the time remaining until the next clock edge, if this one is negative, then the results do not make it possible to validate the correct operation of the circuit.

Moreover, these delays being IR-aware, they become much more accurate and closer to the real values. This induces news reflections. In case of positive critical slack time, it would be possible to use the critical path and its slack time to accurately determine the delay margin of the circuit and to use it to reduce the power consumption or the size of the circuit. As a reminder, the slack time of a logical path is the time remaining between its arrival time and the clock edge. By precisely knowing this time it would be possible to estimate the number of gates that could be reduced in size.

As expected, Fig. 8 also shows that the longer the paths, the greater the impact on delays. We could go further in our thinking. Reducing the size of logic gates means a gain in circuit area but also in power consumption. These gates become slower, but they will generate even less IR Drop. Repeat the iterative analysis of the IR Drop on the circuit with these new gates could be interesting. Obviously, our example of multiplier does not have enough gates for such an optimization to be relevant. However, this reflection leads to interesting optimization perspectives.

In summary, the results of the iterative method show that the performance of this 32bit-multiplier circuit, which didn't seem guaranteed, is actually correct. Moreover, the results on the IR-aware delays open doors for reflection on the optimization of the circuit in terms of area and power consumption. The classical design flow has its limits with the miniaturization of the technologies used in the design of integrated circuits, in particular to manage IR Drop issues. Indeed, there is a double pessimism in this flow :

-the fixed margin of 10% of IR Drop which is not necessary and -the misuse of libraries which leads to an over-estimation of the IR Drop.

In this paper, we have presented a new iterative method to overcome these problems. This method is based on a use of libraries adapted in supply voltage. It consists in conducting several consecutive IR Drop analyses by re-adapting the supply voltage of the libraries used at each iteration and assign a different library to each standard-cell.

By applying this method to the 32bit-multiplier circuit presented, we obtain a better estimation of the IR Drop. The circuit, which had a global IR Drop slightly under 10%, could be considered as operating correctly. However the estimation of real IR Drop leads to IR-aware delays calculation, and IRaware delays showed that the operation of the circuit was not guaranteed.

Computing the IR-aware delays allows to confirm the good operation of the circuit and open some perspectives. Indeed, having a precise margin on the delays could allow to optimize the circuit in terms of area and power consumption. This is one of the main avenues for future works. This method should also be compared to other tools that are emerging such as Tempus PI [3] from Cadence, released in November 2020.

Fig. 1 .

 1 Fig. 1. The different stages of the current design flow.

Fig. 2 .

 2 Fig. 2. Impact of the misuse in libraries on IR Drop calculation

Fig. 3 .Fig. 4 .

 34 Fig. 3. Result of the IR Drop analysis made by Voltus on a test circuit

Fig. 5 .

 5 Fig. 5. Supply voltage of the libraries assigned to each standard-cell on the whole circuit for the first iteration

Fig. 6 .

 6 Fig. 6. Evolution of Vi calculated at each iteration seen by a standar-cell placed in the middle of the circuit

Fig. 8 .

 8 Fig. 8. Final and IR-aware delays obtained by the iterative method, compared to the initial delays