A coupled-mode system with application to nonlinear water waves propagating in finite water depth and in variable bathymetry regions

Coastal Engineering 58 (2011) 337–350

Contents lists available at ScienceDirect

Coastal Engineering j o u r n a l h o m e p a g e : w w w. e l s ev i e r. c o m / l o c a t e / c o a s t a l e n g

A coupled-mode system with application to nonlinear water waves propagating in finite water depth and in variable bathymetry regions K.A. Belibassakis a,⁎, G.A. Athanassoulis b a b

Dept. of Naval Architecture, Technological Educational Institute of Athens, Ag. Spyridonos 12210, Athens, Greece School of Naval Architecture and Marine Engineering, National Technical University of Athens, Zografos 15773, Athens, Greece

a r t i c l e

i n f o

Article history: Received 1 September 2010 Received in revised form 15 November 2010 Accepted 25 November 2010 Available online 6 January 2011 Keywords: Coupled modes Nonlinear water waves Finite depth Variable bathymetry

a b s t r a c t A non-linear coupled-mode system of horizontal equations is presented, modelling the evolution of nonlinear water waves in finite depth over a general bottom topography. The vertical structure of the wave field is represented by means of a local-mode series expansion of the wave potential. This series contains the usual propagating and evanescent modes, plus two additional terms, the free-surface mode and the sloping-bottom mode, enabling to consistently treat the non-vertical end-conditions at the free-surface and the bottom boundaries. The present coupled-mode system fully accounts for the effects of non-linearity and dispersion, and the local-mode series exhibits fast convergence. Thus, a small number of modes (up to 5–6) are usually enough for precise numerical solution. In the present work, the coupled-mode system is applied to the numerical investigation of families of steady travelling wave solutions in constant depth, corresponding to a wide range of water depths, ranging from intermediate depth to shallow-water wave conditions, and its results are compared vs. Stokes and cnoidal wave theories, as well as with fully nonlinear Fourier methods. Furthermore, numerical results are presented for waves propagating over variable bathymetry regions and compared with nonlinear methods based on boundary integral formulation and experimental data, showing good agreement. © 2010 Elsevier B.V. All rights reserved.

1. Introduction The nonlinear water wave–seabed interaction problem is an interesting free-boundary problem for which a broad class of mathematical models and approximation techniques have been developed. An important feature of this problem is that propagation phenomena take place on the horizontal plane, and non-local couplings (wave–wave and seabed–wave) exist through the vertical structure of the flow field. Fully numerical methods, based on FEM (e.g., Wu and Eatock Taylor, 1994, 1995; Turnbull et al., 2003), Finite Difference Methods, in conjunction with vertical transformation (e.g., Bingham and Zhang, 2007; Engsig-Karup et al., 2009), spectral methods (Dommermuth and Yue, 1987), BEM (e.g., Grilli et al., 1989, 2001; Ohyama and Nadaoka, 1991, 1994; Christou et al., 2009), and Smooth Particle Methods (e.g., Dalrymple and Rogers, 2006; Dalrymple et al., 2007) have been developed to treat the problem in general bathymetry domains; see also the reviews by Tsai and Yue (1996) and Dias and Bridges (2006). Some of the above methods, particularly the ones based on BEM and SPH, are able to represent highly nonlinear phenomena like wave folding and breaking, but they are computa⁎ Corresponding author. Tel.: + 30 2105385389; fax: + 30 2105911442. E-mail addresses: [email protected] (K.A. Belibassakis), [email protected] (G.A. Athanassoulis). 0378-3839/$ – see front matter © 2010 Elsevier B.V. All rights reserved. doi:10.1016/j.coastaleng.2010.11.007

tionally intense and thus, their use is more appropriate for shortrange water wave propagation and application to local interaction problems. Various equivalent reformulations of the fully nonlinear-nonlocal water-wave problem, in finite water depth, have been obtained, as e.g., by means of Hamilton's principle and the Dirichlet to Neumann (DtN) map (Craig and Sulem, 1993), by means of Lagrange equations of fluid dynamics and analyticity of the wave potential in the liquid domain (e.g., Craig, 1985, extended by Wu, 1999 in the 3D case, using Clifford analysis), variational principles (e.g., Groves and Toland, 1997) and other approaches. In particular, using as canonical variables the values of the wave potential on the free-surface and the freesurface elevation φ(x, t) = Φ(x, z = η(x, t), t) and η = η(x, t), in conjunction with Hamilton's principle, the following system is obtained in 2D (see also Craig and Sulem, 1993), ηt −GðηÞ∘φ = 0; φt +

  1 2 2
 φx −ðGðηÞ∘φÞ −2ηx φx GðηÞ∘φ + gη = 0; 2 1 + η2x

ð1Þ
1 = 2 ½∂Φ =∂nz = ηðx;t Þ denotes the DtN operawhere GðηÞ∘φ = 1 + η2x tor (∂ Φ/∂ n is the normal derivative), which is defined through the b.v. p on Φ, consisted by the Laplace equation, the bottom (no-entrance) boundary condition, the Dirichlet data on the instantaneous freesurface η = η(x, t) and appropriate lateral conditions. These nonlinear,

338

K.A. Belibassakis, G.A. Athanassoulis / Coastal Engineering 58 (2011) 337–350

nonlocal equations have been extensively studied, in constant depth (see, e.g., Craig and Sulem, 1993; Craig and Nicholls, 2002) and in variable bottom topography, (e.g., Craig et al., 2005). For small (but finite) free-surface elevation, the nonlocal operator G(η) is analytic (Coifman and Meyer, 1985), and thus, it can be expanded in functional Volterra–Taylor series. In constant depth, the successive terms are obtained recursively (Craig and Sulem, 1993). For treating the problem in general bathymetry domains, the DtN map could be based on a boundary integral formulation (see, e.g., Clamond and Grue, 2001; Grue, 2002). On the other hand, various simplified formulations (leading to model equations) of the water-wave propagation problem have been obtained by many authors under various asymptotic assumptions. In this direction, we mention here the weakly nonlinear system with enhanced dispersion characteristics by Trulsen and Dysthe (1996), Trulsen et al. (2001). Another, simple model equation, retaining up to a degree nonlinearity and nonlocality is Whitham's (1967) equation, which for the 1D case (in normalized form) reads as follows: ut + uux + IK ðuÞ = 0;

ð2Þ

where IK(u) is a pseudodifferential operator on the wave field u. The generalised Witham system (see, e.g., Naumkin and Shishmarev, 1994) combines explicit (quadratic) nonlinearity with linear nonlocal terms, modelling to some extent dispersion. From the latter system, Witham and Boussinesq model equations can be obtained as special cases. Applications of the Boussinesq model to water wave propagation and interaction with the seabed in coastal areas are numerous. We mention here specific extensions of this model to the case of mildly varying bathymetry by Peregrine (1967), for which numerical solvers have been presented by various authors; see, e.g., Beji-Battjes (1994). Extended Boussinesq systems fully accounting for nonlinearity and dispersion have been also developed; see, e.g., Liu (1995), Madsen et al. (2002, 2003). Similarly, Green–Naghdi models (e.g., Kim et al., 2001, 2003) have been developed and successfully applied to nonlinear water waves in variable depth. The usual form of Bousinessq models is with respect to free surface elevation η and horizontal velocity u, however, they can be also formulated with respect to the free surface elevation and wave potential φ (Massel, 1989, Chap.4), as follows: ηt + ððh + ηÞφx Þx = 0; φt + gη +

1 h2 2 ðφx Þ − φxxt + … = 0; 2 3

ð3Þ

where h denotes the (possibly varying) depth. In the same direction, extended versions of mild-slope equation have been presented by various authors with applicability also to short wave — bottom interaction problems, e.g., Beji and Nadaoka (1997), Nadaoka et al. (1997), Tang and Ouellet (1997). Also, high-order methods based on the velocity potential formulation have been developed by Jamois et al. (2006), Bingham et al. (2009) and others. In constant, finite water depth nonlinear wave solutions have been obtained, both analytically (see, e.g., Fenton, 1990) and numerically (e.g., Schwartz, 1974; Cokelet, 1977; Rienecker and Fenton, 1981; Drennan et al., 1988 and others). A detailed presentation of nonlinear steady water waves, including the effects of surface tension, and discussing in detail the appearance turning points and bifurcation of highly nonlinear solutions, can be found in Dias and Kharif (1999) and in Okamoto and Shoji (2001). In particular, Schwartz (1974) and Cokelet (1977) extended the Stokes theory to very high order, calculating the speed, momentum, energy and other integral properties of waves in water of arbitrary uniform depth. Parallel to Stokes expansion, nonlinear methods based on the use of truncated Fourier expansions have been also developed, see, e.g., Rienecker and Fenton (1981). The latter approaches are based on stream function formulation, however similar developments based on wave potential

have been presented. In this direction, a homotopy analysis method based on the velocity potential formulation is more recently presented by Tao et al. (2007) to describe the nonlinear progressive waves in water of finite constant depth. Fourier expansion methods have been also proved very useful tools for the description of nonlinear transient wave systems, including application to irregular wave kinematics (see, e.g., Sobey, 1992; Baldock and Swan, 1994; Johannessen and Swan, 1997). In the present work, we consider the problem of non-linear gravity waves propagating over a general bathymetry. An essential feature of this problem is that the wave field is not spatially periodic. Extra difficulties are introduced by the fact that no asymptotic assumptions concerning the free-surface and bottom slope are made. The goal of the present work is to develop an alternative and efficient, from the point of view of the numerical solution, reformulation (in the form of a system of horizontal differential equations) governing the evolution of water waves in variable bathymetry regions. Our formulation is equivalent with the fully nonlinear water–wave problem without introducing any simplifying assumptions concerning the characteristic non-dimensional parameters, and belongs to the same category with high-order models based on velocity potential. Using Luke's (1967) variational principle, in conjunction with an enhanced local-mode series expansion of the wave potential developed by the authors (Athanassoulis and Belibassakis, 2002; Belibassakis and Athanassoulis, 2006), we first present the non-linear coupled-mode system (CMS) on the horizontal plane with respect to the mode amplitudes. The above local-mode series contains the usual propagating and evanescent modes, plus two additional terms, the free-surface mode and the sloping-bottom mode, enabling to consistently treat the non-vertical end-conditions at the free-surface and the bottom boundaries and exhibiting fast convergence. The present coupled-mode system is applied to the numerical investigation of families of steady travelling wave solutions in constant depth regions, corresponding to a wide range of water depths, ranging from intermediate to shallow wave conditions. The derivation of the latter solutions is important for comparison with known theories and validation of the present model, as well as for the consistent formulation of boundary conditions for the CMS in the case of water wave propagation in non-homogeneous environments. Numerical results are compared vs. Stokes and cnoidal wave theories, as well as vs. fully nonlinear Fourier methods, respectively, showing very good agreement. It is illustrated that the present coupled-mode system fully accounts for the effects of non-linearity and dispersion, and the local-mode series exhibits fast convergence. Thus, a small number of modes (up to 5–6) are usually enough for the precise numerical solution, provided that the two new modes, the free-surface and the sloping-bottom ones, are included in the local-mode series (the latter only in the case of non-horizontal bed). Finally, numerical results are presented for waves propagating over variable bathymetry regions and compared with nonlinear methods based on boundary integral formulation and experimental data, showing good agreement.

