A FULLY NONLINEAR BOUSSINESQ MODEL FOR WATER WAVE PROPAGATION

A set of high-order fully nonlinear Boussinesq-type equations is derived from the Laplace equation and the nonlinear boundary conditions. The derived equations include the dissipation terms and fully satisfy the sea bed boundary condition. The equations with the linear dispersion accurate up to [2,2] padé approximation is qualitatively and quantitatively studied in details. A numerical model for wave propagation is developed with the use of iterative Crank-Nicolson scheme, and the two-dimensional fourth-order filter formula is also derived. With two test cases numerically simulated, the modeled results of the fully nonlinear version of the numerical model are compared to those of the weakly nonlinear version.


INTRODUCTION
It is well known that the Boussinesq-type equations are powerful tools to accurately estimate the wave conditions in the near-shore zone.Because they are nonlinear, they can simulate the propagation of large waves with the effects of topography and the nonlinear wave-wave interactions.A large number of researchers (e.g.Peregrine, 1967;Madsen et al. 1991;Nwogu, 1993;Chen and Liu,1995;Wei et al. 1995;Hong,1997;Agnon et al., 1999;Zou, 1999;Gobbi et al. 2000;Madsen et al., 2002) derived or improved different Boussinesq-type equations.It should be particularly mentioned that Hong(1997) presented high-order models with the dissipative terms for non-linear and dispersive wave.Almost at the same time, the different numerical models with the Boussinesq-type equations employed as the governing equations are developed or improved.Most of the numerical models(e.g.Abott et al., 1984;Madsen et al.1991;Nwogu, 1993;Wei and Kirby, 1995;Zhang et al.,2001;Shi et al.2001;Zhang et al.2010) are provided with the use of the finite difference method in the Cartesian coordinates.A review on the various numerical models has been given by Zhang et al.(2010).
In this paper, with the wave velocity potential function expressed as a function defined at an arbitrary water level, a set of fully nonlinear Boussinesq-type equations is derived from the Laplace equation and the nonlinear boundary conditions.The derived set of equations includes the dissipative terms and fully satisfies the sea bottom boundary condition.The equations with the linear dispersion relation accurate up to [2,2] padé approximation is qualitatively and quantitatively studied in details.The nonlinear characteristics are compared between the present fully nonlinear equations and the weakly nonlinear ones.A numerical model for wave propagation is described, and a two-dimension fourth order filter formula is also provided in this paper.With two test cases numerically simulated, the modeled results of the fully nonlinear version of the numerical model are compared to those of the weakly nonlinear version.This indicates that the effects of strong nonlinear terms on wave propagation are studied quantitatively.

Formulations in Non-dimensional Forms and Formal Solutions
A Cartesian coordinate system ( ) , , x y z is adopted, with z * measured upwards from the still water level.Consider a three-dimensional wave field with water surface elevation ( ) where L and d are the characteristic wave length and water depth, respectively; g is the acceleration of gravity ; ρ is the water density; 0 c is the wave celerity; * is the non-dimensional gradient operator; , where a is the characteristic amplitude.Then, the Laplace equation becomes: (2) The boundary condition at the sea bed can be written as: 0, ( , ) (3) The kinematic boundary condition on the free surface can be written as: The dynamic boundary condition on the free surface can be written as: When η is eliminated in these two equations, the kinematic-dynamic boundary condition on the free surface can be combined as: From the Bernoulli equation, the pressure in the wave can be written as: Where 2 2 ( ) Applying the gradient operator to Eq.( 5) yields: ( ) Assuming that the wave velocity potential ϕ can be expressed as the following power series formulations with respect to the parameter of β ： 0 0 , Substituting Eq.( 9) into Eq.( 2) and on the basis of Eq.( 3), one can find readily: 0 0 0 0 0, ( ) 0, ( ) 0, ( , , ) Thus, when The formal solutions of Eq.( 12) with the boundary condition at the sea bed(Eq.13)can be deduced: ( ) Substituting the expressions of n ϕ and its derivatives into Eqs.( 4)and ( 8), one can obtain the Boussinesq-type equations to any order of β : ( 1) 1 1 ( ) 0, and

Governing Equations
Based on Eqs.( 14) and ( 15), the different orders of explicit formulations of wave velocity potential can be obtained.According to Eqs.( 17) and ( 18), the governing equations which include all nonlinear terms that correspond to the dispersion relation accurate up to ( ) where * z is an arbitrary water level and * z sh = .In this paper, 0.5 s = is chosen.

