Ocean Modelling 34 (2010) 16–35

Contents lists available at ScienceDirect

Ocean Modelling
journal homepage: www.elsevier.com/locate/ocemod

Wave–current interaction in an oceanic circulation model with a vortex-force
formalism: Application to the surf zone
Yusuke Uchiyama *, James C. McWilliams, Alexander F. Shchepetkin
Institute of Geophysics and Planetary Physics, University of California, Los Angeles, CA, USAa r t i c l e

i n f oArticle history:
Received 2 October 2009
Received in revised form 12 February 2010
Accepted 12 April 2010
Available online 18 April 2010
Keywords:
Wave–current interaction
Vortex force
ROMS
Littoral currenta b s t r a c t
A vortex-force formalism for the interaction of surface gravity waves and currents is implemented in athree-dimensional (3D), terrain-following, hydrostatic, oceanic circulation model (Regional Oceanic Modeling System: ROMS; Shchepetkin and McWilliams, 2005). Eulerian wave-averaged current equations for
mass, momentum, and tracers are included in ROMS based on an asymptotic theory by McWilliams et al.
(2004) plus non-conservative wave effects due to wave breaking, associated surface roller waves, bottom
streaming, and wave-enhanced vertical mixing and bottom drag especially for coastal and nearshore applications. The currents are coupled with a spectrum-peak WKB wave-refraction model that includes theeffect of currents on waves, or, alternatively, a spectrum-resolving wave model (e.g., SWAN) is used. The
coupled system is applied to the nearshore surf zone during the DUCK94 ﬁeld measurement campaign.
Model results are compared to the observations and effects of parameter choices are investigated withemphasis on simulating and interpreting the vertical proﬁles for alongshore and cross-shore currents.
The model is further compared to another ROMS-based 3D coupled model by Warner et al. (2008) withdepth-dependent radiation stresses on a plane beach. In both tests the present model manifests an onshore
surface ﬂow and compensating offshore near-bed undertow near the shoreline and around the breaking
point. In contrast, the radiation-stress prescription yields signiﬁcantly weaker vertical shear. The currents’cross-shore and vertical structure is signiﬁcantly shaped by the wave effects of near-surface breaker acceleration, vertical component of vortex force, and wave-enhanced pressure force and bottom drag.
Ó 2010 Elsevier Ltd. All rights reserved.

1. Introduction
The effects of wind-driven (primary) surface gravity waves on
oceanic currents and turbulence (hereafter called WEC) have beenrecognized to play a crucial role for scientiﬁc and engineering applications, ranging from wave-induced upper-ocean mixing andcurrent proﬁles to littoral ﬂow, sea level, and sediment transport
relevant to beach management and navigation. An essential featureof most theoretical approaches to WEC is averaging over the fast
oscillations of the primary wind-driven waves, with seminal papers
by Longuet-Higgins (1970), Hasselmann (1971), Craik and Leibovich (1976), and Garrett (1976). Wave averaging is also necessary
for feasible computations of realistic circulations with WEC.
A central arena for WEC is the surf zone where breaking waves
accelerate alongshore and rip currents. The interplay between
waves and currents has been investigated mostly in one- or
two-dimensional, depth-averaged models with fast-wave averag-* Corresponding author. Present address: 3845 Slichter Hall, Los Angeles, CA
90095-1565, USA. Tel.: +1 310 825 5402; fax: +1 310 206 3051.
E-mail addresses: uchiyama@atmos.ucla.edu (Y. Uchiyama), jcm@atmos.ucla.
edu (J.C. McWilliams), alex@atmos.ucla.edu (A.F. Shchepetkin).
1463-5003/$ - see front matter Ó 2010 Elsevier Ltd. All rights reserved.
doi:10.1016/j.ocemod.2010.04.002

ing (recent studies by Ruessink et al. (2001), Yu and Slinn (2003),
Özkan-Haller and Li (2003), Reniers et al. (2004a) and Uchiyama
et al. (2009)). Alternatively, a phase-resolving horizontal twodimensional (2D) approach (e.g., Chen et al., 1999; Terrile et al.,2008) can depict wave–current interaction processes, albeit at aprohibitive computational cost for longer-term, larger-scale current evolution.
Several wave-averaged 3D circulation models have been createdduring the last decade. In Walstra et al. (2000) and Lesser et al.
(2004), the Delft3D-ﬂow code includes WEC by loosely adapting a
set of generalized Lagrangian mean (GLM) equations by Groeneweg(1999), adapted from Andrews and McIntyre (1978a,b). The model
prognostic ﬁeld is Lagrangian mean velocity u‘ and the wave-induced forcing in the ﬂow model is represented by the depth-averaged radiation-stress gradient (e.g., Longuet-Higgins and Stewart,
1962; Hasselmann, 1971; Phillips, 1977), although in practice it isexpressed in terms of breaking and frictional dissipation terms provided by a wave model in accordance with Dingemans et al. (1987)
and imposed in the ﬂow model as surface and bottom stresses. A
simple, geostrophic 3D GLM ocean model was proposed by Perrie
et al. (2003) where Stokes–Coriolis force (Hasselmann, 1971) anda surface-intensiﬁed acceleration due to wave dissipation are taken

Y. Uchiyama et al. / Ocean Modelling 34 (2010) 16–35

into account as WEC. These models neglected the conservative vortex force (VF) and quasi-static set-down (i.e., equivalent to a pressure contribution in the Longuet-Higgins and Stewart (1962)
radiation stress). Another branch of engineering-oriented 3D modeling with WEC is by Xie et al. (2001), following Lewis (1997). It applies
the depth-averaged radiation-stress gradient as a depth-uniformbody force in the Princeton Ocean Model (POM; Blumberg and Mellor (1987)). Later, a depth-dependent form of horizontal radiationstress gradient terms was proposed by Xia et al. (2004) in a theoretically ad hoc way. Warner et al. (2008) (W08) employs a GLM-likevertical mapping approach with a depth-dependent radiation-stressformalism proposed by Mellor (2003, 2005) in the Regional Oceanic
Modeling System (ROMS) code. It is recognized (Ardhuin et al.,2008a) that accurate implementation of this formalism (also the‘‘alternative” GLM equations in Andrews and McIntyre (1978a)) requires knowledge of the wave kinematics to higher order in parameters that deﬁne the large-scale evolution of the wave ﬁeld, such as
the bottom slope; this impracticality is addressed in Mellor (2008),
but it does not yet seem to have been conveyed into the W08 code.
In these 3D models, the WEC are represented as the radiation stressgradient. Based on the Helmholtz decomposition of the advection
terms in the equations of motion, the VF representation comes fromthe identity, u  ru = rjuj2/2 + (r  u)  u, while the radiationstress representation arises from the Reynolds decomposition,u  ru = r  (uu) + u(r  u), together with incompressibilityr  u = 0, where u is the Eulerian velocity. The primary advantageof the wave-averaged VF formalism is its explicit inclusion of a typeof wave–current interaction that few if any available wave models
properly incorporate to allow its complete expression in the radiation stress (e.g., Lane et al., 2007). A conspicuous demonstration of
the utility of a VF representation is Langmuir circulations in theupper ocean (Craik and Leibovich, 1976; Leibovich, 1980; McWilliams et al., 1997). The GLM approach with a VF formalism is taken
in Ardhuin et al. (2008b) and advocated as appropriate for a wide
range of oceanic applications. This formulation is applied to vertical
one-dimensional modeling of the ocean mixed layer by Rascle et al.
(2006, 2009) and tested in a nearshore 2D (cross-shore/vertical) Rascle, 2007. Instead, we utilize an Eulerian reference frame for thewave averaging (e.g., McWilliams and Restrepo, 1999; McWilliams
et al., 2004; Lane et al., 2007), primarily for more direct comparability to most measurements and compatibility with existing circulation models. Within asymptotic approximations the two
approaches are equivalently valid as long as the Lagrangian andEulerian mean velocities are related by u‘ = u + uSt, where uSt is
the 3D Stokes drift velocity. When the model prognostic variable isu‘, care must be taken to retrieve u to estimate horizontal and vertical mixing, bed shear stress, and boundary conditions for realistic
applications.
For littoral currents in the surf zone, Longuet-Higgins (1973) andDingemans et al. (1987) show that conservative and non-conservative contributions in WEC are separable within the radiation stressdivergence when a geometric optics (WKB, ray theory) approximation is applied to the wave ﬁeld; the VF is contained within the radiation stress divergence, but they argued that it is negligible in asurface zone when the breaking-induced acceleration dominates.
However, this assumption has been partially falsiﬁed in several
models with a VF formalism: barotropic models (Smith, 2006;
Uchiyama et al., 2009), a quasi-3D model (Shi et al., 2006), and
3D models (Newberger and Allen, 2007a,b, called ‘‘nearshore POM”
and designated by NA07, and Delft3D-ﬂow). Multiple aspects of
WEC are expected to be important for surf zone currents. In addition to the conservative effects of VF, Bernoulli head, and quasi-static pressure response, there are important non-conservative effects
due to depth-induced breaking (and white capping) near the surface and frictional wave dissipation near the bottom. For 3D conﬁgurations the last two components should be applied at appropriate

17

depths, near the surface for the former and right above the wave
bottom boundary layer for the latter.
In this paper, we develop and test a 3D oceanic circulation model(extending ROMS; Shchepetkin and McWilliams, 2005) withdynamically consistent wave–current interactions suitable to awide range of nearshore, coastal, and open-ocean applications.
We base the model on the Eulerian-averaged, multi-scale, wave–
current asymptotic theory derived in McWilliams et al. (2004)(MRL04). It uses a VF formalism that cleanly separates conservativeand non-conservative WEC mechanisms, unlike the radiation-stressformalism. The conservative part of WEC comprises the VF (Leibovich, 1980), the Stokes–Coriolis force (Hasselmann, 1971), Bernoulli
head, and a quasi-static pressure response known as wave-setup/
down. Non-conservative WEC are included as a surface-concentrated 3D acceleration and wave-enhanced vertical mixing due todepth-induced wave breaking and associated surface rollers; a bottom-conﬁned bottom streaming stress (e.g., Longuet-Higgins,1953) caused by near-bed wave drag; and a wave-enhanced current
bed shear stress (e.g., Soulsby, 1995). The wave-induced verticalmixing is represented as an extension of the KPP model (Largeet al., 1994), a non-local vertical mixing parameterization. The governing current and wave equations are presented in Section 2, and
the WEC implementation in ROMS is detailed in Section 3 with particular attention to non-conservative wave effects adapted fromprevious parameterizations. Section 4 describes the application to
the surf zone during the DUCK94 experiment. In Section 5 a comparison to another ROMS-based 3D wave–current model by W08
and Haas and Warner (2009) (HW09) based on a radiation-stressformalism is made for an idealized plane beach, using an identicalwave ﬁeld from SWAN (Booij et al., 1999). Section 6 provides a summary and an outlook for future applications of the model.
2. Governing equations
The WEC model formulation is built on a sequence of previous
developments. McWilliams et al. (2004) (MRL04) derives a multiscale asymptotic model for the phase-averaged, conservative
dynamical effects of surface gravity waves on currents and infragravity waves with longer space and time scales. MRL04 extends
earlier derivations that feature the central WEC role of VF (Craik
and Leibovich, 1976; Garrett, 1976; McWilliams et al., 1997;
McWilliams and Restrepo, 1999), and this approach is set in thecontext of the larger literature on wave–current interaction in Laneet al. (2007) (LRM07). Uchiyama and McWilliams (2008) present abarotropic ROMS model for the effects of primary wind waves on
infragravity waves, while Uchiyama et al. (2009) presents a vertically averaged (barotropic) ROMS model for surf zone shear instability with WEC. In this and the next section, we describe thegoverning equations and ROMS implementation for 3D WEC andan accompanying surface wave model with effects of currents on
the waves (hereafter called CEW).
2.1. Wave-averaged currents
ROMS is a hydrostatic, incompressible (Boussinesq), free-surface model with non-conservative forcing, diffusion, and bottom
drag. It makes a baroclinic–barotropic mode split, with explicit fasttime-stepping and subsequent conservative averaging of barotropic variables. All present elements in ROMS are retained, but newterms are added to incorporate WEC.1 The ROMS formulation is1
The recent generalization to a non-hydrostatic ROMS model (Kanarska et al.,2007) will include the same wave-averaged effects discussed here. Additional terms
are added to the momentum equations: non-WEC terms for vertical acceleration andfor horizontal Coriolis frequency f y and WEC terms for Stokes–Coriolis force with f y
and for VF with full horizontal vorticity.

18

Y. Uchiyama et al. / Ocean Modelling 34 (2010) 16–35

non-asymptotic in the sense that some additional non-wave terms,
beyond the minimum required for asymptotic consistency as deﬁned
in MRL04, are included for completeness (e.g., the time-derivative of
surface elevation in the kinematic boundary condition and depthintegrated mass balance), along with additional non-conservative
wave effects (e.g., breaker acceleration).
We ﬁrst write the model equations in Cartesian (x,y,z,t) coordinates. The notation is slightly different from MRL04, and the quantities are dimensional. We combine the infragravity wave andcurrent dynamics, which were asymptotically separated in
MRL04. The momentum balance is written in terms of a dynamic
pressure / (normalized by mean density q0) and sea level f after^ and ^fsubtracting the wave-averaged quasi-static components /(n.b., MRL04, Sections 6 and 9.2.3, and LRM07, Eqs. (3.8)–(3.10))
that occur even without currents. All wave quantities are referenced to the local wave-averaged sea level, z ¼ f þ ^f, rather thanthe mean sea level, z = 0. The vertical coordinate range ishðxÞ 6 z 6 f þ ^f. The equations make the particular gauge choicefor the decomposition between VF (J, K) and Bernoulli head K described in MRL04, Section 9.6. The new WEC terms for ROMS are
written on the right side of the equations below. Boldface vectors
are horizontal only, and 3D vectors are designated by (horizontal,
vertical).

ð2Þ

h

h(x) is the resting depth of the ocean; A is the wave amplitude; k is
its wavenumber vector and k is its magnitude;ð3Þ

is its intrinsic frequency; and normalized vertical lengths are
H ¼ k h þ f þ ^f  kD;

and Z ¼ kðz þ hÞ;

ð4Þ

where D ¼ h þ f þ ^f is the wave-averaged thickness of the watercolumn. The horizontal and vertical VF (inclusive of the Stokes–
Coriolis term) and Bernoulli head (after removing quasi-static
terms) are

J ¼ ^z  uSt ðð^z  $?  uÞ þ f Þ  wSt

ou
;
oz

ou
;
oz
Z z 2
1
rA2
o V
0
sinh½2kðz  z0 Þ dz ;K¼
2
4 ksinh ½H h oz0 2

K ¼ uSt 

The quasi-static sea-level component is deﬁned by

A2 k^f ¼  patm :g q0 2 sinh½2H

ð7Þ