2. Variational formulation For simplicity, we restrict ourselves to the 2D problem corresponding to normally incident waves in a variable bathymetry region; see Fig. 1. However, the present analysis can be straightforward generalised to 3D, i.e. the two horizontal dimensions associated with the propagation space and the vertical dimension. The liquid domain is a generally-shaped (non-uniform) strip D, bounded below by the seabed z = − h(x), and above by the free surface z = η(x, t). The function h(x) represents the local depth, measured from the mean water level (z = 0). The functions h(x) and η(x, t) are assumed bounded and smooth (− h(x) b η(x, t)). Moreover, the function η(x, t) is continuously dependent on time t.

K.A. Belibassakis, G.A. Athanassoulis / Coastal Engineering 58 (2011) 337–350

339

On the other hand, the variational equation δηF = 0 models water– wave dynamics,     1 ∂Φ 2 1 ∂Φ 2 ∂Φ + + + gη = 0; x1 bxbx2 ; z = ηðx; t Þ: 2 ∂x 2 ∂z ∂t

Fig. 1. Water waves propagating over a generally-shaped (nonuniform) strip D.

A main feature of the water–wave problem is that the propagation space does not coincide with the physical space. While the latter is the whole liquid domain (an irregularly shaped horizontal strip), the former is only the horizontal plane. This fact, leads to the reformulation of the propagation problem as a non-local wave equation in the propagation (horizontal) space. Under the assumptions of incompressibility and irrotationality, the problem of evolution of water waves, propagating over a variable bathymetry region, can be reformulated as a variational equation by means of Luke's (1967) variational principle. According to this formulation, the admissible fields are free of essential conditions, except, for smoothness and completeness (compatibility) prerequisites. Luke's functional, modelling the nonlinear water–wave problem, in the case of homogeneous liquid, is ( (  )  2 ) t2 x2 z = ηðx;t Þ ∂Φ 1 ∂Φ 2 ∂Φ F ½Φ; η = ∫ ∫ dxdt ∫ + + gz dz; + 2 ∂t ∂x ∂z t x1 z = −hðxÞ 1

ð4Þ where x is the horizontal and z is the vertical (positive upwards) coordinates, Φ = Φ(x, z, t) denotes the velocity potential, η = η(x, t) is the free surface elevation and g is the gravitational acceleration. In this case, the Lagrangian density is the pressure. The nonlinear water-wave problem is then expressed by the variational equation δF ½Φ; η = δΦ F ½Φ; η + δη F ½Φ; η = 0;

ð5Þ

where the first variation of F ½Φ; η can be obtained as the sum of its partial variations with respect to the fields Φ = Φ(x, z, t) and η = η(x, t). On the basis of the above variational equation it is directly seen that the condition of stationarity of functional F ½Φ; η is equivalent to the nonhomogeneous, nonlinear water-wave problem; see, e.g., Whitham (1974). More precisely, the variational equation δΦF = 0 models the water-wave kinematics, 2

Eq. (7) is equivalent to the one obtained by Bernoulli's (energy) integral by setting the Bernoulli constant equal to zero, which is possible in time domain problems by redefining the velocity potential to include this (in general time-dependent) constant. However, in treating the water wave problem in steadily translating coordinate systems, as for example in the derivation of steady travelling wave solutions, this constant comes into play, attaining non-zero values in finite and, especially, shallow water depth. Assuming waves travelling to the positive x-direction, the total wave potential ϕ(ξ, z) in the moving coordinate system is written as follows (see Massel, 1989, Ch. 2) ϕðξ; zÞ = −cξ + Φðξ; zÞ;

∂ Φ ∂ Φ + = 0; x1 bxbx2 ;−hðxÞbzbηðx; t Þ; ∂x2 ∂z2

ð6aÞ

∂Φ ∂η ∂Φ ∂η − + = 0; z = ηðx; t Þ; ∂x ∂x ∂z ∂t

ð6bÞ

∂Φ ∂h ∂Φ + = 0; z = −hðxÞ: ∂x ∂x ∂z

ð6cÞ

where ξ = x−ct;

ð8Þ

and c denotes the phase speed, in accordance with Stokes first definition (see, e.g., Longuet-Higgins, 1975; Fenton, 1990). In this case, under the assumption that the mean value of the free-surface elevation is zero, the dynamic free surface boundary condition (Eq. (7)) becomes     1 ∂Φ 2 1 ∂Φ 2 ∂Φ + −c + gη = K; at z = ηðξÞ; 2 ∂ξ 2 ∂z ∂ξ

ð9Þ

where the Bernoulli constant (K) is given by (see, e.g., LonguetHiggins, 1974; Cokelet, 1977),

K=

u2b 2

ðsince η = 0Þ; with ub =

∂Φðξ; z = −hÞ ; ∂ξ

ð10Þ

the horizontal wave velocity at the bottom and an overbar denotes mean value over a periodic cell. We will come back to this point in Section 4 where the application of the present CMS to the derivation of steady travelling waves will be presented and discussed. 3. Vertical expansion of the wave potential An important feature of the variational Eq. (5) is that it permits us to obtain equivalent reformulations of the water wave problem by introducing valid alternative representations for the involved fields Φ = Φ(x, z, t) and η = η(x, t). In the present work we exploit a localmode series expansion of the wave potential Φ(x, z, t) in variable bathymetry regions, first developed by Athanassoulis and Belibassakis (1999) for the linearised problem and extended by Athanassoulis and Belibassakis (2002) to the non-linear wave problem; see also Belibassakis and Athanassoulis (2006). This expansion has the general form Φðx; z; t Þ =

2

ð7Þ

∞

∑

n = −2

φn ðx; t ÞZn ðz; hðxÞ; ηðx; t ÞÞ;

ð11Þ

where Zn(z ; h(x), η(x, t)) denote vertical functions that are parametrically dependent on the bathymetry and the free-surface elevation, having the property to form a basis in the vertical intervals z∈½−hðxÞ; ηðx; t Þ. In particular, Z−2 ðz; h; ηÞ =

μ0 h0 + 1 μ h +1 2 ð z + hÞ − 0 0 ðηþhÞ + 1; 2ðη + hÞh0 2h0

ð12aÞ

340

K.A. Belibassakis, G.A. Athanassoulis / Coastal Engineering 58 (2011) 337–350

represents the vertical structure of the term φ− 2Z− 2 called the freesurface mode, and Z−1 ðz; h; ηÞ =

μ0 h0 −1 1 2h0 −ðη + hÞðμ0 h0 + 1Þ 2 ðz + hÞ + ; ðz + hÞ + 2h0 ðη + hÞ h0 2h0

ð12bÞ represents the vertical structure of the term φ− 1Z− 1, called the sloping-bottom mode. Furthermore, Z0 ðz; h; ηÞ =

cosh ½k0 ðz + hÞ cos ½kn ðz + hÞ ; Zn ðz; h; ηÞ = ; n = 1; 2; 3; :::; cosh ½k0 ðη + hÞ cos ½kn ðη + hÞ

ð12cÞ are corresponding functions associated with the rest of the terms, called the propagating (φ0Z0) and the evanescent (φnZn, n = 1, 2, …) modes. The quantities k n = k n (h, η), n = 0, 1, 2 …, appearing in Eq. (12c), are obtained as the positive roots of the (dispersion-like) equations μ0 −k0 tanh ½k0 ðh + ηÞ = 0;

μ0 + kn tan ½kn ðh + ηÞ = 0:

extended to second-order Stokes waves (in the frequency domain) in Belibassakis and Athanassoulis (2002), where also the necessity of a free-surface additional mode has been also discussed for the satisfaction of the (second-order) free-surface boundary condition. 4. The coupled-mode system (CMS) The series expansion (Eq. (11)) permits us to obtain corresponding expansion of the variation δ Φ of the wave potential, in terms of the variations of the modal amplitudes δ φn and the free surface elevation δ η. Then, it is shown that the examined hydrodynamic problem in the variable bathymetry region reduces to the following nonlinear coupled-mode system (see Belibassakis and Athanassoulis, 2006), with respect to the mode amplitudes φn(x, t) and the free-surface elevation η(x, t): ∂η + ∂t

∞

∑

n = −2

∂2 φn ∂φ + Bmn ðηÞ n + Cmn ðηÞφn ∂x ∂x2

! = 0;

m = −2; −1; 0; 1; 2…;

ð14aÞ

ð13Þ

