## Abstract

We derive and analyse, in the framework of the mild-slope approximation, a new double-layer Boussinesq-type model that 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.* (Agnon *et al.* 1999*J. Fluid Mech.* **399**, 319–333) and enhanced by Madsen *et al.* (Madsen *et al.* 2003*Proc. R. Soc. Lond. A* **459**, 1075–1104), which consists of constructing infinite-series Taylor solutions to the Laplace equation, to truncate them at a finite order and to use Padé approximants, and the double-layer approach of Lynett & Liu (Lynett & Liu 2004*a**Proc. R. Soc. Lond. A* **460**, 2637–2669), which allows lowering the order of derivatives. We formulate the model in terms of a static Dirichlet–Neumann operator translated from the free surface to the still-water level, and we derive an approximate inverse of this operator that can be built once and for all. The final model consists of only four equations both in one and two horizontal dimensions, and includes only second-order derivatives, which is a major improvement in comparison with so-called high-order Boussinesq models. A linear analysis of the model is performed, and its properties are optimized using a free parameter determining the position of the interface between the two layers. Excellent dispersion and shoaling properties are obtained, allowing the model to be applied up to the deep-water value *k**h*=10. Finally, numerical simulations are performed to quantify the nonlinear behaviour of the model, and the results exhibit a nonlinear range of validity reaching at least *k**h*=3π.

## Introduction

During the past two decades, the original Boussinesq (1872) model for a flat sea-bottom and its uneven sea-bottom version derived by Peregrine (1967) 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 the problem size to be reduced 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 *k**h*, where *k* is the wavenumber and *h* the water depth.

The conventional Boussinesq model for uneven bottoms (Peregrine 1967), 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 *k**h*<0.75, as stated in Madsen *et al.* (2002, 2003), 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 extending the linear and nonlinear range of applicability of Boussinesq-type models. The first historical improvement, proposed by Nwogu (1993), consists of 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 (Bona *et al.* 2002, 2005; Chazel 2007). An effort similar to the one of Nwogu (1993) was pursued by Madsen & Sørensen (1992), which was followed by the works of Liu (1994) and Wei *et al.* (1995), in which the authors efficiently removed the weak nonlinearity assumption, allowing the model to simulate wave propagation with strong nonlinear interactions. According to the reviews proposed by Madsen & Schäffer (1998, 1999), these new Boussinesq-type models allow extending the linear range of applicability to *k**h*=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 Gobbi *et al.* (2000), which uses a fourth-order polynomial approximation: this model exhibits excellent linear properties up to *k**h*=6 for the dispersion relation and up to *k**h*=4 for the vertical profiles of orbital velocities, whereas nonlinear behaviour is fairly well captured up to *k**h*=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 the range of applicability of the model of Gobbi *et al.* (2000) into very deep water without increasing the numerical complexity too much, 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 (Agnon *et al.* 1999). In this work, the authors present a new procedure by which it is 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 Agnon *et al.* (1999) involves velocity variables evaluated at the still-water level and is limited to *k**h*=6, the extended method proposed by Madsen *et al.* (2002, 2003) completely removes the conventional shallow-water limitation, allowing for modelling fully nonlinear and highly dispersive waves up to *k**h*≈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 (Madsen *et al.* 2006). Although a few numerical approaches based on this model have been proposed (Fuhrman & Bingham 2004; Engsig-Karup *et al.* 2006, 2008), the counterpart for this wide range of applicability is the numerical complexity of the underlying model, which includes derivatives up to the fifth order, and consists of more equations (and more unknowns) than the alternative model of Gobbi *et al.* (2000). An interesting alternative approach has been chosen by Jamois *et al.* (2006), where the authors propose using a velocity potential formulation and truncation of 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 *k**h*=10, but entails a much lighter numerical complexity with less equations and with derivatives up to the fourth-order only.