NONLINEAR PROPERTIES OF GOVERNING EQUATIONS
It is assumed that the water depth is constant, namely, ( , ) .h x y const = .The free surface elevation, η , and the velocity in the x-direction, u are expanded with the use of the perturbation method as follows, respectively: (24) Substituting Eqs.( 23) and ( 24) into the one-dimensional forms of Eqs.( 21) and ( 22), one can obtain the following expressions at the first order and second order, respectively: Consider a wave train consisting of two small amplitude periodic waves with amplitudes 1 a and 2 a , and frequencies 1 ω and 2 ω .The water-surface elevation is given by where 1 k and 2 k are the respective wave numbers.The individual waves satisfy the first-order of the Boussinesq equations [(25a) and ( 25b)].Therefore, the horizontal particle velocity can be written as: The second-order wave will consist of a sub harmonic at the difference frequency Where is a quadratic transfer function that relates the amplitude of the second-order wave to the first-order amplitudes.
Differentiating Eq.(26a) with t and differentiating Eq.( 26b) with x firstly , and then manipulating algebraically yields: Substituting Eq.( 30) into Eq.(31), one can obtain: where To illustrate the effects of strong nonlinear terms more clearly, we further ignore the βε , 2 βε and 3 βε terms of Eqs.( 19) and ( 20) and only retain the weakly nonlinear terms , that is, the terms are relevant to ε .Thus, Eqs.( 19) and ( 20) are changed into the weakly nonlinear equations.In the similar way, the corresponding quadratic transfer functions can be obtained., compared to the quadratic transfer functions of the second-order Laplace equation, the set of weakly nonlinear Boussinesq equations underestimates the magnitude of the sum frequency and difference frequency by 42.7% and 62.2%, respectively; In contrast, the set of fully nonlinear Boussinesq equations underestimates the magnitude of the sum frequency by 34.9% and overestimates the magnitude of the difference frequency by 58.5%.Therefore, it is qualitatively illustrated that the precision of the fully nonlinear Boussinesq equations is higher than that of the weakly ones.

Finite Difference Scheme of Governing Equations
The improved Crank-Nicolson method, which consists of three stages, is performed iteratively, with the initial value given by the predictor-corrector method.Firstly, at any given time step n t ∆ , the values of the variables at ( ) ∆ are predicted with the use of the known values at t n t = ∆ .Secondly, the predicted values at ( ) ∆ are then used to compute the values at ( ) stages.Finally, the computed values at ( ) 1 n t + ∆ are then used as an initial estimation in an iterative scheme, and this is repeated until convergence.The approximations of the first-order time and spatial derivatives include the terms of ( ) ( ) ( ) , which involve third-order derivatives.Thus, it is necessary to make correction for the first-order terms.The method of correction is performing back substitution of the truncation terms to obtain the proper calculation results.The correction schemes of the first-order derivatives and the difference discretization of governing equations can be referred to Zhang(2000) and Zhang et al.(2001).

General Boundary Conditions for Open and Fixed Boundaries
Zhang(2000) derived the general boundary conditions for open and fixed boundaries.Assuming the artificial down-wave boundary at M x x = , on the basis of the theory of linear waves, the wave potential at the boundary

M x x =
with arbitrary reflectivity is given by and the free surface elevation is given by 0 where 0 a is the incident wave amplitude; ( , , , ) ( , , , ) where 1/ 2 2 2 1 2 cos 1 2 cos( ) If x ∆ is defined, the expressions for the boundary conditions can be deduced from Eq.( 44): ( , , , ) ( , , , ) The formulations of y A and y τ can be obtained from Eqs.( 45) and ( 49) with x ∆ replaced by y ∆ , respectively.If the waves are outgoing at the down-wave boundary or at the lateral boundaries and the wave directions vary severely, the sponge layer should be set near the corresponding boundaries, respectively, to minimize the numerical errors generated by the waves reflected from the boundaries into the solution for the area of interest.The back ends of the sponge layers are at the corresponding boundaries, respectively.According to the numerical tests (Zhang,2000;Zhang et al, 2009), in the case of regular waves and irregular waves, the general boundary conditions can effectively reflect the influence of different boundaries, and can easily be utilized.

The Improved Specifying Boundary Conditions
If the upstream boundary conditions are prescribed only according to the incident wave only, the time-dependent numerical models can not effectively simulate the wave field when the physical or spurious reflected waves become significant.Zhang et al.(2009) provided a new approach to specifying the incident wave boundary conditions combined with a set sponge layer to damp the reflected waves towards the incident boundary.The improvement method is very simple, and it is described as follows: (1) A sponge layer is set in the neighborhood of the incident boundary.
(2) The analytical solutions are searched in the region where the sponge layer is set.The steps of searching them include: the relevant incident wave factors, such as waveheight, wave period and wavelength, should firstly be provided at the incident wave boundary, according to the relevant wave theory, for example, Stokes wave or cnoidal wave theory, etc.; the analytical solutions of the dependent variables ,which are represented by η and U  in this paper, are then obtained in the range of sponge layer on the basis of the corresponding wave theory.
(3) In the spatial range of the sponge layer, the physical or spurious reflected waves, that is, the "error" , at each calculation step, are corrected as follows: The Formulation of Numerical Filter How to filter numerical noise is very important for developing a numerical model because the numerical noise can make the model collapse.The numerical noise is often encountered in developing the numerical model based on the Boussinesq-type equations, since the governing equations include the high order derivative terms.In this paper, the two-dimensional fourth-order model (Shapiro,1970) is employed: ( 2 ) 2 Where S is the smoothing factor, , i j F represents the original values at point ( , ) i j and * , i j F represents the new values after numerical filtering at the same point.The numerical filtering response function correspondence to Eq.( 53) is: Based on Eqs.( 53) and ( 54), with the use of maple software, The 81-point two-dimensional formulation in the interior of calculation domain is as follows: , , 1 65536