The (numerical) parameters μ0, h0 N 0 are positive constants, which are not subjected to a-priori restrictions. However, the specific values of these parameters influence the accuracy of truncated versions of the local-mode series (Eq. (11)), as it will be illustrated in the sequel in Section 4.2 concerning the dispersion characteristics of the present model. More details concerning the validity of the above expansion are provided in Appendix A. The usefulness of the local-mode representation (Eq. (11)) is that, substituted in the variational Eq. (5), it leads to a non-linear coupled-mode system of differential equations on the horizontal plane, with respect to unknown modal amplitudes φn(x, t) and the unknown free-surface elevation η(x, t). The coupled-mode system, in conjunction with the fast convergence properties of the local mode series, greatly facilitates the numerical solution of the present problem and will be presented in the next section. Similar modal-type series expansion has been previously introduced by Nadaoka et al. (1997) for the development of a fully dispersive, weakly nonlinear, multiterm-coupling model for water waves, with application to slowly varying bottom topography. In that work, the vertical modes have been selected to have the form cosh(kn(z + h))cosh− 1(knh), where the parameters kn N 0 are independent from the free surface elevation η(x, t). This approach can be considered as an extension of Fourier methods above a flat bottom (see, e.g., Baldock and Swan, 1994; Tao et al., 2007) to variable bathymetry regions characterised by mildly sloped bottom. In the present case, the majority of the set of vertical modes fZn ðz; h; ηÞ; n = 0; 1; 2; …g is obtained by solving a Sturm–Liouville problem, formulated at the local vertical interval − h(x) b z b η(x, t), ensuring L2 − completeness; see Appendix A. This set contains both hyperbolic and trigonometric functions, Eq. (12c), dependent on the local depth h(x) and the (instantaneous) free surface elevation η(x, t). However, the boundary conditions satisfied by these local vertical eigenfunctions are not compatible with the boundary conditions of the problem at the bottom surface, if the bottom is not horizontal or mildly sloping, as well as on the free surface. In order to overcome the mild-slope approximation concerning the bottom and to consistently satisfy the free-surface boundary conditions, the present set has been enhanced by including the two additional m o d e s fZ−2 ðz; h; ηÞ; Z−1 ðz; h; ηÞg w i t h u n kn o w n am p l i tu de s fφ−2 ðx; t Þ; φ−1 ðx; t Þg. The latter are the additional degrees of freedom required for the consistent satisfaction of the free-surface and the sloping-bottom boundary conditions, Eqs. (6b) and (6c), respectively. The idea of the sloping-bottom mode has been first presented by Athanassoulis and Belibassakis (1999) for the propagation of linearised waves in general bathymetry regions. The latter work has been

Amn ðηÞ

∞

gη + ∑

n = −2

∞

∑

∞

∑

ℓ = −2 n = −2

  ∂φn ∂η − + ½Wn z = η φn ∂t ∂t ð0;2Þ

aℓn ðηÞφℓ

∂2 φn ∂φ ∂φ ∂φ ð1;1Þ + −aℓn ðηÞ ℓ n + bℓn ðηÞφℓ n + cℓn ðηÞφℓ φn ∂x ∂x ∂x ∂x2

! = 0;

ð14bÞ z;h;ηÞ . (Here and in the sequel, square brackets where Wn ðz; h; ηÞ = ∂Zn ð∂η are used to denote evaluation of involved quantities at the indicated end-points). Since all Zn(z ; h, η) functions are normalized (i.e., Zn (z = η ; h, η) = 1) we easily obtain for the values of these functions on the free-surface

½Wn z = η = −

  ∂Zn ðz = ηÞ −1 = − μo + δ−2n ho ; n = −2; −1; 0; 1; 2…; ∂z ð15Þ

where δmn denotes Kronecker's delta. In the above equations, the 2) (1, 1) matrix-coefficients Amn(η), Bmn(η), Cmn(η), and a(0, mn (η), amn (η), bmn (η), cmn(η), are dependent on the (unknown) free-surface elevation at each horizontal position. The corresponding expressions involve the local vertical modes fZn gn = −2;−1;0;1;… and their derivatives, and can be found in Belibassakis and Athanassoulis (2006), where also details concerning the above derivation can be found. Also, in the latter work second-order version of the CMS (Eqs. (14a) and (14b)) has been developed, and numerical applications are presented demonstrating the rapid decay of the modal amplitudes and the fast convergence of the local modal series (Eq. (11)), in problems characterised by weak non-linearity. The non-linear CMS, Eqs. (14a) and (14b), has been obtained without any assumptions concerning the vertical structure of the wave potential. Thus, this system, being equivalent with the complete formulation, fully accounts for wave non-linearity and dispersion, as it will be demonstrated in the next subsection. Moreover, a distinctive feature of the present model is that no simplifications have been introduced for its derivation. Thus, in principle, various simplified models can be recovered as appropriate limiting forms of Eqs. (14a) and (14b). For example, keeping only the propagating mode Z0(z) in the expansion (Eq. (11)) and linearising the coupled-mode equations, a mild-slope model is obtained, see, e.g., Dingemans (1997). If the evanescent modes Zn(z), n = 1, 2, …, are also retained, an extended mild-slope model is obtained (Massel, 1993). If we keep in the vertical expansion of the wave potential only the vertical mode Z− 2(z), Eq. (12a), which is the quadratic in z and satisfies the flat bottom boundary condition, and retain up to second-order terms in the present CMS, a Boussinesq system (similar to Eq. (3)) is obtained; see Massel (1989, Chap. 4.2). On the other hand, if only the propagating

K.A. Belibassakis, G.A. Athanassoulis / Coastal Engineering 58 (2011) 337–350

mode Z0(z) = cosh(k0(z + h))cosh− 1(k0(η + h)) is kept in the localmode series, and retain up to second-order terms, a two-equation, weakly nonlinear system with respect to φ0(x, t) and η(x, t)is derived, similar to the time-dependent, nonlinear, mild-slope equation, developed by Beji and Nadaoka (1997). 4.1. Reformulation of the CMS as a system of two evolution equations

341

higher-order derivatives can be obtained by termwise differentiation of the local-mode series (Eq. (11)). The above statement can be rigorously proved by multiplying each equation in the system (Eq. (19)) by δφm, m = − 1,0,1,…, and adding by parts. Subsequently, since φ(x, t = t ) is ⁎ given (fixed), Eq. (16) can be used to express the variation of the freesurface mode amplitude δφ− 2 in terms of the rest of the modes, i.e. δφ−2 = − ∑ δφn . Using the latter, in conjunction with the expresn≠−2

A convenient reformulation of the CMS (Eqs. (14a) and (14b)) can be obtained by subtracting by parts (Eq. (14a), for m = − 1, 0, 1, … from Eq. (14a), for m = − 2; see Athanassoulis and Belibassakis, 2007). Thus, it is possible to eliminate the ∂ η/∂ t term from the lefthand side of the former equations. Moreover, by introducing the function φðx; t Þ = ∑ φn ðx; t Þ, which by definition equals to the n = −2

wave potential on the free-surface (since all vertical functions are normalized, i.e. Zn(x, z = η) = 1), i.e., Φðx; z = η; t Þ =

∞

∑

n = −2

φn ðx; t ÞZn ðx; z = ηÞ =

∞

∑

n = −2

φn ðx; t Þ = φðx; t Þ; ð16Þ

the CMS (Eqs. (14a) and (14b)) is equivalently formulated as a set of two evolution equations on the fields f φ; η g, as follows:   ∂η ∂ ∂φ + ðη + hÞ + N1 ðη; φn Þ = 0; ∂t ∂x ∂x

ð17aÞ

  ∂φ 1 ∂φ 2 + N2 ðη; φn Þ = 0; + gη + 2 ∂x ∂t

ð17bÞ

where the nonlinear operators N1, 2(η, φn) are defined as follows: N1 ðη; φn Þ =

∞

∑

n = −2

! 2 ∂ φn ∂φn ^ ^ ^ ð η Þ + B ð η Þ ð η Þφ A + C −2n −2n −2n n ; ∂x ∂x2

∞

− ∑

 ∂η ∞  ∑ ½Wn z = η φn + ∂t n = −2 ∞

∑

ℓ = −2 n = −2

2

^ aℓn ðηÞφℓ

4.2. Dispersion characteristics of the linearised CMS In order to examine the dispersion characteristics of the present CMS, we consider the case of a flat (horizontal) bottom (dh/dx = 0) and thus, the sloping-bottom mode vanishes (φ− 1 = 0) and is not taken into account. By linearising the present CMS (suppressing explicit nonlinearities and expanding various quantities around the mean sea level z = 0 keeping only the first-order terms), and using Eq. (14b) to eliminate the free-surface elevation from Eq. (14a), we finally obtain: ∞

ð18aÞ N2 ðη; φn Þ =

sions of the coefficients (Appendix B) and changing the order of vertical integration and mode summation, we finally obtain a variational equation, which is equivalent to Eqs. (6a), (6c) and Dirichlet data on the free-surface, Eq. (16). In concluding this subsection, we remark here that the present CMS in the form of two evolution equations on the horizontal plane, Eqs. (17a) and (17b), is clearly an extended version of the Boussinesq model, Eq. (3). However, the present system is not local and retains full non-linearity. The non-locality is endowed by the non-linear operators Eqs. (18a) and (18b), which depend on the whole set of modal amplitudes φn that are determined by solving the subsystem Eq. (19), in conjunction with Eq. (16). Based on the fact that the latter subsystem is equivalent to the DtN map, the connection between the present CMS and the Hamiltonian system, Eq. (1), becomes also clear.

∑ n =−2 n≠−1

∂2 φn ∂2 φn + αmn + γmn φn 2 ∂t ∂x2

! = 0;

ð18bÞ Furthermore, the two-equation system (Eqs. (17a) and (17b)) is subjected to the constraints imposed by the (rest of the) equations obtained from Eqs. (14a) and (14b) by subtraction by parts, i.e. ∞

The coefficients αmn and γmn are given by g

−1

z=0

αmn = 〈Zn ; Zm 〉0 =

∫

Zn ðz; h; η = 0ÞZm ðz; h; η = 0Þdz; ð21aÞ

z = −hðxÞ

and g

−1

γmn =

ð19Þ

^ mn ðηÞ; B ^mn ðηÞ appear^mn ðηÞ; C The expressions of the coefficients A ing in the above equations, as well as ^ amn ðηÞ, ^ bmn ðηÞ; ^cmn ðηÞ in Eq. (18b), are given in Appendix B. The subsystem defined by Eq. (19), in conjunction with Eq. (16), is equivalent to the DtN map associated with the calculation of the wave potential in the whole domain D, satisfying the bottom boundary condition; cf. Eq. (1). Given the spatial distribution of the free-surface wave potential φ(x, t = t ), at any time instant (t = t ), and the ⁎ ⁎ corresponding free-surface elevation η(x, t = t ), all the above coeffi⁎ cients are calculated, and the differential system defined by Eqs. (19) and (16) can be solved with respect to the unknown mode amplitudes φn(x, t = t ). The latter substituted in the expansion (Eq. (11)) provides ⁎ us with the wave potential Φ(x, z ; t = t ) all over the water domain D. ⁎ Thus, wave velocities, including the normal one on the free-surface, and

m = −2; 0; 1; 2; …

ð20Þ

!