It contains both an inverse-barometric response to changes inatmospheric pressure patm and a wave-averaged setup/setdown.
With a multi-component wave ﬁeld, A2 is replaced in (2)–(7) by
the sea-level spectrum G(h,r) with integration over wavenumbervector angle h and frequency r. This implies a superposition of the
WEC contributions from different components, consistent with theasymptotic theoretical assumption of small wave slope Ak.
2.2. Boundary conditions
The boundary conditions for ROMS include the usual stress andheat and material ﬂux conditions plus the following kinematic andpressure continuity conditions, again with the additional WEC
terms on the right side:

ð8Þ

with

A2 r

qﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ
gk tanh½H

ð6Þ

of
o^f
 ðujfþ^f  $? Þf ¼ $?  USt þ þ ðujfþ^f  $? Þ^f;ot
ot
gf  /jfþ^f ¼ P;

ð1Þ

r¼


2
1 o A sinh½Z
:2 ot sinh½Hwjfþ^f 

F is the non-wave non-conservative forces, Fw is the wave-induced
non-conservative forces, c is any material tracer concentration (e.g.,
T and S), and C is the non-conservative tracer forcing, where $? isthe horizontal differential operator. The system (1) is completedwith the equation of state.
The 3D Stokes velocity (uSt, wSt) is non-divergent and deﬁned fora monochromatic wave ﬁeld bycosh½2Zk;
22sinh ½HZ z
0wSt ðzÞ ¼ $? uSt dz :E¼

wjh þ ujh  $? h ¼ 0;ououþ ðu  $? Þu þ w
þ f ^z  u þ $? /  F ¼ $? K þ J þ Fw ;otozo/ g q
oK
¼
þ
þ K;
oz q0
oz
ow
$?  u þ
¼ 0;
oz


oc
oc
oc 1 oocþ ðu  $? Þc þ w  C ¼ ðuSt  $? Þc  wSt þ
E
:
ot
oz
oz 2 ozoz

uSt ¼

with V ¼ k  u. The wave-induced tracer diffusivity is deﬁned byð5Þ

(
gA2 tanh½H
oVoVP¼
þ cosh½2Hoz fþ^foz h2r sinh½2H
!
)Z fþ^f 2
o V
0
0
þ

2k
tanh½HVj
cosh½2kz

dz
^
fþf :
oz0 2
h

ð9Þ

In MRL04, Section 9.3, there are additional quasi-static components
in P of higher asymptotic order in the wave slope Ak, but, unlike in ^f
in (8), they have no dynamical coupling with the currents in (1) and(8). So, without a speciﬁc motivation for examining the various deleted quasi-static terms, they are not presently included in ROMS,
although they could easily be added as a diagnostic.
2.3. Barotropic mode
The barotropic mode is derived from (1) as a vertical integral of
the continuity equation and a vertical average of the horizontal
momentum equation. With the WEC terms kept on the right side,
the result isofo^fþ $?  U ¼   $?  USt ;
ot
ot
!
!
 1
ou
U
1
U
þ    ¼  ð$?  USt Þ ujfþ^f 
þ ð$?  UÞ ujfþ^f 
ot D
D
D
D
Z fþ^f1
þ½J  $? K dz þ Fw :
D h
ð10Þ
The dots in the barotropic momentum equation indicate contributions from all the left-side terms in the horizontal momentum
equation in (1) other than the acceleration. HereU¼

Z fþ^fhu dz and USt ¼

Z fþ^fuSt dzð11Þhare the horizontal volume transports by Eulerian and Stokes
 ¼ U=D is the barotropic velocity. (Notecurrents, respectively, and u

19

Y. Uchiyama et al. / Ocean Modelling 34 (2010) 16–35

that the depth integration and averaging violates the strict separation of non-wave and wave-averaged terms on the left and rightsides of the equations above.) We can combine the barotropic continuity and momentum equations in (10) to write the latter in theform used in ROMS:o
U þ    ¼ ujfþ^f $?  U þ USt þotZ fþ^f½J  $? K dz þ DFw :

ð12Þ

h

The free-surface equation in (10) implies the volume conservation
relation,

d
dt

Z Z 
I

h þ f þ ^f dx ¼ 

UþU

St

^ ds:
n

ð13Þ

Thus, mean sea level within a domain is controlled by the boundary
Eulerian and Stokes ﬂuxes.
In a barotropic ROMS model, u f þ ^f ¼ U=D, and the remaining
 in J from (5) and K ¼ 0.
right-side terms are evaluated with u ¼ u
The associated tracer variable is also equated with its depth-averaged value, c ¼ c. The resulting model has been applied to infragravity-wave and nearshore barotropic shear-instability problems(Uchiyama and McWilliams, 2008; Uchiyama et al., 2009).
2.4. Wave dynamics
For a monochromatic wave ﬁeld, a WKB wave model with wave
refraction and conservation of wave action is the following:okkr~ $ u~$? Dþ cg  $? k ¼ k? otsinh 2kDwoAþ $?  ðAcg Þ ¼  ;
ot
r

ð14Þ
ð15Þ

3. Implementation in ROMS
3.1. Wave-averaged 3D currents and tracers
As a prelude to discretization in curvilinear coordinates, we rewrite several of the WEC relations in Section 2 in forms closer tothose used in ROMS, adopting a ﬂux-divergence form of the substantial derivative and deﬁning three new variables,

fc ¼ f þ ^f;/c ¼ / þ K;

ð18Þu‘ ; x‘ ¼ ðu; xÞ þ uSt ; xSt ;where fc is a composite sea level, /c absorbs the Bernoulli head,(u‘, x‘) is the wave-averaged Lagrangian velocity, and x is the vertical velocity in the transformed coordinate system used in ROMS
(see (25) et seq.). The 3D primitive equation (1) become


 o
ou
~ þ Fw ;þ $~?  u~‘ u þ ðw‘ uÞ þ f ^z  u‘ þ $? /c  F ¼ u~St $?  uotozo/c g q
¼ K;
þ
q0
oz
ow‘
$ ?  u‘ þ
¼ 0;
oz


oc
o
1 o
oc
~ ‘ c þ ðw‘ cÞ  C ¼:
þ $~?  u
E
ot
oz
2 oz
oz
ð19Þ
The WEC terms are no longer conﬁned to the right sides of theseequations. The boundary conditions (8) arew‘ jh þ u‘ jh  $? h ¼ 0;
ofc
 u‘ jfc  $? fc ¼ 0;
ot
gfc  /c jfc ¼ P þ g^f  Kjfc :using the tilde convention to identify conjoined horizontal vectors
in a dot product. The wave action is deﬁned by A ¼ E=r where
E ¼ 12 q0 gA2 is the depth-integrated wave energy.

w‘ jfc rds ¼ u  k þ r

The depth-integrated continuity equation in (10) is

ð16Þ is the depthis a Doppler-shifted (CEW) wave frequency, where u
averaged current,2 r is the intrinsic frequency (3), and the associated group velocity isþ
cg ¼ u


1þ
2

r
2k


2kD
k:
sinh 2kD

ð17Þw

 is the depth-integrated rate of wave energy loss (or dissipation).
In the present formulation, we include wave dissipation due todepth-induced breaking and bottom drag, both of which must be
parameterized (Section 3.2).
In some realistic cases, this model is applied with k = kp the pﬃﬃﬃ
pﬃﬃﬃﬃﬃﬃ
spectrum-peak wavenumber and A ¼ 2E=g ¼ Hsig = 2 2 ¼
H =2 the equivalent wave amplitude in terms of the wave energy
E, the so-called signiﬁcant wave height Hsig, or a wave height H
commonly used in breaking parameterizations. The simplest
extension from a monochromatic/spectrum-peak model to a multi-component model is based on superposition of components with
spectrum G. For more general wave dynamics including nonlinear
spectrum evolution and wind generation, a wave simulation modelis used (e.g., SWAN; Booij et al., 1999) to provide G and w. One
may also specify G from available observations, e.g., an offshore
wave buoy.
2
This CEW theory is strictly valid only for depth-uniform currents because strong
vertical current shear invalidates the wave eigenmodes and dispersion relation onwhich the asymptotic WEC theory is based. In practice we use either the depthaverage current in shallow water or an upper-ocean average over a depth / k1 indeep water for the wave model. MRL04 did not include CEW because its scaling
assumptions were that current speed and sea-level elevation are, respectively,
smaller than gravity-wave propagation speed and resting depth.

ð20Þ

ofc
þ $?  U‘ ¼ 0;
ot

ð21Þ

where U‘ is the depth integral of u‘. The associated barotropic horizontal momentum equation (12) isoU þ $~? ot

Z fc~ ‘ u dz ¼ gD$? fc þ F0 þ
u

h

Z fc

R dz;

ð22Þ

h

where F0 is the baroclinic part of the full vertically integrated pressure-gradient force,

F0 ¼ gD$? fc Z fc$? /c dz;

ð23Þ

h

containing the usual terms proportional to gq0 /q0 as well as contributions from wave effects, but excluding the barotropic free-surface
gradient term. Other terms in the 3D momentum equation in (19)have been lumped into a residual horizontal vector,~ þ Fw :R ¼ f ^z  u‘ þ F þ u~St $?  u

ð24Þ

The ROMS equations are expressed in horizontal orthogonal curvilinear and vertical surface- and terrain-following coordinates(n, g, s) (Shchepetkin and McWilliams, 2005, 2008). We use the composite sea-level variable fc and h to deﬁne the s coordinate. m1 andn1 are Lamé metric coefﬁcients, and Hcz is the transformed grid-cell
thickness. The 3D primitive equations for wave-averaged currents
in ROMS are the following:

 




 
o Hcz
o Hcz u‘o Hcz v ‘
o x‘s
þ
þ
þ
¼ 0;
ot mn
on
og m
os mn
n

ð25Þ

20

Y. Uchiyama et al. / Ocean Modelling 34 (2010) 16–35



Du c‘
Hc o/cHcouov
þdFc v St þ z uSt
F v ¼ z
þ v StDton
n on znonHcz  nþF þ F wn ;mn

3.2. Non-conservative wave dissipation and rollers

ð26Þw ¼ b þ wd :



Dv c‘
Hc o/cHcouovd
Fc uSt þ z uSt
þ v StþF u¼ zog
Dtm og z
m
og
þ

Hcz
mn

ðFg þ F wg Þ;

ð27Þ

where D=Dt is the material derivative in conservation form in curvilinear coordinates,









D
o Hcz
o Hcz u‘o Hcz v ‘
o x‘s
¼
 þ
 þ
 þ
 ;
Dt ot mn
on
og m
os mn
n

ð28Þ

c‘ and dwhere FFc are generalized Coriolis frequencies combinedwith the curvilinear metric terms,

 
 
c‘ ¼ Hc f  u‘ o 1 þ v ‘ o 1 ;
F
z
mn
og m
on n

 
 
f
o
1
o
1
c
d
u
þv
:
Fc ¼ Hz
mn
og m
on n



n

ð29Þ



ð31Þ

W ¼

Z s
1

!
oU ‘ oV ‘
1 z þ h ofc
0
ds 


þ
:
mn f þ h ot
on
og


Z 0


gq
/c ¼ g fc  ^f  ðP  KÞjfc þ K Hcz ds:
s

q0

b ¼

ð33Þ

The 3D tracer equation with WEC is  Dc
o E oc¼Cþ
;Dtos Hcz os

where C includes both non-wave and wave-enhanced turbulent
mixing parameterizations (Section 3.4) and E is deﬁned in (6).
The wave model is solved before the predictor stage for thebaroclinic mode (Shchepetkin and McWilliams, 2005), and subsequently the WEC components are computed. For 3D simulations
these are kept unchanged during the barotropic time steps; for2D simulations, however, the WKB wave model is solved at everybarotropic time step, and the WEC terms are updated as fc and u
evolve. All the new terms associated with the conservative WEC
are discretized with the centered ﬁnite-differences in a manner
similar to the other terms in ROMS at the predictor and corrector
stages. The vertical VF in (33) is discretized with the density-Jacobian scheme (Blumberg and Mellor, 1987) to reduce the terrainfollowing-coordinate error. Correction of (u‘, x‘) with the updated
Stokes velocity occurs before the corrector steps. These procedures
enable us to minimize the code changes in ROMS. Notice that the
prognostic variables in ROMS with WEC are composite sea levelfc and Eulerian velocity (u, x).

uw
orb

!0:52
:

ð37Þ

pﬃﬃﬃﬃ
3 p
B3 f
q0 g 4b p5 H7 ;
16
cb D

ð38Þ

where fp is a peak wave frequency (1/Tp = x/2p with Tp a peakwave period), H = 2A (Section 2.4), and Bb and cb are empiricalparameters related to breaker types. An alternative parameterization by Thornton and Whitford (1990) is described by Church andThornton (1993):

b ¼
ð34Þ

rzo

Wave dissipation due to wave breaking b is a combination of deepwater breaking (e.g., white-capping) and depth-induced nearshore
breaking. Deep-water breaking in wind-wave equilibrium is integrally equivalent to the surface wind stress (Sullivan et al., 2007);in this paper for simplicity we will use the latter representation(consistent with ignoring wave generation in the monochromaticwave model in Section 2.4). However, the depth-induced breaking
is essential for surf zone wave–current interaction. A parameterization (Thornton and Guza, 1983) isð32ÞHere U ‘ ¼ Hcz u‘ =n; V ‘ ¼ Hcz v ‘ =m; and W ‘ ¼ x‘s =ðmnÞ are grid-cellvolume ﬂuxes. The geopotential function is evaluated from integration of the vertical momentum equation,

ð36Þ

where uw
orb ¼ rH =ð2 sinh kDÞ is bottom wave orbital velocity, andfw is a wave friction factor (Soulsby, 1997),

fw ¼ 1:39

and the vertical mass ﬂux is calculated as
‘

1
2 p

wd ¼ pﬃﬃﬃﬃ q0 fw uworb 3 ;

g



oz
;
þ u‘  $? z
ot
s

ð35Þ

These wave dissipation processes also imply WEC accelerations(Section 3.3).
Bottom viscous drag on the primary waves causes a dissipation
wd
wd and an associated wave-averaged bottom stress swd
bot ¼  k=r
that induces an Eulerian bottom streaming ﬂow in the directionof wave propagation (Longuet-Higgins, 1953). We use a parameterization by Reniers et al. (2004b) for the realistic regime of a turbulent wave boundary layer, which is based on the Rayleigh wave
height distribution in accordance with Thornton and Guza (1983)
consistent with the present WKB spectrum-peak wave modeling:ð30Þ

where F ¼ ðF ; F Þ is the non-wave body force and parameterized
momentum mixing term; Fw = (Fwn, Fwg) is the non-conservative
wave terms deﬁned later in this section. The vertical motion pasts surfaces isx‘s ¼ w‘ 

The primary wave dissipation rate in (15) is calculated as thesum of the effects of wave breaking b and wave bottom drag wd,

pﬃﬃﬃﬃ


 
3 p
B3 f
H
q0 g b p H3 1 þ tanh 8   1
16
D
cD
2
3
(

2 )52H
5;
 41  1 þ
cb D

ð39Þ

