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

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.

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 𝑘ℎ = 40, where 𝑘 is the wave number and ℎ 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 Section 2. 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: ∇⋅^→𝑢+𝜕𝑤𝜕𝑧=0 (1) and the incompressible Navier-Stokes equations in conservative form lead to Euler’s equations as 𝜕^→𝑢+𝜕𝑡^→𝑢⋅∇^→𝑢𝜕+𝑤^→𝑢1𝜕𝑧=−𝜌∇𝑝,𝜕𝑤+𝜕𝑡^→𝑢⋅∇𝑤+𝑤𝜕𝑤1𝜕𝑧=−𝜌𝜕𝑝𝜕𝑧−𝑔 (2)(3) in the horizontal and vertical directions of a Cartesian coordinates system, respectively. In (1)–(3), ^→𝑢=(𝑢,𝑣) and 𝑤 are the Eulerian velocity components in the horizontal plane and vertical direction, respectively; 𝑝 is the pressure, 𝜌 is the density of the fluid, 𝑔 is the acceleration due to gravity, 𝑡 is the time, and ∇=(𝜕/𝜕𝑥,𝜕/𝜕𝑦) is the horizontal vector operator, where 𝑥, 𝑦, and 𝑧 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 𝑝=0𝑧=𝜂𝑥,𝑦,𝑡,𝑤=𝜕𝜂+𝜕𝑡^→𝑢⋅∇𝜂𝑧=𝜂𝑥,𝑦,𝑡 (4)(5) and the bottom boundary condition on an impermeable seabed is 𝑤=−^→𝑢⋅∇ℎ𝑧=−ℎ𝑥,𝑦 (6) in which 𝜂(𝑥,𝑦,𝑡) is the water surface elevation and ℎ(𝑥,𝑦) 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 𝜕𝜂𝜕𝑡+∇⋅^𝜂→_−ℎ𝑢𝑑𝑧=0. (7) Second, after integrating the momentum equation (3) in the vertical direction, from an arbitrary depth 𝑧 to the free surface, and inserting the boundary condition for the pressure at the free surface, (4), it produces 𝑝𝜌=𝑔𝜂−𝑧+^𝜂_𝑧𝜕𝑤𝜕𝑡𝑑𝑧+^𝜂_𝑧^→𝑢⋅∇𝑤+𝑤𝜕𝑤𝜕𝑧𝑑𝑧. (8)

Substituting (8) into (2) renders the momentum equation for the 𝑥- and 𝑦-axis in the horizontal directions, respectively, as follows: 𝜕^⇀𝑢𝜕𝑡+𝑔∇𝜂+^→𝑢⋅∇^→𝑢𝜕+𝑤^→𝑢𝜕𝑧=−∇^𝜂_𝑧𝜕𝑤𝜕𝑡𝑑𝑧−∇^𝜂_𝑧^→𝑢⋅∇𝑤+𝑤𝜕𝑤𝜕𝑧𝑑𝑧. (9)

It is worth noting that the convective terms, (^→𝑢⋅∇^→𝑢+𝑤𝜕^→𝑢/𝜕𝑧) and (^→𝑢⋅∇𝑤+𝑤𝜕𝑤/𝜕𝑧) in (8) and (9), are transformed to ∇(^→𝑢⋅^→𝑢+𝑤²) and 𝜕(^→𝑢⋅^→𝑢+𝑤²)/𝜕𝑧 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, ^→𝑢_𝑧−∇𝑤≠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 ^⇀𝑢𝑥,𝑦,𝑧,𝑡=^𝑁_𝑚=1cosh𝑘_𝑚(ℎ+𝑧)cosh𝑘_𝑚ℎ^→𝑈_𝑚𝑥,𝑦,𝑡, (10) where 𝑈^→_𝑚(𝑥,𝑦,𝑡) is the 𝑚th component of the horizontal velocity vector at 𝑧=0, and 𝑘_𝑚(𝑥,𝑦) 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, 𝑁=1, it becomes ^⇀𝑢𝑥,𝑦,𝑧,𝑡=𝐹(𝑘ℎ,𝑧)𝑈^→,𝐹()=()𝑥,𝑦,𝑡𝑘ℎ,𝑧cosh𝑘ℎ+𝑧.cosh𝑘ℎ (11) This single-component model has remained popular and useful in practical coastal engineering applications.