∂ φn ∂φ +^ bℓn ðηÞφℓ n + ^cℓn ðηÞφℓ φn : ∂x ∂x2

  ∂2 φ   ∂φ n n ^ ^ ^ ^ ∑ A + B + −2n ðηÞ−Amn ðηÞ −2n ðηÞ−Bmn ðηÞ 2 ∂x ∂x n = −2   ^ ^ + C −2n ðηÞ−Cmn ðηÞ φn = 0; m = −1; 0; 1; …

−

   ∂Zn 〈ΔZn ; Zm 〉0 − : ∂z z = 0

ð21bÞ

In Eqs. (21a) and (21b) the vertical functions Zn(z) are obtained from Eqs. (12a), (12b) and (12c) by setting η = 0, i.e. Zn(z) = (z ; h, η = 0). Thus, in the above linearised case the coefficients and the solution of Eq. (20) are dependent only on the parameters μ0h and μ0h0. In order to investigate the dispersion characteristics of the linear CMS (Eq. (20)), we look for simple harmonic solutions of the form φn(x, t) = fncos(k(x ∓ ct)). Substituting the latter expression in Eq. (20) and truncating the series, we obtain a linear algebraic system with respect to the amplitudes of the modes fn, which is homogeneous, and thus, its determinant should vanish. From this equation we able to calculate the dependence (in non-dimensional form) of the quantity C = c(gh)− 1/2, as a function of the non-dimensional wavenumber (kh) and the parameters μ0 and h0, and compare with the analytical result (from linearised water–wave theory): −1 = 2

cð ghÞ

−1 = 2

= ðkhÞ

tanh

1=2

ðkhÞ:

ð22Þ

342

K.A. Belibassakis, G.A. Athanassoulis / Coastal Engineering 58 (2011) 337–350

Fig. 2 presents such a comparison, obtained by using μ0h = μ0h0 = 0.25 and 2.5, respectively. The results shown are obtained by keeping 1 (only mode 0), and 5 terms (n = − 2,0,1,2,3) in the local-mode series, Eq. (11). We observe that when only mode 0 (the propagating mode) is included in the modal series (which then becomes a single term approximation), the dispersion characteristics of the present approximation matches the analytical curve at the point k h, where ⁎ k h tanh k h = μ0h, indicated in Fig. 2 by using an arrow. (In this ⁎ ⁎ example k h = 0.52 and 2.54, respectively). Thus, by the appropriate ⁎ choice of the numerical parameters μ0 and h0, we are able to obtain a good approximation of the dispersion characteristics of the present system for kh lying in an interval around the point k h, even by using a ⁎ single-term approximation. Furthermore, we observe in Fig. 2 that the inclusion of the free-surface mode (n = − 2), in conjunction with the propagating one (n = 0) and few evanescent modes (n = 1,2,3) in the local-mode series (Eq. (11)), its convergence substantially improves to the exact result, for an extended range of wave frequencies, ranging from very shallow to deep water–wave conditions. In the example shown in Fig. 2a using 5 terms (dashed line), the error

remains less than 1%, for kh up to 10, which is considered enough for usual applications. Extensive numerical investigation of the effects of the numerical parameters μ0 and h0 on the dispersion characteristics of the present CMS has revealed that, if the number of modes retained in the localmode series is equal or greater than 6, the results become practically independent (error less than 0.5%) from the specific choice for the values of the (numerical) parameters μ0 and h0, for all nondimensional wavenumbers in the interval 0 b kh b 24. A good choice for the parameter h0 is the average or the representative depth of the domain, and for the parameter μ0 = ω20/g, where ω0 denotes the peak (angular) frequency associated with the incident waveform. Quite similar results are obtained as concerns the convergence of the wave potential and velocity in the domain. In concluding, a few modes (of the order of 5–6) are sufficient for modelling fully dispersive waves, for an extended range of frequencies, in a constant-depth strip. In the case of general bottom topography, the enhancement of the local-mode series (Eq. (11)) by the inclusion of the sloping-bottom mode (n = − 1) in the representation of the wave potential is of utmost importance, for the consistent satisfaction of the Neumann boundary condition (necessitating zero normal velocity) on the non-horizontal parts of the bottom. This requirement has been discussed in detail in Athanassoulis and Belibassakis (1999) and in Belibassakis et al. (2001), in the case of linear waves in 2D and 3D domains, respectively, and in Belibassakis and Athanassoulis (2002) in the case of second-order Stokes waves in variable bathymetry regions. In the latter works it is also shown that the enhanced series exhibits an improved rate of decay of the modal amplitudes φn(x) of the order O(n− 4). Thus, a small number of modes suffices to obtain a convergent solution, for bottom slopes of the order of 1:1, or higher. 5. Steady travelling waves in constant depth In this section the present nonlinear CMS (Eqs. (17a) and (17b)) is applied to the calculation of steady travelling wave solutions in constant depth regions, corresponding to a wide range of water depths, ranging from intermediate to shallow wave conditions. The derivation of the steady travelling solutions is important for studying non-linear dispersion characteristics of the system, for comparison with other theories and validation of the present model, as well as for the consistent specification of incident (incoming) waves in the case of time evolution problems in non-homogeneous (variable bathymetry) environments. Looking for steady travelling wave solutions in constant, finite but arbitrary depth h, characterised by (unknown) nonlinear wave celerity c, we employ the usual transformation φðx; t Þ = φðξÞ;

φn ðx; t Þ = φn ðξÞ;

ηðx; t Þ = ηðξÞ;

where ξ = x−ct;

ð23Þ and thus, ∂/∂ t = − c ∂/∂ ξ. Without loss of generality, waves propagating in the positive x-direction have been assumed. We also modify Eq. (17b) by including in the right-hand side the (unknown) Bernoulli constant K; cf. Eqs. (9) and (10). Consequently, the present nonlinear CMS (Eqs. (17a) and (17b)) is put in the following equivalent (time independent) form: cL1 ðuÞ−L0 ðuÞ + K = NðuÞ;

ð24Þ

where

u= Fig. 2. Dispersion characteristics of the linearised CMS using the present local-mode expansion, Eq. (11), and keeping 1 (n = 0, shown by dashed lines and denoted singleterm approximation) and 5 terms (n = −2,0,1,2,3, shown also by dashed line) vs. nondimensional wavenumber (kh). The analytical result, Eq. (22), is shown by using solid line. (a) μ0h = μ0h0 = 0.25 and (b) μ0h = μ0h0 = 2.5.

ηðξÞ ; L ðuÞ = φðξÞ 1



∂η = ∂ξ ; L ðuÞ = ∂φ = ∂ξ 0



0 ;N= gη



N1 ;K= N2



0 ; K

ð25Þ and the components of the nonlinear operator NðuÞ are defined by Eqs. (18a) and (18b), without the term n = − 1, since in the present

K.A. Belibassakis, G.A. Athanassoulis / Coastal Engineering 58 (2011) 337–350

case the bottom is flat (dh/dx = 0), and thus, the sloping-bottom mode identically vanishes, i.e. φ− 1 = 0. We stress here that calculation of NðuÞ depends on the solution of the subsystem consisted by Eqs. (19) and (16), which provides us with the modal amplitudes φn(ξ), n = − 2, 0, 1, 2 …., given the distributions of φ(ξ), η(ξ) on the free surface. Also, in the present (steady) case all x-derivatives appearing in Eqs. (18a), (18b), and (19) are replaced by ξ-derivatives. Various formulations of the steady water wave problem in a constant depth strip are possible, depending on the data used (see, e.g., Rienecker and Fenton, 1981, Sec. 2.1). In the present section we consider that the length λ of the periodic cell (the wavelength and thus, also the wavenumber k = 2π/λ), the water depth h (measured from the bottom to the mean water level) and the wave height H are given, and we are seeking the wave potential (Φ), the free-surface elevation (η), the phase speed (c), and the rest properties of the wave, such as the flow rate (Q) under the waves, the wave momentum, energy and the rest integral properties (see, e.g., Longuet-Higgins, 1975). Alternative formulations, based on different set of data are also possible (as, e.g. described in Rienecker and Fenton, 1981). We recall here that in the present work the mean free surface elevation has been assumed to be zero (i.e., η = 0); cf. Section 2 concerning the definition of the origin of the coordinate system. Thus, except of the fields φ(ξ), η(ξ) and the modes φn(ξ), two additional unknowns are introduced in Eq. (24): the wave speed and the Bernoulli constant (c, K). Looking for unimodal wave solutions (only one crest and one trough within the periodic cell, see, e.g., Okamoto and Shoji, 2001), the present non-linear system (Eq. (24)) is supplemented by two additional equations

ξ=λ

∫ ηðξÞdξ = 0; max½ηðξÞ−min½ηðξÞ = H;

ð26Þ

ξ=0

which enable the calculation of fully nonlinear waves in constant depth strip and the corresponding non-linear dispersion relation. We recall here that in this case Garabedian's theorem states that the solutions are symmetric (see, e.g., Okamoto and Shoji, 2001, Sec. 3.9). Non-symmetric periodic waves exist, appearing through symmetrybreaking bifurcations leading to multimodal waves (see, e.g., Zufiria, 1987), which are not examined here. The present fully nonlinear system is iteratively solved by

343

of small wave height (H/h ≤ 0.1) the initial estimate (step i) of the steady travelling solution is based on the linear wave theory, i.e. 2

kc0 = g = tanhðkhÞ; K0 = 0; φ0 ðξÞ = ðgH = 2kc0 Þ sinðkξÞ;

ð29Þ