where Bb and cb are again empirical constants depending on typesof breaking. In the DUCK94 simulation, the latter parameterization
is found to be more successful (Section 4).
To better estimate surf zone currents, an additional wave component is sometimes included (Svendsen, 1984a), i.e., surface rollers, which are onshore-traveling bores of broken primary waves.
The idea is that some fraction ar is converted into rollers that propagate toward the shoreline before dissipating, while the remaining
fraction 1  ar causes local dissipation (hence current acceleration). Here we introduce a surface roller model by primarily following Nairn et al. (1991) and Reniers et al. (2004a) with minor
modiﬁcations. The roller is assumed to have the same k as thebreaking primary wave in (14). The evolution equation for the roller action density Ar is analogous to (15) for A:oArar b  r;
þ r  ðAr cÞ ¼
ot
r

ð40Þ

21

Y. Uchiyama et al. / Ocean Modelling 34 (2010) 16–35

 þ rk2 k iswhere Ar ¼ Er =r; Er is roller energy density; and c ¼ u
the phase speed of the primary wave (not cg; Svendsen, 1984a; Stiveand De Vriend, 1994). While most previous studies, including
Svendsen (1984a), Nairn et al. (1991) and Apotsos et al. (2007), assume that the full primary wave b feeds the roller energy density(i.e., ar = 1), Tajima and Madsen (2006) introduce ar, 0 6 ar 6 1.
We view the latter as useful for correcting b with some ﬂexibility
to depict different breaking wave and beach forms (e.g., spillingor plunging breakers, barred or plane beaches); however, its value
is an ad hoc choice. The roller dissipation rate is then parameterized
as

r ¼

g sin bEr
;
c

ð41Þ

where sin b = 0.1 is an empirical constant (Nairn et al., 1991; Reniers
et al., 2004a). According to Svendsen (1984a), the roller Stokes
transport for a monochromatic primary wave is

Ur ¼

Er

q0 r

k¼

Ar

q0

k;

ð42Þ

ðE þ Er Þ

q0 r

k¼

ðA þ Ar Þ

q0

k:

ð43Þ

We assume that Ur is vertically distributed similarly to the Stokes
drift velocity of the primary waves, henceuSt ðzÞ ¼

r2 cosh 2kðz þ hÞ
2

gsinh kDðA þ Ar Þk:

ð44Þ

The expression for primary-wave Stokes drift in (44) is for nonbreaking, small-slope waves that may not be accurate in the surfzone. We also assume for simplicity that (44) is applicable to theroller Stokes drift, although NA07b represented it with a surfaceintensiﬁed vertical structure. (The same assumption about (44) for
the roller waves is also made in other models such as the Lagrangian-mean radiation-stress model by W08.) There is room for future
investigation.
A coupled wave simulationpmodelsuch as SWAN (Booij et al.,ﬃﬃﬃ1999) provides kp, rp, Hsig ¼ 2H ; b ; wd , and Qs, where Qs is
the fraction of broken waves (0 6 Qs 6 1). As originally given by
Svendsen (1984a), with consideration of Qs, the roller action density is then

Ar ¼

Er

r

¼

q0 gDARQ ;2Lp rp s

ð45Þ

where Lp(= 2p/jkpj) is a primary wave length and AR is a roller area inthe vertical plane estimated from the formulas proposed by Svendsen (1984b) or Okayasu et al. (1986),

AR ¼ asR H2

or AR ¼ aoR H Lp ;ð46Þ

asR ¼ 0:9 and aoR ¼ 0:06.

with empirical constants
The latter isadopted to shift the peak undertow velocity farther shoreward under plunging waves.
3.3. Non-conservative wave accelerations for currents
The wave dissipation processes (b, r, and wd in Section 3.2)
have accompanying WEC accelerations Aw (e.g., Dingemans et al.,1987). We distinguish Aw from the wave-enhanced turbulent vertical mixing Dw and bottom drag stress scd
bot discussed in Section 3.4.
Thus,

Fw ¼ Aw þ Dw ;Aw ¼ Bb þ Bwd ;

Bb ¼

ð47Þ

for brevity Bb contains both the depth-induced breaking and roller
accelerations. We represent the accelerations either as body forcesð1  ar Þb þ r

q0 r

kf b ðzÞ;

ð48Þ

where fb(z) is a vertical distribution function representing vertical
penetration of momentum associated with breaking waves and rollers from the surface. It is normalized as

Z fc

0

f b ðz0 Þdz ¼ 1;

ð49Þ

h

hence the vertical average of Bb (i.e., barotropic acceleration) is

Bb ¼

hence the total Stokes transport is

USt ¼

or as equivalent boundary stresses if the associated turbulent
boundary layers are too thin to be resolved in a particular model
conﬁguration. For simplicity we revert to Cartesian coordinates in
the rest of this section, with implied translation into transformed
coordinates for ROMS along the lines indicated in Section 3.
The breaking and roller accelerations enter as a body force
through Fw in the current momentum equations (26) and (27).
They are expressed asð1  ar Þb þ rk:
q 0 rDð50Þ

We can alternatively incorporate the breaking acceleration as anequivalent surface stress boundary condition for u instead of a bodyforce (e.g., as done in NA07):bssur ¼ swind
sur þ ssur ;ð51Þ

where swind
is the usual oceanic-model representation of surface
surwind stress andsbsur ¼ q0 DBb

ð52Þ

is the stress due to primary wave breaking and rollers. To examine
the sensitivity to the choice of fb, we consider three alternative
shapes:Type I :f b ðzÞ / 1  tanh½kb ðfc  zÞ4 ;Type II :f b ðzÞ / 1  tanh½kb ðfc  zÞ2 ;Type III :

ð53Þ

b

f ðzÞ / cosh½kb ðz þ hÞ

(leaving out the normalization factors for (49)). The vertical length1scale kb controls the penetration depth in each of these shape func1
tions. Usually we represent kb ¼ ab H , where ab is an Oð1Þ constant.
b
The ﬁrst f in (53) is proposed in W08 to account for the depthdependent radiation stresses divergence associated with rollers;
the second one is an alternative that concentrates the breaking effects nearer the surface; and the last one is inspired by the structureof the primary wave and further matches the vertical scale of its
velocity variance with kb = 2k (i.e., a choice based on wavelength
not wave amplitude).
Wave-induced bottom streaming (Longuet-Higgins, 1953; Xu
and Bowen, 1994) is similarly represented as a body force:

Bwd ¼

wd wdkf ðzÞ;q0 r

ð54Þ

where fwd(z) is a vertical distribution function normalized as in (49)for the Reynolds stress divergence associated with the turbulent
wave bottom boundary layer (WBBL). We employ three upward
decaying functions fwd analogous to fb,

Type I :f wd ðzÞ / 1  tanh½kwd ðh þ zÞ4 ;Type II :f wd ðzÞ / 1  tanh½kwd ðh þ zÞ2 ;Type III :

ð55Þ

c

f wd ðzÞ / cosh½kwd ðf  zÞ;1with a decay length kwd ¼ awd dw . awd is a constant, and dw is the
WBBL thickness,

22

Y. Uchiyama et al. / Ocean Modelling 34 (2010) 16–35

dw ¼ 0:09kN

 w 0:82
Aorb
:
kN

ð56Þ

w
Aw
orb ¼ juorb j=r is a semi-orbital excursion of short waves; kN = 30zo
is the Nikuradse roughness; and zo is roughness (Fredsøe and Deigaard, 1995); dw is typically only a few cm. awd = 1 corresponds to thetheoretical turbulent WBBL thickness associated with monochromatic waves, whereas it is known that there exists a signiﬁcant
increase in awd under random waves based on laboratory measurements (Klopman, 1994); e.g., Reniers et al. (2004b) use awd = 3. Its
depth integral has an equivalent effect to a bottom stress,

wd
swd
¼
bot ¼ q0 DB

wd
k:
r

ð57Þ

1

When kwd is too thin to be resolved on the model grid, then thestreaming acceleration is applied only as a stress in the bottom gridcell.
This formula implies bottom streaming occurs in the directionof wave propagation, consistent with the viscous streaming inlaminar (Longuet-Higgins, 1953) and weakly turbulent (LonguetHiggins, 1958) regimes under sinusoidal forcing. In contrast, nonlinear waves with the second-order Stokes theory and asymmetric
forcing in rough turbulent WBBLs reduce the Longuet-Higgins
positive streaming, and under some circumstances even manifest
opposite ﬂow (e.g., Trowbridge and Madsen, 1984; Davies andVillaret, 1999). In the surf zone, near-bed undertow opposite to
the incident waves dominates over streaming (Section 4.6).
There could additionally be a surface streaming ﬂow due
toa wave-viscousboundary layer in a thin layer of thickness,
pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ
ﬃ
2m=r  1 mm, where m is the molecular viscosity (LonguetHiggins, 1953), but we take the view that it is negligibly small,especially in the presence of wave breaking.

based on the Mellor–Yamada level 2.5 closure model (Mellor and
Yamada, 1982; Galperin et al., 1988) by introducing a turbulent kinetic energy input at the surface and a modiﬁed surface roughness
scale with a prescribed bilinear relationship for the turbulent
length scale analogous to the law of the wall. Subsequently this approach has been further developed (Burchard, 2001; Umlauf et al.,
2003; Newberger and Allen, 2007b; Jones and Monismith, 2008).
We take a similar approach by deﬁning K bv ðzÞ consistent with thedepth-averaged eddy viscosity proposed by Battjes (1975) and vertically distributed with a shape function fkV(z) (analogous to fb in
Section 3.3):

K bv ðzÞ ¼ cb


1=3
ð1  ar Þb þ rq0

scd
bot ¼ sc 1:0 þ 1:2

In ROMS the parameterization of vertical mixing is often donein the framework of the non-local, ﬁrst-order turbulent closure
model, K-Proﬁle Parameterization, KPP (Large et al., 1994; Shchepetkin and McWilliams, 2010). Apart from wave effects KPP comprises a shear-convective surface boundary layer (with countergradient ﬂux as well as eddy diffusion), a shear bottom boundary
layer, and interior mixing due to stratiﬁed shear instability, breaking internal waves, and double diffusion. In the KPP formulation foreddy diffusivity Kv, the effects of different mixing processes aremostly superimposed as if they were independent. Exceptions
are that the boundary-layer rule for Kv overrides the interior rule
within the boundary layers and that in shallow locations where
the surface and bottom boundary layers overlap, we choose the local maximum value for the eddy diffusivity, Kv(z) = max[Kvsur(z), Kvbot(z)] (Durski et al., 2004; Blaas et al., 2007). It is likely that thedifferent mixing processes sometimes combine nonlinearly, so
we view the present approach as a preliminary one that has theadvantage of process completeness but should be reconsideredwhen more is known.
In the WEC implementation, Kv is augmented by wave-induced
mixing in both the surface and bottom boundary layers. The incremental wave-enhanced momentum mixing of momentum is due
to surface breaking (b) and bottom current drag (cd):

Kv

ðzÞ;

ð59Þ

where cb is a constant.3 This can be viewed as mixing-length estimate, with the velocity scale based on the breaker dissipation rate
and the length scale based on the local wave height H, modiﬁedby the vertical distribution function fkV(z) and with the depth factor
D to retain the depth-average value of Battjes (1975). Apotsos et al.
(2007) use cb = 1/14 based on deep-water wave dissipation measured by Terray et al. (1996), so we anticipate cb is an Oð0:1Þ parameter. Terray et al. (1996) report that the penetration depth of surface
turbulence is proportional to the wave height, with little reduction
in turbulent intensity to a depth of 0.7H. Hence, the depth-dependency of K bv can be slightly different from that of Bb. To allow distinct
vertical scaling for Bb and K bv , we will deﬁne a vertical scale for K bv
1with kK v ¼ aK v H , where kkV replaces kb in (53) and akV is an Oð1Þparameter.
In the bottom boundary layer due to current shear turbulence,
wave motions enhance the bottom shear; e.g., Soulsby (1995) proposes the drag law,"

3.4. Wave-enhanced vertical mixingH Df

sc ¼ q0





jsw j
jsw j þ jsc j2jlnðzm =zo Þjuju;

3:2 #;
ð60Þ

1
2
jsw j ¼ q0 fw uw
orb ;2

where sc and sw are bottom stresses due to current and waves; j isthe von Kármán constant; zm is a reference depth above the bed,
nominally equivalent to a half bottom-most grid cell height (in abarotropic model zm = D/2; (e.g., Uchiyama et al., 2009)); zo is the
bed roughness length; fw is the wave friction factor given by (37);
and juw
orb j is the bottom wave orbital velocity. As simpler alternatives for sensitivity testing, we deﬁne a linear bottom drag law,scd
bot ¼ q0 lu

ð61Þ(l is a linear drag coefﬁcient (m/s)) and a log-layer drag law,scd
bot ¼ sc ¼ q0



j

2

lnðzm =zo Þ

juju

ð62Þ

(zm and zo are interpreted as in (60)). The magnitude of K cd
v is pro1=2
in a KPP bottom boundary layer scheme.
portional to jscd
bot j
4. DUCK94 experiment

ð58Þ

For both model validation and dynamical interpretation, we
simulate the vertical proﬁle of horizontal velocity measured on anatural sandy beach at Duck, North Carolina, during the DUCK94

plus an equivalent diffusivity for tracers. In the surface region, K bv isadded to Kvsur, by the rationale in the preceding paragraph. In thebottom boundary layer, K cd
v is a generalization for the usual KPP
Kvbot(z) because of the wave-enhanced bottom stress.
Mixing near the surface due to breaking has been modeled withlocal turbulent closures. Craig and Banner (1994) propose a model

3
There is disagreement among local-closure modelers about the shape of K bv near
the surface, primarily because of different assumptions about the length scale proﬁle;e.g., Burchard (2001) has K bv decrease as z ? fc, while Jones and Monismith (2008) has
it increase. Our choice of fkV is monotonically increasing, essentially for proﬁlesimplicity. These distinctions probably matter only on a ﬁner vertical scale (i.e., afraction of H) than we should expect our model to be apt.


 ou
o  b
;
Dw ¼
ðzÞ
K v ðzÞ þ K cd
v
oz
oz

23

Y. Uchiyama et al. / Ocean Modelling 34 (2010) 16–35