The flow velocity 𝑤 in the vertical direction 𝑧 can be obtained upon depth integration of (1) from −ℎ to 𝑧 and the bottom boundary condition equation (6), rendering 𝑤𝑓()𝑥,𝑦,𝑧,𝑡=−∇⋅𝑘ℎ,𝑧𝑘^→𝑈,𝑓()=()𝑥,𝑦,𝑡𝑘ℎ,𝑧sinh𝑘ℎ+𝑧.cosh𝑘ℎ (12)

Inserting (11) into (7) gives the following depth-integrated continuity equation: 𝜕𝜂𝜕𝑡+∇⋅^𝜂_−ℎ𝐹𝑑𝑧𝑈^→=0. (13) Substituting (11) and (12) into (8), and after laborious manipulations, we obtain the pressure expression as 𝑝𝜌=𝑔𝜂−𝑧−𝑝_𝑎𝜕∇⋅^→𝑈−𝑝𝜕𝑡^𝑏𝜕∇ℎ⋅^→𝑈−𝑝𝜕𝑡_𝑐→𝑈⋅∇∇⋅𝑈^→−∇⋅𝑈^→∇⋅𝑈^→−𝑝_𝑑∇⋅𝑈^→∇ℎ⋅𝑈^→+𝑂(∇ℎ)², (14) where 𝑝_𝑎=^𝜂_𝑧𝑓𝑘𝑝𝑑𝑧,_𝑏=^𝜂_𝑧𝜕𝑓𝜕ℎ𝑘𝑝𝑑𝑧,_𝑐=^𝜂_𝑧𝐹𝑓𝑘𝑝𝑑𝑧,_𝑑=^𝜂_𝑧𝐹𝜕𝑓𝜕ℎ𝑘−𝜕𝐹𝑓𝜕ℎ𝑘𝑑𝑧. (15) 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 𝐹 in (11) and then integrating it from 𝑧=−ℎ to 𝜂, as well as after complex mathematical manipulations, produces a depth-integrated expression for the momentum equation in the form of 𝛼₁𝜕𝑈^→+𝛼𝜕𝑡²𝑔∇𝜂+𝛼𝑈₃^→⋅∇𝑈^→+𝛼₄∇⋅𝑈^→𝑈^→+𝛼₅∇ℎ⋅𝑈^→𝑈^→=𝛼₆∇𝜕∇⋅^→𝑈+𝛼𝜕𝑡⁷𝜕∇⋅^→𝑈+𝛼𝜕𝑡∇ℎ₈𝜕∇⋅^→𝑈+𝛼𝜕𝑡∇𝜂₉𝜕∇ℎ⋅∇𝜂^→𝑈+𝛼𝜕𝑡₁₀∇𝑈^→⋅∇∇⋅𝑈^→−∇⋅𝑈^→∇⋅𝑈^→+𝛼₁₁∇⋅𝑈^→(∇ℎ⋅∇)𝑈^→+𝛼𝑈₁₂^→⋅∇∇⋅𝑈^→−∇⋅𝑈^→∇⋅𝑈^→+𝛼∇𝜂₁₃∇⋅𝑈^→∇ℎ⋅𝑈^→∇𝜂+𝑂(∇ℎ)² (16) in which 𝛼₁=^𝜂_−ℎ𝐹²𝑑𝑧,𝛼₂=^𝜂_−ℎ𝐹𝑑𝑧,𝛼₃=^𝜂_−ℎ𝐹³𝛼𝑑𝑧,₄=−^𝜂_−ℎ𝐹𝑓²𝛼𝑑𝑧,₅=^𝜂_−ℎ𝐹²𝜕𝐹𝜕𝜕ℎ−𝑘𝐹𝑓𝑓𝜕ℎ𝑘𝑑𝑧,𝛼₆=^𝜂_−ℎ𝐹𝑝_𝑎𝛼𝑑𝑧,₇=2^𝜂_−ℎ𝐹𝑝_𝑏𝛼𝑑𝑧,₈=^𝜂_−ℎ𝐹𝜕𝑝_𝑎𝛼𝜕𝜂𝑑𝑧,₉=^𝜂_−ℎ𝐹𝜕𝑝_𝑏𝛼𝜕𝜂𝑑𝑧,₁₀=^𝜂_−ℎ𝐹𝑝_𝑐𝛼𝑑𝑧,₁₁=2^𝜂_−ℎ𝐹𝑝_𝑑𝑑𝑧,𝛼₁₂=^𝜂_−ℎ𝐹𝜕𝑝_𝑐𝛼𝜕𝜂𝑑𝑧,₁₃=^𝜂_−ℎ𝐹𝜕𝑝_𝑑𝜕𝜂𝑑𝑧. (17)

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 𝜂=𝜀𝜂₁+𝜀²𝜂₂+⋯,_→𝑈=𝜀𝑈^→₁+𝜀𝑈^2→₂+⋯, (18) where 𝜀 is the dimensionless nonlinearity parameter of wave amplitude, and the subscriptions 𝑖(𝑖=1,2,…) in (18) denote the 𝑖th order approximation.

