# Second-Order Time-Dependent Mild-Slope Equation for Wave Transformation

Authors: Ching-Piao Tsai; Hong-Bin Chen; John R. C. Hsu
Journal: Mathematical Problems in Engineering
 (2014)

---

## Abstract

This study is to propose a wave model with both wave dispersivity and nonlinearity for the wave field without water depth restriction. A narrow-banded sea state centred around a certain dominant wave frequency is considered for applications in coastal engineering. A system of fully nonlinear governing equations is first derived by depth integration of the incompressible Navier-Stokes equation in conservative form. A set of second-order nonlinear time-dependent mild-slope equations is then developed by a perturbation scheme. The present nonlinear equations can be simplified to the linear time-dependent mild-slope equation, nonlinear long wave equation, and traditional Boussinesq wave equation, respectively. A finite volume method with the fourth-order Adams-Moulton predictor-corrector numerical scheme is adopted to directly compute the wave transformation. The validity of the present model is demonstrated by the simulation of the Stokes wave, cnoidal wave, and solitary wave on uniform depth, nonlinear wave shoaling on a sloping beach, and wave propagation over an elliptic shoal. The nearshore wave transformation across the surf zone is simulated for 1D wave on a uniform slope and on a composite bar profile and 2D wave field around a jetty. These computed wave height distributions show very good agreement with the experimental results available.


---

## Body

## 1. Introduction

Accurate estimation of design waves, especially wave heights and their spatial distributions, is essential for the planning and installation of coastal structures, calculation of longshore sediment transport, and prediction of nearshore topographical changes. In practical engineering applications, a narrow-banded sea state centred on a certain dominant wave frequency has remained popular for several decades, despite the use of random wave theories in more recent time. Because wave transformation from deep water into shallow region involves nonlinear process due to shoaling, combined refraction and diffraction, and reflection and breaking, various forms for mild-slope equations have been developed. Previously, wave transformations involving refraction and diffraction had been calculated separately for simple beach topography, until Berkhoff [1] who presented a linear mild-slope equation for the combined effects. Because the original mild-slope equation was in an elliptic type, it was not efficient in performing direct numerical calculation and also in treating the open boundary conditions. The conventional schemes also required large computer storage and computing time. Later, numerical models, based on an elliptic form of the mild-slope equation, were proposed to compute the wave field for a large coastal area [2–5]. These models obtained the solution for the governing equations by iterations using powerful and sophisticated algorithm.The mild-slope equation was also cast in a parabolic form to promote computational efficiency [6] and to handle the problem with weak wave reflection [7]. However, the parabolic approximation method was found inadequate for determining the wave field in the presence of a highly reflective structure. Later, Copeland [8] expressed the mild-slope equation as a pair of first-order hyperbolic equations, without the loss of the reflected wave. Watanabe and Maruyama [9] extended the set of time-dependent mild-slope equations including a wave energy dissipation factor for modelling the nearshore wave field. Tsai et al. [10] incorporated both the nonlinear shoaling factor and the breaking-wave energy dissipation factor into the time-dependent mild-slope equations to improve the accuracy of the calculation of wave transformation across the surf zone. Taking large bottom gradient into the mild-slope equation, the so-called extended mild-slope equation or the modified mild-slope equation were also developed for waves over rapid variation of bottom [11, 12]. On the other hand, the complementary mild-slope equation was derived in terms of a stream function vector rather than in terms of a velocity potential [13, 14]. All of these are in the family of linear mild-slope equations.Based on the second-order Stokes wave theory, the nonlinear mild-slope equation with parabolic type was developed [15, 16]. These theories are adaptable to the relatively deeper water. Belibassakis and Athanassoulis [17] proposed a coupled-mode theory, by using a local mode expansion of the wave potential in conjunction with Luke variational principle, to extend the second-order Stokes theory for the case of a generally shaped bottom profile connecting two half-strips of constant depths; their model was then extended with application to nonlinear water waves propagating in variable bathymetry regions [18]. On the other hand, Boussinesq-type equations developed by Peregrine [19] were extended to cover varying depth adaptable to the shallow water [20]. The improved Boussinesq-type equations can provide the wave dispersion and shoaling characteristics from shallow water to water depth with kh = 40, where k is the wave number and h is the water depth [21]. Recently, Li [22] derived nonlinear wave equations based on combining the second-order Stokes wave theory with the Boussinesq theory.Many investigators have tried to develop nonlinear wave equations adaptable to arbitrary water depth. For instance, a weakly nonlinear and fully dispersive model was introduced by Nadaoka et al. [23], who generated a set of equations that rendered a single time-dependent nonlinear mild-slope equation for water waves [24]. The numerically simulated cnoidal waves and the second-order Stokes waves were found in excellent agreement with theoretical solutions. However, their nonlinear model was only simulated to unidirectional waves, rather than a two-dimensional wave field involving waves combining refraction-diffraction and breaking.In this paper, a new nonlinear wave equation is developed for modelling wave field covering the entire range from deep to shallow water and involving wave breaking. It starts from the depth-integration of the continuity equation and nonlinear momentum equation. A set of fully dispersive and weakly nonlinear wave equations is then developed by introducing a perturbation technique. The theoretical derivation is outlined in Section2. In Section 3, we introduce the finite volume method with a fourth-order Adams-Moulton predictor-corrector technique, rather than the usual finite difference and finite element schemes, for directly computing the nonlinear equations. Finally, the present model is verified against Stokes wave and cnoidal wave propagation on uniform depth and nonlinear wave shoaling on a sloping beach. Examples for wave propagation on an elliptic shoal and wave transformation across the surf zone are also presented.


## 2. Nonlinear Wave Equations

### 2.1. Governing Equations

Consider a three-dimensional wave motion in an incompressible, inviscid, and homogeneous fluid, the principle of conservation of mass gives the following continuity equation:(1)∇·u→+∂w∂z=0
and the incompressible Navier-Stokes equations in conservative form lead to Euler’s equations as
(2)∂u→∂t+u→·∇u→+w∂u→∂z=-1ρ∇p,(3)∂w∂t+u→·∇w+w∂w∂z=-1ρ∂p∂z-g
in the horizontal and vertical directions of a Cartesian coordinates system, respectively. In (1)–(3), u→=(u,v) and w are the Eulerian velocity components in the horizontal plane and vertical direction, respectively; p is the pressure, ρ is the density of the fluid, g is the acceleration due to gravity, t is the time, and ∇=(∂/∂x,∂/∂y) is the horizontal vector operator, where x, y, and z are the Cartesian coordinates.The fluid motion described by (1)–(3) has to satisfy boundary conditions at the free surface and on the seabed. First, the dynamic and kinematic boundary conditions at the free surface are
(4)p=0z=η(x,y,t),(5)w=∂η∂t+u→·∇ηz=η(x,y,t)
and the bottom boundary condition on an impermeable seabed is
(6)w=-u→·∇hz=-h(x,y)
in which η(x,y,t) is the water surface elevation and h(x,y) the local water depth.


### 2.2. Derivation of the Nonlinear Wave Equation

