October 13, 2008

9:18

WSPC/101-CEJ

00187

Coastal Engineering Journal, Vol. 50, No. 4 (2008) 369–395
c World Scientiﬁc Publishing Company and Japan Society of Civil Engineers
A NUMERICAL MODEL ON THE INTERACTION BETWEENNEARSHORE NONLINEAR WAVES AND STRONG CURRENTS

YE TIAN∗ and SHINJI SATO†
Department of Civil Engineering, The University of Tokyo,
7-3-1 Hongo, Bunkyo-ku, Tokyo 113-8656, Japan∗tianye@coastal.t.u-tokyo.ac.jp†
sato@coastal.t.u-tokyo.ac.jp

Received 4 October 2007
Revised 9 June 2008
A new model for nearshore nonlinear water wave is established through the Galerkin
Chebyshev spectral approach on the vertical direction. The model is based on the spatialtemporal separation conception. The Galerkin Chebyshev spectral approach is applied tosolve the Laplace equation, while the Zakharov’s expression of the free surface boundary
condition is considered as the evolution equation of the free surface and integrated temporally. The accuracy and eﬃciency of the model are conﬁrmed by the results of simulationof water waves over even bottom, nonlinear wave shoaling and harmonic generation over a
submerged bar. Finally, the model is used to study the wave blocking phenomenon due tostrong opposing currents. Both the location of the blocking point and the wave structurenear the point are well presented.
Keywords: Chebyshev spectral approach; wave blocking; strong currents.

1. Introduction
The interaction between nearshore waves and currents is a very important andfundamental topic. As the natural consequence of the wave encountering strong opposing current, a striking nonlinear phenomenon called wave blocking often appears
in the region of river mouth. The numerical study of the wave blocking phenomenon
is necessary for a better understanding of the physical processes of nearshore hydrodynamics and morphological changes near the river mouth.
∗ Corresponding author.

369

October 13, 2008

370

9:18

WSPC/101-CEJ

00187

Y. Tian & S. Sato

However, the numerical simulation of the wave blocking phenomenon demands
eﬃcacious numerical model which is capable of dealing with the inherent high nonlinearity and large wave number of this problem. The classical Boussinesq equations(e.g. [Peregrine, 1967; Nwogu, 1993]) do not have a correct form of Doppler shift inconnection with wave-current interaction. Chen et al. [1996] extended the Nwogu[1993]’s Boussinesq equations incorporating correct Doppler shift to simulate thewave transformation with a strong opposing current. However, their dispersion relation was not able to handle the large wave number generated by the ambient current.
Strong bottom friction was also applied in order to stablize the simulation, which
eliminated the wave structures near the blocking point. Madsen et al. [1999] treatedthe case with the same computational conditions through their enhanced high order
Boussinesq type equations and showed a great advancement.
Stimulated by the rapid advancement in the capacities of computer, inspired byMadsen’s enhanced Boussinesq type equations, the nonlinear mild slope equations by
Isobe [1994] and the newly developed High Order Spectral method for the simulationof water waves from Mei et al. [2005], a new model is conceived and implemented forthe 1DH cases in this study. The new model borrows the philosophy of eliminating
the vertical coordinate with approximation of the ﬂow ﬁeld by the polynomials
and the temporal-spatial decoupling-coupling scheme from the High Order Spectral
method. Chebyshev polynomials are adopted to solve the spatial equation. The
nonlinearity is dealt with high order.
It is expected that (1) the new model will greatly enhance the dispersion relationand nonlinear properties along with moderate increment of the need for computational resources, compared with conventional models, and (2) the 1DH model canbe extended to 2DH without more theoretical treatment.
2. Physical description of the problem and governing equations
The physical domain of the nearshore ﬂow ﬁeld is described in Fig. 1. The continuous
water body is enclosed by a free surface, a ﬁxed bottom and two lateral boundaries.
Under the irrotational and inviscid assumption, potential ﬂow theory is convenient for theoretical and numerical treatment. Deﬁne a Cartesian coordinate systemr ≡ (x, y, z), and let (x, y) = x be the horizontal coordinates, z the vertical coordinate, positive upward, and z = 0 the mean free surface. The ﬂow is described by avelocity potential Φ(r, t), which satisﬁes the Laplace equation,∇2 Φ(r, t) = 0(1)

where Φ(r, t) is deﬁned as ∇Φ(r, t) = u. On the free surface, the dynamic boundary
condition can be written as (2),P∂Φ 1+ |∇Φ|2 + gz = −
∂t
2
ρ

(2)

October 13, 2008

9:18

WSPC/101-CEJ

00187

Interaction Between Nearshore Nonlinear Waves and Strong Currents371

z = η(x, y, t)zxy∇2 Φ = 0z = −h(x, y)

Fig. 1. Physical domain of the nearshore ﬂow ﬁeld.

where g is the acceleration of gravity, P is the pressure at the free surface, and ρ is
the water density. Assume the free surface can be represented by z = η(x, t) whereη is a continuous and single-valued function of x, that means that the wave is notoverturned. The kinematic boundary condition is obtained as (3),ηt + ∇x η · ∇x Φ − Φz = 0 ,z = η(x, t)(3)

where ∇x ≡ (∂/∂x, ∂/∂y) denotes the horizontal gradient.
In the context of High-Order Spectral(HOS) Method and the newly developed
high-order Boussinesq type Equations, the Zakharov [1968]’s expression of free surface boundary condition is favored. This re-arrangement of the original Eulerian
form of free surface boundary conditions (2) and (3) is ﬁrst derived by Zakharov[1968] when studying weakly-nonlinear slowly-modulated waves through the nonlinear Schrödinger equation [Zakharov, 1968].
As in Zakharov [1968], we deﬁne the potential on the free surface
Φs (x, t) ≡ Φ(x, η(x, t), t)(4)

Applying chain rules on (4):Φt (x, η, t) = Φst (x, t) − Φz (x, η, t)ηt

(5)

∇x Φ(x, η, t) = ∇x Φs (x, t) − Φz (x, η, t)∇x η(6)

The dynamic and kinematic free surface boundary conditions (2) and (3) can be
rewritten as:s ηt + ∇x η · ∇x Φ − (1 + ∇x η · ∇x η)Φz (x, η, t) = 0
(7)
1
1
P

 Φst + gη + ∇x Φs · ∇x Φs − (1 + ∇x η · ∇x η)Φ2z (x, η, t) = −22
ρAn apparent advantage of using (7) instead of its orginal Eulerian form as freesurface boundary conditions is that given the initial values of η and Φs , (7) can beOctober 13, 2008

372

9:18

WSPC/101-CEJ

00187

Y. Tian & S. Sato