Substituting (18) into (13) and (16) and taking Taylor series expansion on the still water level 𝑧=0, we obtain the continuity equation and the momentum equation, respectively, as follows: 𝜕𝜂𝜕𝑡+∇⋅⁰_−ℎ𝐹𝑑𝑧+𝜂𝑈^→=𝑂𝜀³,_0−ℎ𝐹²𝜕𝑑𝑧^→𝑈+𝜕𝑡⁰_−ℎ+1𝐹𝑑𝑧𝑔∇𝜂2⁰_−ℎ𝐹𝐹²−𝑓²𝑑𝑧−2⁰_−ℎ𝑘𝐹𝑑𝑧²𝑐⁴𝑔²∇𝑈^→⋅𝑈^→=⁰_−ℎ𝐹−𝐹²𝑘²∇𝜕𝑑𝑧∇⋅^→𝑈+𝜕𝑡_0−ℎ𝐹𝜕𝜕ℎ1−𝐹𝑘²𝜕𝑑𝑧∇⋅^→𝑈𝜀𝜕𝑡∇ℎ+𝑂³. (19)(20) These two equations are accurate to the second-order of 𝜀. Substituting the expressions of 𝐹 and 𝑓 given in (11) and (12) into (19) and (20), we obtain the final expression of the continuity equation and momentum equation as follows: 𝜕𝜂𝑐𝜕𝑡+∇⋅²𝑔+𝜂𝑈^→=0,𝑛𝜕𝑈^→1𝜕𝑡+𝑔∇𝜂+2𝑘1−3²𝑐⁴𝑔²∇𝑈^→⋅𝑈^→=1𝑐²∇𝑐²(1−𝑛)𝑘²𝜕∇⋅^→𝑈𝜕𝑡 (21)(22) in which 𝑐²=𝑔𝑘tanh𝑘ℎ,1𝑛=21+2𝑘ℎ().sinh2𝑘ℎ (23)(24)