To develop a nonlinear wave equation, we perform depth-integration on the governing equations (1)–(3) separately and then substitute the boundary conditions (4)–(6) into the appropriate resultant forms. First, upon integrating (1) vertically from seabed to the free surface and substituting the boundary conditions given by (5) and (6), it yields
(7)∂η∂t+∇·∫-hηu→dz=0.
Second, after integrating the momentum equation (3) in the vertical direction, from an arbitrary depth z to the free surface, and inserting the boundary condition for the pressure at the free surface, (4), it produces
(8)pρ=g(η-z)+∫zη(∂w∂t)dz+∫zη(u→·∇w+w∂w∂z)dz.Substituting (8) into (2) renders the momentum equation for the x- and y-axis in the horizontal directions, respectively, as follows:
(9)∂u⃑∂t+g∇η+(u→·∇u→+w∂u→∂z)=-∇∫zη(∂w∂t)dz-∇∫zη(u→·∇w+w∂w∂z)dz.It is worth noting that the convective terms,(u→·∇u→+w∂u→/∂z) and (u→·∇w+w∂w/∂z) in (8) and (9), are transformed to ∇(u→·u→+w2) and ∂(u→·u→+w2)/∂z in Nadaoka et al. [23], based on the irrotational assumption before the integration process described above. However, such transformations do not completely satisfy the irrotationality for variable depth, that is, u→z-∇w≠0, which was mentioned in [23]. Subsequently, the present equation (9) is retained in the original forms of the nonlinear convective terms; that is, we preserve, completely, the irrotational condition. It is noted that Belibassakis and Athanassoulis [18] preserved this flow irrotationality by working with a potential formulation; they derived a nonlinear coupled-mode system of equations on the horizontal plane, with respect to the free surface elevation and the unknown modal amplitudes, in conjunction with Luke variational principle.The horizontal velocity vector can be expressed in terms of a vertical distribution function and a corresponding horizontal velocity vector at the still water level, in the form of(10)u⃑(x,y,z,t)=∑m=1Ncosh⁡km(h+z)cosh⁡kmhU→m(x,y,t),
where U→m(x,y,t) is the mth component of the horizontal velocity vector at z=0, and km(x,y) is the local wave number. The full dispersivity for a broad-banded wave field can be attained by taking only a few terms [23]. However, considering the fact that a narrow-banded sea state centred around a certain dominant wave frequency may be described with sufficient accuracy by the single-component model equations [24], that is, N=1, it becomes
(11)u⃑(x,y,z,t)=F(kh,z)U→(x,y,t),ffffffF(kh,z)=cosh⁡k(h+z)cosh⁡kh.
This single-component model has remained popular and useful in practical coastal engineering applications.The flow velocityw in the vertical direction z can be obtained upon depth integration of (1) from -h to z and the bottom boundary condition equation (6), rendering
(12)w(x,y,z,t)=-∇·[f(kh,z)kU→(x,y,t)],ffffffffffff(kh,z)=sinh⁡k(h+z)cosh⁡kh.Inserting (11) into (7) gives the following depth-integrated continuity equation:
(13)∂η∂t+∇·[(∫-hηFdz)U→]=0.
Substituting (11) and (12) into (8), and after laborious manipulations, we obtain the pressure expression as
(14)pρ=g(η-z)-p~a(∇·∂U→∂t)-p~b(∇h·∂U→∂t)-p~c[(U→·∇)(∇·U→)-(∇·U→)(∇·U→)]-p~d[(∇·U→)(∇h·U→)]+O[(∇h)2],
where
(15)p~a=∫zη(fk)dz,p~b=∫zη∂∂h(fk)dz,p~c=∫zηF(fk)dz,p~d=∫zη[F∂∂h(fk)-∂F∂h(fk)]dz.
The last term on the right-hand-side (RHS) of (14) is the higher order term of the bottom slope which may be omitted under the assumption of mild-slope.Multiplying (9) by the vertical distribution function F in (11) and then integrating it from z=-h to η, as well as after complex mathematical manipulations, produces a depth-integrated expression for the momentum equation in the form of
(16)α~1∂U→∂t+α~2g∇η+α~3(U→·∇)U→+α~4(∇·U→)U→+α~5(∇h·U→)U→=α~6∇(∇·∂U→∂t)+α~7(∇·∂U→∂t)∇h+α~8(∇·∂U→∂t)∇η+α~9(∇h·∇η)∂U→∂t+α~10∇[(U→·∇)(∇·U→)-(∇·U→)(∇·U→)]+α~11[(∇·U→)(∇h·∇)U→]+α~12[(U→·∇)(∇·U→)-(∇·U→)(∇·U→)]∇η+α~13[(∇·U→)(∇h·U→)∇η]+O[(∇h)2]
in which
(17)α~1=∫-hηF2dz,α~2=∫-hηFdz,α~3=∫-hηF3dz,α~4=-∫-hηFf2dz,α~5=∫-hη[F2∂F∂h-kFf∂∂h(fk)]dz,α~6=∫-hηFp~adz,α~7=2∫-hηFp~bdz,α~8=∫-hηF∂p~a∂ηdz,α~9=∫-hηF∂p~b∂ηdz,α~10=∫-hηFp~cdz,α~11=2∫-hηFp~ddz,α~12=∫-hηF∂p~c∂ηdz,α~13=∫-hηF∂p~d∂ηdz.The momentum equation, (16), retains a full nonlinear form under the mild-slope condition. As it is difficult to solve (16) directly by a numerical scheme; a perturbation technique is then proposed. The surface displacement and velocity field are expressed by
(18)η=εη1+ε2η2+⋯,U→=εU→1+ε2U→2+⋯,
where ε is the dimensionless nonlinearity parameter of wave amplitude, and the subscriptions i(i=1,2,…) in (18) denote the ith order approximation.Substituting (18) into (13) and (16) and taking Taylor series expansion on the still water level z=0, we obtain the continuity equation and the momentum equation, respectively, as follows:
(19)∂η∂t+∇·[(∫-h0Fdz+η)U→]=O(ε3),(20)(∫-h0F2dz)∂U→∂t+(∫-h0Fdz)g∇η+12[k2c4g2∫-h0F(F2-f2)dzmmmmcmm-2∫-h0Fdzk2c4g2]∇(U→·U→)=(∫-h0F-F2k2dz)∇(∇·∂U→∂t)+[∫-h0F∂∂h(1-Fk2)dz](∇·∂U→∂t)∇h+O(ε3).
These two equations are accurate to the second-order of ε. Substituting the expressions of F and f given in (11) and (12) into (19) and (20), we obtain the final expression of the continuity equation and momentum equation as follows:
(21)∂η∂t+∇·[(c2g+η)U→]=0,(22)n∂U→∂t+g∇η+12(1-3k2c4g2)∇(U→·U→)=1c2∇[c2(1-n)k2(∇·∂U→∂t)]
in which
(23)c2=gktanhkh,(24)n=12[1+2khsinh⁡(2kh)].The coefficient to the third term∇(U→·U→) on the left-hand-side of (22) consists of two parts. The first part, 1/2, is derived from the nonlinear convective terms of the horizontal velocity; and the second part, -3/2(k2c4/g2), is the nonlinear effect of the vertical velocity. The value of (k2c4/g2) tends to 1 in deep water, implying that the nonlinear variation of vertical velocity is in the same order as the horizontal velocity and it cannot be neglected, whereas it approximates reasonably to zero in the very shallow water. It is worth noting that the c and k can be computed for a prescribed wave period and a local depth based on the dispersion relationship shown in (23). The dispersion characteristics are essentially the same as in [23] because both models are approximated to the second-order of ε.Equations (21) and (22) together become the so-called second-order time-dependent mild-slope equation, because they can be used to represent the wave fields that are weakly nonlinear and fully dispersive and on the basis of the mild-slope assumption. Differing from [23], the present equations (21) and (22) can be directly used to calculate the wave transformation of wave field in plane because the vertical velocity has been embedded in the nonlinear convective terms into a simpler wave model. The consistency of the equation can be simplified into several well-established theories, such as linear mild-slope equation, nonlinear long wave equation, and Boussinesq equation. The process is demonstrated below.


### 2.3. Linear Mild-Slope Equation

Upon dropping the nonlinear term in (21) for the sake of linear (Airy) wave theory, this equation can be rewritten as
(25)∂η∂t+∇·Q⃑=0,
where Q→ is the linear flow rate of the waves defined by
(26)Q→=∫-h0u→dz=c2gU→.
On the other hand, (22) can be rewritten for the linear waves as
(27)n∂U→∂t+g∇η=1c2∇[c2(1-n)k2(∇·∂U→∂t)].
Substituting (26) into (27), then we obtain
(28)∂Q→∂t+c2n∇η=1n∇[(1-n)k2∇·∂Q→∂t-(1-n)c2k2∇c2·∂Q→∂t].
In the right-hand side of (28),
(29)1n∇[(1-n)k2(∇·∂Q→∂t)]=-1n∇[(1-n)k2∂2η∂t2]=1n∇[c2(1-n)η]=(1-n)nη∇c2+c2n∇η-c2n∇(nη),-1n∇[(1-n)c2k2∇c2·∂Q→∂t]=-(1-n)nc2k2∇c2(∇·∂Q→∂t)+O[(∇h)2]=-(1-n)nη∇c2+O[(∇h)2].
The term of O[(∇h)2] can be omitted by the assumption of mild-slope bottom. Upon inserting (29) into (28), the nonlinear equation developed in this study is reduced to
(30)∂Q⃑∂t+c2n∇(nη)=0.
Equations (25) and (30) are identical to the so-called time-dependent mild-slope equation in Copeland [8]. Eliminating flow rate Q→ from (25) and (30) renders the elliptic mild-slope equation in Berkhoff [1].


### 2.4. Nonlinear Long Wave Equation

For long wave theory in shallow water, wave celerity (c), group velocity (cg), and flow rate (Q→) are usually given by
(31)c=cg=gh,Q→=(h+η)U→,
where U→ is the horizontal velocity identical to mean flow velocity in shallow water condition. By dropping the nonlinear terms in the vertical direction in (22) and introducing (31), a system of nonlinear long wave equations can be written in terms of the mean flow velocity U→, given by
(32)∂U→∂t+g∇η+12∇(U→·U→)=0.
For an irrotational flow, the term ∇(U→·U→)/2 can be replaced by (U→·∇)U→, thus rendering the usual expression for the nonlinear long wave equation.


### 2.5. Boussinesq Equations