experiment (e.g., Garcez Faria et al., 1998, 2000; Newberger and Allen, 2007b). The ﬁeld data were obtained on October 12, 1994,
when strong cross-shore currents were present associated with astorm. The vertical proﬁle of the littoral current was measured
with a vertical stack of seven electromagnetic current meters(EMCs) mounted on a mobile sled at elevations of 0.41, 0.68,
1.01, 1.46, 1.79, 2.24 and 2.57 m above the bed. The sled was located at a sequence of sites along a cross-shore transect; each site
sample lasted for 1 h, and seven samples were made across thetransect during the day. Horizontal velocity was also measured
with a spatially ﬁxed cross-shore array of 11 EMCs distributed
around the surf zone (e.g., Feddersen et al., 1998). Directional wave
spectra were measured on an alongshore line of 10 pressure sensors in 8 m depth (Long, 1996), and a cross-shore array of 13 pressure sensors was used to measure wave heights spanning the surfzone (Elgar et al., 1998). Further details of the data acquisition andprocessing are in Gallagher et al. (1996, 1998) and Elgar et al.
(1998).
The vertical current proﬁles in DUCK94 have previously beenmodeled in NA07 using an Eulerian-averaged WEC model with a
VF representation implemented within the POM code (Blumberg
and Mellor, 1987). Their formulation is dynamically consistent
within the shallow-water range, i.e., kh
1. It includes most of
the necessary wave-induced forcing for nearshore applications,
including conservative VF and quasi-static set-down as bodyforces; non-conservative forcing due to wave-breaking and associated surface roller as surface stresses; and wave-enhanced verticalmixing. A limitation of their model arises from the assumption ofkh
1, leading to neglect of the vertical variations in uSt and VF
and distortion of the breaking acceleration proﬁle.
We perform many simulations to expose 3D wave–current modeling sensitivities (Table 2). Run 1 is the baseline numerical experiment. It uses a type III shape function for fb(z) in (53) with ab = 0.2.
The KPP modiﬁcation by breaking relies on a type II shape function
for fkV(z) in (53) with akV = 1.2 and cb = 0.03. These choices represent
breaking acceleration as a shallow, surface-intensiﬁed body force.
For the bottom streaming acceleration, a type III shape is used forfwd(z) in (55) with awd = 3.0. Bottom drag is with the combinedwave–current model by Soulsby (1995) with zo = 0.001 m (Feddersen et al., 1998). The wave ﬁeld is evaluated by the WKB
spectrum-peak model with surface roller (ar = 1.0), the empiricalbreaking parameterization (39), and the wave drag dissipation(36). CEW comprising the frequency Doppler shift and changes inmean water column height due to wave-setup/down are includedin the wave model unless otherwise stated. The wave model istightly coupled with the current model at every ROMS baroclinic
time step. We omit stratiﬁcation (i.e., constant temperature andsalinity with no surface heat and freshwater ﬂuxes). We specify aweak lateral momentum diffusion with a constant coefﬁcient of
0.1 m2/s to obtain smooth solutions. The imposed forcing due towaves, wind, and tides are averaged over the duration of the sledmeasurements and held constant during the simulations. Modelparameters are summarized in Table 1.
The experiments in this section use the 3D code in a vertical (x–z)
2D mode by assuming alongshore y uniformity. The computational
domain is chosen as 768 m in the cross-shore direction (x) with
Dx = 2 m. For all runs 32 vertical s levels are used with grid-height
reﬁnement near the surface and bottom. The ﬁeld-surveyed bottom
topography with a bar around x = 120 m is used without any spatialsmoothing. Alongshore topographic uniformity is assumed, and a
periodic condition is imposed in the alongshore direction. A
Neumann condition is applied at the shoreward boundary to allow
mass and momentum exchange between the interior domain andthe very shallow shoreward region (x < 0). Chapman-type radiation
 are adapted at the offshore open
boundary conditions for fc and u towards ^f þ ftide andboundary with weak nudging for fc and u

Table 1
Common model parameters for DUCK94 simulations.
Variable

Value

Unit

Offshore wave height H
Offshore peak wave period Tp
Offshore incident wave angle ho
CT93 breaking parameter Bb
CT93 breaking parameter cb
Roller dissipation parameter sin b
Offshore tidal elevation ftide

1.6
6.0
193.0
0.64
0.31
0.1
0.7
0.2532

m
s
°
–
–
–
m
Pa

Cross-shore wind stress swind;xsur
Alongshore wind stress swind;ysur
Coriolis frequency f
Lateral momentum diffusion coefﬁcient Kh
Mean water density q0

0.1456Pa

8.5695  1050.1
1027.5

1/s
m2/s
kg/m3

uSt . A Neumann condition is used for all other variables at the offshore boundary. The baroclinic time step is 3 s with a mode-splitting ratio of 30. Simulations are initiated with a resting state andintegrated for 6 h to obtain steady solutions that are y-invariant(i.e., they are stable to littoral shear instability, consistent with the
DUCK94 observations).
4.1. Waves and depth-averaged currents
As background for the 3D simulations, we ﬁrst examine twobarotropic runs with and without the roller model (Runs 2d1 and
2d2). The cross-shore proﬁles of the wave ﬁeld for Run 2d1 andthe bottom topography are in Fig. 1a. The simulated H(x) agrees
reasonably well with the observed wave height (Elgar et al.,1998). The three dissipation terms in the wave model (Fig. 1b) demonstrate that depth-induced breaking b occurs at two locations,around the bar crest and the nearshore region. The roller dissipation
r peaks slightly shoreward of b by design. The frictional bottom
streaming dissipation wd is about one order of magnitude smaller
than the others around the breaking points, but it is dominant in
the offshore region (x > 500 m). Because the WKB ray model is independent of the roller model, b and wd are identical in Run 2d2.
 ¼ ðu; v Þ, and dyWe compare the depth-averaged velocity, unamic sea level, fc  ftide, among Runs 2d1–2d2 and two analogous
3D simulations (i.e., Runs 1 and 2, with and without the roller
model) in Fig. 1c–e. A roller has signiﬁcant effects. The cross-shore is altered by the roller contribution to uSt because it mustvelocity u
be an anti-Stokes ﬂow in alongshore-uniform, steady-state solutions as required by barotropic mass conservation (Section 5.1;
Uchiyama et al., 2009). The differences between 2D and 3D models
 than in u
. v
 generally increases
are appreciable, more so in fc and v
towards the shore, particularly beyond the breaking point aroundthe bar crest, and then it diminishes toward the shore. The roller
 locations shoreward and weakens the crosspushes the peak v
 gradient in the 2D and 3D cases. In the 3D runs, the peak
shore v
v is reduced and the alongshore momentum is distributed further
shoreward than the 2D cases, due to the bottom drag modiﬁcationand vertical momentum imbalance via vertical mixing (Section
4.8). The 3D Run 1 with shallow breaking and roller provides the from Feddersen
best agreement with the observed ﬁxed-array v
 and fc vary in x, the resultant wave ﬁelds
et al. (1998). Although uare nearly the same among the different cases, indicating that CEW
plays a small role in DUCK94 (in contrast to its much larger role inrip currents and littoral shear instabilities Yu and Slinn, 2003;
Özkan-Haller and Kirby, 1999; Uchiyama et al., 2009).
4.2. Vertical structure
To expose Bb depth-dependency and roller contributions in(48), we compare four 3D simulations (Runs 1–4; Table 2). Runs

24

Y. Uchiyama et al. / Ocean Modelling 34 (2010) 16–35

Fig. 1. Cross-shore proﬁles of (a) resting depth h, wave height H (observed and modeled), and modeled wave angle h; (b) q1 times the wave dissipation rates by depth ; and (e) v
 along with the observed alongshore velocity from the ﬁxed-array EMCs (circles) (Feddersen et al.,induced breaking b, roller r, and bottom drag wd; (c) fc; (d) u
1998).

Table 2Computational conﬁgurations for the DUCK94 simulations. ar = 0 means no roller component. Options for bottom drag are the Soulsby model (60) S95, linear (61) LIN, and loglayer (62) LOG. Bottom roughness zo (m) is used for S95 and LOG, while l (m/s) is used for LIN. For fb, fkV, and fwd, roman numericals indicate the shape types deﬁned in (53) and(55). ab, awd, and akV are length scale coefﬁcients for the shape functions. S indicates use of the streaming stress model, either for surface breaking (52) or bottom streaming (57).
SS denotes the Stokes scale, kb = 2k.
Run

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31

Waves

Bottom drag

Bb

Bwd

K bv

b

ab

f

wd

awd

f

kV

akV

cb

uerrorverror

0.2
0.2
SS
SS
0.2
0.2
0.2
0.2
–
0.1
0.5
1.0
–
0.2
0.2
0.2
0.2
0.2
0.2
0.2
0.2
0.2
0.2
0.2
0.2
0.2
0.2
0.2
0.2
0.2
0.2

III
III
III
III
III
III
III
III
–
III
III
III
III
III
III
III
III
III
III
III
III
III
III
III
S
OFF
III
III
III
III
III

3.0
3.0
3.0
3.0
3.0
3.0
3.0
3.0
–
3.0
3.0
3.0
3.0
3.0
3.0
3.0
3.0
3.0
3.0
3.0
3.0
3.0
1.0
5.0
–
–
3.0
3.0
3.0
3.0
3.0

II
II
II
II
II
II
II
II
–
II
II
II
II
I
III
II
II
II
II
II
II
II
II
II
II
II
II
II
II
II
II

1.2
1.2
1.2
1.2
1.2
1.2
1.2
1.2
–
1.2
1.2
1.2
1.2
1.2
1.2
0.4
0.8
1.6
1.2
1.2
1.2
1.2
1.2
1.2
1.2
1.2
1.2
1.2
1.2
1.2
1.2

0.03
0.03
0.03
0.03
0.03
0.03
0.03
0.03
–
0.03
0.03
0.03
0.03
0.03
0.03
0.03
0.03
0.03
0.01
0.05
0.07
0.10
0.03
0.03
0.03
0.03
0.03
0.03
0.03
0.03
0.03

0.42700.75600.81660.87450.65180.56470.4896
0.40960.89520.41300.46170.56580.4336
0.41510.40010.68250.48290.40840.60190.47480.5256
0.58290.43170.42520.43460.47060.45200.55460.44900.43200.43450.09220.42770.2422
0.55820.20720.11130.0806
0.0956
0.6125
0.09650.0816
0.09450.09070.10280.22430.13970.13150.14560.1313
0.09610.12860.18180.09180.0926
0.0916
0.09240.30760.45320.75080.08650.1893

WEC

CEW

ar

Model

zo or l

f

ON
ON
ON
ON
ON
ON
ON
ON
OFF
ON
ON
ON
ON
ON
ON
ON
ON
ON
ON
ON
ON
ON
ON
ON
ON
ON
ON
ON
ON
ON
ON

ON
ON
ON
ON
ON
ON
ON
OFF
OFF
ON
ON
ON
ON
ON
ON
ON
ON
ON
ON
ON
ON
ON
ON
ON
ON
ON
ON
ON
ON
ON
ON

1.0
0.0
1.0
0.0
0.25
0.50
0.75
1.0
–
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0

S95
S95
S95
S95
S95
S95
S95
S95
LOG
S95
S95
S95
S95
S95
S95
S95
S95
S95
S95
S95
S95
S95
S95
S95
S95
S95
S95
S95
LIN
LIN
LOG

0.001
0.001
0.001
0.001
0.001
0.001
0.001
0.001
0.001
0.001
0.001
0.001
0.001
0.001
0.001
0.001
0.001
0.001
0.001
0.001
0.001
0.001
0.001
0.001
0.001
0.001
0.005
0.010
0.004
0.008
0.001

III
III
III
III
III
III
III
III
–
III
III
III
S
III
III
III
III
III
III
III
III
III
III
III
III
III
III
III
III
III
III

1 and 2 have breaking as a shallow body force, while Runs 3 and 4
have a weaker depth variation in Bb by setting kb = 2k (i.e., thesame proﬁle as uSt) in (53); we call the latter runs deep breaking
cases. The roller contribution is included in Runs 1 and 3 and absent in Runs 2 and 4. Figs. 2 and 3 display (u, w), (uSt, wSt) and Kv

Normalized r.m.s. errors

for Run 1 in the x–z-plane (this case has the best match with the ; Fig. 1e). u(x, z) has a surf-zone overturningobserved barotropic v
circulation with a strong onshore ﬂow near the surface and an
opposing, offshore undertow near the bottom. This circulation pattern is most prominent around the bar crest. The largest negative v

25

Y. Uchiyama et al. / Ocean Modelling 34 (2010) 16–35(a) u (m/s)(b) v (m/s)

u (m/s)−0.5−0.25 0 0.25 0.5−0.5v (m/s)−0.8−0.6−0.4−0.2(c) w (m/s)

1
0

z (m)

−1
−2
−3
−4

St

1

(d) u

St

(m/s)

(e) v

0

w (m/s)−0.02 −0.01 0

St

(m/s)

(f) w

0.01

(m/s)

0

z (m)

−1
−2
−3
−4

uSt (m/s)−0.3 −0.2 −0.1
0

50

0

100 150 200 250
distance offshore, x (m)vSt (m/s)−0.03 −0.02 −0.01
300 0

50

0

100 150 200 250
distance offshore, x (m)wSt (m/s)−0.02 −0.01 0
0.01
300 0

50

100 150 200 250
distance offshore, x (m)300

Fig. 2. (a) Cross-shore velocity u, (b) alongshore velocity v, (c) vertical velocity w, (d) cross-shore Stokes velocity uSt, (e) alongshore Stokes velocity vSt, and (f) vertical Stokes
velocity wSt for Run 1. Contour intervals are (a) 0.05 m/s, (b) 0.05 m/s, (c) 0.004 m/s, (d) 0.02 m/s, (e) 0.002 m/s, and (f) 0.002 m/s.

(a)

(b)

(c)

1
0

z (m)

−1
−2

Kcd (m2/s)−3

Kb (m2/s)

v

0

−4
0

50

100 150 200 250
distance offshore, x (m)K (m2/s)v

0.005 0.01 0.015

0
300 0

v

0.015 0.03 0.045
50

100 150 200 250
distance offshore, x (m)

0
300 0

0.015 0.03 0.045
50

100 150 200 250
distance offshore, x (m)

300

b
cd
b
Fig. 3. Distributions of vertical eddy viscosity for Run 1; (a) KPP contribution K cd
v , (b) wave breaking contribution K v , and (c) K v ¼ K v þ K v . Notice that different color scales
are utilized. Contour intervals are (a) 0.001 m2/s, (b) 0.003 m2/s, and (c) 0.003 m2/s.

appears in the trough region shoreward of the bar. An increase of
Kv associated with an increase of K bv is observed around the bar
and near the shoreward boundary (Fig. 3), consistent with theb(x) proﬁle (Fig. 1b). The 3D stokes velocities are strongest at these
two breaking points, with depth-variation even in such a shallow
area (Fig. 2d–f). uSt is much stronger than vSt by an order of magnitude due to the small obliqueness of the incident waves. The divergence implied by uSt induces two wSt dipole circulations, withwave-induced upwelling adjacent to the shoreline and inshore of
the bar, and downwelling offshore of these locations. The verticalvelocity w is comparable in magnitude to wSt, and its primary
upwelling and downwelling centers occur in similar cross-shore
locations, albeit more bottom-concentrated, i.e., up- and downward ﬂows on the inshore and offshore sides of the bar in the offshore-headed near-bed undertow. The vertical variation of uSt
implies a depth-varying VF that leads to a simulation improvementcompared to the NA07 model (Section 4.8). The cross-shore undertow proﬁle is modiﬁed signiﬁcantly by the deep Bb (Fig. 4a–c). The
surf-zone recirculation in u greatly weakens in Runs 3 and 4. Exclusion of the roller shifts the bar recirculation in u towards the shoreward trough region, then it increases the surface onshore u near
the shoreward end. v(x, z) is modiﬁed less than u(x, z) among theseruns (Fig. 4d–f), but the deep breaking acceleration tends to generate the peak v location near the bar crest and deeper than the shallow breaking cases. In turn, the roller acts to shift the peak v
location shoreward, and v is mixed horizontally to reduce its
cross-shore gradient.
Simulated u(x, z) ﬁelds are compared with the observed velocities(Garcez Faria et al., 1998, 2000) in Fig. 5, and the normalized r.m.s.
errors for u and v (as deﬁned in NA07b) at a total of 42 measurement
positions are summarized in Table 2 (last two columns). The errors
for Run 1 in matching the observations are generally the least. The
deep breaking deﬁnitely lacks the recirculation pattern in u, withmuch weaker near-bed offshore undertow and near-surface onshore
ﬂow. All four 3D runs have fairly good agreement in v, while theexclusion of the roller clearly misses the increase of v in the troughregion. As a consequence, both the shallow breaker forcing and theroller shift the peak v location shoreward, and the former acts to generate the recirculating u ﬁeld quite well. The vertical structure of u
and the r.m.s. errors for Run 1 (uerror = 0.43 and verror = 0.092; Table2) are similar to or a bit even better than those in NA07 where uerror
and verror range 0.45–0.70, and 0.12–0.50, respectively; nevertheless, we view NA07 as a generally skillful model of the shallow-water
regime in DUCK94 (kh < 1). The Kv ﬁeld from the modiﬁed KPP model