TEST CASES
To study the effects of nonlinear terms quantitatively, the two versions of the numerical model are employed to simulate the wave propagation.That is, one version is fully nonlinear, its governing equations include all nonlinear terms proportional to powers of ε ; and the other version is weakly nonlinear, its governing equations only include the nonlinear terms proportional to the first order of ε .
The fully nonlinear version is referred to below as the FN model, and the weakly nonlinear version is referred to below as the WN model.The experiments of Ito & Tanimoto (1972) and Berkhoff et al.(1982) are widely used to testify the relevant numerical model.Thus, they are also used in this paper.The dissipative terms are omitted in this paper because the calculation domain is not large (Zhang et al.2005).Ito & Tanimoto(1972) performed an experiment of wave propagation on a submerged shoal with concentric contours.The experimental layout is shown in Fig. 2. The incident wave period T and wave height 0 H are 6.3s and 1.0m , respectively.The lateral boundaries are the fully reflective walls and the down-wave boundary is perfect outflow.The grid steps x ∆ and y ∆ are both chosen as 4.03763m , the time step t ∆ is chosen as / 32 T .The filter is applied every 100 time steps.

COASTAL ENGINEERING 2010
9 It is shown in Fig. 3 that the modeled results along two sections are compared to experiment data, where the solid line, dashed line and symbols represent the results of FN model, those of WN model and experimental data, respectively.It is found that the results of both models correspond to the experimental data.The quantitative mismatch can be measured by the index proposed by Wilmott (1981) as follows: COASTAL ENGINEERING 2010 10 where ( ) x j are experimental data, x is the mean value of ( ) x j , ( ) y j are the modeled results.0 d = means a completely mismatch and 1 d = indicates a perfect agreement.
The calculated index d can be found in Table 1.Although the differences between the WN model and FN model can't be obviously reflected in Fig. 3, it is indicated in Table 1 that the precisions of FN model are indeed higher than those of WN model.Wave Propagation on a Plane Beach with an Elliptic Shoal Berkhoff et al.(1982) conducted an experiment to observe the wave transformation due to reflection, refraction and diffraction.The experimental layout and 8 sections for collecting wave data are shown in Fig. 4. The incident wave period and wave height are1.0s and 0.046m ，respectively.The grid steps x ∆ and y ∆ are chosen as 0.125m and 0.25m , respectively.The time step t ∆ is chosen as / 100 T  The comparisons between the modeled results with experimental data are shown in Fig. 5, where the solid line, dashed line and symbols represent the results of FN model, those of WN model and experimental data, respectively.It can be seen that the modeled results from the FN model are obviously better than those of the WN model.The quantitative comparisons between different models with experimental data are shown in Table 2.It clearly indicates that the accuracy of the FN model is higher than that of the WN model.The average value of index d from the FN model is 0.865.In contrast, it is 0.810 based on the WN model.Thus, if the governing equations include the higher order nonlinear terms, the agreement between the corresponding numerical model and the experiment can be improved.

CONCLUSIONS
A set of high-order fully nonlinear Boussinesq-type equations is derived from the Laplace equation and the nonlinear boundary conditions.The derived equations include the dissipation terms and fully satisfy the sea bed boundary condition.The equations with the linear dispersion accurate up to [2,2] padé approximation is qualitatively and quantitatively studied in details.The difference between the fully nonlinear version of the equations and the weakly version is found by analyzing their secondorder transfer function.A numerical model for wave propagation is developed with the use of iterative Crank-Nicolson scheme, and the two-dimensional fourth-order filter formula is also derived in this paper.With two test cases numerically simulated, the modeled results of the fully nonlinear version of the numerical model are compared to those of the weakly nonlinear version.It is found that the highorder nonlinear terms in the governing equations can indeed improve the precision of the numerical model.

.
It can be expressed as: between the quadratic transfer functions of the Boussinesq equations with those of the second-order Laplace equation are shown in Fig.1.In Fig.1

Fig. 1
Fig.1 Comparisons between quadratic transfer functions of Boussinesq-type equations and those of second order Laplace equation ( : Second order Laplace equation, ：weakly nonlinear Boussinesqtype equations;：fully nonlinear Boussinesq-type equations) coefficient and phase difference, respectively; x k and y k represent the wave number in the x-direction and y-direction, respectively.If the down-wave boundary is an open one, r k should be set as zero; and if it is a fully reflection one, r k should be set as unit and r ε should be set as zero.

( 4 )
are then damped with the sponge layer used.In the end, the calculated values, c η and c U  the modified calculated values of the corresponding variables, respectively.As for the form of ( ) x γ , Zhang et al.(2009) adopted the following formulation:

Table 2 . The calculated index d with the experimental data as and modeled results as
( )