NEW BOUSSINESQ SYSTEM FOR NONLINEAR WATER WAVES

In this paper, a new Boussinesq water wave theory is derived which can simulate highly dispersive nonlinear waves, their depth-varying velocities, and wave-induced currents, from very deep, but still finite, depths through the surf zone to the shoreline.. Boussinesq scaling is employed. We removed the irrotationality assumption by using polynomial basis functions for velocity profile which are inserted into basic equations of motion. Keep terms up to the desired pproximation level and solve the coupled weighted residual system together with vertically integrated mass equation. The computational cost is similar to normal Boussinesq theories although there are more unknown varibles to be solved than that in normal Boussinesq models. Because we can reduce the number of the coupled equations by multiplying some coefficients and subtracting from each other which means the matix to be solved is in similar size as normal Boussinesq models. The models show rapid convergence to exact solutions for linear dispersion, shoaling, and orbital velocities; however, properties may be simultaneously and substantially improved for a given order of approximation using asymptotic rearrangements. This improvement is accomplished using the large numbers of degrees of freedom inherent in the definitions of the polynomial basis functions either to match additional terms in a Taylor series, or to minimize errors over a range. Future work will be focused on rotational performance in 2D model by including viscosity,breaking and turbulence.


INTRODUCTION
Boussinesq-type theories for water waves started from 1960s by Peregrine and have progressed significantly over the past two decades when people started using various methods of asymptotic rearrangment.Researchers have made great efforts to extend the accuracy of linear and nonlinear characteristics to deeper water, moving from the mildly nonlinear, mildly dispersive depth averaged formulation of Peregrine (1967) to the enhanced dispersion and shoaling found in Madsen et al. (1991), Madsen and Sorenson (1992) and the redefinition of the velocity variables in Nwogu (1993), both of which increased accuracy significantly.Further work increased nonlinearity from the mildly nonlinear equations that existed previously to so-called "fully nonlinear" equations with considerably more accurate nonlinear properties (Wei et al., 1995;Madsen and Schaffer, 1998;Kennedy et al., 2001).Formal expansions to higher order increased the accuracy of all properties (Madsen and Schaffer, 1998;Gobbi and Kirby, 1999;Gobbi et al., 2000) but at the cost of much more complex equations.Extensions to include wave breaking and shorelines have made these into true surf zone models able to represent waves and wave-induced currents including wave setup, rip currents, and longshore currents (e.g.Schaffer and Madsen, 1993;Sorensen et al., 1998;Chen et al.,2000;Kennedy et al., 2000a,b;Lynett et al., 2002;Nwogu and Demirbilek, 2010;Bonneton et al., 2011a).However, even though many aspects of Boussinesq theory have progressed significantly, some aspects have reached apparent limits: 1. Higher order equations are extremely complex, so much so that their adoption has been hindered.2. Most Boussinesq equations feature partial irrotationality assumptions which are entirely violated in the surf zone.This means that velocity profiles over the water column will be significantly in error.
3. The most straightforward Taylor series expansions for velocity have a limited radius of covergence, and will diverge strongly in moderate water depths (Kennedy et al., 2002;Madsen and Agnon, 2003).
For these situations, higher order systems could give worse results than lower order systems.Partial solutions have been found: convergent formulations of very high order have been derived and tested (Madsen and Agnon, 2003;Lynett and Liu, 2004;Schaffer, 2009).However, no high-order Boussinesq-type systems exist that are able to simulate accurately highly nonlinear waves, their depth-varying velocities, and wave-induced currents, from very deep, but still finite, depths through the surf zone to the shoreline.Smoothed Particle Hydrodynamics (SPH) and Volume of Fluid (VOF) Navier-Stokes methods can simulate these situations with good accuracy, but have computational costs so large that modeling a 20 wavelength by 20 wavelength coastal region is effectively beyond their capabilities.Here, we derive and test systems of equations for nonlinear water wave transformation.Boussinesq scaling is employed, but without the partial irrotationality assumption of most systems so that rotational surf zone fows may be modeled naturally.
The systems may be extended to arbitrary order and show excellent convergence towards exact solutions for dispersion, shoaling, and orbital velocities.The end results show a resemblance to both Boussinesq and Green-Naghdi systems, and may be recast into different forms that are useful for introducing weakly or moderately nonhydrostatic properties into hydrostatic ocean models.Most of the asymptotic rearrangement techniques used for Boussinesq models may also be employed here SCALING Boussinesq-shallow water scaling for non-dimensional variables is: The dimensionless continuity equation and kinematic boundary conditions are: where Integrating (3) from bottom to surface and applying kinematic boundary conditions give a mass equation in conservation form, Euler's 2-D equations of motion for an inviscid fluid in dimensionless form are, i = 1, 2; j = 1, 2. where µ = k 0 h 0 is a measure of dispersion.Boussinesq equations attempt to go to as high a µ value as possible to increase the range of application.Integrate(9) from z to η , the pressure will be,(P(η) = 0) The pressure profile is similar to most Boussinesq theory and results in mixed space-time derivatives are like typical Boussinesq models.Then plug P(z) back into equation ( 8) with weighting function f m (q): Assume: β n = n when n is even; β n = n + 1 when n is odd.The velocity profile is assumed as a summation of N terms.Each term consists of a µ term, horizontal coefficients u n and basis functions f n .The order of µ stands for the approximation level.The u n depend on x, y, t while the basis functions f n only depend on vertical coordinate q.We removed the irotationality assumption by solving more unknowns( u n ) and seperating the vertical conponent out for each term so that it is easy to get ∂u/∂z.
Plug u and w into ( 7), ( 10) and keep terms up to O(µ N ).The system looks like coupled equations system.
Although the mass equation is explicit, the three coupled momentum equations would seem to need to be solved simultaneously for u 0,t , u 1,t and u 2,t , which would increase the computational cost.However, if it is realized that mixed space-time derivatives only occur for u 0 (i.e.u 0,xxt and related terms), this may easily be reduced to a form where u 0,t and its mixed derivatives are the only unknowns.To do this, u 1,t and u 2,t must be eliminated from one momentum equation, say m = 0.The basic form of the equation will not change, but the integrals will.If we define g0 = (g 0 − d 0 g 1 − e 0 g 2 ), where replacement of g 0 with g and equivalently for all other integrals (e.g.0n is replaced with ˜ 0n ) will result in a replacement equation for m = 0 that does not contain u 1 ;t or u 2 ;t terms.This will be in the standard O(µ 2 ) Boussinesq form and may be arranged into a tridiagonal matrix for 1D to solve for u 0,t using methods that are well known.Velocities u 1,t and u 2,t will in many numerical representations then be solvable purely locally, which is straightforward.
O(µ 4 ) system will look like the fully nonlinear O(µ 2 ) system with additional higher order linear terms.Due to the complexity, the higher order terms will be linearized which makes the system weakly O(µ 4 ) system but should give good nonlinear results.
The mass equation looks similar just extended to O(µ 4 ) The momentum equations will be Similarly to the process for N = 2, it is possible in higher orders to reduce the number of weighted momentum equations that must be solved simultaneously by eliminating u N−1,t and u N,t from weighted momentum equations m = (0; N − 2) through partial Gaussian elimination using momentum equations m = (N − 1; N).The remaining N − 1 equations may be solved for u 0,t to u N−2,t , and the results may then be used to solve for u N−1,t and u N,t using momentum equations m = (N − 1; N).Again, this is possible because mixed space-time derivatives do not appear in the weighted momentum equations for n = [N − 1; N].We will not give details, but the process is very similar to that in (3.6), except that new coeffcients must be found for each equation in m = (0; N − 2).