η0 ðξÞ = ðH = 2Þ cosðkξÞ: In the case of higher waves, the initial estimate is based on an extrapolation (based on simple scaling) of a convergent solution corresponding to lower wave height. Numerical results are presented below in Figs. 4–6, for various representative cases corresponding to intermediate depth and shallow water depth conditions. To distinguish cases, we make use of the diagram presented in Fig. 3 (taken from Fenton, 1990, fig. 2), showing limits to the existence of periodic, steady travelling waves as determined by various theories, for different values of the nonlinearity parameter H/h and the shallowness parameter λ/h. In this diagram the upper limit (Fenton, 1990, Eq. (32)) of wave height to depth (for given shallowness ratio) is indicated by using a solid line, which has been obtained on the basis of results by Williams (1981). In the same diagram, the demarcation line, as estimated by Fenton (1990, Eq. (33)), is shown by using a dashed line, where Stokes theory applicable to short-waves (to the left of this line) and cnoidal theory applicable to long waves (to the right) fail to produce accurate results. To illustrate the convergence of the iterative algorithm, we present in Fig. 4 steady travelling waves of moderate nonlinearity (H/h = 0.4), corresponding to the small open circles of Fig. 3 (marked i and ii), as calculated by the present CMS using 5 modes (n = − 2, 0, 1, 2, 3) and M = 62 points for the discretization of the ξ-interval (0 ≤ ξ ≤ λ), that is found to be sufficient. In particular, in Fig. 4i the solution in the short-wave regime λ/ h = 5 is plotted, and in Fig. 4ii the corresponding one in the longwave regime λ/h = 20. The calculated values of phase speed and Bernoulli constant in the first case are kc2/g = 0.9375, K/gh = 0.0322, and in the second kc2/g = 0.3650, K/gh = 0.0473 and agree perfectly with the predictions by the method of Rienecker and Fenton (1981), using the same discretization. The top subplots in this figure show the calculated nondimensional free surface elevation η˜ ðkξÞ = 2 η˜ ðkξÞ = H and the lower subplots the corresponding non-dimensional free
 ˜ = φðkξÞ = gh3 1 = 2 . In order to derive the surface wave potential φ desired solution we start from a small-amplitude wave H/h = 0.1 in the same periodic cell, denoted in Fig. 4 by Eq. (1), using initial

(i) using initial estimates c0 ; K0 ; u0 (e.g., from the linear waterwave model), (ii) solving Eq. (24),


ck

+ 1 L1 −L0

+ Kk

 + 1

uk +

1

= Nðuk Þ;

in conjunction with Eq. (26), and determining uk and Kk + 1, and (iii) iterating until convergence, i.e.

ð27Þ

+ 1,

‖uk + 1 −uk ‖ + jck + 1 −ck j + jKk + 1 −Kk jbtolerance:

ck + 1

ð28Þ

For the solution of both Eq. (24) and the DtN subsystem, Eqs. (19) and (16), through which the non-linear operator N, Eqs. (18a) and (18b), is calculated, a numerical scheme based on central, secondorder finite differences is used, based on a horizontal grid of M-points (dξ = λ/M), in conjunction with periodicity conditions at the horizontal ends ξ = 0, λ of the domain. The above approximation finally reduces Eq. (27) to a non-linear algebraic one, which is numerically solved by the Newton–Raphson method (Press et al., 1992). In the case

Fig. 3. Limits to the existence of waves as determined by various theories, for different values of the non-linearity parameter H/h and the shallowness parameter λ/h, based on Williams (1981) results. Demarcation line between the validity of Stokes (left) and cnoidal (right side) theories (Fenton, 1990) is shown by using a dashed line.

344

K.A. Belibassakis, G.A. Athanassoulis / Coastal Engineering 58 (2011) 337–350

Fig. 4. Steady travelling solutions of moderate nonlinearity (H/h = 0.4) as calculated by the present CMS using 5 modes: (i) in the short-wave regime λ/h = 5, kc2/g = 0.9375 (left), and (ii) in the long-wave regime λ/h = 20, kc2/g = 0.3650 (right), corresponding to the small open circles of Fig. 3, respectively. Top subplots show the calculated ˜ = 2η = H and lower subplots the nondimensional free surface elevation η
 ˜ = φ = gh3 1 = 2 vs. the noncorresponding non-dimensional free surface potential φ dimensional horizontal coordinate kξ.

estimate from the linear model, Eq. (29). Then, the waves of higher amplitudes H/h = 0.2, 0.3 and finally, H/h = 0.4, denoted in Fig. 4 by Eqs. (2), (3) and (4), respectively, are calculated, based on initial estimates obtained from the previous convergent solution by extrapolation. In all cases 4–5 iterations of Eq. (27), illustrated also in Fig. 4, have been proved sufficient in order to produce convergent results. The final convergent solutions as obtained by the present CMS concerning the free-surface elevation, for H/h = 0.4, λ/h = 5 and λ/ h = 20, respectively, are compared in Fig. 5 vs. the fully non-linear Fourier method of Rienecker and Fenton (1981), showing excellent agreement. Also, in Fig. 5i and ii, respectively, our results (shown by using thin solid lines) are compared with 5th order Stokes solution and nonlinear cnoidal theory (see, e.g., Massel, 1989; Fenton, 1990). We observe in these cases of moderate nonlinearity, corresponding to parameters away from the demarcation line of Fig. 3, where Stokes and cnoidal theories are valid, respectively, that present method results being in perfect agreement with non-linear Fourier method are also consistent with the latter well-known expansions, which are used as standard tools in a variety of engineering applications. For the same two cases (H/h = 0.4, λ/h = 5 and λ/h = 20, respectively), in Fig. 6, the moduli of the modal-amplitude functions jφn ðξÞj; n = −2; 0; 1; 2; 3 are plotted vs. mode-number n, as calcu-

Fig. 5. Same as in Fig. 4 concerning the calculated free-surface elevation η˜ = 2η = H. Comparison of present CMS (thin solid line) with fully non-linear Fourier method (Rienecker and Fenton, 1981), shown by thick solid line only in ½0; π, as well as with 5th order Stokes solution (left) and nonlinear cnoidal theory (right), shown by using dashed lines.

lated by the present nonlinear CMS. We clearly observe the significance of the free-surface (n = − 2) and the propagating (n = 0) modes, in comparison to the evanescent modes (n ≥ 1). In this case the bottom is flat and the sloping-bottom mode (n = − 1) is zero, however, in general bathymetry regions, the latter mode becomes also important, its significance been directly connected with the local value of the bottom slope. We note here that the horizontal axis of Fig. 6 is a multiple replica of the interval ξ a [0, λ], i.e. a sequence of repeated intervals [0, λ], each one associated with a mode, and named after the mode number (n). In the nth replica of interval [0, λ] the modulus of the amplitude of the nth mode ( jφn j) is plotted. In the same figure the curve 0.01n−4 is drawn, which bounds the maxima of the amplitudes of the modal functions (maxξ jφn ðξÞj) as obtained by the present representation, Eq. (11). From this, and many similar numerical results, we obtain that the rate of decay of the mode
 amplitudes is jφn j = O n−4 , quite similarly as in the linearised (Athanassoulis and Belibassakis, 1999) and the weakly non-linear (Belibassakis and Athanassoulis, 2002) problems. The last result ensures the absolute and uniform convergence of the local mode series (Eq. (11)), and its derivatives, up to (and including) the boundaries, and permits us to obtain accurate solutions in the general case keeping a small number of modes (5–6) in the expansions. The same fact also supports efficient application of the present method to

K.A. Belibassakis, G.A. Athanassoulis / Coastal Engineering 58 (2011) 337–350

345

Fenton, 1981), shown by solid lines. Furthermore, for very small waveheight (H → 0) our results are fully consistent with corresponding predictions obtained by the linear dispersion relation (Eq. (29)), which are also indicated in Fig. 7d by using small open circles (on the horizontal axis). 6. Nonlinear waves in variable bathymetry regions In this section, the discrete scheme developed for the numerical solution of the present CMS in the time domain, Eqs. (17a) and (17b), in conjunction with the DtN subsystem Eqs. (16) and (19), is discussed, and numerical results are presented for two characteristic cases of non-linear water wave propagation in variable bathymetry domains. First, the present system is put in an evolution equation form, as follows: ∂u = −RðuÞ; ∂t

ð30Þ

where  8  ∂ ∂φ > > > < ∂x ðη + hÞ ∂x + N1 ðη; φn Þ ηðx; t Þ ; u= ; RðuÞ =   φ ðx; t Þ > > > > 1 ∂φ 2 > > : : gη + + N2 ðη; φn Þ 2 ∂x 8 > > > <

ð31Þ

and the components of the nonlinear operator NðuÞ are defined by Eqs. (18a) and (18b), on the basis of the field uðx; t Þ and the modal amplitudes φn(x, t). The latter, at each time step, are determined as the solution of the subsystem Eqs. (19) and (16). 6.1. The numerical scheme

Fig. 6. Same as in Fig. 4 concerning the moduli jφn ðξÞ j of the modal-amplitude functions vs. mode-number n, as obtained by the present local-mode series expansion, Eq. (11). The bounding curves of the maxima of j φn ðξÞ j are shown by using dashed lines.