Boussinesq equations account for the effect of wave nonlinearity but weak dispersion, that the wave celerity and group velocity may be approximated by the first two terms in the Taylor series,(33)c2=gh(1-13k2h2),n=cgc=(1-13k2h2).
Substituting (33) into (21) and (22) and replacing k2U→ by -∇(∇·U→) for the irrotational flow in shallow water condition, the present nonlinear mild-slope equation gives
(34)∂η∂t+∇·[(h+η)U→]+h33∇2(∇·U→)=0,∂U→∂t+g∇η+12∇(U→·U→)=0.
Equations (34) are the standard form of traditional Boussinesq equations.


## 2.1. Governing Equations

Consider a three-dimensional wave motion in an incompressible, inviscid, and homogeneous fluid, the principle of conservation of mass gives the following continuity equation:(1)∇·u→+∂w∂z=0
and the incompressible Navier-Stokes equations in conservative form lead to Euler’s equations as
(2)∂u→∂t+u→·∇u→+w∂u→∂z=-1ρ∇p,(3)∂w∂t+u→·∇w+w∂w∂z=-1ρ∂p∂z-g
in the horizontal and vertical directions of a Cartesian coordinates system, respectively. In (1)–(3), u→=(u,v) and w are the Eulerian velocity components in the horizontal plane and vertical direction, respectively; p is the pressure, ρ is the density of the fluid, g is the acceleration due to gravity, t is the time, and ∇=(∂/∂x,∂/∂y) is the horizontal vector operator, where x, y, and z are the Cartesian coordinates.The fluid motion described by (1)–(3) has to satisfy boundary conditions at the free surface and on the seabed. First, the dynamic and kinematic boundary conditions at the free surface are
(4)p=0z=η(x,y,t),(5)w=∂η∂t+u→·∇ηz=η(x,y,t)
and the bottom boundary condition on an impermeable seabed is
(6)w=-u→·∇hz=-h(x,y)
in which η(x,y,t) is the water surface elevation and h(x,y) the local water depth.


## 2.2. Derivation of the Nonlinear Wave Equation

To develop a nonlinear wave equation, we perform depth-integration on the governing equations (1)–(3) separately and then substitute the boundary conditions (4)–(6) into the appropriate resultant forms. First, upon integrating (1) vertically from seabed to the free surface and substituting the boundary conditions given by (5) and (6), it yields
(7)∂η∂t+∇·∫-hηu→dz=0.
Second, after integrating the momentum equation (3) in the vertical direction, from an arbitrary depth z to the free surface, and inserting the boundary condition for the pressure at the free surface, (4), it produces
(8)pρ=g(η-z)+∫zη(∂w∂t)dz+∫zη(u→·∇w+w∂w∂z)dz.Substituting (8) into (2) renders the momentum equation for the x- and y-axis in the horizontal directions, respectively, as follows:
(9)∂u⃑∂t+g∇η+(u→·∇u→+w∂u→∂z)=-∇∫zη(∂w∂t)dz-∇∫zη(u→·∇w+w∂w∂z)dz.It is worth noting that the convective terms,(u→·∇u→+w∂u→/∂z) and (u→·∇w+w∂w/∂z) in (8) and (9), are transformed to ∇(u→·u→+w2) and ∂(u→·u→+w2)/∂z in Nadaoka et al. [23], based on the irrotational assumption before the integration process described above. However, such transformations do not completely satisfy the irrotationality for variable depth, that is, u→z-∇w≠0, which was mentioned in [23]. Subsequently, the present equation (9) is retained in the original forms of the nonlinear convective terms; that is, we preserve, completely, the irrotational condition. It is noted that Belibassakis and Athanassoulis [18] preserved this flow irrotationality by working with a potential formulation; they derived a nonlinear coupled-mode system of equations on the horizontal plane, with respect to the free surface elevation and the unknown modal amplitudes, in conjunction with Luke variational principle.The horizontal velocity vector can be expressed in terms of a vertical distribution function and a corresponding horizontal velocity vector at the still water level, in the form of(10)u⃑(x,y,z,t)=∑m=1Ncosh⁡km(h+z)cosh⁡kmhU→m(x,y,t),
where U→m(x,y,t) is the mth component of the horizontal velocity vector at z=0, and km(x,y) is the local wave number. The full dispersivity for a broad-banded wave field can be attained by taking only a few terms [23]. However, considering the fact that a narrow-banded sea state centred around a certain dominant wave frequency may be described with sufficient accuracy by the single-component model equations [24], that is, N=1, it becomes
(11)u⃑(x,y,z,t)=F(kh,z)U→(x,y,t),ffffffF(kh,z)=cosh⁡k(h+z)cosh⁡kh.
This single-component model has remained popular and useful in practical coastal engineering applications.The flow velocityw in the vertical direction z can be obtained upon depth integration of (1) from -h to z and the bottom boundary condition equation (6), rendering
(12)w(x,y,z,t)=-∇·[f(kh,z)kU→(x,y,t)],ffffffffffff(kh,z)=sinh⁡k(h+z)cosh⁡kh.Inserting (11) into (7) gives the following depth-integrated continuity equation:
(13)∂η∂t+∇·[(∫-hηFdz)U→]=0.
Substituting (11) and (12) into (8), and after laborious manipulations, we obtain the pressure expression as
(14)pρ=g(η-z)-p~a(∇·∂U→∂t)-p~b(∇h·∂U→∂t)-p~c[(U→·∇)(∇·U→)-(∇·U→)(∇·U→)]-p~d[(∇·U→)(∇h·U→)]+O[(∇h)2],
where
(15)p~a=∫zη(fk)dz,p~b=∫zη∂∂h(fk)dz,p~c=∫zηF(fk)dz,p~d=∫zη[F∂∂h(fk)-∂F∂h(fk)]dz.
The last term on the right-hand-side (RHS) of (14) is the higher order term of the bottom slope which may be omitted under the assumption of mild-slope.Multiplying (9) by the vertical distribution function F in (11) and then integrating it from z=-h to η, as well as after complex mathematical manipulations, produces a depth-integrated expression for the momentum equation in the form of
(16)α~1∂U→∂t+α~2g∇η+α~3(U→·∇)U→+α~4(∇·U→)U→+α~5(∇h·U→)U→=α~6∇(∇·∂U→∂t)+α~7(∇·∂U→∂t)∇h+α~8(∇·∂U→∂t)∇η+α~9(∇h·∇η)∂U→∂t+α~10∇[(U→·∇)(∇·U→)-(∇·U→)(∇·U→)]+α~11[(∇·U→)(∇h·∇)U→]+α~12[(U→·∇)(∇·U→)-(∇·U→)(∇·U→)]∇η+α~13[(∇·U→)(∇h·U→)∇η]+O[(∇h)2]
in which
(17)α~1=∫-hηF2dz,α~2=∫-hηFdz,α~3=∫-hηF3dz,α~4=-∫-hηFf2dz,α~5=∫-hη[F2∂F∂h-kFf∂∂h(fk)]dz,α~6=∫-hηFp~adz,α~7=2∫-hηFp~bdz,α~8=∫-hηF∂p~a∂ηdz,α~9=∫-hηF∂p~b∂ηdz,α~10=∫-hηFp~cdz,α~11=2∫-hηFp~ddz,α~12=∫-hηF∂p~c∂ηdz,α~13=∫-hηF∂p~d∂ηdz.The momentum equation, (16), retains a full nonlinear form under the mild-slope condition. As it is difficult to solve (16) directly by a numerical scheme; a perturbation technique is then proposed. The surface displacement and velocity field are expressed by
(18)η=εη1+ε2η2+⋯,U→=εU→1+ε2U→2+⋯,
where ε is the dimensionless nonlinearity parameter of wave amplitude, and the subscriptions i(i=1,2,…) in (18) denote the ith order approximation.Substituting (18) into (13) and (16) and taking Taylor series expansion on the still water level z=0, we obtain the continuity equation and the momentum equation, respectively, as follows:
(19)∂η∂t+∇·[(∫-h0Fdz+η)U→]=O(ε3),(20)(∫-h0F2dz)∂U→∂t+(∫-h0Fdz)g∇η+12[k2c4g2∫-h0F(F2-f2)dzmmmmcmm-2∫-h0Fdzk2c4g2]∇(U→·U→)=(∫-h0F-F2k2dz)∇(∇·∂U→∂t)+[∫-h0F∂∂h(1-Fk2)dz](∇·∂U→∂t)∇h+O(ε3).
These two equations are accurate to the second-order of ε. Substituting the expressions of F and f given in (11) and (12) into (19) and (20), we obtain the final expression of the continuity equation and momentum equation as follows:
(21)∂η∂t+∇·[(c2g+η)U→]=0,(22)n∂U→∂t+g∇η+12(1-3k2c4g2)∇(U→·U→)=1c2∇[c2(1-n)k2(∇·∂U→∂t)]
in which
(23)c2=gktanhkh,(24)n=12[1+2khsinh⁡(2kh)].The coefficient to the third term∇(U→·U→) on the left-hand-side of (22) consists of two parts. The first part, 1/2, is derived from the nonlinear convective terms of the horizontal velocity; and the second part, -3/2(k2c4/g2), is the nonlinear effect of the vertical velocity. The value of (k2c4/g2) tends to 1 in deep water, implying that the nonlinear variation of vertical velocity is in the same order as the horizontal velocity and it cannot be neglected, whereas it approximates reasonably to zero in the very shallow water. It is worth noting that the c and k can be computed for a prescribed wave period and a local depth based on the dispersion relationship shown in (23). The dispersion characteristics are essentially the same as in [23] because both models are approximated to the second-order of ε.Equations (21) and (22) together become the so-called second-order time-dependent mild-slope equation, because they can be used to represent the wave fields that are weakly nonlinear and fully dispersive and on the basis of the mild-slope assumption. Differing from [23], the present equations (21) and (22) can be directly used to calculate the wave transformation of wave field in plane because the vertical velocity has been embedded in the nonlinear convective terms into a simpler wave model. The consistency of the equation can be simplified into several well-established theories, such as linear mild-slope equation, nonlinear long wave equation, and Boussinesq equation. The process is demonstrated below.