directly integrated in time for the new values of η and Φs as long as the surface vertical velocity Φz (x, η, t) is known. While the surface vertical velocity Φz (x, η, t) canbe obtained by solving the Laplace equation (1) with the given Φs as the Dirichletboundary condition. At the bottom, the normal velocity of the ﬂow is zero,
∂Φ
=0
(8)
∂n
where n is the normal direction of the boundary. If we deﬁne z = −h(x) to describe
the topographic variation, (8) can be written as:Φ z + ∇x h · ∇x Φ = 0 ,z = −h(x)(9)

3. Galerkin Chebyshev Spectral Approach for the Laplace Equation
3.1. Choice of base functions
The basic idea for solving diﬀerential equations through the spectral approach is to
assume that the unknown function can be approximated by a linear combination ofseveral “base functions”. With strategies to minimize the residual function, all thecoeﬃcients for the “base functions” are settled. It is important to choose a proper
set of base functions.
An appropriate set of base function should be eﬃcient, meaning that the results
converge rapidly as the increase of the number of base functions. Here the hyperbolic
functions cannot be considered to be a good candidate. Indeed the hyperbolic functions are the exact solutions of the Laplace equation, however, the appearance ofthe local wave number k weakens the eﬃciency of approximation. The local wave
number is a function of space and period, and cannot be pre-determined. Therefore
two schemes are adopted in literature. The ﬁrst is using the dispersion relation ofsmall amplitude wave theory directly, as in mild slope equation and Nadaoka [1994]’s
attempt to extend it to the nonlinear cases; the second is the wide band simulationin HOS, which commonly utilizes 1,000 hyperbolic functions with wave number kfrom 2π/L to 1,000 × 2π/L, where L is the length of the computational domain.
The former is not able to represent the phenomena related to nonlinear shoaling
and dispersion, and the latter obviously costs too many computational resources.
Well computational behavior is another requirement for obtaining creditable results. This requirement cannot be satisﬁed by the power functions. As increasing thenumber of base functions (the highest order of the power function), the condition
number of the computational matrices will increase rapidly. In function approximation, the best square polynomial approximation of a given function with powerfunctions to be the base functions generates the notorious Hilbert matrices, which
are canonical examples of ill-conditioned matrices. Large computational error may
be incurred by numercial round-oﬀ.
The shape of convergence domain is disadvantageous for the utilization of powerfunctions as well. For power series, the domain shape is a disk, which means that

October 13, 2008

9:18

WSPC/101-CEJ

00187

Interaction Between Nearshore Nonlinear Waves and Strong Currents

373

a singularity at the imaginary axis will make the model blow up. In fact, it really
happens when deriving formula of cnoidal waves with power series expansion.
In recent development of the spectral methods for partial diﬀerential equations,
Chebyshev orthogonal polynomial functions are frequently adopted. This set of functions is easy to compute, orthogonal, converges rapidly, and has a large convergencedomain.
The Chebyshev polynomials of the ﬁrst kinda are a set of orthogonal polynomials
denoted as Tn (x). The simplest deﬁnition of Chebyshev polynomials is the following
recurrence relation:
T0 (x) = 1
T1 (x) = x ,x ∈ [−1, 1](10)

Tn+1 (x) = 2xTn (x) − Tn−1 (x)The trigonometric deﬁnition (11) is convenient for the calculation with innerproduct.
Tn (cos(θ)) = cos(nθ)
The inner product along with Chebyshev polynomials is deﬁned as:
 1
dx
fg√
(f, g) ≡1 − x2
−1

(11)

(12)

Therefore for a given function f (x) = ai Ti (x), the coeﬃcient ai can be expressed as:
ai =
where


ci =

2
(f, Ti )
πci
2,

i=0

1,

i≤1

(13)

(14)

For the convenience of description, we deﬁneTi , f  ≡ f, Ti  ≡2(f, Ti )πci(15)3.2. Regulation of the computational domain
In order to expand the velocity potential to the Chebyshev series with respect to thevertical coordinate, we map the vertical computational domain z ∈ [−h(x), η(x)]into s ∈ [−1, 1]. Therefore the uneven computational domain is mapped to arectangle.
a Chebyshev polynomials mentioned in this paper are all abbreviated forms for Chebyshev polynomials of
the first kind.

October 13, 2008

374

9:18

WSPC/101-CEJ

00187

Y. Tian & S. Sato

Denoting
h+ (x) = h(x) + η(x)(16)

h− (x) = h(x) − η(x)(17)

2z + h− (x)h+ (x)(18)

s(x) =The velocity potential Φ can be written assh+ (x) − h− (x)= ϕ(x, s)Φ(x, z) = Φ x,
2

(19)

Using chain rules on the Laplace equation (1), we obtainh+ (x)2 ∇2 ϕ(x, s) = mijk si∂ j+k
ϕ(x, s) = 0 ,∂xj ∂ski = 0, 1, 2; j + k ≤ 2

(20)

where
m020 = (h+ )2(21)

m111 = −2h+ ∇x h+ ·(22)m011 = 2h+ ∇x h− ·(23)m101 = 2(∇x h+ )2 − h+ ∇2x h+(24)m001 = −2∇x h− · ∇x h+ + h+ ∇2x h−(25)m202 = (∇x h+ )2(26)m102 = −2∇x h− · ∇x h+(27)

m002 = (∇x h− )2 + 4(28)

The mijk not mentioned from (21) to (28) is equal to 0.
3.3. Chebyshev series expansion for regulated Laplace equation
The target unknown function ϕ(x, s) can be approximated as ϕN (x, s)N

ϕN (x, s) =an (x)Tn (s)

(29)

n=0

Denoting the diﬀerential operator in (20) as
Lijk = si∂ j+k∂xj ∂sk(30)

October 13, 2008

9:18

WSPC/101-CEJ

00187

Interaction Between Nearshore Nonlinear Waves and Strong Currents

375

we obtainN

Lijk ϕN (x, s) =

k
∂j
i ∂
a
(x)s
Tn (s)
n
∂xj
∂sk
n=0

(31)

Therefore the coeﬃcients of the residual function of (20) can be written as:N

rp = Tp (s), mijk Lijk ϕN (x, s) =

mijk
n=0

∂j
∂k
an (x) Tp (s), si k Tn (s)j∂x∂sp = 0, 1, . . . , N(32)

where the last pair of parentheses is a 4-dimensional tensor denoted as
Bpikn =Tp (s), si∂kTn (s)∂sk(33)

The values of Bpikn can be determined analytically.
3.4. Boundary conditions
3.4.1. Free surface boundary conditions
For the free surface, Φs of the present time is already given. Thus, the Dirichlet typeboundary condition can be written as:
NN

an (x) = Φs (x)(34)

dTn ∇x h− − s∇x h+
ds
h+

(35)

an (x)Tn (s)|s=1 =
n=0

n=0