26

Y. Uchiyama et al. / Ocean Modelling 34 (2010) 16–35(a) u: Run 2(b) u: Run 3(c) u: Run 4

u (m/s)−0.5 −0.25 0 0.25 0.5

u (m/s)−0.5 −0.25 0 0.25 0.5

u (m/s)−0.5 −0.25 0 0.25 0.5(d) v: Run 2(e) v: Run 3(f) v: Run 4

1
0

z (m)

−1
−2
−3
−4

1
0

z (m)

−1
−2
−3

v (m/s)−0.8−0.6−0.4−0.2

−4

0

50

v (m/s)−0.8−0.6−0.4−0.2

0100 150 200 250
distance offshore, x (m)

300 0

50

v (m/s)−0.8−0.6−0.4−0.2

0100 150 200 250
distance offshore, x (m)

300 0

50

0

100 150 200 250
distance offshore, x (m)300

Fig. 4. Distributions of u and v for shallow breaking without roller (Run 2; panels a and d); deep breaking with roller (Run 3; panels b and e); and deep breaking without roller(Run 4; panels c and f). Contour intervals are 0.05 m/s.

(a) u (m/s)(b) v (m/s)

0.1
0.2
0.3
0.4

0.2
0.4
0.6
0.8

1

z (m)

0
−1
−2
−3
−4

40

60

80
100
120
140
distance offshore, x (m)

160

180

40

60

obs.
Run 1
Run 2
Run 3
Run 4

80
100
120
140
distance offshore, x (m)

160

180

Fig. 5. Model-data comparison in Runs 1–4 for (a) u and (b) v. The vertical thin axes indicate the measurement locations.

 for different ar values (Runs 1, 2, and 5–7).
Fig. 6. Cross-shore proﬁles of (a) combined breaking and roller dissipation tot and (b) depth-averaged alongshore velocity v(Fig. 3) has wave-enhanced structures both near the surface and at
mid-depth that inﬂuence the currents (Section 3.4). It is qualitatively similar to the Kv ﬁeld in NA07 based on a two-equation local
turbulence closure model with a modiﬁcation due to wave breaking
as in Craig and Banner (1994).
4.3. Effects of waves
The roller model can be viewed as a correction to b (39) and Bb(48). We assess its inﬂuence by varying the value of ar in (48) in 5
cases: Run 1: ar = 1; Run 2: ar = 0; Run 5: ar = 0.25; Run 6: ar = 0.5;

and Run 7: ar = 0.75. As ar increases, the peak tot = (1  ar)b + rlocation moves slightly shoreward, and tot is noticeably intensiﬁed
 proﬁles gainin the trough region, 50 < x < 100 m (Fig. 6a). The v
more alongshore momentum in the trough with larger ar(Fig. 6b). The vertical shears in u are also enhanced with larger ar(Fig. 7a and b). This is due to an increase in the combined breaking
and bottom streaming force Awx (Fig. 7c). The changes in tot alter
the vertical eddy viscosity proﬁles (Fig. 7d), which are partly
responsible for the u proﬁles. Thus, we conﬁrm that a roller component greatly improves the match to DUCK94 (Garcez Faria
et al., 1998, 2000, NA07b).

27

Y. Uchiyama et al. / Ocean Modelling 34 (2010) 16–35

z + h (m)

(a)

(b)

(c)

3

3

3

2

2

2

1

1

1

(d)
3

2

αr = 0 (no roller)
α = 0.25r
α = 0.5r
α = 0.75r
α = 1.0

1r

no CEW
no WECobs.
0−0.3 −0.2 −0.1

0

0.1 0.2 0.3

0
−0.8

−0.6

u (m/s)

−0.4

−0.2

0
−0.03

0

0
−0.02

−0.01

wx

v (m/s)A000.005 0.01 0.015 0.02 0.025
2

2

Kv (m /s)(m/s )

Fig. 7. Model-data comparison in the trough region (x = 82 m) with different ar values (Runs 1, 2, and 5–7) and without CEW (Run 8) or WEC (Run 9): (a) u, (b) v, (c) the crossshore combined non-conservative breaking and bottom streaming forces Awx, and (d) Kv. The horizontal dotted lines are at z + h = D  H or z = fc  H.

z + h (m)

(a)

(b)

(c)

(d)

2.5

2.5

2.5

2.5

2

2

2

2

1.5

1.5

1.5

1

1

1

0.5

0.5

0.5

0
−0.6 −0.4 −0.20

0.2 0.4 0.6

u (m/s)

0
−0.8

−0.6

−0.4

−0.2

v (m/s)

0

0
−0.15

1.5
a = 0.1b
a = 0.2b
ab = 0.5
ab = 1.0

1

0.5

surface stress
Stokes scale
obs.

0
−0.1

−0.05

Bbx (m/s2)

0

0

0.01

0.02

0.03

2

Kv (m /s)Fig. 8. Model-data comparison varying Bb (Runs 1, 3, and 10–13) for (a) u and (b) v around the bar crest (x = 123 m); and modeled (c) Bbx and (d) K bv at x = 123 m. Horizontal
dotted lines are at z + h = D  H or z = fc  H.

In Fig. 7 we show two additional cases that artiﬁcially restrict
the wave–current interaction: Run 8 has no CEW in the WKB model, and Run 9 is entirely without WEC. Ignoring CEW does not havea large effect in this situation because the wave ﬁelds are not very
different among our cases (Section 4.1). However, without WEC,
both the recirculation in u and the ﬂow in v are very weak, while
Kv reverts to the KPP-evaluated K cd
v ; without wave–current interaction, the currents are entirely wind-driven. It is, of course, no surprise that WEC, especially due to Bb, is a primary inﬂuence on
littoral currents.

ﬂows in u are signiﬁcantly reduced, and the r.m.s. errors are muchlarger, while v is overestimated. The Stokes-scale fb in Run 3 gives
the worst agreement, and the shear in u is nearly absent. Verticaleddy viscosity Kv is substantially unaltered for all the cases(Fig. 8d). We conclude (consistent with NA07b) that a surface-concentrated Bb or an equivalent surface stress sbsur is essential toreproduce the surf-zone ﬂow structure. We further test the sensitivity to the types of fb(z) function in (53), but ﬁnd its inﬂuence onu(x, z) to be secondary to the choice of kb.
4.5. Breaking enhancement of vertical mixing

4.4. Depth-dependent breaking acceleration
The vertical scale of Bb(z) set by kb in (53) is crucial in the resultant surf-zone ﬂow structure. Extending the runs in Section 4.2, wetest 6 different settings for the breaking acceleration (Runs 1, 3,
and 10–13 in Table 2). Runs 1 and 10–12 have a shallow breaking
1force with a type III function in (53) with different kb ¼ ab H :ab = 0.1 (Run 10), 0.2 (Run 1), 0.5 (Run 11), and 1.0 (Run 12). Run
3 is for a deep Bb (i.e., kb = 2k). Run 13 speciﬁes the breaking force
as a surface stress as in (52).
For the smaller ab and surface stress cases, u and v around thebar crest (x = 123 m) are rather alike and show a good agreementwith the observations (Fig. 8a and b and the r.m.s. errors in Table2), regardless of the different Bbx proﬁles (Fig. 8c). However, withlarger ab, the surface onshore and near-bed offshore undertow

The wave breaking modiﬁcation to KPP relies on a choice of thevertical shape function fkV(z), its inverse scale akV, and the parameter cb (Section 3.4; (59)). We test their sensitivities by comparisonof Run 1 (type II shape function, akV = 1.2, cb = 0.03) with Runs 14–
15 (types I and III), Runs 16–18 (type II, akV = 0.4, 0.8, and 1.6), andRuns 19–21 and 5 (type II, cb = 0.01, 0.05, 0.07, and 0.1). Unlikewith fb, the diffusivity shape function has a noticeable inﬂuence,by underestimating v when type III is used and K bv exhibits less decay downward to the bed (Fig. 9). Similarly, as akV increases thenear-surface breaking effect deepens to intensify the near-bed
K bv , while the KPP K cd
v changes little among these cases. With smaller aK v the recirculating structure in u and magnitude in v are bothstrengthened. The empirical constant cb in (59) also has moderateinﬂuence on the resultant u ﬁeld (not shown). Its role is somewhat

28

Y. Uchiyama et al. / Ocean Modelling 34 (2010) 16–35

z + h (m)

(a)

(b)

(c)

(d)

2.5

2.5

2.5

2.5

2

2

2

2

1.5

1.5

1.5

1

1

1

0.5

0.5

0.5

a = 1.2 (I)Kv
aKv = 1.2 (II)a = 1.2 (III)Kv
aKv = 0.4 (II)aKv = 0.8 (II)a = 1.6 (II)

1.5

1

0.5

Kv

obs.

0
−0.6 −0.4 −0.20

0.2 0.4 0.6

0−0.8

u (m/s)

0
−0.6

−0.4

−0.2

0

0
0

0.02

v (m/s)

0.04

0.06

Kb (m2/s)

0.08

00.003 0.006 0.009 0.012

Kcd (m2/s)v

v

Fig. 9. Sensitivity to wave-induced mixing (Runs 1 and 14–18) around the bar crest (x = 123 m): model-data comparisons for (a) u and (b) v and modeled (c) K bv and (d) K cdv .
Roman numericals in parentheses in the legend in (c) indicate the shape types for fkV deﬁned in (53). Horizontal dotted lines are at D  H.

(a) x = 123 m
2.5(b) x = 700 m
7
6

2
z + h (m)

5
1.5

4
3

1

2
0.5

z + h (m)

1
0

0

0.5

0.5

0.4

0.4

0.3

0.3

0.2

0.2

0.1

0.1

0
−0.6 −0.4 −0.2

0 0.2 0.4 0.6
u (m/s)awd = 1.0
awd = 3.0
awd = 5.0
bottom stressno streaming
obs.

0−0.08 −0.04

0u (m/s)

0.04

0.08

Fig. 10. Model-data comparisons varying the bottom streaming (Runs 1 and 23–26) for u (a) around the bar crest at x = 123 m and (b) in the outer surf zone at x = 700 m. The
upper panels are the vertical proﬁles for the entire water column, while the lower panels are a zoom near the bed. Horizontal dotted lines indicate the e-folding scale
1
kwd ¼ awd dw in the streaming body force shape function fwd, with the particular value awd = 3.

similar to that by akV: decreasing cb leads to stronger recirculationin u and to speed in v. (Notice that the Run 1 values that best ﬁt
DUCK94 (i.e., akV = 1.2 and cb = 0.03) differ somewhat from the values suggested by Terray et al. (1996) and Apotsos et al. (2007) for
deep-water waves (akV = 0.7 and cb = 0.07).) We infer from the datacomparisons that there is a greater sensitivity to the shape proﬁlein K bv than in Bb, and that the wave-induced mixing scale is significantly larger than the analogous breaking acceleration scale (i.e.,
ab = 0.2).
4.6. Bottom streaming
Wave-induced bottom streaming (Section 3.3) has been lessinvestigated than the other WEC mechanisms because it occurs
within a thin wave bottom boundary layer (WBBL) that makes
measurement and modeling difﬁcult in the natural environment.
We test the sensitivity to the non-dimensional length scale awd

in (54)–(56) with Runs 23, 1 and 24 (awd = 1, 3 and 5, respectively),
as well as Run 25 where the WBBL is unresolved and imposed as abed stress as in (57) and Run 26 where streaming and bottom
wave-drag dissipation are neglected with wd = 0. The r.m.s. errors
in Table 2 demonstrate exclusion of the streaming leads to a modest increase of model error, while the other cases yield approximately the same errors; however, there are few measurements
near the bottom. The model-data comparison around the bar crest(Fig. 10a) shows that u is too strong without streaming because thenear-bottom offshore ﬂow is in the opposite direction to the onshore streaming stress. The streaming inﬂuence is even stronger
in the offshore region with a wind-driven Stokes–Ekman ﬂow(Fig. 10b); the streaming generates an onshore bottom velocity of
0.06 m/s and shifts the proﬁle of u over the entire water column.
Even in the offshore site, the sensitivity to the Bwd(z) (54) proﬁleis modest (including to the function type in fwd; not shown). The
WBBL thickness dw in (56) is estimated as 0.04 ± 0.01 m, which is29

Y. Uchiyama et al. / Ocean Modelling 34 (2010) 16–35

z + h (m)

(a)

(b)

2.5

2.5

2

2

1.5

1.5

1

1

0.5

0.5

S95, zo = 0.001 m
S95, zo = 0.005 m
S95, zo = 0.010 m
LIN, μ = 0.004 m/s
LIN, μ = 0.008 m/s
LOG, z = 0.001 m
o

obs.

0
−0.6 −0.4 −0.2 0 0.2 0.4 0.6
u (m/s)0−1 −0.8 −0.6 −0.4 −0.2
v (m/s)0

Fig. 11. Model-data comparisons for the sensitivity to wave-enhanced bottom drag (Runs 1 and 27–31) around the bar crest (x = 123 m): (a) u and (b) v. Horizontal dottedlines are at D  H.