The coefficient to the third term ∇(𝑈^→⋅𝑈^→) 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(𝑘²𝑐⁴/𝑔²), is the nonlinear effect of the vertical velocity. The value of (𝑘²𝑐⁴/𝑔²) 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 𝑐 and 𝑘 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 𝜕𝜂𝜕𝑡+∇⋅𝑄^⇀=0, (25) where 𝑄^→ is the linear flow rate of the waves defined by 𝑄^→=^0→_−ℎ𝑢𝑐𝑑𝑧=²𝑔^→𝑈. (26) On the other hand, (22) can be rewritten for the linear waves as 𝑛𝜕𝑈^→1𝜕𝑡+𝑔∇𝜂=𝑐²∇𝑐²(1−𝑛)𝑘²𝜕∇⋅^→𝑈.𝜕𝑡 (27) Substituting (26) into (27), then we obtain 𝜕𝑄^→+𝑐𝜕𝑡𝑛²1∇𝜂=𝑛∇(1−𝑛)𝑘²𝜕∇⋅^→𝑄−()𝜕𝑡1−𝑛𝑐²𝑘²∇𝑐²⋅𝜕𝑄^→.𝜕𝑡 (28) In the right-hand side of (28), 1𝑛∇(1−𝑛)𝑘²𝜕∇⋅^→𝑄1𝜕𝑡=−𝑛∇(1−𝑛)𝑘²𝜕²𝜂𝜕𝑡²=1𝑛∇𝑐²()𝜂=()1−𝑛1−𝑛𝑛𝜂∇𝑐²+𝑐²𝑛𝑐∇𝜂−²𝑛∇,−1𝑛𝜂𝑛∇(1−𝑛)𝑐²𝑘²∇𝑐²⋅𝜕𝑄^→()𝜕𝑡=−1−𝑛𝑛𝑐²𝑘²∇𝑐²𝜕∇⋅^→𝑄)𝜕𝑡+𝑂(∇ℎ²()=−1−𝑛𝑛𝜂∇𝑐²+𝑂(∇ℎ)². (29) The term of 𝑂[(∇ℎ)²] 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 𝜕𝑄^⇀+𝑐𝜕𝑡𝑛²∇𝑛𝜂=0. (30) Equations (25) and (30) are identical to the so-called time-dependent mild-slope equation in Copeland [8]. Eliminating flow rate 𝑄^→ 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 (𝑐), group velocity (𝑐_𝑔), and flow rate (𝑄^→) are usually given by 𝑐=𝑐_𝑔=𝑔ℎ,𝑄^→=ℎ+𝜂𝑈^→, (31) where 𝑈^→ 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 𝑈^→, given by 𝜕𝑈^→1𝜕𝑡+𝑔∇𝜂+2∇𝑈^→⋅𝑈^→=0. (32) For an irrotational flow, the term ∇(𝑈^→⋅𝑈^→)/2 can be replaced by (𝑈^→⋅∇)𝑈^→, 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, 𝑐²1=𝑔ℎ1−3𝑘²ℎ²,𝑐𝑛=^𝑔𝑐=11−3𝑘²ℎ². (33) Substituting (33) into (21) and (22) and replacing 𝑘𝑈^2→ by −∇(∇⋅𝑈^→) for the irrotational flow in shallow water condition, the present nonlinear mild-slope equation gives 𝜕𝜂𝜕𝑡+∇⋅ℎ+𝜂𝑈^→+ℎ³3∇²∇⋅𝑈^→𝜕=0,_→𝑈1𝜕𝑡+𝑔∇𝜂+2∇𝑈^→⋅𝑈^→=0. (34) 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: _∀𝜕𝜂𝜕𝑡𝑑∀+_∀𝑐∇⋅²𝑔+𝜂𝑈^→𝑑∀=0,_∀𝑛𝜕𝑈^→𝜕𝑡𝑑∀+_∀1𝑔∇𝜂+21−3𝑘²𝑐⁴𝑔²∇𝑈^→⋅𝑈^→=𝑑∀_∀1𝑐²∇𝑐²(1−𝑛)𝑘²𝜕∇⋅^→𝑈𝜕𝑡𝑑∀. (35) By the divergence theorem, (35) renders 𝜕𝜂𝜕𝑡_𝑚1=−∀_𝑠^→𝑁⋅𝑐²𝑔+𝜂𝑈^→𝜕𝑑𝑆𝑛𝜕𝑡^→𝑈−1∀1𝑐²_𝑚𝑠^→𝑁⋅𝑐²(1−𝑛)𝑘²∇⋅𝑈^→→𝐼𝑑𝑆_𝑚1=−∀𝑔_𝑠𝑁⋅𝜂^→𝐼−1𝑑𝑆∀12𝑘1−3²𝑐⁴𝑔²_𝑚𝑠^→𝑁⋅𝑈^→⋅𝑈^→→𝐼𝑑𝑆 (36) in which subscription 𝑚 denotes the mean quantities referring to the centre of the volume ∀, 𝑆 is the control surface, ^→𝐼 is the unit vector, and 𝑁^→ is the normal vector of 𝑆. By knowing the quantities at a starting time step (denotes as 𝐺_𝑚), the integral form of (36) can be used to calculate the change in 𝐺_𝑚 within a time increment, and thus update the 𝐺_𝑚 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 𝜂^𝑚+1=𝜂^𝑚+Δ𝑡2455Φ^𝑚−59Φ^𝑚−1+37Φ^𝑚−2−9Φ^𝑚−3,[𝐴]𝑈^→𝑚+1=[𝐴]𝑈^→𝑚+Δ𝑡2455Ψ^𝑚−59Ψ^𝑚−1+37Ψ^𝑚−2−9Ψ^𝑚−3. (37)(38) And the equations of the corrector-step are expressed by 𝜂^𝑚+1=𝜂^𝑚+Δ𝑡249Φ^𝑚+1+19Φ^𝑚−5Φ^𝑚−1+Φ^𝑚−2,[𝐴]𝑈^→𝑚+1=[𝐴]𝑈^→𝑚+Δ𝑡249Ψ^𝑚+1+19Ψ^𝑚−5Ψ^𝑚−1+Ψ^𝑚−2, (39)(40) where 1Φ=−∀²_𝐾=1^→𝑆⁺_𝐾⋅𝑐²𝑔+𝜂𝑈^→+_𝐾−^→𝑆⁻_𝐾⋅𝑐²𝑔+𝜂𝑈^→−_𝐾,𝑔Ψ=−∀²_𝐾=1^→𝑆⁺_𝐾⋅𝜂⁺_𝐾^→𝐼−^→𝑆⁻_𝐾⋅𝜂⁻_𝐾^→𝐼−1∀121−3𝑘²𝑐⁴𝑔²×²_𝐾=1^→𝑆⁺_𝐾⋅𝑈^→⋅𝑈^→+_𝐾^→𝐼−^→𝑆⁻_𝐾⋅𝑈^→⋅𝑈^→−_𝐾^→𝐼,[𝐴]𝑈^→=𝑛𝑈^→−1∀1𝑐²²_𝐾=1𝑐²(1−𝑛)𝑘²⁺_𝐾⋅∇⋅𝑈^→+_𝐾^→𝐼−𝑐²(1−𝑛)𝑘²⁻_𝐾⋅∇⋅𝑈^→−_𝐾^→𝐼,∇⋅𝑈^→=1∀²_𝐾=1^→𝑆⁺_𝐾⋅𝑈^→+_𝐾−^→𝑆⁻_𝐾⋅𝑈^→−_𝐾 (41) in which superscript 𝑚 represents the computed values of 𝜂, 𝑈^→, Φ, and Ψ at the 𝑚th time iteration step, and superscript 𝑚+1 denotes the predicted values at the 𝑚+1th time step. The subscript 𝐾 in (41) represents the two directions in 𝑥-𝑦 plane, in which superscript “+” denotes the quantities in the upstream of control-volume and “−” for those in the downstream. The section vectors denoted as ^→𝑆⁺_𝐾 and ^→𝑆⁻_𝐾 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 𝑈^→ 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 Figure 2, the quantities of variables 𝜂 and 𝑈^→ are set in the centre of the control volume. Then the flux terms, Φ and Ψ, in (37) to (40), can be expressed as Φ_𝑖,𝑗1=−2∀_𝑖,𝑗+∀_𝑖+1,𝑗+∀_𝑖,𝑗+1+∀_{𝑖+1,𝑗+1}⋅𝐴₁,Ψ_𝑖,𝑗𝑔=−2∀_𝑖,𝑗+∀_𝑖+1,𝑗+∀_𝑖,𝑗+1+∀_{𝑖+1,𝑗+1}⋅𝐴₂−14∀_𝑖,𝑗+∀_𝑖+1,𝑗+∀_𝑖,𝑗+1+∀_{𝑖+1,𝑗+1}×1−3𝑘²𝑐⁴𝑔²_𝑖,𝑗⋅𝐴₃, (42) where 𝐴₁=^→𝑆_{1𝑖+1,𝑗+1}⋅𝑅^→_{𝑖+1,𝑗+1}+𝑅^→_𝑖+1,𝑗+^→𝑆_1𝑖+1,𝑗⋅𝑅^→_{𝑖+1,𝑗−1}+𝑅^→_𝑖+1,𝑗−^→𝑆_{1𝑖−1,𝑗+1}⋅𝑅^→_{𝑖−1,𝑗+1}+𝑅^→_{𝑖−1,𝑗}+^→𝑆_{1𝑖−1,𝑗}⋅𝑅^→_{𝑖−1,𝑗−1}+𝑅^→_{𝑖−1,𝑗}+^→𝑆_{2𝑖+1,𝑗+1}⋅𝑅^→_{𝑖+1,𝑗+1}+𝑅^→_𝑖,𝑗+1+^→𝑆_2𝑖,𝑗+1⋅𝑅^→_{𝑖−1,𝑗+1}+𝑅^→_𝑖,𝑗+1−^→𝑆_{2𝑖+1,𝑗−1}⋅𝑅^→_{𝑖+1,𝑗−1}+𝑅^→_{𝑖,𝑗−1}+^→𝑆_{2𝑖,𝑗−1}⋅𝑅^→_{𝑖−1,𝑗−1}+𝑅^→_{𝑖,𝑗−1},𝐴₂=^→𝑆_{1𝑖+1,𝑗+1}⋅𝜂_{𝑖+1,𝑗+1}+𝜂_𝑖+1,𝑗^→𝐼+^→𝑆_1𝑖+1,𝑗⋅𝜂_{𝑖+1,𝑗−1}+𝜂_𝑖+1,𝑗^→𝐼−^→𝑆_{1𝑖−1,𝑗+1}⋅𝜂_{𝑖−1,𝑗+1}+𝜂_{𝑖−1,𝑗}^→𝐼+^→𝑆_{1𝑖−1,𝑗}⋅𝜂_{𝑖−1,𝑗−1}+𝜂_{𝑖−1,𝑗}^→𝐼+^→𝑆_{2𝑖+1,𝑗+1}⋅𝜂_{𝑖+1,𝑗+1}+𝜂_𝑖,𝑗+1^→𝐼+^→𝑆_2𝑖,𝑗+1⋅𝜂_{𝑖−1,𝑗+1}+𝜂_𝑖,𝑗+1^→𝐼−^→𝑆_{2𝑖+1,𝑗−1}⋅𝜂_{𝑖+1,𝑗−1}+𝜂_{𝑖,𝑗−1}^→𝐼+^→𝑆_{2𝑖,𝑗−1}⋅𝜂_{𝑖−1,𝑗−1}+𝜂_{𝑖,𝑗−1}^→𝐼,𝐴₃=^→𝑆_{1𝑖+1,𝑗+1}⋅𝑈̈_{𝑖+1,𝑗+1}+𝑈̈_𝑖+1,𝑗^→𝐼+^→𝑆_1𝑖+1,𝑗⋅𝑈̈_{𝑖+1,𝑗−1}+𝑈̈_𝑖+1,𝑗^→𝐼−^→𝑆_{1𝑖−1,𝑗+1}⋅𝑈̈_{𝑖−1,𝑗+1}+𝑈̈_{𝑖−1,𝑗}^→𝐼+^→𝑆_{1𝑖−1,𝑗}⋅𝑈̈_{𝑖−1,𝑗−1}+𝑈̈_{𝑖−1,𝑗}^→𝐼+^→𝑆_{2𝑖+1,𝑗+1}⋅𝑈̈_{𝑖+1,𝑗+1}+𝑈̈_𝑖,𝑗+1^→𝐼+^→𝑆_2𝑖,𝑗+1⋅𝑈̈_{𝑖−1,𝑗+1}+𝑈̈_𝑖,𝑗+1^→𝐼−^→𝑆_{2𝑖+1,𝑗−1}⋅𝑈̈_{𝑖+1,𝑗−1}+𝑈̈_{𝑖,𝑗−1}^→𝐼+^→𝑆_{2𝑖,𝑗−1}⋅𝑈̈_{𝑖−1,𝑗−1}+𝑈̈_{𝑖,𝑗−1}^→𝐼,_→𝑅_𝑖,𝑗=𝑐²_𝑖,𝑗𝑔+𝜂𝑈_𝑖,𝑗^→_𝑖,𝑗,𝑈̈_𝑖,𝑗=𝑈^→_𝑖,𝑗⋅𝑈^→_𝑖,𝑗=𝑈²_𝑖,𝑗+𝑉²_𝑖,𝑗. (43)

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, Δ𝑡≤𝑓_𝑐∀⋅min_→𝑐⋅^→𝑆 (44) in which 𝑓_𝑐 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, 𝑈^→, each time interval in the SOR iteration method needs to satisfy the following converge condition: Δ𝑒^⇀_𝑈=max^→𝑈^𝑚+1−𝑈^→𝑚_→𝑈^𝑚+1≤10⁻⁴ (45) in which Δ𝑒^⇀_𝑈 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 𝜂^𝑡_𝐵′=𝜂_inc^𝑡_𝐵′+𝜂_out𝐵^𝑡−𝜏_→𝑈^𝑡_𝐵′=𝑈^→_inc^𝑡_𝐵′+𝑈^→_out𝐵^𝑡−𝜏, (46) where 𝛼𝜏=Δ𝑟⋅cos_out𝑐,𝜂_out^𝑡_𝐵=𝜂^𝑡_𝐵−𝜂_inc^𝑡_𝐵,_→𝑈_out^𝑡_𝐵=𝑈^→𝑡_𝐵−𝑈^→_inc^𝑡_𝐵, (47) where the subscriptions inc and out represent the quantities of incoming and outgoing waves, respectively, Δ𝑟 is the length of 𝐵𝐵^′, 𝐵^′ is the phantom point of 𝐵, and 𝛼_out is the angle between outgoing wave and 𝐵𝐵^′, 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) 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) where subscriptions St and Sn 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 different 𝑘𝐻 values in intermediate depth (ℎ/𝐿=0.2) are depicted in Figure 4, in which 𝐻 is the wave height, 𝐿 is the wavelength and 𝑎_𝑜=𝐻/2. The temporal and spatial resolutions are Δ𝑡=𝑇/60 and Δ𝑥=𝐿/30, respectively, where 𝑇 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 𝑥/𝐿 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 𝑈^→ is obtained. The maximum nondimensional horizontal velocity calculated from (10) is plotted against dimensionless depth 𝑧/ℎ 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. Figure 6 demonstrates the simulated wave profiles of cnoidal waves in relative depth of ℎ/𝐿=0.05 for two different 𝐻/ℎ values. The temporal and spatial resolutions are Δ𝑡=𝑇/40 and Δ𝑥=𝐿/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 of ℎ/𝐿. Here ℎ/𝐿 is taken as 1/50. Figure 7 illustrates the simulated wave profiles in relative depth of ℎ/𝐿  = 0.02 for 𝐻/ℎ  = 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 𝐻/𝐻_𝑜 versus relative water depth ℎ/𝐿_𝑜 for wave steepness 𝐻_𝑜/𝐿_𝑜 ranging from 0.001 to 0.08, on uniform beach slope of 1/30, in which 𝐻_𝑜 and 𝐿_𝑜 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 ℎ/𝐿_𝑜. 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 ℎ/𝐿_𝑜, the shoaling coefficient increases as wave steepness increases. However, the solution computed by the linear wave shoaling is independent of the wave steepness 𝐻_𝑜/𝐿_𝑜 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 ℎ/𝐿_𝑜 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 Figure 11 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 𝐵 to 𝐸.