3.4.2. Bottom boundary condition
Applying chain rules on (29), we obtain∇x ϕN = ∇x an Tn + an

Substituting (35) into (9) and evaluating at s = −1, we obtain
2an n2 (−1)n−1(∇x h)2 + 1+ (−1)n ∇x an · ∇h = 0
h+

(36)

3.4.3. Incident lateral boundary
The incident waves are formulated through Stream Function Wave theory [Dean,
1965]. Therefore, both the velocity uI and velocity potential ΦI of the incident wave
are represented in the following form:MM

Xi (x, t) cosh ik(z + h) =i=1

Xi (x, t) cosh
i=1

ik(s + 1)h+2(37)

October 13, 2008

376

9:18

WSPC/101-CEJ

00187

Y. Tian & S. Sato

Chebyshev series expansion of (37) arouses the following integration: 1
πIn (ikh+ /2) cosh[ikh+ /2] , n = 0, 2, 4, . . .
cosh[ik(s + 1)h+ /2]√
Tn (s)ds =πIn (ikh+ /2) sinh[ikh+ /2] , n = 1, 3, 5, . . .
1 − s2−1(38)
where In is the nth order modiﬁed Bessel function of the ﬁrst kind.
3.4.4. Far-field conditions
For an impermeable, reﬂective vertical wall, the far-ﬁeld condition is simple as:∇x Φ = 0
Applying Galerkin approach to (35), we obtain


N
∇x h+
∇x h−
an
Bp01n −
Bp11n = 0 ,∇x ap +
h+
h+

(39)

p = 0, 1, . . . , N

(40)

n=0

The open boundary condition for the spectral approach described in this study has
been treated with great care. Although there are many scheme (e.g. [Larsen andDancy, 1983]) to make an open boundary, most of them are designed for equations using the original variables u, η. For the approach using the velocity potential,
those scheme cannot make satisfactory results. Ohyama [1991] proposed a numericalwave absorber for a Mixed Euler-Lagrangian (MEL) method which uses the velocity potential as a variable. Figure 2 shows a schematic diagram of the numericalwave-absorption ﬁlter. An artiﬁcial “sponge layer” with gradually increasing damping coeﬃcient is attached to the computational domain, where µ is the damping
coeﬃcient in the “sponge layer” and xsp0 , xsp1 are the location of the beginning andend of the “sponge layer”. Due to the ineﬃciency of the “sponge layer” to absorb
long wave, behind the “sponge layer”, a radiation boundary condition is set in order
to allow long waves propagating out of the “sponge layer”. Following Ohyama’s idea,
we obtain the free surface boundary condition over the sponge layer:ηt + ∇x η · ∇x Φs − (1 + ∇x η · ∇x η)Φz (x, η, t) = 0






1
1
 s
Φt + gη + ∇x Φs · ∇x Φs − (1 + ∇x η · ∇x η)Φ2z (x, η, t)
(41)
2
2

 xsp1PΦs ∇x µ · dx = −+ µΦs −
ρ
xsp0
Behind the “sponge layer”, the wave is assumed as long wave. The radiative
boundary condition is written as√
η gh, x = xsp1
(42)
∇x a0 =h√where gh is the wave celerity of long wave.

October 13, 2008

9:18

WSPC/101-CEJ

00187

Interaction Between Nearshore Nonlinear Waves and Strong Currents












µmax

µ
















































































Sponge
Layer














































xsp0






Computational
Domain

h



377



xsp1

Fig. 2. Schematic diagram of the numerical wave-absorption ﬁlter.

3.5. Summary of governing equations
For the convenience of description, the equations for the 1DH case is summarized
here. The whole time-dependent problem after the Galerkin spectral approach is
written as follows:
Evolution equations for the free water surfaceηt + ∇x η · ∇x Φs − (1 + ∇x η · ∇x η)Φz (x, η, t) = 0







1
1
P


Φst + gη + ∇x Φs · ∇x Φs − (1 + ∇x η · ∇x η)Φ2z (x, η, t) = − , x < xsp0



2
2
ρ
(43)
1
1


Φst + gη + ∇x Φs · ∇x Φs − (1 + ∇x η · ∇x η)Φ2z (x, η, t)


2
2



 xsp1PΦs ∇x µ · dx = − , x > xsp0
 + µΦs −ρ
xsp0
Continuity equation in the computational domain
N

mijk
n=0

∂j
an (x)Bpikn = 0 ,∂xjp = 0, 1, . . . , N − 2(44)Dirichlet boundary condition for free surfaceNan (x) = Φs (x)

(45)

n=0

Impermeable boundary condition for the bottom
2an n2 (−1)n−1(∇x h)2 + 1+ (−1)n ∇x an · ∇h = 0
h+

(46)

Lateral boundary condition for incident waves
Dirichlet
an = Tn , ϕi  ,n = 0, 1, . . . , N(47)

October 13, 2008

378

9:18

WSPC/101-CEJ

00187

Y. Tian & S. Sato

Neumann boundary condition
N∇x ap +

an
n=0

∇x h−
∇x h+
B
−
Bp11n
p01n
h+
h+


= Tp , ui  ,p = 0, 1, . . . , N(48)

Far field condition
Impermeable vertical wallN

∇x ap +

an
n=0

∇x h+
∇x h−
B
−
Bp11n
p01n
h+
h+


= 0,

x0 = +∞

p = 0, 1, . . . , N(49)Open boundary√
η gh∇x a0 =h(50)

Up to now, the main part of the theoretical derivation is completed. The following
section will discuss the numerical implementation of this model.

4. Numerical Implementation
Through Galerkin approach, the vertical coordinate in the governing equations has
been eliminated. All the diﬀerential operators are applied to the horizontal coordinates. It is then possible to obtain a solvable system with spectral or presudo-spectral
methods on the two horizontal coordinates. However, the beneﬁts of applying spectral methods to the horizontal coordinates is not so signiﬁcant as to the verticalcoordinate. The 10th order Chebyshev polynomials on the vertical coordinate will
satisfy most of the requirements, while the 1,000th order may be necessary for thehorizontal coordinates. Surely a full spectral solver will greatly enhance the accuracy of the results and eﬃciently save the numbers of unknowns, yet the full spectral
approach always induces full matrices, which cost much memory to store and computational time to factorize. Therefore, it is diﬃcult to evaluate the loss and gain.
Although with shortage of accuracy and eﬃciency, the ﬁnite diﬀerence scheme
has several merits. It is easy to program, ﬂexible to carry out local procedures(e.g. to multiply damping coeﬃcient on a certain term) and generates band diagonal
matrices. The 1-dimensional solver is also convenient to extend to 2-dimensional
problem, which is a really tedious work for the full spectral solver.
After weighing those factors mentioned above, in an attempt of this study, we
adopt a semi-spectral solver, which utilizes the ﬁnite diﬀerence scheme along thehorizontal 1 dimension cases. A schematic diagram for solving (43) to (50) is shownin Fig. 3.