PROPERTIES AND ASYMPTOTIC REARRANGEMENT OPTIMIZATION
In order to understand the range of applicability for different orders of approximation, the magnitude and behavior of the error, and to compare to previous work, it is useful to obtain theoretical representations for linear dispersion and shoaling.These are found from the first two orders of a multiple scales expansion, with the linear dispserion found at first order and the shoaling properties at second order.
The multiple scales expansion in space (assuming that the system is steady in time) then has fast and slow spatial derivatives with similar expressions for η.
Defining the water depth to be only slowly varying and thus have only slow X 1 derivatives, the horizontal derivative of the pressure equation, to O( h ,X 1 ), is then A multiple scales expansion for the mass equation gives Now expand each component in a perturbation series: η = η (0) + η (1) , u n = u (0) n + u (1) n .Insert these into ( 7) and ( 11) , and collect all terms into the various orders to get perturbation equations: At O( 1) At O( ), where Taking each component (at any order) where ψ ,x ≡ k(X 1 ), ψ ,t ≡ −σ all necessary derivatives are then u n,xxt = iσk 2 ũn e iψ u n,xX 1 t = σk ,X 1 ũn e iψ + σk ũn,X 1 e iψ u n,X 1 xt = σk ũn,X 1 e iψ (29) At first order, the system is closed and the linear dispersion relations may be found by setting the determinant of (24) to zero.They could be compared to exact linear dispersion At second order, the system is still not closed because there are more unknowns than equations, and more constraints must be developed.First, we specify η(2) 0 = 0, so that the wave height on the slope is the same as on a flat bed.Next we relate k ,X 1 to h ,X 1 .If we write the dispersion relation as where the derivative of Q is with respect to kh.This works for all dispersion relations, when the appropriate expressions are used for Green-Naghdi or exact quantities.Finally, we must relate ũ(0) n,X 1 to η(0) ,X 1 .As shown in section ??, we may relate these quantities through the dispersion matrix.All of these relationships will have the form which, once it is noted that η(0) ,X 1 = η(0) ,h h ,X 1 , closes the system.The revised equations are then All η(0) ,h terms are then moved from the RHS to the LHS (as they are unknowns) and a linear matrix is then solved for η(0) ,h and ũ(1) n , n = 0, 1, . . ., N. When written in the form they may be compared to the infinite order linear Airy solution Obviously we have a lot of choice of basis functions like monomials or polynomials.Here we use Shifted Legendre Polynomial basis function which is orthognal with good properties but still can be improved by asymptotic rearrangement.
Figure (1) shows linear dispersion and shoaling for orders of approximation O(µ 2 , µ 4 , µ 6 , µ 8 ) using Gauss-Legendgre basis functions f n ≡ T n (q).A clear increase in accuracy is seen for increasing level of approximation.At O(µ 2 ), phase speeds have several percent error by kh = 1.5, while O(µ 4 ) remains good until at least kh = 6.By O(µ 6 ), accuracy extends to kh = 15, while the highest level of approximation, O(µ 8 ), has accuracy extending past kh = 20.Since the nominal deep water limit for water waves is kh = π, these higher levels of approximation are very accurate.
For lower order solutions, analytical representations may be found for phase speeds.At O(µ 2 ), the phase speed using Gauss-Legendre basis functions is which is the same as is found for Peregrine's depth-averaged Boussinesq equations, and for Green-Naghdi level I theory (Shields and Webster, 1986)  which is identical to Green-Naghdi theory III.At levels greater than O(µ 4 ), dispersion becomes too complex for closed form representations and values at each wavenumber must be solved numerically.
For lower order systems, it is possible to arrive at dispersion results for generalized basis functions.If we define then the general dispserion relation for an O(µ 2 ) system with any choice of (a, b, c) will be Thus, although there appear to be three free coefficients, only one combination has any influence on dispersion at O(µ 2 ).For the shifted Legendre polynomial basis functions (suitably normalized so that the coefficient of the highest degree polynomial in each basis function is unity), we find b − ac = 1/3, while to arrive at the Pade [2,2] approximant as seen in figure(2) (Madsen et al., 1991) we set b − ac = −1/5.For the simple basis functions f n = q n , we find that b − ac = 0 which does not yield accurate dispersion.If we define For the shifted Legendre polynomial basis functions we used, the dispersion relation at O(µ 4 ) will be The dispersion relation at O(µ 6 ) will be The analytical solutions for shoaling using Shifted Legendre basis functions will not be shown here while they are plot in figure(1).
Second harmonic could be optimized simultaneously.Figure (3) shows the ratio of η 2 to the Stokes solution for O(µ 2 )andO(µ 4 ) system using shifted legendre polynomials and optimized basis functions.