## 2.3. Linear Mild-Slope Equation

Upon dropping the nonlinear term in (21) for the sake of linear (Airy) wave theory, this equation can be rewritten as
(25)∂η∂t+∇·Q⃑=0,
where Q→ is the linear flow rate of the waves defined by
(26)Q→=∫-h0u→dz=c2gU→.
On the other hand, (22) can be rewritten for the linear waves as
(27)n∂U→∂t+g∇η=1c2∇[c2(1-n)k2(∇·∂U→∂t)].
Substituting (26) into (27), then we obtain
(28)∂Q→∂t+c2n∇η=1n∇[(1-n)k2∇·∂Q→∂t-(1-n)c2k2∇c2·∂Q→∂t].
In the right-hand side of (28),
(29)1n∇[(1-n)k2(∇·∂Q→∂t)]=-1n∇[(1-n)k2∂2η∂t2]=1n∇[c2(1-n)η]=(1-n)nη∇c2+c2n∇η-c2n∇(nη),-1n∇[(1-n)c2k2∇c2·∂Q→∂t]=-(1-n)nc2k2∇c2(∇·∂Q→∂t)+O[(∇h)2]=-(1-n)nη∇c2+O[(∇h)2].
The term of O[(∇h)2] can be omitted by the assumption of mild-slope bottom. Upon inserting (29) into (28), the nonlinear equation developed in this study is reduced to
(30)∂Q⃑∂t+c2n∇(nη)=0.
Equations (25) and (30) are identical to the so-called time-dependent mild-slope equation in Copeland [8]. Eliminating flow rate Q→ from (25) and (30) renders the elliptic mild-slope equation in Berkhoff [1].


## 2.4. Nonlinear Long Wave Equation

For long wave theory in shallow water, wave celerity (c), group velocity (cg), and flow rate (Q→) are usually given by
(31)c=cg=gh,Q→=(h+η)U→,
where U→ is the horizontal velocity identical to mean flow velocity in shallow water condition. By dropping the nonlinear terms in the vertical direction in (22) and introducing (31), a system of nonlinear long wave equations can be written in terms of the mean flow velocity U→, given by
(32)∂U→∂t+g∇η+12∇(U→·U→)=0.
For an irrotational flow, the term ∇(U→·U→)/2 can be replaced by (U→·∇)U→, thus rendering the usual expression for the nonlinear long wave equation.


## 2.5. Boussinesq Equations

Boussinesq equations account for the effect of wave nonlinearity but weak dispersion, that the wave celerity and group velocity may be approximated by the first two terms in the Taylor series,(33)c2=gh(1-13k2h2),n=cgc=(1-13k2h2).
Substituting (33) into (21) and (22) and replacing k2U→ by -∇(∇·U→) for the irrotational flow in shallow water condition, the present nonlinear mild-slope equation gives
(34)∂η∂t+∇·[(h+η)U→]+h33∇2(∇·U→)=0,∂U→∂t+g∇η+12∇(U→·U→)=0.
Equations (34) are the standard form of traditional Boussinesq equations.


## 3. Numerical Methods

### 3.1. Finite Volume Method

A finite volume method is used to calculate the wave heights based on the governing equations given by (21) and (22). These two equations are first proceeded using a volume integration over a specific control volume ∀, as follows:
(35)∫∀(∂η∂t)d∀+∫∀[∇·(c2g+η)U→]d∀=0,∫∀(n∂U→∂t)d∀+∫∀[g∇η+12(1-3k2c4g2)∇(U→·U→)]d∀=∫∀1c2[∇(c2(1-n)k2∇·∂U→∂t)]d∀.
By the divergence theorem, (35) renders
(36)(∂η∂t)m=-1∀∫s[N→·(c2g+η)U→]dS∂∂t{nU→-1∀(1c2)m∫sN→·[c2(1-n)k2(∇·U→)I→]dS}m=-1∀g∫sN·(ηI→)dS-1∀[12(1-3k2c4g2)]m∫sN→·[(U→·U→)I→]dS
in which subscription m denotes the mean quantities referring to the centre of the volume ∀, S is the control surface, I→ is the unit vector, and N→ is the normal vector of S. By knowing the quantities at a starting time step (denotes as Gm), the integral form of (36) can be used to calculate the change in Gm within a time increment, and thus update the Gm values of the next time step.Regarding the computational mesh, it is envisaged that the grid size should be fine enough to resolve important characteristic details yet coarse enough to keep storage and computational time requirements reasonable. For instance, a fine mesh is used in the area where the wave transformation is significant, such as the region around the structure or near wave breaking. For the efficient and accurate computation for an irregular geometric boundary, the boundary-fitted mesh system is obtained in the present computation according to a method based on Poisson equations [25].


### 3.2. Adams-Moulton Predictor-Corrector Schemes

The fourth-order Adams-Moulton predictor-corrector technique is introduced to compute (36). The discrete equations of the predictor-step are given by
(37)ηm+1=ηm+Δt24(55Φm-59Φm-1+37Φm-2-9Φm-3),(38){[A]U→}m+1={[A]U→}m+Δt24(55Ψm-59Ψm-1+37Ψm-2-9Ψm-3).
And the equations of the corrector-step are expressed by
(39)ηm+1=ηm+Δt24(9Φm+1+19Φm-5Φm-1+Φm-2),(40){[A]U→}m+1={[A]U→}m+Δt24(9Ψm+1+19Ψm-5Ψm-1+Ψm-2),
where
(41)Φ=-1∀∑K=12{S→K+·[(c2g+η)U→]K+-S→K-·[(c2g+η)U→]K-},Ψ=-g∀∑K=12[S→K+·(η)K+I→-S→K-·(η)K-I→]-1∀[12(1-3k2c4g2)]×∑K=12[S→K+·(U→·U→)K+I→-S→K-·(U→·U→)K-I→],[A]U→=nU→-1∀(1c2)∑K=12[(c2(1-n)k2)K+·(∇·U→)K+I→mmmmmmmmmmm-(c2(1-n)k2)K-·(∇·U→)K-I→(c2(1-n)k2)K+],∇·U→=1∀∑K=12[S→K+·(U→)K+-S→K-·(U→)K-]
in which superscript m represents the computed values of η, U→, Φ, and Ψ at the mth time iteration step, and superscript m+1 denotes the predicted values at the m+1th time step. The subscript K in (41) represents the two directions in x-y plane, in which superscript “+” denotes the quantities in the upstream of control-volume and “−” for those in the downstream. The section vectors denoted as S→K+ and S→K- represent opposing faces of the middle volume. The definition of variables in the finite volume method is illustrated in Figure 1. In the discrete continuity equations, that is, (37) and (39), the new values can be computed explicitly. But the implicit process is required to obtain the new values for the momentum equations, (38) and (40), because the dispersive term involves the time and spatial derivatives simultaneously. That is, the surface elevation η can be solved explicitly, but the velocity U→ is solved by an implicit numerical scheme. The successive over-relaxation iteration method (SOR) is adopted for the implicit scheme.Figure 1
Definition sketch of the variables in finite volume method.As shown in Figure2, the quantities of variables η and U→ are set in the centre of the control volume. Then the flux terms, Φ and Ψ, in (37) to (40), can be expressed as
(42)Φi,j=-12(∀i,j+∀i+1,j+∀i,j+1+∀i+1,j+1)·A1,Ψi,j=-g2(∀i,j+∀i+1,j+∀i,j+1+∀i+1,j+1)·A2-14(∀i,j+∀i+1,j+∀i,j+1+∀i+1,j+1)×(1-3k2c4g2)i,j·A3,
where
(43)A1={[(S→1)i+1,j+1·(R→i+1,j+1+R→i+1,j)mmc+(S→1)i+1,j·(R→i+1,j-1+R→i+1,j)]mm-[(S→1)i-1,j+1·(R→i-1,j+1+R→i-1,j)mmmc+(S→1)i-1,j·(R→i-1,j-1+R→i-1,j)]mm+[(S→2)i+1,j+1·(R→i+1,j+1+R→i,j+1)mmmc+(S→2)i,j+1·(R→i-1,j+1+R→i,j+1)]mm-[(S→2)i+1,j-1·(R→i+1,j-1+R→i,j-1)mmmc+(S→2)i,j-1·(R→i-1,j-1+R→i,j-1)]},A2={[(S→1)i+1,j+1·(ηi+1,j+1+ηi+1,j)I→mmc+(S→1)i+1,j·(ηi+1,j-1+ηi+1,j)I→]mm-[(S→1)i-1,j+1·(ηi-1,j+1+ηi-1,j)I→mmcm+(S→1)i-1,j·(ηi-1,j-1+ηi-1,j)I→]mm+[(S→2)i+1,j+1·(ηi+1,j+1+ηi,j+1)I→mcmm+(S→2)i,j+1·(ηi-1,j+1+ηi,j+1)I→]mm-[(S→2)i+1,j-1·(ηi+1,j-1+ηi,j-1)I→mmmc+(S→2)i,j-1·(ηi-1,j-1+ηi,j-1)I→]},A3={[(S→1)i+1,j+1·(U¨i+1,j+1+U¨i+1,j)I→mmc+(S→1)i+1,j·(U¨i+1,j-1+U¨i+1,j)I→]mm-[(S→1)i-1,j+1·(U¨i-1,j+1+U¨i-1,j)I→mmmc+(S→1)i-1,j·(U¨i-1,j-1+U¨i-1,j)I→]mm+[(S→2)i+1,j+1·(U¨i+1,j+1+U¨i,j+1)I→mmmc+(S→2)i,j+1·(U¨i-1,j+1+U¨i,j+1)I→]mm-[(S→2)i+1,j-1·(U¨i+1,j-1+U¨i,j-1)I→mmmc+(S→2)i,j-1·(U¨i-1,j-1+U¨i,j-1)I→]},R→i,j=(ci,j2g+ηi,j)U→i,j,U¨i,j=U→i,j·U→i,j=Ui,j2+Vi,j2.Figure 2
Mesh arrangement and corresponding variables around the centre of a control volume.