October 13, 2008

9:18

WSPC/101-CEJ

00187

Interaction Between Nearshore Nonlinear Waves and Strong Currents

379

Solving(44)∼ (50)##s , η, tΦ:

Evaluation of Bijkn

Φ, tRenewing Data}Φs , η, t +c ∆tΦz z=η , t(43)

Temporal Integration

Fig. 3. Schematic diagram for sloving (43) to (50).

4.1. Evaluation of Bijkn
As shown in Fig. 3, the tensor Bijkn is evaluated before the main loop of temporal
integration. All the evaluations are based on three fundamental values:Bi00n = Ti , Tn  =

0,

i=n

1,

i=n


2
Bi01n = Ti , dTn /ds =ci
1
Bi02n = Ti , d Tn /ds  =
ci
2

2

(51)

n,

i ∈ {n − 1, n − 3, n − 5, . . .}0,otherwisen(n2 − i2 ) ,i ∈ {n − 2, n − 4, n − 6, . . .}0,otherwise

(52)

(53)

Thus, dTn /ds and d2 Tn /ds2 can be expressed as

dTn
=
ds
d2 Tn
=
ds2

n−1

Bi01n Ti

(54)

Bi02n Ti

(55)

i=0
n−2
i=0

From the recursive deﬁnition of Chebyshev Polynomials (10), we obtain Tk−1 + Tk+1,
sTk (s) =
2
T
1

k≥1
k=0

(56)

October 13, 2008

380

9:18

WSPC/101-CEJ

00187

Y. Tian & S. Sato

and



 Tk−2 + 2Tk + Tk+2 ,


4




3T1 + T3
s2 Tk (s) =
,

4






 T0 + T2 ,
2

k≥2
k=1

(57)

k=0

Substituting (56) and (57) into (54) and (55), and applying (51)–(53) we can evaluate
all the other Bijkn’s.
Take Bi11n as an example:
n−1
 Tr−1 + Tr+1
dTn
, r≥1
(58)
=
Br01n
s
2
T ,
ds
r=0
r=0

1

Therefore
Bi11n =

Ti , sdTnds Tr−1 + Tr+1 , r ≥ 1 \
2
i
r01n

r=0/
T1 ,
\ r=0

n−1
 Ti , Tr−1  + Ti , Tr+1  , r ≥ 1
2
Br01n
=

r=0
r=0
Ti , T1  ,
n−1 Bi00(r−1) + Bi00(r+1) , r ≥ 1
2
Br01n
=

r=0
r=0
Bi001 ,
/ n−1

B
= T,

(59)

By the similar virtue, all the Bijkn are able to be evaluated.
4.2. Spatial discretization
By knowing the value of Bijkn , spatial discretization can be carried out. The unknowns are the coeﬃcients of the Chebyshev polynomials, an (x). A non-staggeredmesh is adopted.
Equation (44) can be written as∂ 2 an
+
∂x2

N
m=0

∂am
+
C(x)nm∂x

ND(x)nm am = 0 ,m=0where C(x) and D(x) are coeﬃcient matrices.

n = 1, 2, . . . , N − 2(60)October 13, 2008

9:18

WSPC/101-CEJ

00187

Interaction Between Nearshore Nonlinear Waves and Strong Currents

381

For the ﬁrst term in (44), the 2nd order central diﬀerence scheme is applied: 2 i
i
i−1
∂ an
ai+1
n − 2an + an=
, O(∆x2 )
(61)
∂x2
∆x2
For the second term in (44), the diﬀerence scheme changes with the m subscript.
If m = n, we adopt the second order upwind diﬀerence scheme:

1
i
i−1
i−2

i 

O(∆x2 ) , Cnm < 0 2∆x (3am − 4am + am ) ,
∂am
(62)
=

∂x

i+2
2
 1 (−3aim + 4ai+1O(∆x ) , Cnm > 0
m − am ) ,2∆x
If m = n, the second order central diﬀerence scheme is applied:


1
∂am i
(ai+1 − ai−1=
O(∆x2 )
(63)
m ),
∂x
2∆x m
4.3. Incident boundary conditions
In the ﬁrst few periods, a moderation coeﬃcient α(t) is multiplied to the incident
conditions.


t
(64)
α = 1 − exp −βTwhere β is usually set to 1. Thus a tentative incident boundary condition is adopted:
ΦstI = α(t)ΦI

(65)

ust
I = α(t)uI(66)

ηIst = α(t)ηI(67)

4.3.1. Dirichlet boundary condition
For the Dirichlet boundary condition, another factor has to be considered. From
(7) and (1) we can ﬁnd that those equations still hold after adding a constant Con Φ. That means there may be a spatial constant diﬀerence between ΦI and Φcal
due to diﬀerent datums (The datum for the velocity potential from stream function
theory is a moving point x − ct = 0, while the datum for the velocity potential in
the computational domain is a ﬁxed point located at x = +∞). Fortunately the 0thorder Chebyshev polynomials T0 is a constant value 1, meaning that the verticalderivative Φz will not be inﬂuenced by the coeﬃcent before T0 . Therefore we use
the value (T1 , Φi ), (T2 , Φi ), . . . , only leaving the T0 for a derivative-type relation inorder to eliminate the inconsistency due to the diﬀerence between datums.
−1.5a00 − 2a10 − 0.5a20 = T0 , u(x0 )I 
−1/6a00 + a10 − 0.5a20 − 1/3a30 = T0 , u(x1 )I 

(68)

October 13, 2008

382

9:18

WSPC/101-CEJ

00187

Y. Tian & S. Satoa0n = Tn , Φ(x0 )I  ,n = 1, 2, . . . , N

a1n = Tn , Φ(x1 )I  ,n = 1, 2, . . . , N(69)

4.3.2. Neumann boundary condition
The Neumann boundary condition is simple as follows:N−1.5a0p − 2a1p − 0.5a2p +

n=0

a0n

∇x h+
∇x h−
B
−
Bp11n
p01n
h+
h+


= Tp , uI  ,p = 0, 1, . . . , NN−1/6a0p + a1p − 0.5a2p − 1/3a3p

n=0

a0n

∇x h+
∇x h−
B
−
Bp11n
p01n
h+
h+

(70)


= Tp , uI  ,p = 0, 1, . . . , N

Impermeable reflective boundary condition
From (49), we simply obtain−1 − 2aN −2 + 0.5aN −3 +
1.5aN
p
p
p

N
n=0

−1
aN
n




∇x h+
∇x h−
Bp01n − + Bp11n = 0 ,h+hp = 0, 1, . . . , N(71)Open boundary condition
−1
−2
−3
− 2aN
+ 0.5aN
=
1.5aN
0
0
0