only marginally resolved with a bottom-most grid height that varies from 0.03 m at the shoreward boundary to 0.1 m at the offshore
boundary. In other tests we have seen that marginal resolution, oreven an unresolved bottom stress, is sufﬁcient to capture the bulk
effect of streaming on u(z). In offshore regions wd typically dominates over b, and its effect is known to be signiﬁcant (Xu and Bowen, 1994; Lentz et al., 2008). It leads to a cross-shore bottom
velocity convergence offshore of the bar crest in our simulations,and its associated sediment transport may help to maintain thebar structure in the surf zone.
4.7. Wave-enhanced bottom dragBottom drag is well known to be the most important factor that
 (e.g., Uchiyama et al.,determines barotropic alongshore velocity v
2009). We examine its sensitivity on 3D u with alternative dragformulations: the combined wave–current model (60) withzo = 0.001 m (Run 1), 0.005 m (Run 27) and 0.01 m (Run 28); thelinear drag model (61) with l = 0.004 m/s (Run 29) and 0.008 m/
s (Run 30); and the log-layer model (62) with z0 = 0.001 m, omitting bed shear stress enhancement due to waves (Run 31). The vertical proﬁles of u around the bar crest (Fig. 11) vary substantiallyamong these cases, as expected from barotropic modeling experience. Runs 1 and 30 lead to approximately equivalent ﬂow ﬁelds
with good model skills (Table 2). Increasing zo with the Soulsby
model increases model error by weakening the recirculation in u
and reducing v. The linear drag model with small l = 0.004 andthe log-layer model yield good u proﬁles, but they have larger v
than measured. Other sensitivities related to bottom drag are that
Kv(z) is slightly altered through CEW on K bv and bed shear stress on
the KPP K cdv (not shown). Among the runs presented here, the Soulsby model with zo = 0.001 m provides the best overall agreementwith the data; the same conclusion was drawn with barotropic
 for DUCK94 (Feddersen et al., 1998). In summary,simulations of v
bottom drag also plays a crucial role in the 3D structure of u.
4.8. Horizontal momentum balance
To diagnose the inﬂuential mechanisms in our 3D solutions, weanalyze the momentum balances in Runs 1–3 (i.e., the case thatbest matches Duck94, a case removing the roller, and a caseincreasing the vertical scale of the breaking acceleration). The
advection and horizontal VF terms in (26) and (27) are rearranged
into non-ﬂux forms to separate the Eulerian advection (u  r)u andthe horizontal VF J in (5), which includes the Stokes–Coriolis force.
The vertical VF K in (5) is extracted from the second term in the
integral in (33) and added to the horizontal VF terms to combine

them into the total VF terms. The horizontal gradient of the restof (33) is the total pressure gradient force (PGF) Ptot. We
decompose the PGF into the non-WEC current contribution Pc;
the quasi-static response Pqs; the Bernoulli head contribution Pbh
from the interaction between vertical current shear and depthdependent Stokes drift; and the WEC surface pressure boundary
correction Ppc:

Ptot ¼ Pc þ Pwec ¼ Pc þ Pqs þ Pbh þ Ppc


Z z
gq
¼ r? gfc þ
dz þ g r?^f þ r? Kjfc þ r? Pjfc :
h

q0

ð63Þ

where Pwec = Pqs + Pbh + Ppc denotes the combined WEC contribution. All the momentum terms shown here are volume-averaged
and placed on the right side of (26) and (27).
The barotropic cross-shore momentum balance in Fig. 12 is led
by the pressure gradient (Ptot) and breaking forces. This is consistent with the classical view of surf-zone momentum balance
between wave-setup and breaking acceleration (cf., Bowen et al.,
1968; Raubenheimer et al., 2001). In the shallow breaking cases(Runs 1 and 2), the advection and the VF terms provide partiallycanceling cross-shore transports, whereas they are quite small inthe deep breaking case (Run 3). Note that the cross-shore horizontal
VF in Runs 1 and 2 is dominated by the vertical VF contribution. The
alongshore momentum balance is led by the breaking accelerationand the bottom drag, again consistent with the classical view. Thereis a secondary balance between the advection and VF as required by ﬂow
the barotropic mass balance that results in the anti-Stokes u
for an alongshore-uniform, steady circulation (Uchiyama et al.,2009). These secondary terms are again cross-shore transport,
and they have a larger reach between the bar and the nearshore
when the roller is present (Runs 1 and 3). The alongshore VF generally opposes y-advection where the latter is strong, but they are notcompletely canceled out, apparently because of different verticalstructures of u and uSt. The VF is also opposite to Bby < 0 aroundthe bar region and near the shoreline, while accelerating v < 0 in
the trough. The trough acceleration occurs in a more shoreward
location with shallow breaking and roller (Run 1) than otherwise(Runs 2–3). In addition to the broader tot in the trough with roller, occurs
the alongshore VF 3D shape is the reason why the peak vwell inshore of the bar in Run 1 with better agreement with themeasurements (Fig. 1e; Sections 4.1 and 4.2).
The depth-averaged cross-shore PGF decomposition (63) showsthat WEC contribution Pwec x is a signiﬁcant modiﬁcation of Ptot x,
although shallow breaking (Runs 1 or 2) has a larger Pwec x thandeep breaking (Fig. 13). Ptot x for Run 3 is primarily led by Pcx andthe quasi-static P qsx, corresponding to the classical barotropic

30

Y. Uchiyama et al. / Ocean Modelling 34 (2010) 16–35

× 10−23(b) Run 2: cross−shore(a) Run 1: cross−shore(c) Run 3: cross−shore
advectionvortex forcepressure
breakingbottom drag(m/s2)

2
1
0
−1
−2

× 10−3(d) Run 1: alongshore(e) Run 2: alongshore(f) Run 3: alongshore
advectionvortex forcepressure
breakingbottom drag2

(m/s2)

1
0
−1
−2
0

50

100 150 200 250
distance offshore, x (m)

300 0

50

100 150 200 250
distance offshore, x (m)

300 0

50

100 150 200 250
distance offshore, x (m)300

Fig. 12. Depth-averaged (top) cross-shore and (bottom) alongshore momentum tendencies for Runs 1 (baseline case), 2 (no roller), and 3 (deep B).

0.04(a) Run 1(b) Run 2(c) Run 3
bh x

P
qs x
P
Ppc x
cx
P
tot x
P

0.03

2

(m/s )

0.02
0.01
0
−0.01
0

50

100 150 200 250
distance offshore, x (m)

300 0

50

100 150 200 250
distance offshore, x (m)

300 0

50

100 150 200 250
distance offshore, x (m)300

Fig. 13. Depth-averaged cross-shore PGF decomposed as in (63) for Runs 1–3.

wx

1(a) x: A2

2(m/s )wx

wx(b) x: advection (m/s )(c) x: D−A−2 −1 0 1 2−2 −1 0 1 2

2(m/s )

2(d) x: vortex force (m/s )

z (m)

0
−1
−2
−3 −2 −1 0 1 2−4
× 10−22

1(e) x: pressure (m/s )× 10−22−2 −1 0 1 2
× 10−2wy(f) y: advection (m/s )(g) y: D−2 −1 0 1 2−2 −1 0 1 2

wy−A

× 10−22(m/s )2(h) y: vortex force (m/s )

z (m)

0
−1
−2
−3 −2 −1 0 1 2
−2
−4
× 10
0

50 100 150 200 250 300 0
distance offshore, x (m)

−3

× 10

50 100 150 200 250 300 0
distance offshore, x (m)−2 −1 0 1 2
−3

× 10

50 100 150 200 250 300 0
distance offshore, x (m)

−3

× 10

50 100 150 200 250 300
distance offshore, x (m)Fig. 14. Leading terms in the 3D horizontal momentum balance for Run 1: (a) x wave acceleration Awx, (b) x Eulerian advection, (c) x vertical mixing Dwx, (d) x total VF, (e) x
PGF, (f) y Eulerian advection, (g) y vertical mixing Dwy, and (h) y total VF.

31

Y. Uchiyama et al. / Ocean Modelling 34 (2010) 16–35
wx

1(a) x: A2

2(m/s )wx

wx(b) x: advection (m/s )(c) x: D−A−2 −1 0 1 2−2 −1 0 1 2

2(m/s )

2(d) x: vortex force (m/s )

z (m)

0
−1
−2
−3 −2 −1 0 1 2−4
× 10−22

1(e) x: pressure (m/s )× 10−22−2 −1 0 1 2
× 10−2wy(f) y: advection (m/s )(g) y: D−2 −1 0 1 2−2 −1 0 1 2

wy−A

× 10−22(m/s )2(h) y: vortex force (m/s )

z (m)

0
−1
−2
−3 −2 −1 0 1 2−4
× 10−2050 100 150 200 250 300 0
distance offshore, x (m)× 10−3−2 −1 0 1 2
× 10−350 100 150 200 250 300 0
distance offshore, x (m)50 100 150 200 250 300 0
distance offshore, x (m)× 10−350 100 150 200 250 300
distance offshore, x (m)Fig. 15. Same as Fig. 14, but for Run 3.

cross-shore momentum balance (cf., Bowen et al., 1968; Uchiyama
et al., 2009). For Runs 1 and 2, Pwec x is dominated by Pbhx and Pqsx;
i.e., Ptot x is modiﬁed signiﬁcantly by the vertical current shear
through Pbhx. This mechanism causes the difference in the wave-induced sea-level rise (setup) between the 2D and 3D cases in Fig. 1c.
The WEC surface pressure boundary correction Ppcx is small for allcases here.
Leading terms in the 3D cross-shore momentum balance forboth Runs 1 and 3 (Figs. 14a–e and 15a–e) show the central roleof vertical mixing D wx in the bar and shore regions. It verticallyconnects the surface-intensiﬁed breaking to the bottom drag inthe cross-shore force balance (in combination with a nearly barotropic PGF). The increase of D wx < 0 with height and the bottomintensiﬁed structure of cross-shore advection and VF induce theu recirculation in Fig. 2a; however, the features are much weakerwith deep breaking (Run 3), as is the recirculation itself (Fig. 4b).
The WEC role in the alongshore momentum balance (Fig. 14f–h)
is even more striking because VF has a signiﬁcant decrease withdepth in balancing the baroclinic distributions of both advection
and vertical mixing (whereas VF merely balances advection in thedepth-averaged budget; Fig. 12). With shallow breaking (Run 1),
advection and vertical mixing primarily redistribute v vertically,
while VF is mainly a horizontal redistribution. Even with deep
breaking (Run 3), VF has signiﬁcant vertical variation (Fig. 15h). Thisimplies that in shallow littoral regions like DUCK94, it is necessary tohave a fully 3D structure for Stokes drift and VF, unlike in the NA07
model. The primary alongshore momentum balance occurs between
advection and vertical mixing in Run 1, particularly over the bar andnear the shoreward boundary, where the recirculating u is mostprominent, and VF enters at the same order of magnitude. In Run 3
VF contributes relatively more to balance the advection, while thevertical mixing plays a secondary role (Fig. 15f–h). The alongshore
non-conservative forcing Awy (not shown) is proportional to Awx by
deﬁnition, hence always negative to drive the negative v.
In summary, WEC is quantitatively important for u(x, z) in thesurf zone momentum balance through both VF and especially Pwec.
A similar conclusion is in NA07 regarding VF, although its verticalstructure is neglected there.
5. Model comparisons for a plane beach
Another test case is a littoral ﬂow driven by obliquely incident
waves on a plane beach with a uniform slope of 1:80. This problem
was posed in HW09 to comparatively assess a quasi-3D model(SHORECIRC) and a fully-3D model (developed in W08 using

ROMS), both based upon a radiation-stress formalism. The W08
model is based on a set of wave-averaged equations derived by
Mellor (2003, 2005) using a GLM-like vertical mapping approachwith a depth-varying radiation stresses. Here we compare the latter solution with ones using the present ROMS model with identical wave ﬁelds obtained using SWAN without CEW. As described in
HW09, the offshore spectral waves are speciﬁed by the JONSWAP
spectrum for 2 m signiﬁcant wave height Hsig with a peak periodof 10 s at an angle of 10° off the shore-normal direction. A quasibarotropic quadratic bottom drag model is used:scd
bot ¼ q0 c D juju;

ð64Þwith cD = 0.0015. For 3D runs u in (64) is the horizontal velocity at the in a 2D barotropic run. Therebottom-most grid cells, whereas it is u
are minor functional differences between the two models – the W08
model has a wetting/drying capability (whereas our model imposes asolid wall boundary with the minimum water depth, hmin = 0.01 m,
at the shoreward boundary), and it relies on a GLS k   turbulence
closure model (Warner et al., 2005) (whereas our model uses either
the modiﬁed KPP or the analytical model in (65); however, these do
not dominate the solution differences we present).
The horizontal extent of the domain is 1180 m in x (crossshore)  140 m in y (alongshore) with grid spacings of Dx = Dy =20 m. The model coordinates have a west-coast orientation, withthe offshore open boundary located at x = 0. The resting depth hvaries from 12 m to 0.01 m (hmin), and is discretized with 20 uniform vertical s levels. Boundary conditions are alongshore periodicity, zero normal ﬂux and tangential Neumann conditions at theshoreline boundary, and Chapman-type radiation with weak nudg and fc (as in Section 4) at the offshore boundary. Rotation
ing for u
is excluded with f = 0. There is no lateral momentum diffusion,
stratiﬁcation, nor surface wind/heat/freshwater ﬂuxes. Roller
waves and bottom streaming are excluded. We conduct four simulations with the present model: a 2D barotropic case (Run a), andthree 3D cases with breaking acceleration Bb (48) using a type III
shape function (53) with either kb = (0.2H)1 (Run b), kb = 2k
(Run c), or the same as Run c but with a parabolic proﬁle forvertical mixing (Run d). The last run mimics the eddy viscosity(Kv) distribution in the W08 model, viz.,hþzK v ðzÞ ¼ 0:011ðh þ zÞ 1:0 
:
D

ð65Þ

In our modiﬁed KPP, we use a type II shape function with kkV =(1.2H)1 and cb = 0.03 for Runs b and c. For each run the time integration is continued until a steady solution occurs.

32

Y. Uchiyama et al. / Ocean Modelling 34 (2010) 16–35

 ; and (d) v
 . Runs a, b, and d are compared
Fig. 16. Plane-beach cross-shore proﬁles of (a) depth h, signiﬁcant wave height Hsig, and breaking dissipation rate b/q0; (b) fc; (c) uwith the HW09 solution. Plus marks depict the analytical solutions (66)–(68).

5.1. Barotropic ﬁelds
Uchiyama et al. (2009) show for steady, alongshore-uniform,
unstratiﬁed, non-rotating solutions that the barotropic continuity
balance can be integrated in x from the shoreline to yield the offshore-directed, anti-Stokes cross-shore transport velocity, ¼ uSt :uð66Þ

 is obtained from the mean alongshore
An analytical solution for the v
momentum balance, where (66) implies separate sub-balances ox v ¼ uSt
between the VF and alongshore advection (i.e., u Þ) and between bottom drag and breaking acceleration:ð^z  r?  ub k
q0 cD ju jv ¼ y ;r

ð67Þ

pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ
j ¼ u 2 þ v 2 . The cross-shore barotropic momentum balwhere ju
ance is dominated by PGF and breaking acceleration:



b kx
gox fc  ^f :
q0 Dr

with DUCK94 results (Sections. 4.4 and 4.8). In the breaking zone in coincides with that for the W08 model. The W08 modelRun d, v
 opposing the incident wave direction
produces a small positive v
in the offshore (x < 500 m). Run c has a sufﬁciently similar answer
to Run d that we do not show it here, indicating an insensitivity tothe particular Kv(z). All the cases shown here generally agree wellwith the analytical depth-averaged solutions, although HW09 has
a slightly larger mismatch and a non-zero value offshore (in the
direction opposite to the incident waves), again probably a boundary condition artifact.
In the present model, the vertically integrated WEC accelerations by breaking and the quasi-static PGF may be combined as