### 3.3. Initial and Converge Conditions

The condition of zero flow or cold start has been set as the initial condition for most numerical simulations. In this case, an unsteady fluctuation will commonly be produced in the initial computing stage, but it will become negligible when the computation time is long enough, and finally a steady state can be obtained. In the explicit iteration for solving the surface elevation, the Courant stability condition has to be satisfied at each time step, that is,(44)Δt≤fc·min⁡(∀|c→·S→|)
in which fc is the safe coefficient (commonly taken as 0.6). In the computations, several iteration cycles are usually required to achieve a steady result. The computed results are regarded as reaching a steady state when the iterative errors of the computed wave height are less than 2%. As for the velocity in the horizontal plane, U→, each time interval in the SOR iteration method needs to satisfy the following converge condition:
(45)ΔeU⃑=max⁡(|U→m+1-U→mU→m+1|)≤10-4
in which ΔeU⃑ denotes the absolute maximum error.


### 3.4. Boundary Specifications

Three types of boundary conditions are specified for computing the wave field. These are applied using phantom point exterior from the computational mesh. The first one is an offshore boundary which is the incident wave conditions. In general, there exist not only incident waves but also outgoing wave components across the offshore boundary due to the reflection from structures or the shore. Referring to [9, 26], such outgoing waves should be transmitted freely through the boundary. Thus the offshore boundary is set as
(46)(η)B′t=(ηinc)B′t+(ηout)Bt-τ(U→)B′t=(U→inc)B′t+(U→out)Bt-τ,
where
(47)τ=Δr·cos⁡⁡(αout)c,(ηout)Bt=(η)Bt-(ηinc)Bt,(U→out)Bt=(U→)Bt-(U→inc)Bt,
where the subscriptions inc and out represent the quantities of incoming and outgoing waves, respectively, Δr is the length of BB′, B′ is the phantom point of B, and αout is the angle between outgoing wave and BB′, as shown in Figure 3.Figure 3
Sketch of the offshore open boundary.The second type is for the side or lateral boundary of the computation domain, which is generally a virtual open boundary. It is treated as having no wave reflection at this boundary, hence, allowing waves to pass freely through them. This can be represented by(48)(η)B′t=(η)Bt-τ,(U→)B′t=(U→)Bt-τ.
And the third is the internal boundary. If an impermeable structure exists, where full reflection from the structure surface is assumed for the internal boundary associated with wave reflection, it can be set as
(49)(η)B′t=(η)Bt,(U→St)B′t=(U→St)Bt,(U→Sn)B′t=-(U→Sn)Bt,
where subscriptionsSt andSn represent the tangent and normal components at the boundary, respectively.


## 3.1. Finite Volume Method

A finite volume method is used to calculate the wave heights based on the governing equations given by (21) and (22). These two equations are first proceeded using a volume integration over a specific control volume ∀, as follows:
(35)∫∀(∂η∂t)d∀+∫∀[∇·(c2g+η)U→]d∀=0,∫∀(n∂U→∂t)d∀+∫∀[g∇η+12(1-3k2c4g2)∇(U→·U→)]d∀=∫∀1c2[∇(c2(1-n)k2∇·∂U→∂t)]d∀.
By the divergence theorem, (35) renders
(36)(∂η∂t)m=-1∀∫s[N→·(c2g+η)U→]dS∂∂t{nU→-1∀(1c2)m∫sN→·[c2(1-n)k2(∇·U→)I→]dS}m=-1∀g∫sN·(ηI→)dS-1∀[12(1-3k2c4g2)]m∫sN→·[(U→·U→)I→]dS
in which subscription m denotes the mean quantities referring to the centre of the volume ∀, S is the control surface, I→ is the unit vector, and N→ is the normal vector of S. By knowing the quantities at a starting time step (denotes as Gm), the integral form of (36) can be used to calculate the change in Gm within a time increment, and thus update the Gm values of the next time step.Regarding the computational mesh, it is envisaged that the grid size should be fine enough to resolve important characteristic details yet coarse enough to keep storage and computational time requirements reasonable. For instance, a fine mesh is used in the area where the wave transformation is significant, such as the region around the structure or near wave breaking. For the efficient and accurate computation for an irregular geometric boundary, the boundary-fitted mesh system is obtained in the present computation according to a method based on Poisson equations [25].


## 3.2. Adams-Moulton Predictor-Corrector Schemes