√
η gh
h

(72)

4.4. Analysis of the linear algebra system
The linear algebra system generated from the above discretization approach is a
blocked sparse matrix. The distribution of the non-zero element is shown in Fig. 4.
We can ﬁnd that the matrix consists of several (N + 1) × (N + 1) block. Exactly
speaking, it is a band diagonal matrix with band width 2 × (N + 1). The LAPACK
routines for solving band diagonal linear algebra system is adopted. Denoting thenumber of total points is Np , the routine for LU factorization in LAPACK named
DGBTRF, carries out approximately 2Np (2N + 1)(2N + 2) ∼ 2Np (2N + 1)(4N +3) times of ﬂoating-point operations, depending on the interchanges. Routine forforward and backward substitution to solve the equation named DGBTRS procedes
2Np (6N + 3) times ﬂoating-point operations. Therefore, the temporal compexity
for the spatial solving approach can be estimated as O(Np N 2 ). Compared with full

October 13, 2008

9:18

WSPC/101-CEJ

00187

Interaction Between Nearshore Nonlinear Waves and Strong Currents

383

N +1

_ _ _  _ o _ _/ O ••••••◦◦◦◦•◦◦◦◦• • • • •  ◦ • ◦ ◦ ◦ ◦ • ◦ ◦ ◦N +1 • • • • • ◦ ◦ • ◦ ◦ ◦ ◦ • ◦ ◦• • • •  ◦ ◦ ◦ • ◦ ◦ ◦ ◦ • ◦•• • • • _◦ ◦ _◦ ◦ _• ◦ ◦ ◦ ◦ •_•• ◦ _◦ ◦ _◦ • • • • • _• ◦ _◦ ◦ _◦ _• ◦ _◦ ◦ _◦ ◦ • ◦ ◦ ◦ • • • • • ◦ • ◦ ◦ ◦ ◦ • ◦ ◦ ◦◦◦•◦◦•••••◦◦•◦◦◦◦•◦◦ ◦ ◦ ◦ • ◦ • • • • • ◦ ◦ ◦ • ◦ ◦ ◦ ◦ • ◦◦ ◦ ◦ • • • • • • _◦ ◦ _◦ ◦ _• ◦ ◦ ◦ ◦ •_◦• ◦ _◦ ◦ _◦ _• • _• • _• • • • • • _• • _• • _• _• ◦ _◦ ◦ _◦ ◦ • ◦ ◦ ◦  • • • • • • • • • • • • • • • ◦ • ◦ ◦ ◦◦◦•◦◦•••••••••••••••◦◦•◦◦ ◦ ◦ ◦ • ◦  • • • • • • • • • • • • • • • ◦ ◦ ◦ • ◦_◦ ◦ _◦ ◦ _• _• • _• • _• _• • • • • _• • _• • _• ◦ ◦ ◦ ◦ •• ◦ ◦ ◦ ◦  • • _• • _• • • • • • _• • _• • _• _• ◦ _◦ ◦ _◦ ◦ • ◦ ◦ ◦  • • • • • • • • • • • • • • • ◦ • ◦ ◦ ◦◦◦•◦◦•••••••••••••••◦◦•◦◦ ◦ ◦ ◦ • ◦  • • • • • • • • • • • • • • • ◦ ◦ ◦ • ◦_◦ ◦ _◦ ◦ _• _• • _• • _• • • • • • _• • _• • _• _◦ ◦ _◦ ◦ _••◦◦◦◦••••••••••••••••◦◦◦◦◦•◦◦◦•••••••••••••••◦•◦◦◦◦◦•◦◦•••••••••••••••◦◦•◦◦◦◦◦•◦•••••••••••••••◦◦◦•◦◦◦◦◦••••••••••••••••◦◦◦◦•···•◦◦◦◦••••••••••••••••◦◦◦◦◦•◦◦◦•••••••••••••••◦•◦◦◦◦◦•◦◦•••••••••••••••◦◦•◦◦◦◦◦•◦•••••••••••••••◦◦◦•◦◦◦◦◦••••••••••••••••◦◦◦◦••◦◦◦◦•◦◦◦◦••••••◦◦◦◦◦•◦◦◦◦•◦◦◦•••••◦•◦◦◦◦◦•◦◦◦◦•◦◦•••••◦◦•◦◦◦◦◦•◦◦◦◦•◦•••••◦◦◦•◦◦◦◦◦•◦◦◦◦••••••◦◦◦◦••◦◦◦◦•◦◦◦◦•••••◦•◦◦◦◦•◦◦◦•••••◦◦•◦◦◦◦•◦◦•••••◦◦◦•◦◦◦◦•◦•••••◦◦◦◦•◦◦◦◦••••••
Fig. 4. Distribution of the non-zero element of the LAS for 1DH.

matrix which arouses a temporal compexity of O(Np3 ), the present scheme saves thecomputational time signiﬁcantly at least for 1DH cases.
It is necessary to point out that although a compact band matrix can be obtained
for the 1DH case, the matrix elements must be scattered for the 2DH, which may
reduce computational eﬃciency. The Alternative-Direction-Implicit(ADI) scheme isexpected to be applied to the 2DH cases in order to enhance the computationaleﬃciency.
4.5. Temporal integration and calculation of vertical derivatives
The temporal integration of (43) is carried out through the 4th order Runge-Kutta
method without too much diﬃculty.
While a key value in (43) is Φz , which is also easy to be calculated as follows.
2
ds
= +
dz
h

(73)

Applying the chain rule, we obtainN2
dTn (s)anΦz = +
h n=0
ds
N

=

2
an n2
h+ n=1(74)

October 13, 2008

384

9:18

WSPC/101-CEJ

00187

Y. Tian & S. Sato

5. Model Verification
In order to verify the model, we use the model to simulate three phenomena: thepropagation of ﬁnite amplitude waves over an even bottom with diﬀerent relative
water depth and wave number; nonlinear wave shoaling; harmonic generation forwaves propagating over a submerged breakwater. The ﬁrst two cases have theoretical solutions which can be used as the criteria and in the third case the resultof simulation will be compared with laboratory data, and with the results fromBoussinesq-type models. The 12th order Chebyshev polynomial is adopted in all
the simulations.
5.1. Waves on horizontal bottom
For the waves on a horizontal bottom, two cases of incident waves have been selected.
The ﬁrst case is the test of short wave and the second case is the test of nonlinear
wave in the intermediate water depth. The characteristic parameters of the incident
waves are listed in Table 1.
Table 1. Characteristic parameters of theincident waves.

6

x 10

T (s)

h(m)

H(m)

∆t(s)

case 1

0.6

1.0

0.001

0.01

case 2

5.0

1.0

0.3

0.1

−4

Simulated
10th SFM

4

η (m)

2

0

−2

−4

−6

0