DISCUSSION AND CONCLUSIONS
For extension into the surf zone and through to the shoreline, additional viscous/turbulent stresses need to be specified corresponding to breaking wave dissipation and bottom stresses; these are under development and will be detailed in a future publication.The present formulation has the great advantage that because it is derived without any irrotationality conditions, it will be able to make use of more standard turbulence closures and thus needs to make fewer ad-hoc breaking assumptions than many other Boussinesq-type breaking models.
The generalization of velocity basis functions combined with Boussinesq scaling allows for asymptotic rearrangements similar in character to those of Nwogu (1993); Madsen and Schaffer (1998) and others, and with similar improvements in accuracy.For linear properties of O(µ 2 ) and higher, generalized analysis using arbitrary basis functions becomes highly complex but the use of shifted Legendre polynomials, which have excellent orthogonality properties, provides accuracy that is considerably higher than the formal level of approximation for both linear dispersion and shoaling (O(µ 2N−2 )).This is particularly obvious when compared with simple monomial basis functions.
Fully nonlinear systems up to O(µ 2 ) are straightforward to derive and code by hand.Linear components up to arbitrary order may also be explicitly written and coded without great difficulty, but the nonlinear components become extremely complex.As such, O(µ 2 ) may represent a practical limit to nonlinearity for these types of systems.It is quite possible to develop codes that automatically sum and integrate the nonlinearity without ever writing down the system on paper.However, related developments using many of the same concepts as the present paper, including Boussinesq and rotational shallow water scalings and asymptotic rearrangement, may prove to be better candidates for very high order representation of nonlinearity.In any case, as mentioned in the introduction, the assumption of a single-valued free surface η(x; y; t) will impose an upper limit on surf zone accuracy, particularly in regions of strong plungers.

Figure 1 :
Figure 1: Shoaling and Dispersion for different approximation levels using Shifted Legendre Polynomial basis functions