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

1 Department of Civil Engineering, National Chung Hsing University, Taichung 402, Taiwan 2Department of Leisure and Recreation Management, Chihlee Institute of Technology, New Taipei 220, Taiwan 3Department of Marine Environment & Engineering, National Sun Yat-Sen University, Kaohsiung 804, Taiwan 4 School of Civil & Resources Engineering, University of Western Australia, Crawley, WA 6009, Australia


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 mildslope 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][3][4][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 2 Mathematical Problems in Engineering 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 halfstrips of constant depths; their model was then extended with application to nonlinear water waves propagating in variable bathymetry regions [18].On the other hand, Boussinesqtype 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.

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: and the incompressible Navier-Stokes equations in conservative form lead to Euler's equations as in the horizontal and vertical directions of a Cartesian coordinates system, respectively.In ( 1)-( 3), ⃗  = (, V) 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 and the bottom boundary condition on an impermeable seabed is in which (, , ) is the water surface elevation and ℎ(, ) the local water depth.

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 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 Substituting ( 8) into (2) renders the momentum equation for the and -axis in the horizontal directions, respectively, as follows: It is worth noting that the convective terms, ( ⃗  ⋅ ∇ ⃗  +  ⃗ /) and ( ⃗ ⋅∇+/) in ( 8) and ( 9), are transformed to ∇( ⃗  ⋅ ⃗  +  2 ) and ( ⃗  ⋅ ⃗  +  2 )/ 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 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 ℎ . ( 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) Inserting ( 11) into (7) gives the following depth-integrated continuity equation: Substituting (11) and ( 12) into (8), and after laborious manipulations, we obtain the pressure expression as 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 in which 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 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: 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: in which The coefficient to the third term ∇( ⃗  ⋅ ⃗ ) on the lefthand-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( 2  4 / 2 ), is the nonlinear effect of the vertical velocity.The value of ( 2  4 / 2 ) 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 mildslope 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..Upon dropping the nonlinear term in (21) for the sake of linear (Airy) wave theory, this equation can be rewritten as

Linear Mild-Slope Equation
where ⃗  is the linear flow rate of the waves defined by On the other hand, ( 22) can be rewritten for the linear waves as Substituting ( 26) into ( 27), then we obtain In the right-hand side of (28), The term of [(∇ℎ) 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 Equations ( 25) and ( 30) are identical to the so-called timedependent mild-slope equation in Copeland [8].Eliminating flow rate ⃗  from ( 25) and ( 30) renders the elliptic mild-slope equation in Berkhoff [1]..For long wave theory in shallow water, wave celerity (), group velocity (  ), and flow rate ( ⃗ ) are usually given by

Nonlinear Long Wave Equation
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 For an irrotational flow, the term ∇( ⃗ ⋅ ⃗ )/2 can be replaced by ( ⃗ ⋅∇) ⃗ , thus rendering the usual expression for the nonlinear long wave equation.

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, 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 Equations ( 34) are the standard form of traditional Boussinesq equations.

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: By the divergence theorem, (35) renders 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].

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 And the equations of the corrector-step are expressed by where 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.
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 where

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, 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: ) ≤ 10 −4 (45) in which Δ ⇀  denotes the absolute maximum error.

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 where 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.
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 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 where subscriptions St and Sn represent the tangent and normal components at the boundary, respectively.

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 twodimensional 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.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].

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.
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.
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.[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.

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.
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 .

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 Following Kabiling and Sato [31], the mixing term due to energy dissipation can be expressed by where ]  is the eddy viscosity related to an energy dissipation factor,   , due to breaking,  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: 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 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: 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.
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 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.

Conclusions
In this paper, the development of a second-order timedependent 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.

Figure 2 :
Figure 2: Mesh arrangement and corresponding variables around the centre of a control volume.

Figure 3 :
Figure 3: Sketch of the offshore open boundary.

Figure 4 :Figure 5 :
Figure 4: Simulated Stokes waves in intermediate water depth compared with analytical solution (circles: present solution, line: 2nd-order Stokes wave theory).

Figure 7 :
Figure 7: Simulated wave profile for solitary waves in very shallow water compared with analytical solution (circles: present solution, line: solitary wave theory).

Figure 8 :
Figure 8: Comparison on the wave height from the present nonlinear wave shoaling solution with other solutions.

5 Figure 10 :
Figure 10: Simulation results on wave height distribution behind an elliptic shoal.

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

13 HFigure 12 :Figure 13 :
Figure 12: Comparisons among numerical solutions and laboratory experiment for wave transformation on uniform slope.

Figure 14 :
Figure 14: The computed wave field in the vicinity of a jetty.

Figure 15 :
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.
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.
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 Figure4, in which  is the wave height,