(a)

(b)

(c)

(d)

(e)

(f)

(g)

(h)

Figure 11 

Simulation results of wave height distribution along eight transects shown in Figure 9 and comparison with experimental data and other numerical models.

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, 𝑀_𝐷. In this way, (22) becomes 𝜕𝑈^→+𝑔𝜕𝑡𝑛1∇𝜂+𝑘2𝑛1−3²𝑐⁴𝑔²∇𝑈^→⋅𝑈^→=1𝑛𝑐²∇𝑐²(1−𝑛)𝑘²𝜕∇⋅^→𝑈𝜕𝑡+𝑀_𝐷. (50) Following Kabiling and Sato [31], the mixing term due to energy dissipation can be expressed by 𝑀_𝐷=𝜈_𝑒∇𝑈^2→, (51) where 𝜈_𝑒 is the eddy viscosity related to an energy dissipation factor, 𝑓_𝐷, due to breaking, 𝜈_𝑒=𝑔ℎ𝜔²𝑓_𝐷;𝑓_𝐷=𝛼_𝐷tan𝛽𝑔ℎ𝑄_𝑚−𝑄_𝑟𝑄_𝑠−𝑄_𝑟. (52) In (52), 𝑓_𝐷 is the energy dissipation factor, 𝜔(=2𝜋/𝑇) 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]), 𝑄_𝑚 is the amplitude of flow rate, and 𝑄_𝑟 and 𝑄_𝑠 are determined from the following equations: 𝑄_𝑠=0.40.57+5.3tan𝛽𝑔ℎ³,𝑄_𝑟=0.135𝑔ℎ³. (53) In (52), the dissipation factor 𝑓_𝐷 is set equal to zero outside the surf zone and in any region where 𝑄_𝑚<𝑄_𝑟.