1

2

3

4

5
x/L

6

7

8

9

10

Fig. 5. Spatial distribution for calculated surface elevation, kh = 11.

October 13, 2008

9:18

WSPC/101-CEJ

00187

Interaction Between Nearshore Nonlinear Waves and Strong Currents

385

0.25
Simulated
10th SFM
0.2

0.15

η (m)

0.1

0.05

0

−0.05

−0.1

0

1

2

3
x/L

4

5

6

Fig. 6. Spatial distribution for calculated surface elevation, H/h = 0.3.
Table 2. Computational conditions for diﬀerent cases of shoaling study.
Case

H0 /L0h1 [m]

h2 [m]KS1

H1 /2[m]Slope Length[m]

Slope

1

0.001

1.0

0.5

1.6

0.08

50

0.01

2

0.005

6.0

1.37

1.00710.2518463

0.01

Figures 5 and 6 demonstrate the spatial distribution of the water surface elevation for cases 1 and 2, comparing with the analytical results from the 10th order
Stream Function wave theory. In both ﬁgures, the diﬀerence between the numericalsimulation and the theory cannot be easily recognized. Here it is pointed out thatfor the short wave in case 1, the results from most of the Boussinesq type equations
will have signiﬁcant phase shift referring to the theoretical solution.
5.2. Nonlinear wave shoaling
In order to examine the ability of this model to predict wave shoaling over slowly
varying depth, two incident waves and their corresponding shoaling conditions arelisted in Table 2, where h1 is the oﬀshore depth, h2 the onshore depth, H1 theoﬀshore wave height calculated from the deepwater wave height H0 , using the linear
wave theory shoaling coeﬃcient Ks . Figure 7 shows a sketch for the computational
domain used in this study.
The shoaling coeﬃcents are calculated, and compared with the results of Stiassnie and Peregrine [1980] and Shuto [1974] in Figs. 8 and 9. It is clear from these

October 13, 2008

386

9:18

WSPC/101-CEJ

00187

Y. Tian & S. Sato

Slope length
H1

h2
h1

Fig. 7. Spatial distribution for calculated surface elevation, wave height, and bottom depth.
2.8
Simulation
Shuto (1974)
Stiassnie and Peregrine (1980)

2.6
H 0 /L0 = 0.001
2.4

H/H 0

2.2
2
1.8
1.6
1.4
0.005

0.006

0.007
h/L0

0.008

0.009

0.01

Fig. 8. Shoaling coeﬃcient obtained by the present model, Stiassnie and Peregrine [1980] and Shuto[1974].
2.6
Present model
Shuto (1974)
Stiassnie and Peregrine (1980)

2.4
2.2

H 0 / L 0 = 0.005

H/H0

2
1.8
1.6
1.4
1.2
1
0.8
10

−1.9

10

−1.8

10

−1.7

10

−1.6

−1.5

10
h / L0

10

−1.4

10

−1.3

10

−1.2

10

−1.1

Fig. 9. Shoaling coeﬃcient obtained by the present model, Stiassnie and Peregrine [1980] and Shuto[1974].

October 13, 2008

9:18

WSPC/101-CEJ

00187

Interaction Between Nearshore Nonlinear Waves and Strong Currents

387

0Depth (m)

−0.1
−0.2

1/10

1/20

−0.3
−0.4
−0.5

0

5

10

15

20
Distance (m)

25

30

35

40

Fig. 10. Bathymetry of Luth’s experiment.

ﬁgures that the model developed through Chebyshev Spectral approach has a good
capability to simulate water wave shoaling over slowly varying topography.
5.3. Harmonic generation for waves propagating over a submerged
breakwater
Wave transformation over submerged bars in shallow water is an intriguing problem,
which involves the phenomena of wave-wave interaction and harmonic generation.
Here we simulate the experiment by Luth [1994]. The experiments were equivalent
to the set up of Beji and Battjes [1993] except for a factor 2 scaling. In the reduced
scale the bathymetry is deﬁned as shown in Fig. 10, a horizontal ﬂume with a depthof 40 cm, a trapezoidal bar with an upward slope of 1/20, a downward slope of 1/10,
and a depth of 10 cm on top of the bar. Active absorption was applied at both endsof the ﬂume. Their test case A is used for the veriﬁcation in the following. Case A
is based on regular incident waves with a wave period of 2.02 s and a wave heightof 2.0 cm.
Figures 11 to 13 show the the simulated time series of surface elevations at threelocations, together with the simulation by Madsen [1999] for the same condition
using Nwogu [1993] and Wei and Kirby [1995]’s equation. From Figs. 11 to 13, theagreement between the simulation and measurement is fair. The slight overshooting of the wave height for the ﬁrst order harmonic wave is related to the lack ofdissipation of the model, while either Wei [1995]’s or Nwogu [1993]’s model fails tosimulate the harmonic generation at the location where x = 19 m and x = 21 m.
The main reason for the success of the model in this study is believed to be due to its
fully-nonlinear characteristics and high accuracy to describe shoaling and dispersion.

6. Numerical Study of Nonlinear Wave/Current Interaction
After the model veriﬁcation via diﬀerent computational conditions, the reliability
and accuracy of the model have been conﬁrmed. The direct application of the modelis the numerical study of nonlinear wave-current interaction. This study focuses inparticular on the numerical simulation for wave blocking phenomenon.

October 13, 2008

388

9:18

WSPC/101-CEJ

00187

Y. Tian & S. Sato
0.03
Present model
Luth(1994)

0.025
0.02
0.015

η (m)

0.01
0.005
0
−0.005−0.01−0.015
−0.02
136

136.5

137

137.5

138
Time (s)

138.5

139

139.5

140

Fig. 11. Harmonic generation over a submerged bar, x = 14.5 m.

0.02
0.015
0.01

η (m)

0.005
0
−0.005−0.01−0.015Present model
Luth(1994)
Nwogu(1993)
Wei(1995)−0.02−0.025
136

136.5

137

137.5

138
Time (s)

138.5

139

139.5

140

Fig. 12. Harmonic generation over a submerged bar, x = 19 m.

6.1. Simulation for wave blocking phenomenon
Here it is necessary to eluciate the condition of the wave blocking phenomenion
again: when propagating waves encounter an opposing current of increasing magnitude, wave blocking phenomenon occurs at a place where the current speed matches
the group velocity of the wave. In the blocking region, the wave becomes shorterand steeper, and eventually breaks as it is reﬂected and travels back.
In the small amplitude wave theory, an explanation from the perpsective ofenergy is that: when a wave train encounters a current in which the surface velocity
varies, the excess momentum ﬂux results in an interchange of energy between waves
and current. When the convection velocity of energy is equal and opposite to the
local group velocity of the waves, the energy can no longer be propagated against
the stream and accumulated before the breaking point. The accumulation of theOctober 13, 2008