The fourth-order Adams-Moulton predictor-corrector technique is introduced to compute (36). The discrete equations of the predictor-step are given by
(37)ηm+1=ηm+Δt24(55Φm-59Φm-1+37Φm-2-9Φm-3),(38){[A]U→}m+1={[A]U→}m+Δt24(55Ψm-59Ψm-1+37Ψm-2-9Ψm-3).
And the equations of the corrector-step are expressed by
(39)ηm+1=ηm+Δt24(9Φm+1+19Φm-5Φm-1+Φm-2),(40){[A]U→}m+1={[A]U→}m+Δt24(9Ψm+1+19Ψm-5Ψm-1+Ψm-2),
where
(41)Φ=-1∀∑K=12{S→K+·[(c2g+η)U→]K+-S→K-·[(c2g+η)U→]K-},Ψ=-g∀∑K=12[S→K+·(η)K+I→-S→K-·(η)K-I→]-1∀[12(1-3k2c4g2)]×∑K=12[S→K+·(U→·U→)K+I→-S→K-·(U→·U→)K-I→],[A]U→=nU→-1∀(1c2)∑K=12[(c2(1-n)k2)K+·(∇·U→)K+I→mmmmmmmmmmm-(c2(1-n)k2)K-·(∇·U→)K-I→(c2(1-n)k2)K+],∇·U→=1∀∑K=12[S→K+·(U→)K+-S→K-·(U→)K-]
in which superscript m represents the computed values of η, U→, Φ, and Ψ at the mth time iteration step, and superscript m+1 denotes the predicted values at the m+1th time step. The subscript K in (41) represents the two directions in x-y plane, in which superscript “+” denotes the quantities in the upstream of control-volume and “−” for those in the downstream. The section vectors denoted as S→K+ and S→K- represent opposing faces of the middle volume. The definition of variables in the finite volume method is illustrated in Figure 1. In the discrete continuity equations, that is, (37) and (39), the new values can be computed explicitly. But the implicit process is required to obtain the new values for the momentum equations, (38) and (40), because the dispersive term involves the time and spatial derivatives simultaneously. That is, the surface elevation η can be solved explicitly, but the velocity U→ is solved by an implicit numerical scheme. The successive over-relaxation iteration method (SOR) is adopted for the implicit scheme.Figure 1
Definition sketch of the variables in finite volume method.As shown in Figure2, the quantities of variables η and U→ are set in the centre of the control volume. Then the flux terms, Φ and Ψ, in (37) to (40), can be expressed as
(42)Φi,j=-12(∀i,j+∀i+1,j+∀i,j+1+∀i+1,j+1)·A1,Ψi,j=-g2(∀i,j+∀i+1,j+∀i,j+1+∀i+1,j+1)·A2-14(∀i,j+∀i+1,j+∀i,j+1+∀i+1,j+1)×(1-3k2c4g2)i,j·A3,
where
(43)A1={[(S→1)i+1,j+1·(R→i+1,j+1+R→i+1,j)mmc+(S→1)i+1,j·(R→i+1,j-1+R→i+1,j)]mm-[(S→1)i-1,j+1·(R→i-1,j+1+R→i-1,j)mmmc+(S→1)i-1,j·(R→i-1,j-1+R→i-1,j)]mm+[(S→2)i+1,j+1·(R→i+1,j+1+R→i,j+1)mmmc+(S→2)i,j+1·(R→i-1,j+1+R→i,j+1)]mm-[(S→2)i+1,j-1·(R→i+1,j-1+R→i,j-1)mmmc+(S→2)i,j-1·(R→i-1,j-1+R→i,j-1)]},A2={[(S→1)i+1,j+1·(ηi+1,j+1+ηi+1,j)I→mmc+(S→1)i+1,j·(ηi+1,j-1+ηi+1,j)I→]mm-[(S→1)i-1,j+1·(ηi-1,j+1+ηi-1,j)I→mmcm+(S→1)i-1,j·(ηi-1,j-1+ηi-1,j)I→]mm+[(S→2)i+1,j+1·(ηi+1,j+1+ηi,j+1)I→mcmm+(S→2)i,j+1·(ηi-1,j+1+ηi,j+1)I→]mm-[(S→2)i+1,j-1·(ηi+1,j-1+ηi,j-1)I→mmmc+(S→2)i,j-1·(ηi-1,j-1+ηi,j-1)I→]},A3={[(S→1)i+1,j+1·(U¨i+1,j+1+U¨i+1,j)I→mmc+(S→1)i+1,j·(U¨i+1,j-1+U¨i+1,j)I→]mm-[(S→1)i-1,j+1·(U¨i-1,j+1+U¨i-1,j)I→mmmc+(S→1)i-1,j·(U¨i-1,j-1+U¨i-1,j)I→]mm+[(S→2)i+1,j+1·(U¨i+1,j+1+U¨i,j+1)I→mmmc+(S→2)i,j+1·(U¨i-1,j+1+U¨i,j+1)I→]mm-[(S→2)i+1,j-1·(U¨i+1,j-1+U¨i,j-1)I→mmmc+(S→2)i,j-1·(U¨i-1,j-1+U¨i,j-1)I→]},R→i,j=(ci,j2g+ηi,j)U→i,j,U¨i,j=U→i,j·U→i,j=Ui,j2+Vi,j2.Figure 2
Mesh arrangement and corresponding variables around the centre of a control volume.


## 3.3. Initial and Converge Conditions

The condition of zero flow or cold start has been set as the initial condition for most numerical simulations. In this case, an unsteady fluctuation will commonly be produced in the initial computing stage, but it will become negligible when the computation time is long enough, and finally a steady state can be obtained. In the explicit iteration for solving the surface elevation, the Courant stability condition has to be satisfied at each time step, that is,(44)Δt≤fc·min⁡(∀|c→·S→|)
in which fc is the safe coefficient (commonly taken as 0.6). In the computations, several iteration cycles are usually required to achieve a steady result. The computed results are regarded as reaching a steady state when the iterative errors of the computed wave height are less than 2%. As for the velocity in the horizontal plane, U→, each time interval in the SOR iteration method needs to satisfy the following converge condition:
(45)ΔeU⃑=max⁡(|U→m+1-U→mU→m+1|)≤10-4
in which ΔeU⃑ denotes the absolute maximum error.


## 3.4. Boundary Specifications

Three types of boundary conditions are specified for computing the wave field. These are applied using phantom point exterior from the computational mesh. The first one is an offshore boundary which is the incident wave conditions. In general, there exist not only incident waves but also outgoing wave components across the offshore boundary due to the reflection from structures or the shore. Referring to [9, 26], such outgoing waves should be transmitted freely through the boundary. Thus the offshore boundary is set as
(46)(η)B′t=(ηinc)B′t+(ηout)Bt-τ(U→)B′t=(U→inc)B′t+(U→out)Bt-τ,
where
(47)τ=Δr·cos⁡⁡(αout)c,(ηout)Bt=(η)Bt-(ηinc)Bt,(U→out)Bt=(U→)Bt-(U→inc)Bt,
where the subscriptions inc and out represent the quantities of incoming and outgoing waves, respectively, Δr is the length of BB′, B′ is the phantom point of B, and αout is the angle between outgoing wave and BB′, as shown in Figure 3.Figure 3
Sketch of the offshore open boundary.The second type is for the side or lateral boundary of the computation domain, which is generally a virtual open boundary. It is treated as having no wave reflection at this boundary, hence, allowing waves to pass freely through them. This can be represented by(48)(η)B′t=(η)Bt-τ,(U→)B′t=(U→)Bt-τ.
And the third is the internal boundary. If an impermeable structure exists, where full reflection from the structure surface is assumed for the internal boundary associated with wave reflection, it can be set as
(49)(η)B′t=(η)Bt,(U→St)B′t=(U→St)Bt,(U→Sn)B′t=-(U→Sn)Bt,
where subscriptionsSt andSn represent the tangent and normal components at the boundary, respectively.


## 4. Numerical Results

Validation of the present model is made by comparing the numerical results with analytical solutions and/or experimental data for Stokes wave and cnoidal wave propagation on uniform depth, one-dimensional wave shoaling, and two-dimensional wave transformation across an elliptic shoal. Wave transformation involving energy dissipation due to wave breaking is finally simulated for wave propagation on a uniform slope beach and a bar-type beach and wave around a jetty.


### 4.1. Stokes Waves

The first case of validation is the numerical simulation of Stokes waves to the second-order. The computed results for two differentkH values in intermediate depth (h/L=0.2) are depicted in Figure 4, in which H is the wave height, L is the wavelength and ao=H/2. The temporal and spatial resolutions are Δt=T/60 and Δx=L/30, respectively, where T is the wave period. It can be seen that the numerical results are in good agreement with the theoretical solution of the Stokes waves to a second-order. The minor discrepancy in phase for large x/L is possibly caused by the truncation error in the discrete numerical process which affects the solution of the phase velocity. Nevertheless, the verification indicates the present model is capable of producing identical results to the Stokes waves in water of finite depth.Figure 4
Simulated Stokes waves in intermediate water depth compared with analytical solution (circles: present solution, line: 2nd-order Stokes wave theory).The vertical distribution of the horizontal velocity can be calculated from (10) once the solution of U→ is obtained. The maximum nondimensional horizontal velocity calculated from (10) is plotted against dimensionless depth z/h in Figure 5. Good agreement is found upon comparing the results with the experimental data of Kim et al. [27].Figure 5
Vertical distribution of the computed maximum horizontal velocity compared with experimental data.


### 4.2. Cnoidal Waves

It is well known that cnoidal wave theory is preferred than Stokes wave theory in shallow water. Figure6 demonstrates the simulated wave profiles of cnoidal waves in relative depth of h/L=0.05 for two different H/h values. The temporal and spatial resolutions are Δt=T/40 and Δx=L/20 in this computation. The present numerical solution is also found in good agreement with the theoretical solution of cnoidal waves.Figure 6
Simulated wave profiles for cnoidal waves compared favourably with analytical solution (circles: present solution, line: Cnoidal wave theory).In very shallow water region, cnoidal waves approach the limit of a solitary wave. Since it is not possible to define a finite wavelength for a solitary wave, as it is infinite in theory, the simulation of the solitary wave is often made by taking a very small value ofh/L. Here h/L is taken as 1/50. Figure 7 illustrates the simulated wave profiles in relative depth of h/L  = 0.02 for H/h  = 0.20 and 0.30, respectively.Figure 7
Simulated wave profile for solitary waves in very shallow water compared with analytical solution (circles: present solution, line: solitary wave theory).From the comparisons above for simulations of nonlinear waves on uniform depth, it reveals that the present model is adequate for the entire range of water depth from deep to very shallow. The following simulations will cover nonlinear wave shoaling and two-dimensional wave field involving combined wave refraction-diffraction, even wave breaking.


### 4.3. Wave Shoaling on Uniform Slope