Btot ¼

Z fc

Bb dz þ gDr?^f:

In W08 the 3D horizontal radiation-stress divergence S is

S ¼ ðSx ; Sy Þ ¼




ð68Þ

By integrating (68) from the offshore boundary with fc jx¼0 ¼ ^f(assuming bjx=0 = 0), the fc(x) proﬁle is approximately retrieved(e.g., Longuet-Higgins and Stewart, 1962; Bowen et al., 1968). With(66)–(68) and (7) plus the SWAN-produced wave data for b, k, and ðxÞ and fc(x).
A, we can readily calculate analytical solutions for u
The shoreward decrease in Hsig occurs because of b, whilerefractive wave shoaling has a weak inﬂuence (Fig. 16a). Wave
breaking occurs around 500 < x < 1000 m with the peak nearx = 700 m. The ROMS fc in Fig. 16b generally agrees well with theanalytical fc from (68). A minor deviation arises in the case witha shallow breaking force (Run b) near the shoreline by Bernoulli
head effect (as in Section 4.8). A slightly larger deviation in fc isalso found with the W08 model; because it extends far offshore,
we suspect it is an artifact from the open boundary condition,
but it causes only a modest discrepancy compared to the 3D differ in Fig. 16c is in almost perfectences emphasized in Section 5.2. u also corresponds to theagreement with (66) for all the cases. v
analytical solution (67) fairly well, particularly for Run a (2D)
(Fig. 16d) while the 3D cases show deviations because of the 3D peaks at almost the same
current in the bottom drag force. v
cross-shore location for all the cases, except for Run b with a more
 j maximum because of a cross-shore momentum imbalonshore jv
ance associated with the Aw and vertical mixing Dw (47), consistent

ð69Þ

h

oSxx oSyx oSpx
oSxy oSyy oSpy

þ
;

þ
ox
oy
oz
ox
oy
oz


ð70Þ

(Mellor, c 2003; Warner et al., 2008), with an associatedRf
Btot ¼ h S dz. Fig. 17 shows that the Btotx proﬁles are essentially
identical between the two models and unchanged among our different cases. There is some small difference in Btoty(x) (Fig. 17b), where
the HW09 proﬁle has a slight positive bias associated with the artiﬁcial positive v in the offshore region (x < 500 m; Fig. 16d).
5.2. Vertical shear
In the surf zone, the shallow breaking case (Run b) has thestrongest recirculation in u(x, z), with near-surface onshore ﬂow
and near-bed offshore undertow (Fig. 18). This pattern is similar
to that in the DUCK94 breaking regions (Fig. 4). The other 3D cases(Run d and HW09) induce much weaker recirculation in association with deeper breaker acceleration proﬁles. In contrast, Run b
has the weakest v(x, z), while Run d and HW09 have similar v.
The modiﬁed KPP scheme concentrates Kv(x, z) near the surface inthe breaking zone (Run b), while the other cases have a mid-depth
maximum for Kv that increases with depth offshore (Fig. 18: lowermiddle). Because Runs b and c (not shown) are more similar in u
and v than either is with Run d and HW09, we conclude that themost important distinction is the vertical structure of B, with Kv
providing a lesser distinction (cf., Section 4.5). All models yield ðxÞ (Fig. 16), but the W08 model generates
the same anti-Stokes u

33

Y. Uchiyama et al. / Ocean Modelling 34 (2010) 16–35(a) Btot x

−2

× 10

0

1.5
(m2/s2)(b) Btot y

−3

× 10

1

−0.5

0.5

Run a (2D)
Run b
Run d
Haas & Warner (2009)

−1

0
−0.5
0

200

400

600

800

−1.5
0

1000

200

x (m)

400

600

800

1000

x (m)

Fig. 17. Plane-beach proﬁles of (a) cross-shore and (b) alongshore components of Btot for all four runs. In addition, the gray dashed line in (a) is the quasi-static PGF, gDox ^f,for Run a.

(a) Run b(b) Run d(c) Haas & Warner (2009)

0

z (m)

−2
−4
−6
u (m/s)−8−0.3 −0.2 −0.1

−10

u (m/s)0u (m/s)0.1

0.2−0.3 −0.2 −0.1

0

0

0.3

−1.2 −0.9 −0.6 −0.3

0.1

0.2

−0.3 −0.2 −0.1

0

0

0.3

−1.2 −0.9 −0.6 −0.3

0.1

0.2

0

0.3

−12
0

z (m)

−2
−4
−6
−8

v (m/s)−10−1.2 −0.9 −0.6 −0.3

v (m/s)v (m/s)

−12
0

z (m)

−2
−4
−6
2

2

Kv (m /s)

−8
0

−10

1

2

3

4

−2

0

1

Kv (m /s)

2

3

4

−2

× 10

−12

2

Kv (m /s)
0

1

2

3

4

y

2

−0.4

−0.1

0.2

400
600
x (m)

800

−2

× 10

× 10

0

z (m)

−2
−4
−6
by

B

−8
−1

−10

−0.7

2

by

(m/s )

−0.4

−0.1

B
0.2

× 10−3

−12
0

200

−1

−0.7

2

(m/s )−0.4−0.1

S (m/s )0.2

× 10−3
400
600
x (m)

800

1000

200

−1

−0.7

× 10−3
400
600
x (m)

800

1000

200

1000

Fig. 18. Model comparisons among Runs b and d for the present model and the W08 model (in HW09): (top) u; (upper-middle) v; (lower-middle) Kv; and (bottom) Bb (Runs b
and d) or S (HW09).

an onshore surface u(x, z) even in the offshore region outside thebreaking zone, which is not seen with the present model. (HW09
shows that the SHORECIRC model also yields a less sheared u proﬁle in the offshore region than the W08 model.)In the W08 model, the radiation stress tensor (70), e.g., Sxy, canbe rearranged as

Sxy ¼

E kx ky cosh 2kðh þ zÞ þ 1:sinh 2kD

q0 k

ð71Þ

This 3D radiation stress has a depth-dependency of cosh[2k(h + z)]
function, consistent with our type III vertical shape function fb(z)for Bb with kb = 2k in (53). This becomes nearly constant in shallow

water (kh ? 0), including over much of the surf zone. Thus, S in the
W08 is vertically homogeneous (Fig. 18: bottom panels), in contrast
to the shallow Bb in Run b. In Run d, Bb resembles S. Therefore Runsb and d (and HW09) are closely related to Runs 1 and 3 in the DUCK
experiment (Section 4, Fig. 5) in terms of the depth-dependency
parameter choice. Accordingly, it is anticipated that the 3D radiation-stress model (e.g., Mellor, 2003) could have noticeable deﬁciency for surf-zone applications. In addition, a more verticallysheared velocity ﬁeld with the appropriate Bb forcing is essential
to more correct VF representation that leads to signiﬁcant modiﬁcation in the momentum balance as seen in Section 4.8, particularly
through the horizontal and vertical VF and the Bernoulli-headpressure force K. The radiation-stress divergence contains multiple

34

Y. Uchiyama et al. / Ocean Modelling 34 (2010) 16–35

aspects of WEC: the conservative VF (if the accompanying wave
model takes CEW into account appropriately; Lane et al., 2007),the conservative gradient of the quasi-static PGF g r?^f (as part of
Sxx and Syy; (e.g., Longuet-Higgins and Stewart, 1964)), and thenon-conservative acceleration Bb. When S is evaluated with the vertical structure of the leading-order primary wave solution (e.g.,
Mellor, 2003), then it causes an underestimation of vertical shear
in u. The present model cleanly separates the different WEC inﬂuences and allows them to have different spatial distributions.
6. Summary and prospects
Wave–current interaction with a vortex-force (VF) formalism isimplemented in a fully-3D oceanic circulation model (ROMS) intended for use in a wide range of conditions. The Eulerian waveaveraged, multi-scale asymptotic theory by MRL04 is adapted tobe appropriate for ROMS. Conservative wave effects include theVF, the Stokes–Coriolis force, Bernoulli head, and quasi-static pressure and sea-level responses. Non-conservative accelerations areincluded through parameterizations of depth-induced breaking,associated surface rollers, and wave-induced streaming dissipation
at the bottom. Wave-enhanced mixing from surface breakers is included (adapted from Battjes (1975)), as is wave-enhanced bottom
stress and bottom-boundary layer mixing in a KPP parameterization (Large et al., 1994). Here the model is applied to the surf zone
off Duck, NC and demonstrates a good agreement with in situ
velocity data for appropriate choices of several model parameters.
The model is further compared to another ROMS-based 3D modelwith a depth-dependent radiation-stress formalism (Warner
et al., 2008) for a littoral current on a gently-sloping plane beach.
Littoral currents are caused by depth-induced wave breaking in ðxÞ
the problems solved here. In the DUCK94 problem, alongshore v
is largest near the topographic bar and has a momentum balance of
breaking acceleration against deceleration by bottom drag and vor ðxÞ is an offshore, anti-Stokes transtex force, while cross-shore u
port locally enhanced near the bar and shoreline by breaking andshifted shoreward by VF. Surface-intensiﬁed breaking (on ascale < H) by the primary and roller waves Bb is essential to reproduce the measured current proﬁles. v(x, z) peaks shoreward of thebar and has a modest degree of vertical shear, while u(x, z) has twostrong recirculations (onshore at surface, offshore at depth) near
the bar and shoreline. Wave effects of Bb, vertical mixing, PGF,
and vortex force all contribute to the maintenance of the current
proﬁles. Offshore of the breaking region, the wave-induced bottom
streaming stress shifts the maximum of the anti-Stokes u(z) > 0 upward to mid-depth. Similar conclusions obtain in the plane beachproblem, where in particular a previously proposed deep radiation-stress representation greatly underestimates the recirculationin u(x, z) compared to WEC with shallow Bb, K bv , and vortex force.
The two applications presented here are of limited generality
due to the weakness of the effect of currents on the waves (CEW),
the absence of strong alongshore variation (e.g., rip currents), density stratiﬁcation, interaction with eddying currents, suspended
sediments, and even non-hydrostatic current dynamics. We anticipate that signiﬁcant additional wave–current interaction phenomena will be abundant in these variously more general regimes.
Acknowledgments
We thank John Warner for his comments and for sharing themodel data used in Section 5. Thanks are due to those who carriedout the DUCK94 ﬁeld experiments and provided online data access.
We thank members of the NOPP-CSTMS project and three anonymous reviewers for their constructive suggestions. This work is
ﬁnancially supported by the National Science Foundation (DMS-

0723757) and Ofﬁce of Naval Research
N00014-06-1-0945, N00014-08-1-0597).