9:18

WSPC/101-CEJ

00187

Interaction Between Nearshore Nonlinear Waves and Strong Currents

389

0.03
Present model
Luth(1994)
Nwogu(1993)
Wei(1995)

0.025
0.02
0.015

η (m)

0.01
0.005
0
−0.005−0.01−0.015
−0.02
136

136.5

137

137.5

138
Time (s)

138.5

139

139.5

140

Fig. 13. Harmonic generation over a submerged bar, x = 21 m.

energy leads to inﬁnite wave amplitude which is unrealistic and appears in the formof wave breaking (e.g. [Longuet-Higgins and Stewart, 1961]; [Phillips, 1977]).
Analytical and semi-analytical studies of the nonlinear wave blocking phenomenon point out that the breaking of wave before the blocking point is in theform of reﬂection and backward travel of waves with high wave number (e.g. [Smith,1975]; [Shyu and Phillips, 1990]; [Shyu and Tung, 1999]). However, theoretical treatment for wave blocking by Peregrine and Thomas [1979] does not include the eﬀect
of wave reﬂection and fails in the blocking region where the assumption of slowly
varying waves breaks down.
Numerical approaches for the wave blocking phenomenon have been carried out
by Chen [1996] through a set of modiﬁed Boussinesq equations. A current-wave
seperation skill was applied to Nwogu [1993]’s Boussinesq equation with length scaleparameters. However, this approach requires a very strong bottom friction in orderto stablize simulation. The strong dissipation introduced by the bottom friction
remove all the key structures near the blocking point.

6.2. Computational conditions
In this study, similar computational conditions as in Chen [1996] are adopted. The
spatial variation of current is generated by spatial change of topography. This caseconsiders a wave train propagating against a current in a channel with a subermged
bar. A sketch of the bathymetry is shown in Fig. 14.
The channel is 60 m long, 0.8 m deep on both sides of the bar and 0.2 m deep on
top of the bar. The bar is located from 8.0 m to 55.0 m. The slopes of the bar are
1/50 on the incoming side and 1/20 on the transmission side. The two corners of thebar near x = 40 m are rounded for the convenience to obtain the initial boundary
conditions.

October 13, 2008

390

9:18

WSPC/101-CEJ

00187

Y. Tian & S. Sato

Depth (m)

0
−0.4

1/20

1/50
−0.8
0

10

20

30
Distance (m)

40

50

60

Fig. 14. Bathymetry of the channel.
0.02

η (m)

0
−0.02
−0.04

0

10

20

30
Distance (m)

40

50

60

0

10

20

30
Distance (m)

40

50

60

Φs (m 2 /s)

10
0
−10
−20

Fig. 15. Initial spatial distribution of water elevation, and velocity potential.

The depth uniform steady current has a discharge of −0.136 m2 /s (from rightto left). As shown in Fig. 15, the steady solution of the spatial distribution of thecurrent and its corresponding velocity potential on the free surface are obtained
through a program using boundary element methods. On the bar crest the current
velocity is −0.8090 m/s, the setdown of surface elevation is −0.0319 m. According
to linear theory only waves with a period larger than 3.35 s can be expected to pass
the bar without being blocked by the current.
A sinusoidal wave train propagation from left to right on top of the steady current
has been imposed. The incoming wave has a period of 2.0 s and an initial heightof 0.002 m. According to linear theory this wave are expected to be blocked where√
the Froudes number Fc = U/ gd is −0.475, where d is the local water depth. Thislocation is found at x = 36.78 m, where the water depth of the initial steady current
from the free surface is 0.203 m.
6.3. Results of simulation
The computed surface elevations are shown in Figs. 16–22 for diﬀerent instants
during the simulation.

October 13, 2008

9:18

WSPC/101-CEJ

00187

Interaction Between Nearshore Nonlinear Waves and Strong Currents

391

K (m)

0.02
0
Ŧ0.02
Ŧ0.04

0
x 10

5

10

15

20
25
Dis tance (m) (a)

30

35

40

5

10

15

20
25
Dis tance (m) (b)

30

35

40

Ŧ3

K Ŧ K0 (m)

4
2
0
Ŧ2
Ŧ4

0

Fig. 16. Spatial distribution of surface elevation with and without setdown from current, t = 35 s.

K (m)

0.02
0
Ŧ0.02
Ŧ0.04

0
x 10

5

10

15

20
25
Dis tance (m) (a)

30

35

40

5

10

15

20
25
Dis tance (m) (b)

30

35

40

Ŧ3

K Ŧ K0 (m)

4
2
0
Ŧ2
Ŧ4

0

Fig. 17. Spatial distribution of surface elevation with and without setdown from current, t = 38 s.

In Fig. 16, the incoming wave train has just reached the area with strong currents
and there is yet no sign of blocking. In Fig. 17, the wave amplitude increases rapidly,
implying that the wave energy begins to accumulate. In Fig. 18, the reﬂected wave
has been generated and can be observed in a narrow zone.
From Figs. 19 to 22, the wave amplitude increases to 10 times of that of the incident wave, while the reﬂected waves are kept in traveling backward. For the reﬂected
waves, the eﬀects of the opposing “original wave” together with the accumulated
wave energy will make them shorter and steeper. In reality, those short wave will
be rapidly dissipated by surface tension and viscosity ([Trulsen and Mei, 1993]). In
the simulation, the steep and short reﬂected waves demand very ﬁne spatial andtemporal resolution, if no dissipation is applied. No dissipation is imposed in thepresent computation.

October 13, 2008

392

9:18

WSPC/101-CEJ

00187

Y. Tian & S. Sato

η (m)

0.02
0
−0.02
−0.04

0

5

10

15

20
25
Distance (m) (a)

30

35

40

0

5

10

15

20
25
Distance (m) (b)

30

35

40

η − η0 (m)

0.01
0.005
0
−0.005−0.01

Fig. 18. Spatial distribution of surface elevation with and without setdown from current, t = 42.5 s.

η (m)

0.02
0
−0.02
−0.04

0

5

10

15

20
25
Distance (m) (a)

30

35

40

0

5

10

15

20
25
Distance (m) (b)

30

35

40

η − η0 (m)

0.01
0.005
0
−0.005−0.01

Fig. 19. Spatial distribution of surface elevation with and without setdown from current, t = 50 s.

The results of simulation shown here indicate that the model not only accurately
predicts the location of the blocking point, but also gives a lucid description of the
wave structure in the blocking zone (the zone near the blocking point). It is believedthat the stability of the model for simulating the wave blocking phenomenon can be
improved by introducing a certain wave breaking model.

7. Conclusion
In this study, a new set of equations based on the spatial-temporal separation conception is established through the Galerkin Chebyshev spectral approach on thevertical direction. The 1DH version of this set of equation has been implemented