propagation of general wave systems in 3D domains with complex bathymetric variations. First results in this direction concerning the linearised CMS have been presented in Belibassakis et al. (2001) and in Gerostathis et al. (2008). In Fig. 7, we present families of non-linear travelling waves as calculated by the present method (using as before 5 modes in the local-mode series and M = 62 gridpoints for the spatial discretization of the CMS), in intermediate water depth (λ/h = 4), around the border between short and long waves (λ/h = 9), and in very shallow water depth (λ/h = 25). In these plots only the final convergent solutions are shown concerning the non-dimensional free surface elevation ˜ = 2η = H vs. the non-dimensional horizontal coordinate kξ, for η various waveheights, from H/h = 0.1 up to the limit for each given shallowness ratio, as indicated in Fig. 3 by using crosses. Finally, in Fig. 7d plots of the non-linear dispersion relation are presented, showing the dependence of the non-dimensional phase speed (kc2/g) on the corresponding wave height (H/h), for the above values of shallowness parameter, respectively. Present method results, shown in Fig. 7d by using crosses, are again found to be in excellent agreement with the non-linear Fourier method (Rienecker and

The discrete version of Eqs. (30) and (31) is obtained by truncating the local-mode series (Eq. (11)) to a finite number Nm = Ne + 2 of terms (modes), where Ne is the number of evanescent modes retained (the upper index of truncation of the series). Moreover, second-order, central finite differences are used to approximate the horizontal xderivatives appearing in the system. On the basis of the above considerations, the discrete system in the case of two dimensional problems (one horizontal and one vertical dimension) is finally reduced to a linear algebraic system. At present the system is numerically solved by applying direct method in Matlab®. The coefficient matrix of the system is block structured, each block being a sparse 3-diagonal submatrix, and has a total dimension Nd = NmMx, where Mx is the number of gridpoints subdividing the horizontal domain. In the case of three-dimensional problems the total dimension of the discrete system will be Nd = NmMxMy where My is the number of gridpoints subdividing the transverse horizontal direction, and the coefficient matrix will be again block structured, each block being a band-diagonal submatrix with 5 non-zero elements per line. The problem is treated in a finite horizontal interval ranging from x = a to x = b, containing strictly into its interior the depth inhomogeneity. Thus, the whole domain is decomposed into a finite variable bathymetry part (x1 ≤ x ≤ x2, where the depth h(x) is variable), joining two subregions of constant but, possibly, different depth: the region of incidence (a b x b x1, where the depth is constant end equal to h1) and the region of transmission (x2 b x b b, where the depth is also constant end equal to h2). Assuming waves propagating to the positive x-direction, the wave field (free-surface elevation and potential, as well as the amplitudes of the modes) are specified at the left end (x = a) on the basis of the steady travelling wave solution which is calculated as described in the previous section in the incidence, constant-depth subregion.

346

K.A. Belibassakis, G.A. Athanassoulis / Coastal Engineering 58 (2011) 337–350

˜ = 2η = H of steady travelling waves as calculated by the present CMS (a) in inter-mediate water depth (λ/h = 4), (b) around the border between short Fig. 7. Free-surface elevation η and long waves (λ/h = 9), and (c) in very shallow water depth (λ/h = 25). Derivation of the whole family of solutions from relatively small waveheight (H ≈ 0) up to the limit of wave height to depth (for each given shallowness ratio), as also indicated in Fig. 3 by using crosses. (d) Plot of non-dimensional phase speed (kc2/g) vs. wave height (H/h) for the above values of shallowness parameter, respectively.

A perfectly matched layer (PML) model, as e.g., described in Berenger (1994) and Turkel and Yefet (1998), is applied to absorb wave energy reaching the vicinity of the right end (x = b), in conjunction with homogeneous Neumann condition exactly at x = b (also modelled by second-order, backward finite difference scheme). Following previous analysis presented by Diaz and Joly (2006) for the scalar wave equation in the time domain (see also Hu, 2001 for similar analysis concerning Euler equations) the time-derivative operator in the left hand-side of Eq. (30) is substituted by the mixed operator ∂u + σ ðxÞu, where the PML-parameter σ(x) is a positive absorption ∂t coefficient with support in the subdomain b − ℓ b x b b, where ℓ denotes the thickness of the layer. In the present work we exploit results by Collino and Monk (1998), applied by Belibassakis et al. (2001) to the case of linear harmonic water waves in variable bathymetry, concerning the optimum distribution of the absorption coefficient within the layer, in the form 3

σ ðxÞ = σ0 ðx =ℓÞ ; in

b−ℓbxbb;

ð32Þ

and zero elsewhere. The selection of the parameter σ0 is based on minimization of the reflection coefficient associated with the amplitude the back-propagating waves in flat subdomains (as, e.g. in the region of transmission). The optimum value is dependent from the depth h2, the central frequency ω of the waveform, as well as from numerical parameters as the grid resolution (see Belibassakis et al., 2001, Sec. 5). In the numerical examples presented and discussed below, the region of incidence has been artificially extended and the CMS has been numerically integrated up to a time ensuring that the reflected/ diffracted wave due to depth inhomogeneity (back-propagating towards x = a) does not affect the incoming wave condition. Future extensions of the present method will be enhanced by including a second PML (in the vicinity of x = a) designed to absorb the reflected wave energy. On the basis of all the above, standard Runge–Kutta methods are applicable for integration of the evolution equations, starting from rest, formulated as an initial-boundary value problem. In the examples shown below computations were based on 2nd/3 rd order

K.A. Belibassakis, G.A. Athanassoulis / Coastal Engineering 58 (2011) 337–350

Runge–Kutta method, using Matlab®-ode23 (see, e.g., Shampine and Reichelt, 1997). The scaling of the solution effort by the present method is estimated by examining the rate of increase of computation time (on a 2.27 GHz Pentium-P8400 with 4 Gb of RAM memory) with Nd = NmMx. In particular, the initial development for 100 time-steps (starting from rest) of a non-linear wave solution in constant depth, corresponding to λ/h= 4 and H/h= 0.4 (i.e. about 80% of the limiting steepness for this value of the shallowness parameter) has been examined, retaining Nm = 5 modes (n = −2,0,1,2,3) in the present local-mode series. It is found that the present method computation time increases linearly with the discrete dimension of the system, as CPU ≈ 0.4Nd (or log10CPU ≈ log10Nd − 0.4). This result is found to be compatible (although somewhat more expensive computationally) with similar results from high-order methods based on finite differences in conjunction with vertical transformation; see, e.g. Bingham and Zhang (2007, fig. 1). We note however, that iterative methods (as also discussed in the latter work) are applicable and could significantly reduce the computational effort, especially at large values of Nd. This is expected to be quite important for three-dimensional applications and is left to the examined in a future work. 6.2. Presentation of numerical results in general bathymetry domains To demonstrate the applicability of our non-linear model, in this section we discuss two particular examples, dealing with monochromatic waves propagating above a smooth but very steep shoal and above a trapezoidal bar. In both examples, a total number of 6 modes (n = − 2, − 1, 0, 1, 2, 3) have been retained in the local-mode series (Eq. (11)). In all cases a uniform grid has been used to approximate horizontal derivatives by central finite differences, and the grid resolution varies from about 70 (min) to 140 (max) gridpoints per (local) wavelength. The first example presented in Figs.7 and 8, concerns the case of harmonic incident waves, propagating over an infinite underwater step, characterised by shoaling ratio h2/h1 = 0.3, taken from Ohyama and Nadaoka (1994). In order treat the vertical wall of the underwater step in the framework of the present continuous CMS, a smooth but very steep bottom profile has been used, hðxÞ = h1 = 0:65−0:35 tanhðαðx−x Þ = h1 Þ; with α = 10:

ð33Þ

Fig. 8. Propagation of harmonic waves over a smooth but very steep shoal simulating an underwater step (characterised by h2/h1 = 0.3), as obtained by the present CMS. The frequency parameter in the deeper region is ω2h1/g = 0.8, and the corresponding nonlinearity parameter H/h1 = 0.1.

347

In this case, the frequency parameter in the deeper water region is ω2h1/g = 0.8 and the corresponding non-linearity parameter H/ h1 = 0.1. Although, in the region of incidence the wave non-linearity remains at low levels, it substantially increases (above H/h3 = 0.3) in the transmission region. The corresponding values of the shallowness parameter (estimated on the basis of linear dispersion relation) are λ1/h1 = 6.08, and λ1/h1 = 12.31, and the wave transformation is schematically depicted in Fig. 3 by using an arrow connecting two open circles. This example is a characteristic one where harmonic generation takes place over the shoal transferring significant energy to higher-harmonics. The transformation of the incident wave in a subregion around the position of the underwater step, for a whole wave period, is shown in Fig. 8, as calculated by the present CMS. This plot concerns the spatialtemporal evolution of the free-surface elevation, some periods after the wave front has passed the step. We observe in this figure the significant steepening of the wave-crests as they pass above the step and the appearance secondary crests, which is evidence of the generation of higher harmonics in the region of transmission. Comparison with the results by Ohyama and Nadaoka (1994, fig. 5) is presented in Fig. 9, concerning the spatial evolution of the first, second and third harmonics in the region, which are obtained by Fourier transforming the calculated free-surface elevation. Noticeable agreement is observed between the present method and the nonlinear, boundary integral equation method (see also Ohyama and Nadaoka, 1991). As a second example, numerical results obtained by the present CMS are presented in Fig. 10, concerning the propagation of harmonic waves of period T = 2 s and height H = 2 cm over the trapezoidal bar of Beji and Battjes (1994), for which experimental data are available. In this case, the depth ranges from a maximum h = 0.4 m, in the regions of incidence and transmission, to a minimum h = 0.1 m above the bar. The slope of the face of the trapezoidal shoal is 5% and of the back (downwave side) 10%; see fig. 1 in Beji and Battjes (1994). Comparisons between the results obtained by the present CMS and the experimental data are presented in Fig. 10a and b, for two positions: Station 2 (just before the top of the shoal) and Station 3 (just after the top of the shoal), and the agreement is considered to be very good. Finally, in Fig. 10c a snapshot of the calculated free-surface elevation is shown, at an instant some wave periods after the

Fig. 9. Same as in Fig. 8 concerning the spatial distribution of the first (crosses), second (circles) and third (dots) harmonics of the wave elevation over the underwater step, as calculated by Ohyama and Nadaoka (1994). The corresponding results obtained by the present CMS are shown by using solid line, dashed–dot line and dashed line, respectively.

348

K.A. Belibassakis, G.A. Athanassoulis / Coastal Engineering 58 (2011) 337–350


 f ðzÞ∈ C 2 ð J Þ∩C 1 J , where J = J∪∂J. Let us now define the following boundary operators Bη ðφÞ =

∂φðx; z = ηðx; t Þ; t Þ −μ0 φðx; z = ηðx; t Þ; t Þ; ∂z

ðA:1Þ

Bh ðφÞ =

∂φðx; z = −hðxÞ; t Þ ; ∂z

ðA:2Þ

involving the mixed derivative of the wave potential at the upper end z = η(x, t) and its vertical derivative at the bottom z = − h(x), respectively. Except of the case of linearised (infinitesimal amplitude) monochromatic waves of frequency ω = ω0, Bη(φ) is generally non-zero. Moreover, Bh(φ) is also non-zero, except for waves propagating in a uniform-depth strip (h(x) = h = const). These two quantities Bη(φ) and Bh(φ) are unknown, in the general case of nonlinear waves propagating in a variable bathymetry region. Let us define the freesurface and the sloping-bottom mode amplitudes (φn, n = − 2, − 1) as follows: Fig. 10. Evolution of incident harmonic waves over a trapezoidal bar. The period of the incoming waves is T = 2 s and waveheight H = 2 cm. Comparison of the present numerical solution with experimental data from Beji and Battjes (1994), shown in the upper two subplots (in Stations 2 and 3, respectively) by using dots. For clarity, the freesurface elevation shown in the last subplot has been 5-times exaggerated.

wavefront passes the shoal, showing the left-incident harmonic regular waves, under the combined effects of nonlinearity and dispersion, to transform over the face of the shoal to nonlinear shallow-water waves. 7. Conclusions The present non-linear CMS has been obtained without assumptions concerning the degree of retained nonlineariry and the vertical structure of the wave potential, being thus, equivalent with the complete water– wave formulation in the framework of potential flow. The theoretical value and practical effectiveness of the present model is that a small number of modes is found to be enough for numerical convergence, even in cases of very steep free-surface elevation and for general, but finite water depth, characterized by possibly large bottom slope. Extensive numerical evidence suggests that the rate of decay of the mode amplitudes is fast, φn = O(n− 4), and thus, in most applications, truncation of the modal series to its first few terms (retaining up to 5–6 modes) is sufficient for an accurate numerical solution. The present system is applied to the derivation of steady travelling waves over horizontal bottom, corresponding to various water depths, and its results are found to be consistent with fully nonlinear Fourier method, as well as with Stokes 5th-order and nonlinear cnoidal wave theory (in cases of limited non-linearity). Moreover, numerical results obtained by the present model in the case of waves propagating over general bottom topography (including shoals, underwater steps and bars) are compared with nonlinear methods based on boundary integral formulation and experimental data, showing good agreement. Finally, the analytical structure of the present model facilitates its extension to various directions, such as 3D waves and wave-body-seabed interaction problems in general bathymetry regions. Appendix A. The local-mode series expansion Consider the restriction f(z)of the wave potential Φ(x, z, t), at any vertical section x = const, and for any time instant (see, e.g., Fig. 1). From the properties of the solution of Laplace's equation in the present case, we expect that this function, defined on the vertical interval J = f−hðxÞbzbηðx; t Þg, is a smooth (continuously differentiable) one,

φ−2 ðx; t Þ = h0 Bη ðφÞ;

φ−1 ðx; t Þ = h0 Bh ðφÞ;

ðA:3Þ

where h0 is a scaling parameter that can be arbitrarily selected. Using Eqs. (A.1) and (A.2) in Eq. (A.3), we obtain that φ− 2(x, t), φ− 1(x, t) are continuously differentiable functions with respect to both x and t. The vertical structure of these modes, Eqs. (12a) and (12b), see Fig. A.1, 1 has been appropriately selected so that Bη(Z− 2) = h− 0 , Bh(Z− 2) = 0 1 and Bη(Z− 1) = 0, Bh(Z− 2) = h− , so that the reduced wave potential 0 φR ðx; z; t Þ = φðx; z; t Þ−φ−2 ðx; t ÞZ−2 ðz; η; hÞ−φ−1 ðx; t ÞZ−1 ðz; η; hÞ; ðA:4Þ is a smooth function, twice continuously differentiable in D, which, for all x and t, satisfies at the free-surface and at the bottom surface Bη ðφR Þ = 0;

Bh ðφR Þ = 0:

ðA:5Þ

The above conditions are sufficient to ensure that φR(x, z, t), at any vertical section and any time, can be represented in the form of the vertical eigenfunction expansion, ∞

φR ðx; z; t Þ = ∑ φn ðx; t ÞZn ðz; hðxÞ; ηðx; t ÞÞ;

ðA:6Þ

n=0

where the set of vertical functions fZn ðz; hðxÞ; ηðx; t ÞÞ; n = 0; 1; 2; 3:::g, see Fig. A.1, and the set of numbers fkn ; n = 0; 1; 2; 3…g, as given by Eqs. (12c) and (13), respectively, are obtained as the solution of the local (for each x and t) Sturm–Liouville problem: ∂2 Zn 2 −kn Zn = 0; −hðxÞbzbηðx; t Þ; ∂z2 ∂Zn ðz = ηðx; t ÞÞ Bη ðZn Þ = −μ0 Zn ðz = ηðx; t ÞÞ = 0; ∂z Bh ðZn Þ =

ðA:7Þ

∂Zn ðz = −hðxÞÞ = 0: ∂z

From the properties of regular eigenvalue problems (see e.g., Coddington and Levinson, 1955), the set of eigenfunctions fZn ðz; h; ηÞ; n = 0; 1; 2; 3; …g constitutes an L2 − basis in the z-interval − h(x) b z b η(x, t), for each x and t. Moreover, since the function φR(x, z ; t) fulfils the same boundary conditions, Eq. (A.5), as the eigenfunctions do, Eqs. (A.8b), (A.8c), the series (A.6) converges uniformly to the function φR(x, z ; t) in − h(x) ≤ z ≤ η(x, t). Then, by substituting Eq. (A.6) to Eq. (A.4), and passing the newly introduced modes φnZn, n = − 2,− 1, to the left-hand side, we obtain the present enhanced local-mode representation, as given by Eq. (11).

K.A. Belibassakis, G.A. Athanassoulis / Coastal Engineering 58 (2011) 337–350

We remark here that from the theory of spectral methods (see, e.g., Gottlieb and Orszag, 1977; Boyd, 2001) the local mode series defined by the right-hand side of Eq. (A.6), without the additional terms φnZn, n = − 2, − 1, exhibits slow-convergence characterised by the decay 2 rate φn = O(k− n ), in the general case, since the wave potential does not satisfy the same boundary conditions as the vertical eigenfunctions fZn ; n = 0; 1; 2; 3…g do. On the basis of Eq. (13), it holds kn (η + h) ≈ nπ, for large n, and thus φn = O(n− 2). This is remedied by the enhanced local-mode series (Eq. (11)) containing the additional terms, which converges much faster, and the decay rate of the mode amplitudes becomes φn = O(n− 4). The last result is obtained by using integration by parts in the vertical interval, in conjunction with the definition of the spectral coefficients and the properties of the solution of the present problem. This fact implies also convergence of the corresponding series obtained by termwise differentiation of the local-mode series (Eq. (11)) up to and including the (free surface and bottom) boundaries. Appendix B. Coefficients of the nonlinear coupled-mode system ^mn ðηÞ; B ^mn ðηÞ involved in Eqs. (18a) and ^mn ðηÞ; C The coefficients A (19) are defined as follows:

^ ðηÞ = A mn

usual inner product in the vertical interval − h(x) b z b η(x, t), z = ηðx;t Þ

〈f ; g〉 =

∫

ð−1 + Zn ðz; h; ηÞZm ðz; h; ηÞÞdz;

〈

〉


 ^ ðηÞ = 2 ∂Zn ; Z + ∂h −1 + ½Z Z  B mn n m z = −h ; ∂x m ∂x

^ ðηÞ = 〈ΔZ ; Z 〉 + C mn n m

z = ηðx;t Þ

^ aℓn ðηÞ = 〈Zn ; Wℓ 〉 =

Zn ðz; h; ηÞWℓ ðz; h; ηÞdz;

ðB:5Þ

 ∂Zn ∂h ∂Zn ^ bℓn ðηÞ = 2 ½Zn Wℓ z = −h − ; ; Wℓ + ∂x ∂x z = η ∂x

ðB:6Þ

h

  ∂h ∂Zn ∂Zn 1 ∂Zn ∂Zℓ − ; + Wℓ 2 ∂z ∂z z = η ∂z ∂x ∂x z = −h ðB:7Þ

where the functions Wℓ are calculated from the vertical basis Zn (defined by Eqs. (12a), (12b), and (12c)), as follows:

References

Fig. A.1. Vertical structure of the modes n = − 2, − 1, 0, 1, 2, 3, for μ0h = μ0h0 = 1 and H/ h = 0.4 (under the crest). Observe that all vertical modes satisfy zero vertical derivative at the bottom, except the sloping-bottom mode, which for this particular selection of parameters is linear in z.

i



ðB:2Þ

where ΔZn = ∂ 2Zn/∂ x2 + ∂ 2Zn/∂ z2, and the angle brackets denote the

∫ z = −hðxÞ

^cℓn ðηÞ = 〈ΔZn ; Wℓ 〉 +

Wn ðz; h; ηÞ =

ðB:3Þ

ðB:4Þ

bmn ðηÞ; ^cmn ðηÞ, involved Moreover, the matrix-coefficients ^ amn ðηÞ;^ in Eq. (18b), are also dependent on the free-surface elevation (η), and are defined as follows:

ðB:1Þ

    ∂h ∂Zn ∂Zn ∂η ∂Zn ∂Zn + Zm − Zm + ; ∂x ∂x ∂x ∂x ∂z ∂z z = −h z=η

f ðzÞg ðzÞdz:

∫ z = −hðxÞ

z = ηðx;t Þ z = −hðxÞ

349

∂Zn ðz; h; ηÞ : ∂η

ðB:8Þ

Athanassoulis, G.A., Belibassakis, K.A., 1999. A consistent coupled-mode theory for the propagation of small-amplitude water waves over variable bathymetry regions. J. Fluid Mech. 389, 275–301. Athanassoulis, G.A., Belibassakis, K.A., 2002. A Nonlinear Coupled-mode Model for Water Waves over a General Bathymetry. Proc. 21st Int. Conf. on Offshore Mechanics and Arctic Engineering. OMAE2002, Oslo, Norway 2002. Athanassoulis, G.A., Belibassakis, K.A., 2007. A coupled-mode method for non-linear water waves in general bathymetry with application to steady travelling solutions in constant, but arbitrary, depth. J. Discrete Cont. Dyn. Syst. DCDS B 75–84 special volume of selected papers from 6th Int. Conference on Dynamical Systems and Differential Equations, Poitiers Meeting, June 25–28, 2006. Baldock, T.E., Swan, C., 1994. Numerical calculations of large transient water waves. Appl. Ocean Res. 16, 101–112. Beji, S., Battjes, J.A., 1994. Numerical simulation of wave propagation over a bar. Coast. Eng. 23, 1–16. Beji, S., Nadaoka, K., 1997. A time-dependent nonlinear mild-slope equation for water waves. Proc. R. Soc. Lond. A 453, 319–332. Belibassakis, K.A., Athanassoulis, G.A., 2002. Extension of second-order Stokes theory to variable bathymetry. J. Fluid Mech. 464, 35–80. Belibassakis, K.A., Athanassoulis, G.A., 2006. A coupled-mode technique for weakly nonlinear wave interaction with large floating structures lying over variable bathymetry regions. Appl. Ocean Res. 28, 59–76. Belibassakis, K.A., Athanassoulis, G.A., Gerostathis, Th., 2001. A coupled-mode model for the refraction-diffraction of linear waves over steep three-dimensional bathymetry. Appl. Ocean Res. 23 (6), 319–336. Berenger, J.-P., 1994. A perfectly matched layer for the absorption of electromagnetic waves. J. Comput. Phys. 114, 185–200. Bingham, H., Zhang, H., 2007. On the accuracy of finite-difference solutions for nonlinear water waves. J. Eng. Math. 58, 211–228. Bingham, H.B., Madsen, P.A., Fuhrman, D.R., 2009. Velocity potential formulations of Boussinesq-type models. Coastal Eng. 56, 467–478. Boyd, J.P., 2001. Chebyshev and Fourier Spectral Methods. Dover, New York. Christou, M., Hague, C.H., Swan, C., 2009. The reflection of nonlinear irregular surface water waves. Eng. Anal. Bound. Elem. 33 (5), 644–653. Clamond, D., Grue, J., 2001. A fast method for fully nonlinear water wave computations. J. Fluid Mech. 447, 337–355. Coddington, E.A., Levinson, N., 1955. Theory of Ordinary Differential Equations. McGraw Hill, New York. Coifman, R., Meyer, Y., 1985. Nonlinear harmonic analysis and analytic dependence. Pseudodifferential operators and applications. Proc. Symp. Pure Math. Am. Math. Soc. Providence RI 43, 71–78. Cokelet, E.D., 1977. Steep gravity waves in water of arbitrary uniform depth. Philos. Trans. R. Soc. Lond. A Math. Phys. Sci. 286, 183–230. Collino, F., Monk, P., 1998. Optimizing the perfectly matched layer. Comput. Meth. Appl. Mech. Eng. 164, 157–171. Craig, W., 1985. An existence theory for water waves and the Bousinesq and the Korteweg-de Vries scaling limits. Comm. Partial Differ. Equ. 10, 787–1003.

350

K.A. Belibassakis, G.A. Athanassoulis / Coastal Engineering 58 (2011) 337–350

Craig, W., Nicholls, D.P., 2002. Traveling gravity water waves in two and three dimensions. Eur. J. Mech. B. Fluids 21, 615–641. Craig, W., Sulem, C., 1993. Numerical simulation of gravity waves. J. Comput. Phys. 108, 73–83. Craig, W., Guyenne, P., Nicholls, D., Sulem, C., 2005. Hamiltonian long-wave expansions for water waves over a rough bottom. Proc. R. Soc. Lond. A 461, 839–873. Dalrymple, R., Rogers, B.D., 2006. Numerical modelling of water waves with the SPH method. Coastal Eng. 53 (2–3), 141–147. Dalrymple, R., Rogers, B.D., et al., 2007. Smoothed Particle Hydrodynamics for Water Waves. Proc. 26th Int. Conf. on Offshore Mechanics and Arctic Engineering. OMAE2007, San Diego California 2007. Dias, F., Bridges, Th., 2006. The numerical computation of freely propagating timedependent irrotational water waves. Fluid Dyn. Res. 38, 803–830. Dias, F., Kharif, C., 1999. Nonlinear gravity and capillary-gravity waves. Annu. Rev. Fluid Mech. 31, 301–346. Diaz, J., Joly, P., 2006. A time domain analysis of PML models in acoustics. Comput. Meth. Appl. Mech. Eng. 195, 3820–3853. Dingemans, M., 1997. Water Wave Propagation over Uneven Bottoms. World Scientific. Dommermuth, D.G., Yue, D.K.P., 1987. A high-order spectral method for the study of nonlinear gravity waves. J. Fluid Mech. 184, 267–288. Drennan, W.M., Hui, W.H., Tenti, G., 1988. Accurate calculations of Stokes water waves of large amplitude. Z. Angew. Math. Phys. ZAMP 43, 367–384. Engsig-Karup, A.P., Bingham, H.B., Lindberg, O., 2009. An efficient flexible-order model for 3D nonlinear water waves. J. Comput. Phys. 228, 2100–2118. Fenton, J.D., 1990. Nonlinear wave theories. In: Mehaute, B., Hanes, D. (Eds.), The Sea, Vol. 9A. J. Wiley & Sons. Gerostathis, T., Belibassakis, K.A., Athanassoulis, G.A., 2008. A coupled-mode model for the transformation of wave spectrum over steep 3d topography. A parallelarchitecture implementation. J. Offshore Mech. Arct. Eng. 130. Gottlieb, D., Orszag, S., 1977. Numerical Analysis of Spectral Methods: Theory and Applications. Soc. Ind. Appl. Mathematics (SIAM). Grilli, S.T., Skourup, J., Svendsen, I.A., 1989. An efficient boundary element method for nonlinear water waves. Eng. Anal. Bound. Elem. 6 (2), 97–107. Grilli, S.T., Guyenne, P., Dias, F., 2001. A fully non-linear model for three-dimensional overturning waves over an arbitrary bottom. Int. J. Numer. Methods Fluids 35, 829–867. Groves, J., Toland, M., 1997. On variational formulations of steady water waves. Arch. Ration. Mech. Anal. 137, 203–226. Grue, J., 2002. On four highly nonlinear phenomena in wave theory and marine hydrodynamics. Appl. Ocean Res. 24, 261–274. Hu, F.Q., 2001. A stable, perfectly matched layer for linearized Euler equations in unsplit physical variables. J. Comput. Phys. 173 (2), 455–480. Jamois, E., Fuhrman, D.R., Bingham, H.B., Molin, B., 2006. A numerical study of nonlinear wave run-up on a vertical plate. Coastal Eng. 53, 929–945. Johannessen, T.B., Swan, C., 1997. Nonlinear transient water waves — part I. A numerical method of computation with comparisons to 2D laboratory data. Appl. Ocean Res. 19, 293–308. Kim, J.W., Bai, K.J., Ertekin, R.C., Webster, W.C., 2001. A derivation of the Green–Naghdi equations for irrotational flows. J. Eng. Math. 40 (1), 17–34. Kim, J.W., Bai, K.J., Ertekin, R.C., Webster, W.C., 2003. A strongly-nonlinear model for water waves in water of variable depth: the irrotational Green–Naghdi model. J. Offshore Mech. Arct. Eng. 125 (1), 25–32. Liu, P.L.-F., 1995. Model equations for wave propagation from deep to shallow water. In: Liu, P.L.-F. (Ed.), Advances in Coastal and Ocean Engineering. World Scientific. Longuet-Higgins, M.S., 1975. Integral properties of periodic gravity waves of finite amplitude. Proc. R. Soc. Lond. A 342, 157–174. Luke, J.C., 1967. A variational principle for a fluid with a free surface. J. Fluid Mech. 27, 395–397.

Madsen, P.A., Bingham, H.B., Liu, H., 2002. A new Boussinesq method for fully nonlinear waves from shallow to deep water. J. Fluid Mech. 462, 1–30. Madsen, P.A., Bingham, H.B., Schafer, H.A., 2003. Boussinesq-type formulations for fully nonlinear and extremely dispersive water waves: derivation and analysis. Proc. R. Soc. Lond. A 459, 1075–1104. Massel, S., 1989. Hydrodynamics of Coastal Zones. Elsevier. Massel, S., 1993. Extended refraction–diffraction equations for surface waves. Coast. Eng. 19, 97–126. Nadaoka, K., Beji, S., Nakagawa, Y., 1997. A fully dispersive weakly nonlinear model for water waves. Proc. R. Soc. Lond. A 453, 303–318. Naumkin, P.I., Shishmarev, I.A., 1994. Nonlinear nonlocal equations in the theory of waves. Mathematical Monographs, Vol. 133. American Mathematical Society. Ohyama, T., Nadaoka, K., 1991. Development of a numerical wave tank for analysis of nonlinear and irregular wave field. Fluid Dyn. Res. 8, 231–251. Ohyama, T., Nadaoka, K., 1994. Transformation of a nonlinear wave train passing over a submerged shelf without breaking. Coast. Eng. 24, 1–22. Okamoto, H., Shoji, M., 2001. The mathematical theory of permanent progressive water waves. Advanced Series in Nonlinear Dynamics, Vol.20. World Scientific. Peregrine, D.M., 1967. Long waves on a beach. J. Fluid Mech. 27, 815–827. Press, W.H., Flannery, B.P., Teukolsky, S.A., Vetterling, W.T., 1992. Numerical Recipes in Fortran. Cambridge University Press. Rienecker, M.M., Fenton, J.D., 1981. A Fourier approximation method for steady water waves. J. Fluid Mech. 104, 119–137. Schwartz, L.W., 1974. Computer extension and analytic continuation of Stokes expansion for gravity waves. J. Fluid Mech. 62, 553–578. Shampine, L.F., Reichelt, M.W., 1997. The MATLAB ODE Suite. SIAM J. Sci. Comput. 18, 1–22. Sobey, R.J., 1992. A local Fourier approximation method for irregular wave kinematics. Appl. Ocean Res. 14, 93–105. Tang, Y., Ouellet, Y., 1997. A new kind of nonlinear mild-slope equation for combined refraction–diffraction of multifrequency water waves. Coast. Eng. 31, 3–36. Tao, L., Song, H., Chakrabarti, S., 2007. Nonlinear progressive waves in water of finite depth — an analytic approximation. Coastal Eng. 54, 825–834. Trulsen, K., Dysthe, K.B., 1996. A modified nonlinear Schrödinger equation for broader bandwidth gravity waves on deep water. Wave Motion 24, 281–289. Trulsen, K., Gudmestad, O.T., Velarde, M.G., 2001. The nonlinear Schrödinger method for water wave kinematics on finite depth. Wave Motion 33, 379–395. Tsai, W., Yue, D.K., 1996. Computation of nonlinear free-surface flows. Annu. Rev. Fluid Mech. 28, 249–278. Turkel, E., Yefet, A., 1998. Absorbing PML boundary layers for wave-like equations. Appl. Numer. Math. 27 (4), 533–557. Turnbull, M.S., Borthwick, A., Eatock, Taylor R., 2003. Wave–structure interaction using coupled structured–unstructured finite element meshes. Appl. Ocean Res. 25, 63–77. Whitham, G.B., 1967. Variational method and applications to water wave. Proc. R. Soc. Lond. 229, 6–25. Whitham, G.B., 1974. Linear and Nonlinear Waves. Wiley, NY. Williams, J.M., 1981. Limiting gravity waves in water of finite depth. Proc. R. Soc. Lond. A 302, 139–188. Wu, S., 1999. Well-posedness in Sobolev spaces of the full water wave problem in 3D. J. Am. Math. Soc. 12, 445–495. Wu, G.X., Eatock Taylor, R., 1994. Finite element analysis of two-dimensional non-linear transient water waves. Appl. Ocean Res. 16, 363–372. Wu, G.X., Eatock Taylor, R., 1995. Time stepping solutions of the two dimensional nonlinear wave radiation problem. Ocean Eng. 22 (8), 785–798. Zufiria, J.A., 1987. Weakly non-linear non-symetric gravity waves on water of finite depth. J. Fluid Mech. 180, 371–385.