(N00014-04-1-0166,

References
Andrews, D.G., McIntyre, M.E., 1978a. An exact theory of nonlinear waves on a
Lagrangian-mean ﬂow. J. Fluid Mech. 89, 609–646.
Andrews, D.G., McIntyre, M.E., 1978b. On wave action and its relatives. J. Fluid
Mech. 89, 647–664.
Apotsos, A., Raubenheimer, B., Elgar, S., Guza, R.T., Smith, J.A., 2007. Effects of wave
rollers and bottom stress on wave setup. J. Geophys. Res. 112, C02003.
doi:10.1029/2006JC003549.
Ardhuin, F., Jenkins, A.D., Belibassakis, K.A., 2008a. Comments on ‘‘The threedimensional current and surface wave equations”. J. Phys. Oceanogr. 38, 1340–
1350.
Ardhuin, F., Rascle, N., Belibassakis, K.A., 2008b. Explicit wave-averaged primitive
equations using a generalized Lagrangian mean. Ocean Modell. 20, 35–60.
Battjes, J.A., 1975. Modeling of turbulence in the surfzone. In: Proceedings of the
Symposium on Modeling Techniques 1975 San Francisco, CA, American Societyof Civil Engineers, New York, pp. 1050–1063.
Blaas, M., Dong, C., Marchesiello, P., McWilliams, J., Stolzenbach, K., 2007. Sedimenttransport modeling on Southern Californian shelves: a ROMS case study. Contin.
Shelf Res. 27, 832–853.
Blumberg, A.F., Mellor, G.L., 1987. A description of a three-dimensional coastal
ocean circulation model. In: Three-Dimensional Coastal Ocean Models,American Geophysical Union, pp. 1–16.
Booij, N., Ris, R.C., Holthuijsen, L.H., 1999. A third generation wave model for coastal
regions 1. Model description and validation. J. Geophys. Res. 104 (C4), 7649–
7666.
Bowen, A.J., Inman, D.L., Simmons, V.P., 1968. Wave ‘set-down’ and wave setup. J.
Geophys. Res. 73, 2569–2577.
Burchard, H., 2001. Simulating the wave-enhanced layer under breaking surface
waves with two-equation turbulence models. J. Phys. Oceanogr. 31, 3133–3145.
Chen, Q., Dalrymple, R.A., Kirby, J.T., Kennedy, A.B., Haller, M.C., 1999. Boussinesq
modeling of a rip current system. J. Geophys. Res. 104, 20617–20637.
Church, J.C., Thornton, E.B., 1993. Effects of breaking wave induced turbulence
within a longshore current model. Coastal Eng. 20, 1–28.
Craig, P.D., Banner, M.L., 1994. Modeling wave-enhanced turbulence in the ocean
surface layer. J. Phys. Oceanogr. 24, 2546–2559.
Craik, A.D.D., Leibovich, S., 1976. A rational model for Langmuir circulations. J. Fluid
Mech. 73, 401–426.
Davies, A.G., Villaret, C., 1999. Eulerian drift induced by progressive waves above
rippled and very rough beds. J. Geophys. Res. 104 (C1), 1465–1488.
Dingemans, M.W., Radder, A.C., Vriend, H.J.D., 1987. Computation of the driving
forces of wave-induced currents. Coastal Eng. 11, 539–563.
Durski, S.M., Glenn, S.M., Haidvogel, D., 2004. Vertical mixing schemes in the coastal
ocean: comparision of the level 2.5 Mellor–Yamada scheme with an enhanced
version of the K proﬁle parameterization. J. Geophys. Res. 109, C01015.
doi:10.1029/2002JC001702.
Elgar, S., Guza, R.T., Raubenheimer, B., Herbers, T.H.C., Gallagher, E.L., 1998. Spectral
evolution of shoaling and breaking waves on a barred beach. J. Geophys. Res.
103, 15797–15805.
Feddersen, F., Guza, R.T., Elgar, S., Herbers, T.H.C., 1998. Alongshore momentum
balances in the nearshore. J. Geophys. Res. 103, 15667–15676.
Fredsøe, J., Deigaard, R., 1995. Mechanics of Coastal Sediment Transport. World
Scientiﬁc, Singapore.
Gallagher, E.L., Boyd, W., Elgar, S., Guza, R.T., Woodward, B., 1996. Perfomance of asonar altimeter in the nearshore. Mar. Geol. 133, 241–248.
Gallagher, E.L., Elgar, S., Guza, R.T., 1998. Observations of sand bar evolution on anatural beach. J. Geophys. Res. 103, 3203–3215.
Galperin, B., Kantha, L.H., Hassid, S., Rossati, A., 1988. A quasi-equilibrium turbulent
energy model for geophysical ﬂows. J. Atmos. Sci. 45, 55–62.
Garcez Faria, A.F., Thornton, E.B., Stanton, T.P., Soares, C.V., Lippmann, T.C., 1998.
Vertical proﬁles of longshore currents and related bed stress and bottomroughness. J. Geophys. Res. 103, 15667–15676.
Garcez Faria, A.F., Thornton, E.B., Lippmann, T.C., Stanton, T.P., 2000. Undertow overa barred beach. J. Geophys. Res. 105, 16999–17010.
Garrett, C., 1976. Generation of Langmuir circulations by surface waves – a feedback
mechanism. J. Mar. Res. 34, 116–130.
Groeneweg, J., 1999. Wave–current interaction in a generalized Lagrangian mean
formulation. Ph.D. Thesis, Delft University of Technology, Delft, The Netherland.
Haas, K.A., Warner, J.C., 2009. Comparing a quasi-3D to a full 3D nearshore
circulation model: SHORECIRC and ROMS. Ocean Modell. 26, 91–103.
Hasselmann, K., 1971. On the mass and momentum transfer between short gravity
waves and larger-scale motions. J. Fluid Mech. 50, 189–201.
Jones, N.L., Monismith, S.G., 2008. Modeling the inﬂuence of wave-enhanced
turbulence in a shallow tide- and wind-driven water column. J. Geophys. Res.
113, C03009. doi:10.1029/2007JC004246.
Kanarska, Y., Shchepetkin, A.F., McWilliams, J.C., 2007. Algorithm for nonhydrostatic dynamics in the Regional Oceanic Modeling System. Ocean
Modell. 18, 143–174.
Klopman, G., 1994. Stokes transport. WL-Delft Hydraulics Report H840.30, Part II.
Lane, E.M., Restrepo, J.M., McWilliams, J.C., 2007. Wave–current interaction: acomparison of radiation-stress and vortex-force representations. J. Phys.
Oceanogr. 37, 1122–1141.

Y. Uchiyama et al. / Ocean Modelling 34 (2010) 16–35
Large, W.G., McWilliams, J.C., Doney, S.C., 1994. Oceanic vertical mixing: a review
and a model with nonlocal boundary layer parameterisation. Rev. Geophys. 32,
363–403.
Leibovich, S., 1980. On wave–current interaction theories of Langmuir circulations.
J. Fluid Mech. 99, 715–724.
Lentz, S.J., Fewings, M., Howd, P., Fredericks, J., Hathaway, K., 2008. Observations
and a model of undertow over the inner continental shelf. J. Phys. Oceanogr. 38,
2341–2357.
Lesser, G.R., Roelvink, J.A., van Kester, J.A.T.M., Stelling, G.S., 2004. Development andvalidation of a three-dimensional morphological model. Coastal Eng. 51, 883–
915.
Lewis, J.K., 1997. A three-dimensional ocean circulation model with wave effects.
In: Spaulding, M.L., Blumberg, A.F. (Eds.), Proceedings of the Fifth International
Conference on Estuarine and Coastal Modeling, Alexandria, VA. American
Society of Civil Engineers, Reston, VA, pp. 584–600.
Long, C.E., 1996. Index and bulk parameters for frequency-direction spectra
measured at CERC Field Research Facility, June 1994 to August 1995. US Army
Corps of Engineers Waterway Experiment Station, Vicksburg, MI, USA.
Longuet-Higgins, M.S., 1953. Mass transport in water waves. Philos. Trans. R. Soc.
Lond., Ser. A 245, 535–581.
Longuet-Higgins, M.S., 1958. The mechanics of the boundary later near the bottom
in a progressive wave. Appendix to ‘‘An experimental investigation of drift
proﬁles ub a closed channel” by R.C.H Russel and J.D.C. Osorio. In: Proceedingsof the Sixth Coastal Engineering International Conference 1957, Florida,
American Society of Civil Engineers, New York, pp. 184–193.
Longuet-Higgins, M.S., 1970. Longshore currents generated by obliquely incident
sea waves, 1 & 2. J. Geophys. Res. 75 (33), 6778–6801.
Longuet-Higgins, M.S., 1973. The mechanics of the surfzone. In: Becker, E.,
Mikhailov, G.K. (Eds.), Proceedings of the 13th International Congress ofTheoretical and Applied Mathematics 1972, Moscow, USSR. Springer, Berlin,
pp. 213–228.
Longuet-Higgins, M.S., Stewart, R.W., 1962. Radiation stress and mass transport ingravity waves, with application to ‘surf beats’. J. Fluid Mech. 13, 481–504.
Longuet-Higgins, M.S., Stewart, R.W., 1964. Radiation stresses in water waves: aphysical discussion, with applications. Deep-Sea Res. 11, 529–562.
McWilliams, J.C., Restrepo, J.M., 1999. The wave-driven ocean circulation. J. Phys.
Oceanogr. 29, 2523–2540.
McWilliams, J.C., Sullivan, P.P., Moeng, C.-H., 1997. Langmuir turbulence in theocean. J. Fluid Mech. 334, 1–30.
McWilliams, J.C., Restrepo, J.M., Lane, E.M., 2004. An asymptotic theory for theinteraction of waves and currents in coastal waters. J. Fluid Mech. 511, 135–
178.
Mellor, G.L., 2003. The three-dimensional current and surface wave equations. J.
Phys. Oceanogr. 33, 1978–1989.
Mellor, G.L., 2005. Some consequences of the three-dimensional currents andsurface wave equations. J. Phys. Oceanogr. 35, 2291–2298.
Mellor, G.L., 2008. The depth-dependent current and wave interaction equations: arevision. J. Phys. Oceanogr. 38, 2587–2596.
Mellor, G.L., Yamada, T., 1982. Development of a turbulent closure model for
geophysical ﬂuid problems. Rev. Geophys. 20, 851–875.
Nairn, R.B., Roelvink, J.A., Southgate, H.N., 1991. Transition zone width andimplications for modelling surfzone hydrodynamics. In: Edge, B.L. (Ed.),
Proceedings of the 22nd Coastal Engineering International Conference 1990,
Delft, The Netherlands. American Society of Civil Engineers, New York, pp. 68–
82.
Newberger, P.A., Allen, J.S., 2007a. Forcing a three-dimensional, hydrostatic,
primitive-equation model for application in the surf zone: 1. Formulation. J.
Geophys. Res. 112, C08018. doi:10.1029/2006JC003472.
Newberger, P.A., Allen, J.S., 2007b. Forcing a three-dimensional, hydrostatic,
primitive-equation model for application in the surf zone: 2. Application to
DUCK94. J. Geophys. Res. 112, C08019. doi:10.1029/2006JC003474.
Okayasu, A., Shibayama, T., Mimura, N., 1986. Velocity ﬁeld under plunging waves.
In: Edge, B.L. (Ed.), Proceedings of the 20th Coastal Engineering International
Conference 1986, Taipei, Taiwan. American Society of Civil Engineers, New York,
pp. 660–674.
Özkan-Haller, H.T., Kirby, J.T., 1999. Nonlinear evolution of shear instabilities of thelongshore current: a comparison of observations and computations. J. Geophys.
Res. 104, 25953–25984.
Özkan-Haller, H.T., Li, Y., 2003. Effects of wave–current interaction on shear
instabilities of longshore currents. J. Geophys. Res. 108 (C5), 3139. doi:10.1029/
2001JC001287.
Perrie, W., Tang, C.L., Hu, Y., DeTracy, B.M., 2003. The impact of waves on surface
currents. J. Phys. Oceanogr. 33, 2126–2140.
Phillips, O.M., 1977. The Dynamics of the Upper Ocean. Cambridge University Press,
Cambridge.
Rascle, N., 2007. Impact of waves on the ocean circulation. Ph.D. Thesis, Universitéde Bretagne Occidentale – Brest, Brest, France.
Rascle, N., Ardhuin, F., Terray, E.A., 2006. Drift and mixing under the ocean surface: acoherent one-dimensional description with application to unstratiﬁed
conditions. J. Geophys. Res. 111, C03016. doi:10.1029/2005JC003004.
Rascle, N., Ardhuin, F., Terray, E.A., 2009. Drift and mixing under the ocean surface
revisited: stratiﬁed conditions and model-data comparisons. J. Geophys. Res.
114, C02016. doi:10.1029/2007JC004466.
Raubenheimer, B., Guza, R.T., Elgar, S., 2001. Field observation of wave-driven
setdown and setup. J. Geophyis. Res. 106 (C3), 4629–4638.

35

Reniers, A.J.H.M., Roelvink, J.A., Thornton, E.B., 2004a. Morphodynamic modeling ofan embayed beach under wave group forcing. J. Geophys. Res. 109, C01030.
doi:10.1029/2002JC001586.
Reniers, A.J.H.M., Thornton, E.B., Stanton, T.P., Roelvink, J.A., 2004b. Vertical ﬂow
structure during sandy duck: observations and modeling. Coastal Eng. 51, 237–
260.
Ruessink, B., Miles, J., Feddersen, F., Guza, R., Elgar, S., 2001. Modeling thealongshore current on barred beach. J. Geophyis. Res. 106 (C10), 22451–22463.
Shchepetkin, A.F., McWilliams, J.C., 2005. The Regional Oceanic Modeling System: asplit-explicit, free-surface, topography-following-coordinate oceanic model.
Ocean Modell. 9, 347–404.
Shchepetkin, A.F., McWilliams, J.C., 2008. Computational kernel algorithms for ﬁnescale, multiprocess, longtime oceanic simulations. In: Temam, R., Tribbia, J.
(Eds.), Handbook of Numerical Analysis: Computational Methods for the Ocean
and the Atmosphere. Elsevier, Amsterdam, pp. 119–181.
Shchepetkin, A.F., McWilliams, J.C., 2010. KPP revisited, in preparation.
Shi, F., Kirby, J.T., Haas, K., 2006. Quasi-3D nearshore circulation equations: a CLvortex force formulation. In: Proceedings of the 30th International Conference
on Coastal Engineering, American Society of Civil Engineering, pp. 1028–1039.
Shi, F., Kirby, J.T., Haas, K., 2006. Quasi-3D nearshore circulation equations: a CLvortex force formulation. In: Proceedings of the 30th International Conference
on Coastal Engineering, American Society of Civil Engineering, pp. 1028–1039.
Smith, J., 2006. Wave–current interactions in ﬁnite depth. J. Phys. Oceanogr. 36,
1403–1419.
Soulsby, R.L., 1995. Bed shear-stresses due to combined waves and currents. In:
Stive, M., Fredsøe, J., Hamm, L., Soulsby, R., Teisson, C., Winterwerp, J. (Eds.),
Advances in Coastal Morphodynamics. Delft Hydraulics, Delft, The Netherlands.
pp. 420–423.
Soulsby, R.L., 1997. Dynamics of Marine Sands: A Manual for Practical Applications.
Thomas Telford, London.
Stive, M.J.F., De Vriend, H.J., 1994. Shear stresses and mean ﬂow in shoaling andbreaking waves. In: Edge, B.L. (Ed.), Proceedings of the 24th Coastal Engineering
International Conference 1994, Kobe, Japan. American Society of Civil Engineers,
New York, pp. 594–608.
Sullivan, P.P., McWilliams, J.C., Melville, W.K., 2007. Surface gravity wave effects inthe oceanic boundary layer: large eddy simulation with vortex force andstochastic breakers. J. Fluid Mech. 593, 405–452.
Svendsen, I.A., 1984a. Mass ﬂux and undertow in a surf zone. Coastal Eng. 8, 347–
365.
Svendsen, I.A., 1984b. Wave height and set-up in a surf zone. Coastal Eng. 8, 303–
329.
Tajima, Y., Madsen, O.S., 2006. Modeling near-shore waves, surface rollers, andundertow velocity proﬁles. J. Waterway Port Coastal Ocean Eng. 132, 429–438.
Terray, E.A., Donelan, M.A., Agrawal, Y.C., Drennan, W.M., Kahma, K.K., Williams, A.J.,
Hwang, P.A., Kitaigorodskii, S.A., 1996. Estimates of kinetic energy dissipation
under breaking waves. J. Phys. Oceanogr. 26, 792–807.
Terrile, E., Brocchini, M., Christensen, K.H., Kirby, J.T., 2008. Dispersive effects onwave–current interaction and vorticity transport in nearshore ﬂow. Phys. Fluid
20, 036602. doi:10.1063/1.2888973.
Thornton, E.B., Guza, R.T., 1983. Transformation of wave height distribution. J.
Geophys. Res. 88 (C10), 5925–5938.
Thornton, E.B., Whitford, D.J., 1990. Longshore currents over a barred beach – part
II: model, Naval Postgraduate School, Monterey, California, pp. 1–30.
Trowbridge, J., Madsen, O.S., 1984. Turbulent wave boundary layers, 2. Secondorder theory and mass transport. J. Geophys. Res. 89 (C5), 7999–8007.
Uchiyama, Y., McWilliams, J.C., 2008. Infragravity waves in the deep ocean:generation, propagation, and seismic hum excitation. J. Geophys. Res. 113,
C07029. doi:10.1029/2007JC004562.
Uchiyama, Y., McWilliams, J.C., Restrepo, J.M., 2009. Wave–current interaction innearshore shear instability analyzed with a vortex-force formalism. J. Geophys.
Res. 114, C06021. doi:10.1029/2008JC005135.
Umlauf, L., Burchard, H., Hutter, K., 2003. Extending the k–x turbulence model
towards oceanic applications. Ocean Modell. 5, 195–218.
Walstra, D.J.R., Roelvink, J.A., Groeneweg, J., 2000. Calculation of wave-driven
currents in a 3D mean ﬂow model. In: Edge, B.L. (Ed.), Proceedings of the 27th
Coastal Engineering International Conference 2000, Sydney, Australia. American
Society of Civil Engineers, New York, pp. 1050–1063.
Warner, J.C., Sherwood, C.R., Arango, H.G., Signell, R.P., 2005. Performance of four
turbulence closure models implemented using a generic length scale method.
Ocean Modell. 8, 81–113.
Warner, J.C., Sherwood, C.R., Signell, R.P., Harris, C.K., Arango, H.G., 2008.
Development of a three-dimensional, regional, coupled wave, current, andsediment transport model. Comput. Geosci. 34, 1284–1306.
Xia, H., Xia, Z., Zhu, L., 2004. Vertical variation in radiation stress and wave-induced
current. Coastal Eng. 51, 309–321.
Xie, L., Wu, K., Pietrafesa, L.J., Zhang, C., 2001. A numerical study of wave–current
interaction through surface and bottom stresses: wind-driven circulation in theSouth Atlantic Bight under uniform winds. J. Geophys. Res. 106 (C8), 16841–
16855.
Xu, Z., Bowen, A.J., 1994. Wave- and wind-driven ﬂow in water of ﬁnite depth. J.
Phys. Oceanogr. 24, 1850–1866.
Yu, J., Slinn, D.N., 2003. Effect of wave–current interaction on rip currents. J.
Geophys. Res. 108 (C3), 3088–3106.