October 13, 2008

9:18

WSPC/101-CEJ

00187

Interaction Between Nearshore Nonlinear Waves and Strong Currents

393

0

η (m)

−0.01
−0.02
−0.03
−0.04
30

31

32

33

34
35
36
Distance (m) (a)

37

38

39

40

31

32

33

34
35
36
Distance (m) (b)

37

38

39

40

0

η − η (m)

0.01
0.005
0
−0.005−0.01
30

Fig. 20. Spatial distribution of surface elevation with and without setdown from current, t = 50 s.
0

η (m)

−0.01
−0.02
−0.03
−0.04
30

31

32

33

34
35
36
Distance (m) (a)

37

38

39

40

31

32

33

34
35
36
Distance (m) (b)

37

38

39

40

0

η − η (m)

0.01
0.005
0
−0.005−0.01
30

Fig. 21. Spatial distribution of surface elevation with and without setdown from current, t = 55 s.

and veriﬁed with theoretical and experimental data. The model is proved to be ofhigh accuracy with moderate cost of computational resources and time.
The model can describe most of the wave and current motions of the potential
ﬂow, before the wave overturn and/or breaking. All the nonlinear characteristics
have been retained in this model. The dispersion relation of the model is also veriﬁedby the simulation of harmonic generation of waves propagating over a submergedbar.
With this model, the wave blocking phenomenon in wave-current interaction has
been studied. The results of simulation indicate that the location of blocking point
determined by the linear wave theory is accurate for nonlinear phenomenon, while
the nonlinear wave structure in the blocking zone is also observed. Besides the above
success, several issues emerged during the development of the model.

October 13, 2008

394

9:18

WSPC/101-CEJ

00187

Y. Tian & S. Sato
0

η (m)

−0.01
−0.02
−0.03
−0.04
30

31

32

33

34
35
36
Distance (m) (a)

37

38

39

40

31

32

33

34
35
36
Distance (m) (b)

37

38

39

40

η − η0 (m)

0.01
0.005
0
−0.005−0.01
30

Fig. 22. Spatial distribution of surface elevation with and without setdown from current, t = 65 s.

(1) An incident wave boundary condition which is able to deal with the reﬂectedwave is neccesary in order to simulate the ﬂow pattern around coastal structures.
(2) 2DH version of the numerical model is expected in order to be applied to real
problems.
(3) A breaking wave model is needed for further simulation.
References
Beji, S. & Battjes, J. A. [1965] “Experimental investigation of wave-propagation over a bar,” Coastal
Engineering 19, 151–162.
Chen, Q., Madsen, P. A., Sørensen, O. R. & Basco, D. R. [1996] “Boussinesq equations with improved doppler shift and dispersion for wave/current interaction,” in Proc. 25th Int. Conf.
on Coastal Eng., ASCE, pp. 1060–2073.
Dean, R. G. [1965] “Stream function representation of nonlinear ocean waves,” J. Geophys. Res.
70(18), 4561–4572.
Isobe, M. [1994] “Time dependent mild slope equations for random waves,” in Proc. 24th Int. Conf.
on Coastal Eng., ASCE, pp. 305–325.
Larsen, J. & Dancy, H. [1983] “Open boundaries in short-wave simulations — A new approach,”
Coastal Engineering 7(3), 285–297.
Longuet-Higgins, M. S. & Stewart, R. W. [1961] “The changes in amplitude of short gravity waves
on steady non-uniform currents,” J. Fluid Mech. 10, 529–549.
Luth, H. R., Klopman, G. & Kitou, N. [1961] “Kinematics of waves breaking partially on an oﬀshorebar; LDV measurements for waves with and without a net onshore current,” Delft Hydraulics
Report, H1573, Delft Hydraulics, The Netherlands.
Madsen, P. A. [1999] “A review of boussinesq-type equations for surface gravity waves,” Advancesin Coastal and Ocean Engineering, Vol. 5, Chapter 1 (World Scientiﬁc), pp. 1–89.
Mei, C. C., Stiassnie, M. & Yue, D. K. P. [2005] Theory and Application of Ocean Surface Wave
Part 2: Nonliear Aspects (World Scientiﬁc).
Nadaoka, K., Beji, S. & Nakagawa, Y. [1994] “A fully dispersive nonlinear wave model and its
numerical solution,” in Proc. 24th Int. Conf. on Coastal Engg., ASCE, pp. 427–441.
Nwogu, O. [1993] “Alternative form of boussinesq equations for nearshore wave propagation,” J.
Waterways, Port, Coastal, Ocean Engrg., ASCE 119(6), 618–638.

October 13, 2008

9:18

WSPC/101-CEJ

00187

Interaction Between Nearshore Nonlinear Waves and Strong Currents

395

Ohyama, T. & Nadaoka, K. [1991] “Development of a numerical wave tank for analysis of nonlinear
and irregular wave ﬁeld,” Fluid Dynamics Res. 8(5–6) (November) 231–251.
Phillips, O. M. [1977] The Dynamics of the Upper Ocean (Cambridge Univ. Press).
Peregrine, D. H. [1967] “Long waves on a beach,” J. Fluid Mech. 27, 815–820.
Peregrine, D. H. & Thomas, G. P. [1979] “Finite-amplitude deep-water waves on currents,” Phi.
Trans. Roy. Soc. London 292, 371–390.
Shuto, N. [1974] “Nonlinear long waves in a channel of variable section,” Coastal Engineering inJapan 17, 1–12.
Shyu, J. H. & Phillips, O. M. [1990] “The blockage of gravity and capillary waves by longer waves
and currents,” J. Fluid Mech. 217, 115–141.
Shyu, J. H. & Tung, C. C. [1999] “Refraction of oblique waves by currents: analytical solutions andtheir application to numerical computations,” J. Fluid Mech. 396, 143–182.
Smith, R. [1975] “The reﬂection of short gravity waves on a non-uniform current,” Math. Proc.
Camb. Phil. Soc. 78, 517–525.
Stiassnie, M. & Peregrine, D. H. [1980] “Shoaling of ﬁnite-amplitude surface-waves on water ofslowly-varying depth,” J. Fluid Mech. 97, 783–805.
Trulsen, K. & Mei, C. C. [1993] “Double reﬂection of capillary/gravity waves by a non-uniform
current: A boundary-layer theory,” J. Fluid Mech. 251, 239–271.
Wei, G. & Kirby, J. T. [1995] “Time-dependent numerical code for extended boussinesq equations,”
J. Waterways, Port, Coastal, Ocean Engrg., ASCE 121, 251–261.
Zakharov, V. E. [1968] “Stability of periodic waves of ﬁnite amplitude on the surface of a deep
ﬂuid,” J. Appl. Mech. Tech. Phys. (Eng. transl.) 2, 190–194.