The second approach that has been studied is the double-layer approach, as proposed, among others, by Lynett & Liu (2004*a*,*b*) and Audusse (2005). 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 (2004*a*,*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 wave propagations up to *k**h*=6 to be modelled accurately, 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 Madsen *et al.* (2002, 2003) and the one of Lynett & Liu (2004*a*,*b*). The primary goal of this paper is to offer an efficient alternative to the models of Madsen *et al.* (2002, 2003) by mixing their procedure, the simplifications of Jamois *et al.* (2006) and the double-layer concept of Lynett & Liu (2004*a*,*b*). We aim at deriving a model that is (i) applicable to complex domains, such as coastal areas, islands or estuaries, and (ii) 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 *k**h*=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 (i.e. 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 relations, can be found in the work of Dommermuth & Yue (1987), Craig & Sulem (1993), Smith (1998), and in the work of Matsuno (1993) and its extension to very general bathymetries by Artiles & Nachbin (2004*a*,*b*). However, the application of the above approaches to complex 2DH domains has not been reported yet. Thus, with an eye towards coastal engineering applications, we prefer to base our approach on approximating the static Dirichlet–Neumann operator. Furthermore, we mention the new promising approach of Ablowitz *et al.* (2006) 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 to reduce 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

### Governing equations and boundary conditions

We aim at formulating a double-layer Boussinesq-type model for the three-dimensional 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 owing 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 the horizontal coordinates by *X*=(*x*,*y*) and the vertical one by *z*, with 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 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 , where σ is an arbitrary parameter in ]0,1[. Thus, the thickness of the two layers are constant fractions of the still-water 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 and . 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 still-water depth *h* verifies |∇*h*|≪1, which corresponds to the classical mild-slope approximation. This approximation consists of neglecting all the quadratic (and higher order) terms in ∇*h*, along with the derivatives of *h* of order greater than 1. 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 (2004*a*,*b*). For clarity, all the equations derived using this mild-slope approximation are indicated in this work by the symbol 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 ∇ is 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.1a
2.1a
where *P*_{i} is the reduced pressure field in each layer *g* is the gravitational acceleration and *R* is 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 to be 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:
2.2a
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 between the layers, the following natural continuity conditions are enforced:
2.3
Observe that ϕ_{1}=ϕ_{2} and *w*_{1}=*w*_{ 2} at imply ∇ϕ_{1}=∇ϕ_{2} and thus *u*_{1}=*u*_{ 2} at , hence recovering the continuity conditions on the horizontal and vertical velocities enforced by Lynett & Liu (2004*a*).

Finally, the system of equations is closed by the usual kinematic condition at the sea bottom ,
2.4
which expresses that the sea bottom is a bounding surface. We now introduce
Following Zakharov (1968), we can reformulate the equations as
2.5a
2.5b
where (2.5a) is the Euler equation expressed at the free surface and (2.5b) is the kinematic condition at the free surface,
2.5c
2.5d
where (2.5c) and (2.5d) are the Laplace equations in each layer, and
2.5e
2.5f
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 mild-slope hypothesis to neglect the |∇*h*|^{2} term. We point out that we work with a velocity potential formulation, unlike Madsen *et al.* (2002, 2003) 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 us to consider two less equations than with a velocity formulation.

In the following, 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 vertical velocity at the free surface in terms of the velocity potential 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.

### A translated Dirichlet–Neumann operator

The Dirichlet–Neumann operator associated with the problems (2.5c)–(2.5g) is denoted by and is defined by 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 , and η as
2.6
This Dirichlet–Neumann operator is at the heart of the derivation of Boussinesq-type 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 from the equations by plugging (2.6) into (2.5a) and (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 Madsen *et al.* (2002, 2003) and consists of keeping in the equations, which entails solving (2.5a) and (2.5b) and then computing through the use of the Dirichlet–Neumann operator at each time step. This is the method we have chosen to use to keep the equation’s 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 at *z*=η, on a time-dependent domain bounded from above by the free surface *z*=η. Keeping in the equations involves constructing an approximation to at each time step, which can be a serious computational problem as it increases the numerical cost. An interesting solution to this issue consists of 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. Therefore, we introduce ϕ_{0}(*t*,*X*)=ϕ_{1}(*t*,*X*,*z*=0) and *w*_{0}(*t*,*X*)=*w*_{1}(*t*,*X*,*z*=0) and denote by the Dirichlet–Neumann operator corresponding to η=0, i.e.
2.7
The translated operator is such that , where (ϕ_{1},ϕ_{2}) solves the boundary value problem
2.8
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 .

### 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 us to artificially define ϕ_{0} and *w*_{0} when η<0. This tool has been implicitly used by many authors, such as Nwogu (1993), Wei *et al.* (1995), Gobbi *et al.* (2000) and Madsen *et al.* (2002, 2003), who derived models based on a horizontal velocity (or potential) variable taken at free surface-independent levels, thus allowing these levels to exceed the bounding value *z*=η for η negative enough.

### Closure relation and model formulation

Our goal is twofold. First, to construct an approximation to the translated Dirichlet–Neumann operator and second, to look for closure relations between the unknowns , , ϕ_{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 the fourth and third order, respectively, (see remark 2.2) and the Laplace equation Δϕ_{1}=−∂_{z}*w*_{1} at *z*=0 yields the desired closure relations, namely
2.9
where , and .

We can now state our double-layer model as follows: 2.10 where is an approximation to the static Dirichlet–Neumann operator , which is detailed in the following section.

The main advantages of this model are that (i) it consists of only four equations both in 1DH and 2DH, (ii) it can be used in complex domains, and (iii) the approximate Dirichlet–Neumann operator can be computed once and for all at *t*=0 since this operator is static. Furthermore, we will see in §3 that it includes at most only second-order horizontal derivatives. This is a major improvement in comparison with high-order Boussinesq-type models, such as those of Jamois *et al.* (2006) and Madsen *et al.* (2002, 2003), which contain fourth- and fifth-order derivatives, respectively, 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 orders *O*(ϵ^{4}μ^{2},ϵ^{5}μ^{2}) and *O*(ϵ^{4}μ^{2},ϵ^{3}μ^{2}), where ϵ=*a*/*h*_{0} and 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 , whose maximum value is admittedly for a stable wave (Williams 1981). 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

### Theoretical solutions to the Laplace equations

The first step in the derivation of the approximate static Dirichlet–Neumann operator is to look for solutions to the Laplace equations
3.1
where we have redefined the upper-layer domain as . To this end, we follow the generalized Boussinesq procedure introduced by Madsen *et al.* (2002, 2003), which consists of 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 Madsen *et al.* (2002, 2003) allows finding 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 *w*_{i}, such that ϕ_{i}=ϕ_{i}(*X*,*z*=*z*_{i}(*X*)) and *w*_{i}=*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 (*z*−*z*_{i}),
3.2
where the choice of the vertical variable (*z*−*z*_{i}) instead of *z* actually allows us to save one step compared with the procedure of Madsen *et al.* (2002, 2003). Substituting (3.2) into the Laplace equations (3.1) and using the mild-slope assumption leads to the recurrence relation
3.3
Observing that for all *k* yields the expression of and in terms of ϕ_{i} and *w*_{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 *w*_{i} for all (*X*,*z*) belonging to Ω_{i},
3.4
where and are infinite-series pseudo-differential operators defined by
3.5
where the slope terms Γ_{ϕi} and Γ_{wi} are given by
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.

### 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 Madsen *et al.* (2002, 2003), and even of Jamois *et al.* (2006). We take advantage of this possibility and truncate the operators and by retaining only the first two terms of the series, which leads to the following approximations
3.7
We first plug these approximations into (3.4), which leads to the following truncated expressions of ϕ_{i} and *w*_{i}
3.8
where
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},
3.10
The advantage of this new formulation will become clear in the following.

### 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 *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 corresponds to the dispersion parameter in each layer. Supposing that all the derivatives are of order *O*(1), we can analyse the order of the third terms of each series in deep water, for instance *k**h*_{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: ; ; and . 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).

### 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 Madsen *et al.* (2002, 2003) and expand the variables ϕ_{i} and *w*_{i} in terms of auxiliary variables and through the relations
3.11
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:
3.12
where we have again kept the first two terms in each modified truncated series, except the fifth derivative of 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 Madsen *et al.* (2002, 2003), 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 and *q*_{i}+4*p*_{i} (in front of the fourth-order and third-order derivatives, respectively) 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 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}, for the interface and for the sea bottom. Consequently, requiring that the quantity vanishes at the still-water level and at the interface, and that the quantity vanishes at the interface and at the sea bottom eliminates all the fourth-order derivatives from the model. Using (3.9), this yields
Then, taking *q*_{i}=−4*p*_{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}, defined by *z*_{1}=−σ_{1}*h* and *z*_{2}=−σ_{2}*h*, respectively, are thus taken at the middle of the upper layer at rest and at the middle of the lower one, respectively.

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, and . We first obtain the following expressions at the still-water level:
3.13
where
3.14
At the interface , we use the mild-slope approximation to obtain
3.15
where
3.16
Finally, the expressions at the sea bottom are
3.17

### Formulation of the approximate static Dirichlet–Neumann operator

The final step in deriving our approximate static Dirichlet–Neumann operator is to reformulate the boundary conditions at and 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 and on one side, and and on the other side, differ from each other. A simple reformulation of the continuity condition is
Consequently, the continuity conditions (2.5e) and (2.5f ) take the form
3.18
Finally, applying the operator to (2.5g) leads to the following reformulation of the kinematic condition at the bottom:
3.19

Gathering all the previous results, we are able to construct the system of five equations on the six unknowns ϕ_{0}, *w*_{0}, , , and that defines the approximate operator linking *w*_{0} to ϕ_{0}. The first 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 and are obtained by plugging (3.15) into (3.18). The fourth is the condition at the sea bottom derived by plugging (3.17) into (3.19), and the last 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
3.20
with the differential operators
3.21
Incidentally, we observe that it is possible to eliminate all the third-order derivatives from the operators , , and in these equations. Indeed, applying the operator ∇*h*⋅∇ to the fourth, first and second equation of (3.20), using the mild-slope approximation and combining the results yields
3.22

Plugging these results into (3.20), we formulate our approximate static Dirichlet–Neumann operator as
3.23
where the differential operators *P*_{0}, , , and are given by (3.21) and the new operators are defined by
3.24

We denote by the matrix differential operator linking , , and to ϕ_{0} in (3.23). We can write each of the operators as a sum of a first-order (in *h*) operator and a mild-slope operator , yielding
3.25
where
We denote by ** U** the vector and by

**the right-hand side (**

*F**P*

_{0}ϕ

_{0},

*Q*

_{1}ϕ

_{0},

*Q*

_{2}ϕ

_{0},

*Q*

_{3}ϕ

_{0})

^{T}. The differential system in (3.23) then takes the form 3.26 A Fourier analysis of the differential operator shows that it can be inverted in the case of a flat bottom. Considering an uneven bottom, we can write , where ,

*h*

_{0}is the mean depth and is of order

*O*(∇

*h*). This ensures that the differential operator is invertible when ∇

*h*is small enough. Finally, the differential operator is defined by 3.27 where Id is the identity operator, and is thus an approximate inverse of the operator up to

*O*(|∇

*h*|

^{2}) terms. Hence, 3.28 which yields the explicit expressions of , , and in terms of ϕ

_{0}.

The very last step consists of introducing the operators
3.29
and plugging the expressions of and 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
3.30
This expression completes the formulation (2.10) of our double-layer Boussinesq-type 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 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 , and using (3.28). Therefore, we can recover the vertical profiles of the velocity potentials ϕ_{1} and ϕ_{2} and the vertical velocities *w*_{1} and *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
3.32
for if *i*=1 and for if *i*=2, where
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 Madsen *et al.* (2002, 2003). 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 analyse 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.

### Linearization of the governing equations

In order to investigate these linear properties, we restrict the analysis to the one-dimensional case. We linearize the governing equations (2.10) around steady state, which yields and , and leads to the linearized model
4.1a
4.1b
4.1c
or equivalently, up to terms,
4.2a
4.2b
4.2c
4.2d
where *h*_{x} is the bottom slope and the differential operators are left unchanged, except that they are one-dimensional operators here. 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) and (4.2b) to obtain
4.3a
4.3b
4.3c
4.3d
where ϕ^{0}=*P*_{0}ϕ_{0},*N*=*P*_{0}η and *W*_{0}=*P*_{0}*w*_{0}, and where we have plugged the relation into to obtain for *i*∈{1,3}, since *Q*_{i} are differential operators of order *O*(*h*_{x}).

Plugging the expressions of ϕ^{0} and *W*^{0} in terms of and (given by the first line of the differential system (4.3c) and by (4.3d), respectively) into equations (4.3a) and (4.3b) leads to the final reformulation of the linearized model
4.4

We now look for solutions of the classical form
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 Madsen *et al.* (2002, 2003), 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.

### 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:
4.6
where *c* is the wave celerity and where the (*a*_{i}) and (*b*_{i}) coefficients are given by
where *S*=σ(1−σ)/12. This dispersion relation is compared, in §4*e*, to the exact linear dispersion relation given by the Stokes linear theory, namely
4.7
and to its Padé [6,8] approximation, which has the same rational form as (4.6).

### 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) and (3.33) since
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*) and *w*_{1}(*z*) in the upper layer and ϕ_{2}(*z*) and *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
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
4.10

### Linear shoaling

We now aim at determining the linear shoaling gradient γ_{0} of the double-layer model defined by
4.11
In order to determine this shoaling gradient, we use the method proposed by Madsen *et al.* (2002, 2003). Coming back to the linear formulation (4.4) together with the 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 expressing *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}, *k**h* 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 Madsen & Sørensen (1992) using energy flux conservation combined with Stokes linear theory, namely
4.12

### 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. Therefore, we tune the free parameter σ that defines the interface level .

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
4.13
where α∈{*c*,ϕ,*w*,γ}, *K* is a reference relative water depth, and
4.14
where *c*, *c*_{s}, (ϕ,*w*), (ϕ_{s},*w*_{s}), γ_{0} and γ_{s} come, respectively, from equations (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/*k**h* helps to keep the errors to a minimum for low wave numbers, as in Madsen *et al.* (2002, 2003). 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/*k**h* is especially well suited for the shoaling gradient errors, which are far more critical in shallow water than in deep water.

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 *k**h* is limited. We thus choose to minimize and simultaneously to infer the optimal value for σ, and then optimize the shoaling gradient error differently. We start with the minimization of the errors and through the average error , 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 per cent within the range *k**h*∈[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 *k**h*=3. Indeed, we can see that a 2 per cent error is reached at the very-deep-water value *k**h*=24 for σ_{opt}=0.473 with a very small 0.01 per cent error peak at *k**h*≈3, whereas the same error is reached at *k**h*=28 for σ_{opt}=0.314, but with a 0.04 per cent error peak at *k**h*≈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 Jamois *et al.* (2006) and the one of Madsen *et al.* (2002, 2003). 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 Jamois *et al.* (2006); a 2 per cent error is reached at the very-deep-water value *k**h*≈28 for the double-layer model, whereas the same error is already reached at *k**h*≈18 for the Padé [6,8] approximation and *k**h*≈12 for the model of Jamois *et al.* (2006). In comparison with these results, Lynett & Liu (2004*a*,*b*) showed that their double-layer model reaches the same 2 per cent error (not plotted here) at *k**h*=8. As for the model derived by Madsen *et al.* (2002, 2003), a 2 per cent error is reached at *k**h*=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*_{ϕ} (*a*) and *E*_{w} (*b*) 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 per cent error is reached at *k**h*≈4 and a 2 per cent error at *k**h*≈8. As far as the velocity potential is concerned, a 1 per cent error is reached at *k*≈6.5 and a 2 per cent error is reached at *k**h*≈10. By comparison, the model derived by Madsen *et al.* (2002, 2003) yields a 2 per cent error at *k**h*=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 sacrificing 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 *k**h*=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 the order 0.2 per cent in shallow and intermediate waters. 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.

### Improved model with tightened shoaling properties

We now consider the linear optimization of the shoaling gradient properties. As shown in figure 4, the model linear shoaling gradient γ_{0} (dashed line) only matches the exact linear shoaling gradient γ_{s} (solid line) up to *k**h*=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 *k**h*=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 (3.23), we introduce a new constant parameter *r* and apply the differential operator 1+*r**h*∇*h*⋅∇ to the last equation of (3.23), which yields
4.15
where .

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 in the slope terms. By doing this, we greatly improve the shoaling gradient properties, as we will see in the following. For the moment, we rewrite (4.15) as
4.16
We then take advantage of (3.22) to eliminate the third-order derivative from the previous equation, and obtain
4.17
This implies redefining the approximate operator as follows:
4.18
where
4.19
Our modified model is thus identical to (2.10), but with as redefined above. Starting from it, the new linearized model remains essentially the same, except that the second equation of (4.4) now is
with . 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 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 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 *k**h*≈12. This accuracy is the same as that reached by the model of Madsen *et al.* (2002, 2003).

To summarise, our final double-layer Boussinesq-type model consists of (2.10) and (4.18)–(4.19) with σ=0.314 and *r*=0.0076. This model exhibits excellent linear properties: the phase celerity is accurate up to *k**h*=28; the vertical velocity profiles are accurate up to *k**h*=10 for the velocity potential and up to *k**h*=8 for the vertical velocity component; and the shoaling gradient is accurate up to *k**h*=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 Madsen *et al.* (2002, 2003), 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 (Dean 1965; Rienecker & Fenton 1981). 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 software called Stream_HT by one of the authors (Benoit *et al.* 2002). For the selected application, the domain of interest covers one wavelength (*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 *k**h*=3π≈9.425, which corresponds to deep-water conditions. The wave height is chosen as *H*=6.4 m, so that the steepness is *H*/*L*=0.1 or *k**H*/2=π/10≈0.314, i.e. approximately 70 per cent of the theoretical maximum value of the steepness for a stable wave (Williams 1981). 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 eighth-order Savitsky–Golay smoothing filter is applied twice after each time step to η, and . 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 wavelength) 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*. Computations have been performed with the deep-water (*k**h*=10) optimal value σ=0.314. Some simulations (that are not reported here for brevity) have shown that any value in [0.28,0.36] leads to very similar results, while it quickly deteriorates outside this range. Results obtained after durations of 10 and 25 *T* are plotted in figure 5 for the free surface elevation and the free surface velocity potential. The results are compared with the reference solution obtained with the stream function approach (which propagates at constant 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 with 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 approximately 0.08 per cent, 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 *k**h*=10, at least for flat-bottom conditions.

## Footnotes

- Received December 5, 2008.
- Accepted April 7, 2009.

- © 2009 The Royal Society