Shuto [28] proposed an equation for nonlinear wave shoaling prior to breaking along a uniform slope based on a K-dV formulation. His results were found to agree reasonably well with the experimental data for a wide range of wave steepness [29]. The validity of the present nonlinear model is further tested against the results of Shuto [28] by comparing the shoaling coefficients H/Ho versus relative water depth h/Lo for wave steepness Ho/Lo ranging from 0.001 to 0.08, on uniform beach slope of 1/30, in which Ho and Lo are the wave height and wavelength in deep water. As seen in Figure 8, the results computed from the present model (in solid circles) are in good agreement with the nonlinear shoaling theory of Shuto [28], despite small fluctuations about the simulated curves within some range of h/Lo. Tsai et al. [10] incorporated a nonlinear shoaling factor into the linear time-dependent mild-slope equations; their results (circles in Figure 8) are also in good agreement with Shuto [28]. At a particular h/Lo, the shoaling coefficient increases as wave steepness increases. However, the solution computed by the linear wave shoaling is independent of the wave steepness Ho/Lo which apparently underestimates any shoaling wave in shallow water compared with the present nonlinear shoaling model, but has the same value between the nonlinear and linear wave shoaling as h/Lo greater than 0.15.Figure 8
Comparison on the wave height from the present nonlinear wave shoaling solution with other solutions.


### 4.4. Wave Transformation across an Elliptic Shoal

The hydraulic model test [30] for wave propagation over a 3D elliptical shoal on a plane beach with uniform slope 1 : 50 has usually been employed for wave model verification. The wave height distribution was measured along eight cross sections indicated in Figure 9. An incident wave with a period of 1.0 sec and wave height of 2.32 cm was used in the simulation. The finer rectangular grids are set in the shoal region to better resolve the wave transformation in the computation. The simulation results for wave height distribution behind the shoal (seen in Figure 10) exhibit significant contribution of the combined effect of refraction-diffraction.Figure 9
Geometry showing wave propagation over an elliptic shoal.Figure 10
Simulation results on wave height distribution behind an elliptic shoal.The simulation results for wave height distribution along the eight transects indicated in Figure11 are compared with the experimental data and the numerical model of Li and Anastasiou [3], who derived an improved solution using an elliptic form of the linear mild-slope equation. The comparison in Figure 11 demonstrates that the present nonlinear model can yield the wave height oscillations spatially distributed due to combined refraction-diffraction effect behind an elliptical shoal, particularly on transects B to E.Simulation results of wave height distribution along eight transects shown in Figure9 and comparison with experimental data and other numerical models.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)


### 4.5. 1D and 2D Wave Transformation across Surf Zone

Moreover, the momentum equation used in the present model can be applied to simulate wave transformation across the surf zone, upon incorporating a mixing term,MD. In this way, (22) becomes
(50)∂U→∂t+gn∇η+12n(1-3k2c4g2)∇(U→·U→)=1nc2∇[c2(1-n)k2(∇·∂U→∂t)]+MD.
Following Kabiling and Sato [31], the mixing term due to energy dissipation can be expressed by
(51)MD=νe∇2U→,
where νe is the eddy viscosity related to an energy dissipation factor, fD, due to breaking,
(52)νe=ghω2fD;fD=αDtanβgh(Qm-QrQs-Qr).
In (52), fD is the energy dissipation factor, ω(=2π/T) is the wave frequency, tanβ is an average bottom slope from shoreline to wave breaking point, αD is a modification coefficient (usually taken as 2.5 as suggested in [32]), Qm is the amplitude of flow rate, and Qr and Qs are determined from the following equations:
(53)Qs=0.4(0.57+5.3tanβ)gh3,Qr=0.135gh3.
In (52), the dissipation factor fD is set equal to zero outside the surf zone and in any region where Qm<Qr.On the basis of shallow water assumption for wave breaking, we can simplify∇2U→=-k2U→ and render MD=-fDU→, thus (50) becomes
(54)∂U→∂t+gn∇η+12n(1-3k2c4g2)∇(U→·U→)+fDU→=1nc2∇[c2(1-n)k2(∇·∂U→∂t)].
This modified momentum equation can now be used for the wave transformation after wave breaking. To estimate the location of wave breaking, a breaking criterion proposed by Goda [33] is used as follows:
(55)HbLo=0.17{1-exp⁡[-1.5πhbLo(1+15tan4/3β)]}
in which Hb and hb are the wave height and water depth of the breaking wave, respectively, and Lo is the deep water wavelength. The computed wave height at any location based on (21) and (22) is always compared with the value calculated from the right-hand side of (55). If the computed wave height reaches the breaking criterion, the momentum equation (54), instead of (22), is then used for calculating the wave height transformation after wave breaking.Numerical verification on one-dimensional wave transformation across the surf zone is now considered for two different bottom configurations, a uniform slope in 1/20 and a composite profile with a bar in 1/20 slope, as shown in Figures12 and 13, respectively. The numerical results of the present nonlinear model are compared with the linear solution of wave transformation [9] and the experimental data of Nagayama [34]. The results of Tsai et al. [10] are also included in these figures. The local wave heights calculated using the present nonlinear model are in better agreement with the experimental data than the linear model.Figure 12
Comparisons among numerical solutions and laboratory experiment for wave transformation on uniform slope.Figure 13
Comparisons among numerical solutions and laboratory experiment for wave transformation on a composite profile of bar type.In Figure12, the computed wave heights increases as it shoals, which ultimately approaches a maximum at breaking, then decreases rapidly until reaching the shoreline. For the resultant wave height transformation on a composite slope, as in Figure 13, wave height reaches a maximum prior to the bar crest where first breaking occurs on the leading uniform slope of 1/20. Wave height recovers soon after it passes the bar trough connected to the trailing uniform slope on which the second breaking appears. Finally, the surviving wave diminishes at the shoreline. For the region before the breaking on a uniform slope (Figure 12) and prior to the first breaking on a composite slope (Figure 13), the present nonlinear solution and linear model agree with each other, otherwise the linear model is found to deviate from the experimental results. This confirms that the present nonlinear model is adequate for wave decay, recovery, and secondary breaking on a composite beach profile.The example of two-dimensional wave transformation was conducted for wave field in the vicinity of a jetty, for which experimental and numerical results by Watanabe and Maruyama [9] were available. The jetty in the experimental tank was placed perpendicular to the shoreline on a beach with a slope of 1/50. The model jetty was extended to an offshore distance of 4 m, where the water depth was 8 cm. The waves were incident with angle 60° to the contour lines, and the wave height and period in deep water were 2.0 cm and 1.2 s, respectively.The computed wave height distribution around the jetty is shown in Figure14. Aclapotis gaufré wave system is visible in front of the jetty by the reflection from it. On the other hand, the wave diffraction was observed in the lee of the jetty. The comparisons of the wave height distributions at cross-shore sections y=3.0 m, 4.8 m, and 5.2 m are shown in Figure 15, which indicates the present nonlinear solutions are in better agreement with the experimental results.Figure 14
The computed wave field in the vicinity of a jetty.Figure 15
Comparisons among numerical solutions and experimental results for distributions of wave height at sectionsy=3.0 m, 4.8 m, and 5.2 m shown in Figure 14.


## 4.1. Stokes Waves

The first case of validation is the numerical simulation of Stokes waves to the second-order. The computed results for two differentkH values in intermediate depth (h/L=0.2) are depicted in Figure 4, in which H is the wave height, L is the wavelength and ao=H/2. The temporal and spatial resolutions are Δt=T/60 and Δx=L/30, respectively, where T is the wave period. It can be seen that the numerical results are in good agreement with the theoretical solution of the Stokes waves to a second-order. The minor discrepancy in phase for large x/L is possibly caused by the truncation error in the discrete numerical process which affects the solution of the phase velocity. Nevertheless, the verification indicates the present model is capable of producing identical results to the Stokes waves in water of finite depth.Figure 4
Simulated Stokes waves in intermediate water depth compared with analytical solution (circles: present solution, line: 2nd-order Stokes wave theory).The vertical distribution of the horizontal velocity can be calculated from (10) once the solution of U→ is obtained. The maximum nondimensional horizontal velocity calculated from (10) is plotted against dimensionless depth z/h in Figure 5. Good agreement is found upon comparing the results with the experimental data of Kim et al. [27].Figure 5
Vertical distribution of the computed maximum horizontal velocity compared with experimental data.


## 4.2. Cnoidal Waves

It is well known that cnoidal wave theory is preferred than Stokes wave theory in shallow water. Figure6 demonstrates the simulated wave profiles of cnoidal waves in relative depth of h/L=0.05 for two different H/h values. The temporal and spatial resolutions are Δt=T/40 and Δx=L/20 in this computation. The present numerical solution is also found in good agreement with the theoretical solution of cnoidal waves.Figure 6
Simulated wave profiles for cnoidal waves compared favourably with analytical solution (circles: present solution, line: Cnoidal wave theory).In very shallow water region, cnoidal waves approach the limit of a solitary wave. Since it is not possible to define a finite wavelength for a solitary wave, as it is infinite in theory, the simulation of the solitary wave is often made by taking a very small value ofh/L. Here h/L is taken as 1/50. Figure 7 illustrates the simulated wave profiles in relative depth of h/L  = 0.02 for H/h  = 0.20 and 0.30, respectively.Figure 7
Simulated wave profile for solitary waves in very shallow water compared with analytical solution (circles: present solution, line: solitary wave theory).From the comparisons above for simulations of nonlinear waves on uniform depth, it reveals that the present model is adequate for the entire range of water depth from deep to very shallow. The following simulations will cover nonlinear wave shoaling and two-dimensional wave field involving combined wave refraction-diffraction, even wave breaking.