On the basis of shallow water assumption for wave breaking, we can simplify ∇𝑈^2→=−𝑘𝑈^2→ and render 𝑀_𝐷=−𝑓𝑈_𝐷^→, thus (50) becomes 𝜕𝑈^→+𝑔𝜕𝑡𝑛1∇𝜂+𝑘2𝑛1−3²𝑐⁴𝑔²∇𝑈^→⋅𝑈^→+𝑓𝑈_𝐷^→=1𝑛𝑐²∇𝑐²(1−𝑛)𝑘²𝜕∇⋅^→𝑈.𝜕𝑡 (54) 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: 𝐻_𝑏𝐿_𝑜ℎ=0.171−exp−1.5𝜋^𝑏𝐿_𝑜1+15tan^4/3𝛽 (55) in which 𝐻_𝑏 and ℎ_𝑏 are the wave height and water depth of the breaking wave, respectively, and 𝐿_𝑜 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 Figures 12 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 Figure 12, 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 Figure 14. A clapotis 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 𝑦=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 sections 𝑦=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.

Conflict of Interests

The authors declare that there is no conflict of interests regarding the publication of this paper.

Acknowledgment

This material was based upon the work supported by the Ministry of Science and Technology, Taiwan, under Grant no. NSC 99-2221-E-005-116-MY3.