## 4.3. Wave Shoaling on Uniform Slope

Shuto [28] proposed an equation for nonlinear wave shoaling prior to breaking along a uniform slope based on a K-dV formulation. His results were found to agree reasonably well with the experimental data for a wide range of wave steepness [29]. The validity of the present nonlinear model is further tested against the results of Shuto [28] by comparing the shoaling coefficients H/Ho versus relative water depth h/Lo for wave steepness Ho/Lo ranging from 0.001 to 0.08, on uniform beach slope of 1/30, in which Ho and Lo are the wave height and wavelength in deep water. As seen in Figure 8, the results computed from the present model (in solid circles) are in good agreement with the nonlinear shoaling theory of Shuto [28], despite small fluctuations about the simulated curves within some range of h/Lo. Tsai et al. [10] incorporated a nonlinear shoaling factor into the linear time-dependent mild-slope equations; their results (circles in Figure 8) are also in good agreement with Shuto [28]. At a particular h/Lo, the shoaling coefficient increases as wave steepness increases. However, the solution computed by the linear wave shoaling is independent of the wave steepness Ho/Lo which apparently underestimates any shoaling wave in shallow water compared with the present nonlinear shoaling model, but has the same value between the nonlinear and linear wave shoaling as h/Lo greater than 0.15.Figure 8
Comparison on the wave height from the present nonlinear wave shoaling solution with other solutions.


## 4.4. Wave Transformation across an Elliptic Shoal

The hydraulic model test [30] for wave propagation over a 3D elliptical shoal on a plane beach with uniform slope 1 : 50 has usually been employed for wave model verification. The wave height distribution was measured along eight cross sections indicated in Figure 9. An incident wave with a period of 1.0 sec and wave height of 2.32 cm was used in the simulation. The finer rectangular grids are set in the shoal region to better resolve the wave transformation in the computation. The simulation results for wave height distribution behind the shoal (seen in Figure 10) exhibit significant contribution of the combined effect of refraction-diffraction.Figure 9
Geometry showing wave propagation over an elliptic shoal.Figure 10
Simulation results on wave height distribution behind an elliptic shoal.The simulation results for wave height distribution along the eight transects indicated in Figure11 are compared with the experimental data and the numerical model of Li and Anastasiou [3], who derived an improved solution using an elliptic form of the linear mild-slope equation. The comparison in Figure 11 demonstrates that the present nonlinear model can yield the wave height oscillations spatially distributed due to combined refraction-diffraction effect behind an elliptical shoal, particularly on transects B to E.Simulation results of wave height distribution along eight transects shown in Figure9 and comparison with experimental data and other numerical models.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)


## 4.5. 1D and 2D Wave Transformation across Surf Zone

Moreover, the momentum equation used in the present model can be applied to simulate wave transformation across the surf zone, upon incorporating a mixing term,MD. In this way, (22) becomes
(50)∂U→∂t+gn∇η+12n(1-3k2c4g2)∇(U→·U→)=1nc2∇[c2(1-n)k2(∇·∂U→∂t)]+MD.
Following Kabiling and Sato [31], the mixing term due to energy dissipation can be expressed by
(51)MD=νe∇2U→,
where νe is the eddy viscosity related to an energy dissipation factor, fD, due to breaking,
(52)νe=ghω2fD;fD=αDtanβgh(Qm-QrQs-Qr).
In (52), fD is the energy dissipation factor, ω(=2π/T) is the wave frequency, tanβ is an average bottom slope from shoreline to wave breaking point, αD is a modification coefficient (usually taken as 2.5 as suggested in [32]), Qm is the amplitude of flow rate, and Qr and Qs are determined from the following equations:
(53)Qs=0.4(0.57+5.3tanβ)gh3,Qr=0.135gh3.
In (52), the dissipation factor fD is set equal to zero outside the surf zone and in any region where Qm<Qr.On the basis of shallow water assumption for wave breaking, we can simplify∇2U→=-k2U→ and render MD=-fDU→, thus (50) becomes
(54)∂U→∂t+gn∇η+12n(1-3k2c4g2)∇(U→·U→)+fDU→=1nc2∇[c2(1-n)k2(∇·∂U→∂t)].
This modified momentum equation can now be used for the wave transformation after wave breaking. To estimate the location of wave breaking, a breaking criterion proposed by Goda [33] is used as follows:
(55)HbLo=0.17{1-exp⁡[-1.5πhbLo(1+15tan4/3β)]}
in which Hb and hb are the wave height and water depth of the breaking wave, respectively, and Lo is the deep water wavelength. The computed wave height at any location based on (21) and (22) is always compared with the value calculated from the right-hand side of (55). If the computed wave height reaches the breaking criterion, the momentum equation (54), instead of (22), is then used for calculating the wave height transformation after wave breaking.Numerical verification on one-dimensional wave transformation across the surf zone is now considered for two different bottom configurations, a uniform slope in 1/20 and a composite profile with a bar in 1/20 slope, as shown in Figures12 and 13, respectively. The numerical results of the present nonlinear model are compared with the linear solution of wave transformation [9] and the experimental data of Nagayama [34]. The results of Tsai et al. [10] are also included in these figures. The local wave heights calculated using the present nonlinear model are in better agreement with the experimental data than the linear model.Figure 12
Comparisons among numerical solutions and laboratory experiment for wave transformation on uniform slope.Figure 13
Comparisons among numerical solutions and laboratory experiment for wave transformation on a composite profile of bar type.In Figure12, the computed wave heights increases as it shoals, which ultimately approaches a maximum at breaking, then decreases rapidly until reaching the shoreline. For the resultant wave height transformation on a composite slope, as in Figure 13, wave height reaches a maximum prior to the bar crest where first breaking occurs on the leading uniform slope of 1/20. Wave height recovers soon after it passes the bar trough connected to the trailing uniform slope on which the second breaking appears. Finally, the surviving wave diminishes at the shoreline. For the region before the breaking on a uniform slope (Figure 12) and prior to the first breaking on a composite slope (Figure 13), the present nonlinear solution and linear model agree with each other, otherwise the linear model is found to deviate from the experimental results. This confirms that the present nonlinear model is adequate for wave decay, recovery, and secondary breaking on a composite beach profile.The example of two-dimensional wave transformation was conducted for wave field in the vicinity of a jetty, for which experimental and numerical results by Watanabe and Maruyama [9] were available. The jetty in the experimental tank was placed perpendicular to the shoreline on a beach with a slope of 1/50. The model jetty was extended to an offshore distance of 4 m, where the water depth was 8 cm. The waves were incident with angle 60° to the contour lines, and the wave height and period in deep water were 2.0 cm and 1.2 s, respectively.The computed wave height distribution around the jetty is shown in Figure14. Aclapotis gaufré wave system is visible in front of the jetty by the reflection from it. On the other hand, the wave diffraction was observed in the lee of the jetty. The comparisons of the wave height distributions at cross-shore sections y=3.0 m, 4.8 m, and 5.2 m are shown in Figure 15, which indicates the present nonlinear solutions are in better agreement with the experimental results.Figure 14
The computed wave field in the vicinity of a jetty.Figure 15
Comparisons among numerical solutions and experimental results for distributions of wave height at sectionsy=3.0 m, 4.8 m, and 5.2 m shown in Figure 14.


## 5. Conclusions

In this paper, the development of a second-order time-dependent mild-slope equation is reported, based on the fundamental principles of conservation of mass and momentum, together with a depth-integration procedure. After some appropriate mathematical simplifications, the present nonlinear equation can yield the usual form of linear mild-slope equation, nonlinear long wave equation, and Boussinesq equation, respectively. We have also demonstrated that the present nonlinear model is applicable to the entire range of water depth from deep to very shallow, even wave breaking.The resultant nonlinear mild-slope equation is computed by a finite volume method with a fourth-order Adams-Moulton predictor-corrector technique, rather than the usual finite difference or finite element scheme. The validity of the present nonlinear model is also satisfactorily verified with examples that cover a wide range of practical problems, such as Stokes wave to a second-order, cnoidal wave propagation on uniform depth, nonlinear wave shoaling on a sloping beach, and wave transformation on an elliptic shoal. By introducing an extra term representing energy dissipation due to wave breaking to the momentum equation, the present nonlinear model is suitable for the wave transformation within the surf zone on a uniform slope and on a composite profile with a bar on which second wave breaking occurs after wave recovering. The comparisons of computed wave height distributions around a jetty indicate the present nonlinear solutions are in very good agreement with the experimental results.



---

*Source: 341385-2014-06-04.xml*