PUBLICATIONSJournal of Geophysical Research: Oceans
RESEARCH ARTICLE
10.1002/2017JC013005
Key Points:
 The ROMS-WEC model is applied toanalyze 3-D turbulent rip currents
containing signiﬁcant very low
frequency motions (VLFs)
 The intrinsic VLFs are excited onrealistic bathymetries under steady,uniform incident offshore wave
conditions
 Three-dimensionality, the current
effects on waves (CEW), andbathymetric and forcing conditions
on the intrinsic VLFs are investigatedCorrespondence to:
Y. Uchiyama,uchiyama@harbor.kobe-u.ac.jp

Citation:
Uchiyama, Y., J. C. McWilliams, andC. Akan (2017), Three-dimensionaltransient rip currents: Bathymetric
excitation of low-frequency intrinsic
variability, J. Geophys. Res. Oceans, 122,
5826–5849, doi:10.1002/
2017JC013005.
Received 24 APR 2017
Accepted 19 JUN 2017
Accepted article online 5 JUL 2017
Published online 25 JUL 2017

C 2017. American Geophysical Union.
V

All Rights Reserved.

UCHIYAMA ET AL.

Three-dimensional transient rip currents: Bathymetric
excitation of low-frequency intrinsic variability
Yusuke Uchiyama1,2,3, James C. McWilliams2, and Cigdem Akan2

1

Department of Civil Engineering, Kobe University, Kobe, Japan, 2Department of Atmospheric and Oceanic Sciences,
University of California, Los Angeles, California, USA, 3Coastal and Estuarine Environmental Department, Port and Airport
Research Institute, Yokosuka, Japan

Abstract The ROMS-WEC model [Uchiyama et al., 2010] based on an Eulerian wave-averaged vortex-force
asymptotic theory of McWilliams et al. (2004) is applied to analyze 3-D transient wave-driven rip currents andassociated intrinsic very low-frequency (VLF) variability in the surf zone on a surveyed bathymetry under spatiotemporally uniform offshore incident waves. The 3-D rip currents are substantially depth-dependent due tothe vertical recirculation, composed of pairs of counter-rotating longitudinal overturning roll cells that promote surface convergence. The vortex force plays an important role in vorticity budget, preconditioning overall vorticity reduction. These rip currents are intrinsically unstable and contribute about 70% to kinetic energy(KE) as eddy kinetic energy (EKE), consistent with observations. The dominant ﬂuctuation period ﬁts the VLF
band, at about 18 min. The current effect on waves (CEW) alters not only the mean rip structure, but also theassociated turbulence as the modiﬁed cross-shore EKE proﬁle with considerable accentuation in the inner surfzone. Increased alongshore bathymetric variability proportionally intensiﬁes KE and intrinsic EKE, whereas it
reduces the VLF period. With a guide of a pseudo 2D model, we reveal that vortex tilting effect due to the horizontal vorticity inherent in the 3-D rip currents promotes collapse of the 3-D eddies through an enhanced forward kinetic energy cascade, leading to short-lived, laterally-stretched 3-D eddies resulting in elongatedﬁlaments that decay more quickly than coherent, long-lived, circular 2-D eddies.

1. Introduction
Surf zones act as a barrier between inner-shelf and shoreline by wave-driven littoral currents. However, thesurf zone retention is occasionally collapsed by offshore-headed rip currents [Reniers et al., 2009], resulting
in horizontal circulation referred to as rip cells [e.g., Bowen, 1969]. Rip currents are associated with vortex
shedding due to lateral shear and thus are readily unstable [Kennedy, 2003; Haller and Dalrymple, 2001] ashave been observed in numerical models, laboratories and ﬁeld experiments [e.g., Sonu, 1972; Aagaard
et al., 1997; Brander and Short, 2001; Haas and Svendsen, 2002; Yu and Slinn, 2003; MacMahan et al., 2005;
Reniers et al., 2009]. Rip current pulsations under near-normal wave incidence and alongshore-propagating
shear waves [e.g., Oltman-Shay et al., 1989] under oblique wave incidence occur intermittently or quasiperiodically based on both observations [MacMahan et al., 2004a,2004b] and theoretical works [Haller andDalrymple, 2001; Yu, 2006]. The vortical motions in the littoral currents are expected to occur at signiﬁcantly
lower frequencies (f < 0.004 Hz) than infragravity waves (0.004 Hz < f < 0.04 Hz). They are often referred toas very low frequency motions (VLFs), because of the vortical motions at time scales of O (1–10) minutes€
induced by (1) shear instabilities [e.g., Ozkan-Haller
and Kirby, 1999; Haller and Dalrymple, 2001; Noyes et al.,€2004], (2) wave groups [Sonu, 1972; MacMahan et al., 2004a; Long and Ozkan-Haller,
2009], and (3) vorticalforcing due to individual breaking waves [e.g., Spydell and Feddersen, 2009; Feddersen, 2014]. The shear
instability is well explained and represented as an inﬂection point, normal mode instability of a speciﬁed
current proﬁle [cf. Bowen and Holman, 1989; Dodd et al., 1992, 1994, 2000; Allen et al., 1996; Reniers et al.,1997]. In turn, the inﬂuences of wave groups have typically been interpreted in relation to changes in mass
ﬂux associated with the incident groups themselves [Munk, 1949] or infragravity waves generated by thewave groups [Sonu, 1972; MacMahan et al., 2004a]. Reniers et al. [2007] suggested from a two-way coupled,time-dependent phase-averaged wave and barotropic circulation model on a rip-channeled beach that thegroup scale forcing is substantial to excitation of VLFs. These VLFs are most energized by directionally broad
random waves in the surf zone, and they rapidly decay in the offshore direction just outside the surf zone.

3-D TRANSIENT RIP CURRENTS

5826

Journal of Geophysical Research: Oceans

10.1002/2017JC013005€
Long and Ozkan-Haller[2009] reported with a similar model that random wave group forcing signiﬁcantly
broadens frequency-wavenumber velocity spectra of shear wave for alongshore-uniform bathymetry. The
consequent wave-forcing ﬁeld appears as patches with alongshore length scales of O (100) m, leading totransient vortical lateral torque on the scale of the rip-channel spacing. The similar conclusion has beendrawn from Spydell and Feddersen [2009], Feddersen et al. [2011], and Feddersen [2014], who demonstrated
that VLFs are excited as intensely as observations, as a consequence of wave-breaking eddy generation in aphase-resolving, barotropic Boussinesq model with directionally spread, random incident waves on analongshore-uniform beach. Such a stochastic wave forcing also stimulates transient rip currents in a phaseaveraged 3-D model [Kumar and Feddersen, 2017a]. A consensus is being achieved such that the transient
irrotational forcing due to group waves and/or randomness of ﬁnite-crest length wave breaking is responsible for a part of VLFs (hereinafter, ‘‘forced’’ VLFs). Nevertheless, alongshore-uniform bathymetry is not the
optimal test bed to measure spontaneous VLFs (‘‘intrinsic’’ VLF) because it does not provoke alongshorevariable vortical torque due to wave breaking, which is a main source of the intrinsic VLFs [Geiman and
Kirby, 2013].
The feedback mechanisms of the current effect on waves (CEW) [cf. Uchiyama et al., 2009, hereinafter
UMR09], consisting of the frequency Doppler shifting and refraction of primary wave on ambient phaseaveraged currents and surface elevation (viz., wave set-up/down), is anticipated to alter VLFs because it is€
known to attenuate instabilities of shear waves [Ozkan-Hallerand Li, 2003; UMR09] and rip currents [Yu andSlinn, 2003; Haas et al., 2003; Weir et al., 2011]. Bottom drag is another important factor determining the stability limits [e.g., Bowen and Holman, 1989; Haas et al., 2003]. Wave-driven currents are generally accompanied by strong vertical shear due to vertically recirculating ﬂows caused by surface-intensiﬁed, shoreward
breaker acceleration that primarily balances with cross-shore pressure gradient to promote seaward undertow currents conﬁned near the bed [cf. Uchiyama et al., 2010, hereinafter UMS10]. While bottom boundary
layer streaming acceleration is of second importance to surf-zone dynamics, it is essential to near-bed offshore transport of negatively buoyant larvae [Fujimura et al., 2014]. Rip currents are also observed to havesigniﬁcant vertical structure [Haas and Svendsen, 2002; MacMahan et al., 2005]. The vorticity around the surfzone promotes a nontrivial vortex force particularly in the alongshore direction interplayed with Stokes drift[Kumar et al., 2012; Kumar and Feddersen, 2017a]. Therefore, 3-D wave-averaged circulation models basedon a vortex force formalism [McWilliams et al., 2004, hereinafter MRL04] with CEW are inevitable for conducting this type of study, whereas most of the numerical studies have been conducted with a depthaveraged (i.e., barotropic), 2-D wave-averaged model, based on the classical radiation stress formalism [e.g.,
Longuet-Higgins and Stewart, 1962; Hasselmann, 1971; Phillips, 1977].
Among a limited number of the numerical studies on the 3-D transient littoral currents, Reniers et al. [2009]
extended their earlier work based on a barotropic framework [Reniers et al., 2007] by exploiting the Delft3D
model [Walstra et al., 2000; Lesser et al., 2004], forced by a spectral peak WKB wave model operated on the
wave group timescale. They successfully reproduced the combined forced and intrinsic VLFs for the periodsof 4–60 min in the rip current during the RCEX ﬁeld experiment on the surveyed variable bathymetry[MacMahan et al., 2010a] and further described lateral mixing efﬁciency as a Lagrangian coherent structure
depicted by the ﬁnite-time Lyapunov exponent [Reniers et al., 2010]. Geiman et al. [2011] compared the resultof Reniers et al. [2009] with that from the phase-resolving 2-D Boussinesq model for the RCEX data andreported that both the models yielded comparable Lagrangian mixing efﬁciency such as absolute and pair
separations, consistent with the observations. Marchesiello et al. [2015] applied the two-way coupled WKB
spectral peak wave model with the 3-D circulation model developed in UMS10 to transient rip currents on
Biscarrosse Beach, France. They successfully reproduced energetic intrinsic VLFs at O (10) minutes, mainly
driven by slowly modulating waves on timescales of O (1–2) hours due to tides and wind. The aforementioned
study of Kumar and Feddersen [2017a] utilized a wave-averaged 3-D circulation model forced by stochasticallybreaking directional waves on an alongshore-uniform planar beach to evaluate cross-shelf transition of themixing efﬁciency due to the forced VLFs in transient rip currents. Newberger and Allen [2007] examined the
intrinsic VLFs in 3-D shear waves during the DUCK94 ﬁeld experiment. They pointed out that 3-D shear waves
are more stable than 2-D, consistent with the result from a quasi-3-D model [Zhao et al., 2003], and that omission of the vortex force term and CEW increase the errors in the alongshore velocity.
However, 3-D numerical studies on the transient rip currents in the surf zone have not been extensively
conducted yet to fully answer questions such as: What is the 3-D nature of the transient rip currents? Is

UCHIYAMA ET AL.

3-D TRANSIENT RIP CURRENTS

5827

Journal of Geophysical Research: Oceans

10.1002/2017JC013005

bathymetric variability essential in promoting the 3-D intrinsic VLFs? What controls the excitation of the
intrinsic VLFs? How and why are the 3-D intrinsic VLFs different from the 2-D ones? In the present study we
aim to answer those questions by using a 3-D numerical framework, ROMS-WEC (UMS10), which is a waveaveraged version of ROMS (Regional Oceanic Modeling System) [Shchepetkin and McWilliams, 2005, 2008]
based on a vortex force formalism of MRL04, tightly coupled with a spectral peak WKB wave refraction
model. Three-dimensional numerical experiments are conducted for turbulent rip currents driven by normal
incident waves. The ROMS-WEC model has been adapted to 3-D coastal applications with prevailing wavedriven currents, showing quite good agreements with measurements [cf. Marchesiello et al., 2015; Kumar
et al., 2015]. We show here that time-invariant and spatially homogeneous offshore incident waves, which
are not expected to induce the forced VLFs, excite 3-D turbulent rip currents on timescales of O (1–10)
minutes on an immobile beach with realistic spatial bathymetric variability ubiquitous in nature. Our focus
is on the intrinsic part of the VLFs, yet the forced VLFs usually contain equivalent turbulence intensity [e.g.,€
Reniers et al., 2007] or enstrophy [e.g., Long and Ozkan-Haller,
2009] with that in the intrinsic VLFs. We further demonstrate rip currents are spontaneous for spatially variable bathymetry with alongshore undulations, resulting in unstable motions that contain signiﬁcantly energetic intrinsic VLFs, comparable to thosemeasured during the SandyDuck ﬁeld campaign [MacMahan et al., 2010b; de Schipper et al., 2012] without
wave groups and random wave breaking. Supplemental simulations are carried out to examine inﬂuencesfrom three-dimensionality, alongshore variability in the bathymetry, wave-current interaction (i.e., CEW),
and different forcing conditions on the intrinsic VLFs.
The remainder of the paper is organized as follows. The numerical model and experimental conﬁguration
are described in section 2, supplemented by Appendices. Section 3 demonstrates 3-D transient rip currents
containing signiﬁcant intrinsic VLFs are driven by steady, normal incident waves for the surveyed beachbathymetry. In section 4, the eddy generation mechanism is discussed through the eddy statistic analysis
and the barotropic vorticity budget analysis, followed by a series of model experiments to demonstrate
modiﬁcation of the intrinsic VLFs in response to CEW, amplitudes of bathymetric variability, and forcing conditions. A comparative spin-down experiment for decaying processes of the initially forced eddies in rip currents is presented in section 5 to further examine 3-D effects on the fate of eddies. Section 6 gives
conclusions.

2. Methods
2.1. The ROMS-WEC Model
The ROMS-WEC model [Uchiyama et al., 2010; UMS10] is a 3-D wave-averaged ocean circulation model
based on a terrain-following, hydrostatic ocean model, ROMS [e.g., Shchepetkin and McWilliams, 2005, 2008],
which accounts for the interaction of surface gravity waves and currents with a vortex force formalism.
Eulerian wave-averaged conservation equations for mass, momentum, and tracers are included in the
ROMS based on an asymptotic theory of MRL04, where conservative WEC (wave effects on currents) include
the vortex force, the Stokes-Coriolis force, Bernoulli head, and quasi-static pressure and sea-level responses.
The ROMS-WEC further considers nonconservative WEC due to wave breaking, associated surface roller
waves, bottom streaming, wave-enhanced vertical mixing, and wave-enhanced bottom drag suitable tocoastal and nearshore applications (Appendix A). The ROMS-WEC model was applied to 3-D transient rip
currents, and validated against in situ wave ﬁeld observed by a video imaging technique [Holman and Haller, 2013] and the ADCP data in Marchesiello et al. [2015].
The depth-averaged, barotropic 2-D model is conﬁgured and operated as the stand-alone model, based on
the barotropic mode of the ROMS-WEC by omitting the barotropic-baroclinic coupling requisite for fully 3-D
simulations (Appendix B). The boundary conditions and our choice of the submodels such as the bottom
drag prescription are provided in Appendix C. Both of the 3-D and barotropic 2-D ROMS-WEC models aretightly coupled with the built-in spectrum-peak WKB wave-refraction model that takes into account theeffect of currents on waves (CEW) as described in Appendix D.
2.2. The Experimental Setup
The ROMS-WEC model is applied to transient rip currents on a realistic beach bathymetry driven by steady,time-invariant, alongshore-uniform normal incident offshore waves. We choose a surveyed bathymetry near
the Field Research Facility (FRF) in Duck, NC, on 22 February 2010, when a pronounced rip channel emerged

UCHIYAMA ET AL.

3-D TRANSIENT RIP CURRENTS

5828

Journal of Geophysical Research: Oceans

10.1002/2017JC013005

near the 613 m-long FRF pierat around y 5 460 m (Figure1). The beach is double-barredwith the bathymetric peaks
approximately located for
220  x  280 m (the innerbar) and 360  x  420 (theouter bar). The model domain
extends laterally to 1,024 m
offshore 3 1,152 m northwardwith uniform grid spacing and
20 stretched vertical s-layers(Table 1). We also carry out asupplemental barotropic 2-D
simulation for the same bathymetric and forcing conditions.
We show the results forced by
the time-invariant offshore
incident waves with the RMS
wave height Hrms of 1.2 m, thepeak wave period Tp of 10 s,
and the peak wave direction
hp of 08 as the offshore
clamped boundary conditionof the WKB wave model(Appendix D). The feedback
mechanism of CEW is always
Figure 1. Model bathymetry showing the entire model domain used in the present study, reptaken into account unless othresenting a sandy beach at Duck, NC, based on the surveyed bathymetry on 22 February 2010.
erwise stated. Although the
The horizontal coordinate follows the deﬁnition relative to Field Research Facility (FRF) of U.S.
Army Corps of Engineers. Contours with colors show depth h (m). The black thick rectangle is
ROMS-WEC accounts forthe blow-up region used in Figure 3. Open triangles are four monitoring stations A–D in Figure
buoyancy effects by solving
6. The dash-dot lines represent the boundaries outward from which bathymetries are blended
the wave-averaged heat andfor the alongshore periodicity.
salinity conservation equations and the equation of state, we assume constant density in the present study. The imposed external forcing
only comes from steady and uniform incident waves. Other forcing such as surface wind stress, heat, and freshwater ﬂuxes are all omitted. The Coriolis and Stokes-Coriolis terms are also neglected because of the small
domains and the short time scales. The surface roller model (Appendix D) is not activated to avoid further complexity. We leave these terms and equations merely for the completeness. No additional lateral mixing is introduced except for small hyper diffusion inherent in the third order upstream-biased advection scheme of ROMS.
The resting ocean is assumed initially and offshore waves are imposed instantaneously at t 5 0. Because instability
intrinsically takes place immediately after the initialization from the resting state, neither random noise nor initialramping of the wave forcing is imposed as was done in the shear wave experiment in UMR09. The baroclinic andbarotropic time steps are typically 1 sec and 0.15 s. The model output occurs every minute until t 5 600 min. The
analyses are mostly performed for 100  t  600 min to exclude the initial adjustment and spin-up effects. Notethat the results shown below are quite robust qualitatively even with different incident wave and bathymetric conditions unless the incident angle |hp| > 38.
Table 1. Model Conﬁguration
Grid cells
256 3 288 3 vertical 20 layers
Grid spacing
Dx 5 Dy 5 4 m
s-coordinate parameters
hs 5 hb 5 5, hc52.0 m
1.2 m
Incident wave height, Hrms
10.0 sIncident wave period, Tp
08 (to W)
Incident wave angle, hp
Non-wave surface elevation, f0m
CT breaking parameters
Bb 5 1.3, cb 5 0.38

UCHIYAMA ET AL.

3-D TRANSIENT RIP CURRENTS

3. Transient 3-D Rip Currents
3.1. Mean Cross-Shore Profiles
Figure 2 shows cross-shore proﬁles of the time- andalongshore-averaged (denoted by <>) waves, barotropic current ubt , and barotropic Stokes drift uSt
bt ,
5829

Journal of Geophysical Research: Oceans

10.1002/2017JC013005

Figure 2. Cross-shore proﬁles of time- and alongshore-averaged (denoted by the <> operator) model results forced by steady offshore
waves with the RMS wave height Hrms of 1.2 m, the peak wave period Tp of 10 s, and the peak wave direction hp of 08 (normal incidence)for the fully 3-D (solid lines) and the barotropic 2-D (dotted lines) cases: (a) Hrms (blue lines) and Tp (red lines), (b) hp (blue) relative to theshore-normal direction and the kinematic breaker dissipation rate b =q (red), and (c) the cross-shore (offshore) and (d) the alongshore(northward) components of the barotropic velocity ubt (blue) and barotropic Stokes drift uSt
bt (red). Notice the line color corresponds to thatof the y axis and its labels. The thick gray line in Figure 2c represents the alongshore-averaged bathymetry, -h (m), which shares the labelson the right red axis. The two vertical dashed lines denote the approximate offshore extents of the inner (x < 340 m) and outer(340  x < 500 m) surf zones.

from the 3-D and 2-D models. The alongshore-averaged bathymetry is overlaid with the thick gray curve in(c). The mildly sloping beach with a mean slope of about 1/100 has two plateaus segmented by two longshore bars peaked at around x 5 240 (the inner bar) and 400 m (the outer bar). The wave ﬁelds are almost
identical in the 2-D and 3-D models. Depth-induced breaking b peaks at three locations with the maximum
appearing in the vicinity of the shoreline and the others on the seaward faces of the bars. We thus deﬁne
the inner surf zone for x < 340 m and the outer surf zone for 340  x < 500 m. The wave amplitudedecreases toward the shore in a stepwise manner corresponding to the occurrences of wave breaking.
The wave period and direction are slightly modulated within the inner surf zone, while the waves propagate
westward approximately normal to the shore (|hp| < 0.168). Likewise, the resultant mean barotropic Eulerian
offshore velocity ubt , whose time average is approximately anti-Stokes (cf. UMR09) because the steadycross-shore mass balance requires the anti-Stokes proﬁle, @x ulbt 50, where uSt
bt < 0. The small deviation fromthe normal direction in hp invokes the weak time-averaged northward current of vbt < 0.11 m/s for the 3-D
case, whereas vbt < 0.05 m/s with the seaward peak shift by 30–40 m manifests in the 2-D model. The assoStSt
ciated vbt
is negligibly small with jvbt
j < 1 mm/s. The quantitative consistency in the rip currents (ubt ) andthe perpendicular component (vbt ) ensures that the 3-D and 2-D models represent the equivalent results oftransient rip currents, which allow us feasible inter-model comparisons.
3.2. 3-D Structures
Instantaneous 3-D ﬂow ﬁelds of unstable, transient rip currents at t 5 299 and 327 min are depicted with
the near-surface and near-bed velocity vectors on the associated barotropic relative vorticity (Figure 3).
Snapshots of the 8 cross-sectional distributions of 3-D velocity at t 5 299 min are plotted in Figure 4. At
t 5 299 min, the rip current bifurcates beyond the inner surf zone by x 5 340 m, and the eruption of amushroom-like, coherent vortex dipole with paired positive and negative vorticity is evident. In turn, ﬁner
vortices and elongated ﬁlaments dominate at t 5 327 min, with strong vorticity aligning alongshore with asign change at intervals of about 60–80 m near the shore (x < 160 m) that feeds initial development of rip
currents. The torque of the breaker acceleration (see also section 4.2) yields positive and negative relative
vorticity in the surf zone, leading to alternate offshore and onshore currents aligned in the alongshore
direction that can readily be unstable [e.g., Haller and Dalrymple, 2001; Yu, 2006]. The developed rip cells at
t 5 299 min are not quite persistent with time, but short-lived particularly in the 3-D model. The coherent
eddies subsequently collapse into much ﬁner vortices and elongated ﬁlaments such as those at t 5 327
min, followed by reorganization to evolve into coherent eddies and perceptible rip currents (not shown).
This process occurs quasi-periodically at a period of O (1–10) min, which ﬁts the VLF band (section 3.4).

UCHIYAMA ET AL.

3-D TRANSIENT RIP CURRENTS

5830

Journal of Geophysical Research: Oceans

10.1002/2017JC013005

Figure 3. Plan view plots of instantaneous 3-D rip cells in the blow-up region shown in Figure 1 at (left) (a) t 5 299 and (b) (right) 327 minafter the initialization. The near-surface, top-most s-layer horizontal velocity (black) and the near-bed, bottom-most s-layer horizontal
velocity (green) and the barotropic relative vorticity (color) are depicted. The vertical dashed lines denote the approximate offshore
extents of the inner surf zone (x < 340 m). The thin contours are isobaths at 1 m intervals. At t 5 299 min, eruption of a mushroom-like,coherent vortex dipole with paired positive and negative vorticity is evident, while ﬁner vortices and elongated ﬁlaments dominate at
t 5 327 min. The velocity vectors are subsampled every 4 grid-points.

Figure 4. Snapshots of the 3-D structure of rip currents at t 5 299 min (see Figure 3a for the plan view). (left) Offshore velocity u, (middle)
northward alongshore velocity v, and (right) vertical upward velocity w are plotted along 8 alongshore cross-sections at 32 m intervals in y.
The vertical dashed lines denote the approximate offshore extents of the inner surf zone (x < 340 m). The main offshore rip current isdeveloped around y 5 450–482 m with the peak u at 1–3 m depths along with strong up- and down-welling w up to 610 cm/s indicatingprominent convergence and divergence.

UCHIYAMA ET AL.

3-D TRANSIENT RIP CURRENTS

5831

Journal of Geophysical Research: Oceans

10.1002/2017JC013005

Vigorous offshore eruption of the rip current is demonstrated at t 5 299 min around the rip channel formednear y 5 460 m (Figure 1). It gradually veers northward, followed by a further left turn, resulting in a laterally
circulating cyclonic rip cell with a diameter of about 250 m. The seaward rip current velocity extends up tox 5 450 m, about 1.7 times as wide as the surf zone width. Affected by the near-surface breaker-driven ﬂow
and the near-bed undertow, the maxima of the offshore-directed rip current in u appear at depth, typicallyat 1–3 m below the mean sea surface, rather than right at the surface. In contrast, the onshore return ﬂow(u < 0) has a surface conﬁned structure (Figure 4a: 514  y  610 m, x < 300 m) due to the coincident
onshore breaker acceleration near the surface. Depth-dependency is more pronounced both in u and v
especially where velocity magnitude is small, while the main body of the rip current has less prominent vertical structure (Figure 3). Associated vertical velocity w reaches several cm/s particularly near the main rip.
At y 5 450, 482 and 514 m in Figure 4c, upwelling occurs on the both sides of the main rip, whereas downward ﬂow is formed at the rip center where vertical vorticity changes its sign (Figure 3). These results indicate that the 3-D rip currents are composed of pairs of the counter-rotating longitudinal (i.e., following therip current) overturning roll cells that lead to convergence of the surface ﬂow at the rip center while divergence at the bottom.
3.3. Spatiotemporal Development
The 3-D transient rip currents are compared with the 2-D counterpart. Figure 5 shows sequential blow-up
snapshots of the barotropic relative vorticity q5rH 3 ubt for 290  t  330 min at 10 min intervals from the2-D and 3-D models. Time-varying transient rip currents and associated vortices are evident with the quasipersistent strong rip current formed in the rip channel around y 5 460 m. In addition, offshore eruption of
rip currents occurs at other places, including those originating from the shore at around y 5 250, 650 and
750 m. The approximate maximum offshore extent of the rip currents is slightly beyond y 5 500 m. The 2-D
model clearly illustrates more organized, coherent vortices that roll up into circular-shaped eddies to form
positive and negative vortex dipoles. By contrast, such organized structures are immediately destroyed inthe 3-D model. The 3-D rip-induced vortices tend to be stretched seaward and to evolve into ﬁlaments andstreaks. As a consequence, the horizontal aspect ratio of the 3-D vortices is distorted with longer crossshore length scales as compared to the 2-D case.
3.4. Intrinsic VLFs
The rip currents for the surveyed bathymetry under the near-normal incident steady waves are intrinsicallyunstable with low frequency variability (Figure 5). We deﬁne the four sampling stations, Sta. A–D (Figure 1),

Figure 5. Temporal evolution of the barotropic relative vorticity q5rH 3 ubt in color for the 2-D (top row) and 3-D (bottom row) cases at
10-minute intervals. The two vertical dashed lines denote the approximate offshore extents of the inner (x < 340 m) and outer(340  x < 500 m) surf zones. The thin contours are isobaths at 1 m intervals. The main rip is formed and ﬂuctuates in the middle of thedomain with a pair of counter-rotating vortex dipole. This structure tends to roll up into circular shapes in the 2-D particularly at t 5 310
and 330 min, while the 3-D vortex pair is more stretched longitudinally and elongated as ﬁlaments.

UCHIYAMA ET AL.

3-D TRANSIENT RIP CURRENTS

5832

Journal of Geophysical Research: Oceans

10.1002/2017JC013005

Figure 6. Time series of (a) cross-shore barotropic velocity ubt, (b) alongshore barotropic velocity vbt, and (c) frequency power spectra of
the barotropic kinetic energy, Sbt, from the fully 3-D model, at the 4 sampling stations (Figure 1) chosen for the prevailing seaward rip
current (Sta. A), the offshore quiescent rip current (Sta. B), the bifurcation on the shoal (Sta. C), and the shoreward return ﬂow (Sta. D). The
shaded area in Figure 6c represents the VLF band for the periods between 5 and 100 min.

to monitor temporal ﬂuctuations of the transient rip currents and their frequency spectra. The choice ismade intentionally for the prevailing seaward rip current (Sta. A), the offshore quiescent rip current (Sta. B),
the bifurcation on the shoal (Sta. C), and the shoreward return ﬂow (Sta. D). The barotropic velocity timeseries plots (Figures 6a and 6b) clarify that the rip current, the return ﬂow, and the bifurcating ﬂow ﬂuctuatewith time and are composed of many frequency components as depicted in their kinetic energy frequency
spectra (Figure 6c). The predominant ﬂuctuations occur for the periods of 30–80 min with higher frequency
variability superposed. The higher frequency components are signiﬁcantly weak at the offshore Sta. B,
whereas the lower frequency energy is still apparent. The energy is extensively distributed in broadband
and thus the spectra are ‘‘red’’ to some extent, which may require more elaborate methods such as withfrequency-wavenumber spectra to selectively extract the vortical component of the VLFs from that irrotational part involving infragravity waves [e.g., MacMahan et al., 2010b; Feddersen, 2014]. Nevertheless, with avisual guide from Figures 6a and 6b, spectral peaks are detected at about 5 3 1024 Hz (33.3 min) and 2 31023 Hz (8.3 min) at all the 4 stations (Figure 6c), ﬁtting the VLF band of f < 0.004 Hz [e.g., Reniers et al.,2007]. Therefore, we deﬁne our VLF band for the periods between 5 and 100 min and extract the intrinsic
VLF variability from the model results using a frequency band-pass ﬁlter.

4. Intrinsic VLFs in 3-D Transient Rip Currents
4.1. Eddy Statistics
To extract 3-D effects on turbulent rip currents and associated vortices, we ﬁrst examine cross-shore proﬁlesof depth- and alongshore-averaged eddy statistics (Figure 7). For the 3-D case, we introduce decompositions into the depth-invariant barotropic component u0 bt , the depth-dependent baroclinic component u0 bc ,and their sum as the total, where the prime and the subscrirpt bc respectively stand for the deviation fromthe ensemble (time) average denoted by the over bar, anddepth-dependent component (i.e., the u5u1u0 ; u5ubt 1ubc ). Based on these notations, the EKE (5 u0 2 1v 0 2 =2) and kinematic lateral Reynolds
stress (52u0 v 0 ) are decomposed into the barotropic component with u0bt and the baroclinic component
with u0bc for the 3-D case. Furthermore, the barotropic EKE production rate, denoted by KmKe, in an EKE conservation equation is
Km Ke 52u0 2bt  @x ubt 2u0 bt v 0 bt  @y ubt 2u0 bt v 0 bt  @x vbt 2v 0 2bt  @y vbt ;(1)which accounts for the mean to eddy kinetic energy conversion rate responsible for feeding shear instability
if it is positive.
In Figure 7a, the barotropic and total EKE of the 3-D model are approximately the same as the barotropic
EKE of the 2-D model, although the 3-D EKEs are slightly larger. The baroclinic EKE contributes only 5.0% to
the total EKE. The depth-averaged horizontal Reynolds stress (Figure 7b) shows the largest difference
between the two models among the other statistics. The peak magnitudes appear near the shore

UCHIYAMA ET AL.

3-D TRANSIENT RIP CURRENTS

5833

Journal of Geophysical Research: Oceans

10.1002/2017JC013005

Figure 7. Cross-shore proﬁles of depth- and alongshore-averaged (a) eddy kinetic energy (EKE), (b) lateral kinematic Reynolds stress, (c)time-averaged enstrophy of barotropic vorticity (q2/2), and (d) barotropic EKE production rate (equation (1)). The RMSs in Figures 7c and7d are the deviations from the alongshore average. The two vertical dashed lines denote the approximate offshore extents of the inner(x < 340 m) and outer (340  x < 500 m) surf zones.

(x 5 140 m) and shoreward of the inner bar (x 5 240 m), both of which are close to the peaks of the breakeracceleration b (Figure 2b) where the total EKE increases (Figure 7a). The third peak from the shore emerges
further offshore at x 5 320 m in both the models, whereas the 3-D Reynolds stress is more widely distributed
in x with the baroclinic contribution having 6.2% of the total. The time- and alongshore-averaged enstrophyof the barotropic vorticity q signiﬁcantly increases near the shore (feeder region; x < 200 m) in the 3-D case
(Figure 7c), consistent with the barotropic EKE (Figure 7a), while their cross-shore averages for x  600 m arealmost identical. In Figure 7d, the alongshore-averaged Km Ke is always positive with the 3 peaks approximately
coinciding with those of b and EKE, which evidently demonstrates that shear instability intensiﬁes EKE as the
major source of transient rip currents [e.g., Haller and Dalrymple, 2001] in particular in the feeder region andthe inner surf zone (x < 340 m). The 3-D Km Ke is somewhat larger than that of the 2-D, consistent with the
EKE. The difference is greatest around the peaks, in particular in the vicinity of the shore (x < 140 m) and theinner bar crest (x  270 m). In summary, the 3-D model enhances overall eddy activity in transient rip currents
due to the increased production rate around the breaking points, particularly on the sloping shore face andthe inner bar crest. The former accentuates enstrophy signiﬁcantly right at the shore, while the latter correlates with modiﬁcation in the cross-shore proﬁle of the horizontal Reynolds stress with enhanced KmKe, leading to more pronounced shear instability. Nevertheless, the depth-dependent transient velocity u0bc makes
moderate contributions, by about 5–6%, to the generation of the EKE and the Reynolds stress.
4.2. Barotropic Vorticity Budget
The conservation equation for the barotropic vorticity q5rH 3 ubt is derived by taking the curl of the barotropic momentum equation in equation (A5).
 wd 
scd
ssur
2
St
k 2rH 3 bot@t q1 rH 3½ðq1fc Þ^z 3 ubt  52rH 3ðfc z^ 3 uSt
Þ
2r
3ðq^
z
3
u
Þ
1
r
3
2r
3
H
H
H
bt
bt
|ﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄ{zﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄ}|ﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄ{zﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄ}
q0 D
q0 rD
q0 D
|ﬄﬄﬄﬄﬄﬄ{zﬄﬄﬄﬄﬄ
ﬄ} |ﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄ{zﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄ}|ﬄﬄﬄﬄﬄﬄﬄﬄ{zﬄﬄﬄﬄﬄﬄﬄ
ﬄ}
advection ðadvÞ
curl J
curl Bcurl S

curl D(2)

Notice that the left hand side (LHS) terms are the nonwave contributions comprising the tendency andadvection effects combined with vorticity stretching (adv), while all the right hand side (RHS) terms reﬂect
WEC. Because we assume no rotation with fc 5 0, the Stokes-Coriolis (the ﬁrst RHS) term vanishes. The other
RHS terms are the vortex force torque (curl J), the breaker acceleration torque (curl B), the bottom boundary
streaming torque (curl S), and the bottom drag torque (curl D). Such barotropic vorticity analyses havebeen conducted in Bonneton et al. [2010] and Bruneau et al. [2011] with 2-D models based on a waveaveraged shallow water equation considering the vortex force [Smith, 2006; Ardhuin et al., 2008]. Both concluded that the curl B is the main source of vorticity, yet the curl J was not discussed and curl S wasomitted.

UCHIYAMA ET AL.

3-D TRANSIENT RIP CURRENTS

5834

Journal of Geophysical Research: Oceans

10.1002/2017JC013005

Figure 8 shows snapshots of the barotropic vorticity q, the Lagrangian transport Dulbt , the RMS wave heightHrms, the vorticity advection term (adv), and the four WEC torques for the 3-D case at t 5 300 min when theenergetic rip current and the associated counterclockwise rotating transport develop on the rip channel aroundy 5 460 m. The advection is the leading term that approximately balances with the sum of the tendency term(not shown) and the four RHS WEC terms. In the feeder region for x < 200 m, q is forced by the intense curl B,
where wave breaking occurs as partially depicted by shoreward decreasing Hrms, but compensated for by thecurl D. Thus, alongshore-varying vorticity sources for the intrinsic VLFs are generated with the strong bathymetric constraint as argued in Geiman and Kirby [2013]. Note that random directional incident waves lead to ﬁnite
breaking crest-length vorticity generation, which also provokes alongshore-variable q source for the forced VLF
excitation [Spydell and Feddersen, 2009]. In the large portion of the offshore area (x > 200 m), the curl D substantially absorbs q as exhibited by the similar spatial distribution with the opposite sign. In the neck region(200 < x < 340 m), q is fed again by curl B where the EKE production rate is considerably intensiﬁed (Figure 7d).
The curl S resembles the curl B in the patterns and is generally negative to decelerate offshore-directed ﬂow
including rip currents, while the curl S is much smaller than the curl B by only a few percent.
The vortex force torque curl J plays a dominant role in the vorticity budget other than the breaker torquecurl B as the primary vorticity source. The contribution of curl J reaches about 60–70% of curl B, which isquite signiﬁcant. In the 2-D momentum balance, the vortex force is known to counteract the breaker acceleration [Bruneau et al., 2011]. To examine the role played by the vortex force in the 2-D vorticity balance, we
further rearrange curl J into:

Figure 8. Snapshots of (a) barotropic vertical vorticity q (color) and the barotropic velocity ubt (vectors), (b) Lagrangian transport Dulbt(vectors) and its magnitude (color), (c) RMS wave height and barotropic Stokes drift uSt
bt (vectors), (d) vorticity advection term, (e) breakeracceleration torque curl B, (f) bottom-streaming torque curl S, (g) bottom drag torque curl D, and (h) vortex force torque curl J, at t 5 300
min for the 3-D case. Notice that different color scales are used. The two vertical dashed lines denote the approximate offshore extents of
the inner (x < 340 m) and outer (340  x < 500 m) surf zones. The thin contours are isobaths at 1 m intervals. The main rip vorticity andassociated offshore Lagrangian transport appear around the middle of the domain, forced by the curl B that determines the initial widthof the rip jet near the shore (i.e., the feeder region). Nevertheless, the vortex dipole develops with longer alongshore wavelength in andbeyond the feeder region (200 < x < 340 m; the neck region) where the predominant torque is curl J, while curl D is secondarily important.
In contrast, curl S is generally negative to decelerate the offshore-headed ﬂow but is much less than the other torques by an order ofmagnitude. Note that the 2D vorticity balance is qualitatively the same and thus is not shown.

UCHIYAMA ET AL.

3-D TRANSIENT RIP CURRENTS

5835

Journal of Geophysical Research: Oceans

10.1002/2017JC013005





St
St
St
curl J52rH 3 q^z 3 uSt
bt 52rH  ubt q 52 @x ubt q 1@y vbt q :(3)St
Because hp  0 (Figure 2b) and vbt
 0 (Figure 2d), curl J  2@x uSt
bt g . On the present coordinate system,St
ubt generally increases from the shore except for the decreasing regions shoreward of the breaking points
near the two longshore bars (Figure 2c). In regions of the seaward increasing uSt
bt , positive (negative) q
results in negative (positive) curl J, and thus q always reduces its magnitude. The opposite is true where uSt
bt
decreases seaward. Consequently, the rip currents and associated vorticity are extensively suppressed
under the near-normal onshore uSt
bt condition, whereas they recharge vorticity in the shoreward of the
breaking points before going out to the offshore rip-head region (340 < x < 500 m; the outer surf zone). For
instance, the positive q of the main rip eddy indicated by the black boxes in Figures 8a, h is increased by
the positive curl J in the shoreward of x 5 340 m, whereas decreased by the negative curl J beyond that.

Therefore, in the transient rip currents, the four wave-induced forcing terms substantially modify the usual
nonwave vorticity balance. The breaker torque (curl B) is the primary vorticity source particularly in the feeder
region near the shore and is generally reduced by the bottom drag torque (curl D) and the bottom streaming
torque (curl S). The vortex force torque (curl J) has magnitudes comparable to curl B and works not only as an
overall vorticity sink, but also as a source in the shoreward of the breaking points, where uSt
bt decreases offshore, coincidently with curl B. In the rip-head region (x > 340 m), as the curl B almost vanishes, the nonwave
advection term takes over the vorticity balance to compete with the tendency term. Note that the barotropic
vorticity budget of the 2-D case is qualitatively the same as the 3-D case, and thus is not shown here. In addition, the above argument is most striking under conditions of strong rip currents, although the temporal
RMSs of the terms provide less pronounced curl J effects reduced by 44–71% of those in Figure 8.
4.3. CEW and Three-Dimensionality on Intrinsic VLFs
In previous studies, attention has not been paid to inﬂuences of CEW and three-dimensionality on transient
processes in the VLF band deﬁned in section 3.4. Because CEW has been well known to reduce the offshore
extent of the rip currents [Haas et al., 2003; Yu and Slinn, 2003; Weir et al., 2011], we further conduct two
more simulations using the same 3-D and 2-D models without CEW (denoted by NC). Figure 9 shows the
cross-shore proﬁles of the time- and alongshore-averaged total kinetic energy KEbt, EKE for the VLF band
EKEvlf, the ratio of EKEvlf to the total EKE (EKEtot), and a mean VLF period Tm2 for the 4 cases. The last item is
evaluated as the square-root of the ratio of two spectral moments,
ð1
pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ
Tm2 5 m0 =m2 ; where mn 5
fbtn Sbt ðfbt Þdfbt
(4)
21

fbt: frequencies of the 2-D transient rip currents, and Sbt (fbt): the barotropic kinetic energy frequency spectrum. The volume averages of the quantities in Figure 9 for x  600 m are summarized in Table 2. Note thatthe barotropic velocity ubt is exploited in the present analysis.

Figure 9. Cross-shore proﬁles of time- and alongshore-averaged (a) total barotropic kinetic energy KEbt, (b) barotropic eddy kinetic energy
in the VLF band EKEvlf, (c) ratios of EKEvlf/EKEtot, where EKEtot is the total barotropic EKE, and (d) the spectral mean VLF period Tm2, plotted
for the 3-D (blue) and 2-D (magenta) cases. ‘‘NC’’ denotes the cases without CEW (broken curves). The two vertical dashed lines denote theapproximate offshore extents of the inner (x < 340 m) and outer (340  x < 500 m) surf zones. Differences are evident in KEbt and EKEvlf
between the cases with CEW and the NC cases.

UCHIYAMA ET AL.

3-D TRANSIENT RIP CURRENTS

5836

Journal of Geophysical Research: Oceans

10.1002/2017JC013005

The 2-D and 3-D differences are subtle as expected for our choice of thebottomdragparameterization
3-D
3D-NC
2-D
2D-NC
Unit
(AppendixC).
The
inclusionof CEW
22
2 22
6.25
8.49
5.75
8.54
10 m sKEbt
EKEtot
4.28
3.45
3.83
3.67
1022 m2s22
prominently alters not only <KEbt>,
EKEvlf
2.80
2.53
2.62
2.58
1022 m2s22
but also <EKEvlf>. As the effect of
Tm2
17.9
17.7
18.9
15.4
min
reducing the offshore extent of the
71.9
46.1
69.5
49.2
%
EKEtot/KEbt
46.1
32.2
46.5
34.4
%
EKEvlf/KEbt
rip current, <KEbt> decreases remark64.3
71.1
67.1
70.1
%
EKEvlf/EKEtot
ably due to CEW in the neck and riphead regions (i.e., 160 < x < 350 m)
where the rip currents are most active (e.g., Figure 8). The overall reduction of KEbt due to CEW reaches 26–
33%. However, intensiﬁcation of <EKEvlf> is caused by the CEW with the most inﬂuential region in the vicinity of the shore (i.e., the feeder region). In addition, CEW enhances <EKEvlf> around x 5 260 m, where intensive depth-induced breaking occurs with the large b (Figure 2b). The <EKEtot> proﬁles are more or lesssimilar to those of <EKEvlf> (not shown). Thus, the EKEvlf to EKEtot ratios are mostly unaltered in x with themodest seaward decreasing trend in the offshore region (x > 500 m), where rip currents barely reach. The
EKEvlf to EKEtot ratio ranges from 64.3 to 71.1% with a moderate decrease by the CEW. The mean ﬂuctuation
period Tm2 monotonically increases seaward with the spatial averages ranging from 15.4 to 18.9 min, thus
the predominant variability is in the VLF band. In contrast, the EKEvlf to KEbt ratio considerably increases
from 32.2 and 34.3% without CEW (NC) to 46.1 and 46.5% with CEW.
Table 2. Volume-Averages of the Variables Plotted in Figure 9 Averaged forx  600 m

The above result manifests that CEW intensiﬁes overall variability as well as that in the VLF band in the transient rip currents on the realistic bathymetry either with the 3-D or 2-D model. Furthermore, the impact of
CEW on EKEvlf emerges as a modiﬁcation of the cross-shore proﬁle of EKEvlf. The intensiﬁcation by CEW is
most evident where b is magniﬁed. As argued in UMR09, in the barotropic shear waves on an idealized
bathymetry, CEW decelerate the intrinsic VLFs associated with shear waves due to the wave ray bending
and the resultant magniﬁcation of the shoreward breaker momentum that counteracts with the seaward
momentum transport of the background alongshore currents. The plausible underlying mechanism for thetransient rip currents is that the shoreward contraction effect due to CEW [cf. Yu and Slinn, 2003] hinders
seaward transport of EKE to accumulate EKE within the surf zone (Figure 9b). Several other mechanisms are
possible causes of the difference from the shear wave case. For instance, the broadband nature of the turbulent rip currents generated on the complex, realistic bathymetry (Figure 6) may alter the CEW responses
because waves feel current ﬂuctuations comparable to or greater than their own wavelength according to
the ray theory. Alignments of the principal direction of the background rip currents and incoming waves
may provoke their mutual interaction through the vortex force that can lead to the vorticity reduction andmagniﬁcation in the neck region (section 4.2).
4.4. Bathymetric Variability on Intrinsic VLFs
In literature, intrinsically unstable rip currents have hardly been explicit in the modeling studies on relatively simpliﬁed bathymetries forced by steady, uniform incident waves [e.g., Yu and Slinn, 2003; Weiret al., 2011; Kumar et al., 2012], except for Geiman and Kirby [2013] as argued in section 4.2. In this subsection, a quantitative assessment is performed to examine inﬂuences of the bathymetric variability on the
intrinsic VLFs. A set of synthetic bathymetries is introduced based on the surveyed bottom proﬁle hðx; y Þ
as depicted in Figure 1. We deﬁne a parameter a that controls the alongshore variability in the following
form:hðx; y Þ5hhiðx Þ1ah0 ðx; y Þ; h0 ðx; y Þ5hðx; y Þ2hhiðx Þ;(5)where hhiðx Þ: the alongshore-averaged bed proﬁle, and h0 ðx; y Þ: the perturbation relative to hhiðx Þ. The
original, surveyed proﬁle is retrieved with a 5 1, while the alongshore-uniformity is achieved by a 5 0. We
carry out the full 3-D simulations with CEW for a 5 0, 0.01, 0.1, 0.2, 0.5, 1, 1.2, and 1.5. Note that the sea
bed partially pokes out to the sea surface during the simulations when a  1.7. The other conditions,
including the offshore incident waves, are unchanged as conﬁgured in the previous sections. This aparameter is contexurized with the following bulk metric that measures alogshore bathymetric
nonuniformity:

UCHIYAMA ET AL.

3-D TRANSIENT RIP CURRENTS

5837

Journal of Geophysical Research: Oceans

10.1002/2017JC013005

Figure 10. Same as Figure 9, but for the 3-D runs with CEW on several bed proﬁles, controlled by the perturbation parameter, a (equation(5)). Proportionality to a is pronounced in KEbt, EKEvlf, and EKEvlf/EKEtot although Tm2 are inversely proportional to a.

v2 5

1
Lx Ly

ð Lx ð Ly00hðx; y Þ2hhiðx Þ 2
dy dx;
hhiðx Þ(6)

where Lx and Ly are the cross-shore and alongshore lengths of the domain. For the Duck beach, Ruessink
et al. [2001] and Feddersen and Guza [2003] reported v2 ranges 0.002–0.021 with the average of
0.0038 6 0.008. Our original bathymetry (a51) gives us v250.01346, from which we can evaluate
v2 5 0.01346 a2 for any a. Thus, a  0.53 corresponds to the averaged Duck condition, while a 51 ﬁts the
typical nonumiform regime that has frequently been observed.
Propotionality to a is manifested in the crosshore proﬁles of <KEbt> (Figure 10a) and <EKEvlf> (Figure 10b).
They both are positively correlated with a (Figures 11a and 11b) and well approximated by cubic polynomials while retaining the similarity in the cross-shore proﬁles. Unless a is too small (viz., a  0.1), the ratios ofEKEvlf/EKEtot are about 60% on average (Figure 11c) and do not vary that much with moderate shoreward ampliﬁcations with increased a (Figure 10c). Similarly, the ratios of EKEvlf to KEbt are about
40%. Note that a 5 0 does not generate any rip currents and merely is in the steady cross-shore momentum
balance between the breaker acceleration and the pressure gradient force associated with wave set-up g@x
f  sxsur =ðq0 DÞ5b kx =ðq0 rDÞ and the cross-shore anti-Stokes mass balance, @x ulbt 50 (cf. UMR09, UMS10)
with vbt 50, leading to a neglegibly small <KEbt> and zero EKEs.
In contrast, the spectral mean intrinsic VLF period <Tm2> is inversely proportional to a with a seaward
increasing trend (Figure 10d), ranging from 13 to 36 min on average (Figure 11d). The linear stability analysis of Haller and Dalrymple [2001] indicates the transient rip currents are viewed as the classical Bickely jet[e.g., Marchesiello et al., 2015]. The corresponding oscillation frequency fHD may be determined by the wavenumber k 50.639 and frequency r 50:255 for the fastest growing sinuous mode as:

Figure 11. Scatter plots of the time- and volume-average (denoted by  ) of (a) KEbt, (b) EKEvlf, (c) EKEvlf/EKEtot, and (d) Tm2 as functionsof a (open circles) with the estimated cubic polynomial regressions (dashed curves) and the error bars. The ratios of EKEvlf/KEbt aresuperposed with cross marks in Figure 11c. The cross-shore averages are taken for x  400 m.

UCHIYAMA ET AL.

3-D TRANSIENT RIP CURRENTS

5838

Journal of Geophysical Research: Oceans
fHD 5

r 2U
U
 0:081  ;2p dd10.1002/2017JC013005(7)

where d: width of the jet, and U: jet velocity. In addition, a universal Strouhal number St 5fB d=U  0:181[Bearman, 1967] may alternatively be diagnosed because it is known to explain well the oscillating frequency fB due to wakes and shear turbulence even in geophysical ﬂuid for a wide range of Reynolds number. If we assume d 5 50–200 m from the alongshore bathymetric variability (Figure 1) and jet widths(Figures 3 and 5) and U 5 0.2–0.8 m/s (Figure 6), then fHD gives the period of 13–204 min, while fB gives 5–
92 min. Although the latter agrees moderately better with the model result, both the frequency estimatessuggest that the transient rip currents and their intrinsic VLFs are initiated by shear instability (see also section 4.1; Figure 7) and developed as shear turbulence. Consistent with Figure 10d, fHD and fB exhibit that theoscillation period becomes longer as U decreases, which is the case for the seaward decaying rip currents
and for the attenuated KEbt due to the decreased bathymetric perturbation, a (Figure 11a). Therefore,bathymetric variability inherent in nature is viewed as essential to the onset and maintenance of the intrinsic VLFs in the transient rip currents even under time-independent, alongshore-uniform incident offshore
waves, while the group waves and/or random waves introduce additional turbulent motion as the forced
VLFs. In particular, the alongshore variability is responsible for the intrinsic VLFs through segmentation of
rip cells due to the distributed curl B ﬁeld (section 4.2).
4.5. Steady Forcing Magnitude on Intrinsic VLFs
To examine responses of the intrinsic VLFs to forcing magnitudes, the 3-D model is run for several steadyincident wave height Hrms ranging from 0.1 to 2 m at the offshore open boundary with the other conditions
kept unaltered. In the inner surf zone (x < 340 m), <KEbt>, <EKEvlf>, and <EKEvlf/EKEtot> are evidently proportional to Hrms for Hrms < 1.2 m, whereas they decrease for Hrms > 1.2 m but are enhanced in the offshore(Figures 12 and 13). This outcome is qualitatively consistent with the 2-D model results on much simpler
rip-channeled beaches in Yu and Slinn [2003]. The transition of the intrinsic VLF response to the incident
Hrms is related to the wave breaking around the inner bar crest at x 5 240 m (Figure 2c) that hinders Hrms
and Eb/q to develop <KEbt> and <EKEvlf> in the inner surf zone for Hrms > 1.2 m. Instead, the offshore
breaking on the seaward face of the outer bar at x 5 400 m is invoked in the higher wave conditions, which
eventually takes over the breaking on the inner bar for Hrms 5 2.0 m. The breaking around the inner- andouter-bars and associated shoreward momentum are known to suppress the offshore momentum transport(e.g., UMR09), resulting in the reduced KEbt and EKEvlf for higher waves (Figure 13). As expectedfrom the discussion in section 4.4, the mean intrinsic VLF period is inversely proportional to Hrms for smaller
waves, while it is approximately constant at 18 min for Hrms > 0.8 m (Figure 13d). Note that the intrinsic
VLF responses to the incident waves are anticipated to moderately vary with a choice of the wave breaking
and bottom drag parameterizations [e.g., Yu and Slinn, 2003; Weir et al., 2011].

Figure 12. Cross-shore proﬁles of alongshore-averaged (a) total barotropic kinetic energy KEbt, (b) barotropic eddy kinetic energy in the
VLF band EKEvlf, (c) RMS wave height Hrms, and (d) kinematic breaking acceleration Eb/q for the 3-D CEW runs forced by the different incident Hrms in m indicated in the legend in Figure 12b. The two vertical dashed lines denote the approximate offshore extents of the inner(x < 340 m) and outer (340  x < 500 m) surf zones. KEbt and EKEvlf are proportional to Hrms for Hrms < 1.2 m, although the offshore breaking
in Figure 12d predominates to rather suppress the VLFs in the surf zones.

UCHIYAMA ET AL.

3-D TRANSIENT RIP CURRENTS

5839

Journal of Geophysical Research: Oceans

10.1002/2017JC013005

Figure 13. Same as Figure 11, but as a function of the imposed offshore incident wave height, Hrms.

Because the Duck bathymetry is known to be quite variable [e.g., Lippmann and Holman, 1990] due tomorphodynamic processes primarily induced by waves and littoral currents, a qualitative comparison
between the present model and the observation is a nontrivial task. In the present study, we take advantage of the estimation conducted by MacMahan et al. [2010b], who evaluated in situ VLFs on the same
beach at Duck, NC, using the data set for 13–15 November 1997, during the SandyDuck ﬁeld experiment.
The energetic VLFs of rip currents were intentionally extracted by the conditional sampling for the offshore Hrms > 0.5 m and hp ranging between 29 and 138 to avoid strong alongshore currents. The VLFs aredecomposed into the irrotational and rotational parts: the forced VLFs due to infragravity waves arisingfrom group waves and the intrinsic VLFs relevant to the present study. The resultant intrinsic VLFs are proportional to Hrms with pronounced magniﬁcation in the surf zone. The estimated intrinsic EKEvlf in the surfzone ranges 0.04–1.81 3 1022 m2s22 for the local Hrms between 0.2 and 1.36 m, interpreted from equation (1) and Figure 4 of MacMahan et al. [2010b]. The forced VLFs due to infragravity waves were at thesame magnitude as the intrinsic EKEvlf. A quantitatively similar result was reported by de Schipper et al.
[2012], who extended the analysis of MacMahan et al. [2010b] for the much longer period between
August and December 1997. In turn, the ROMS-WEC model yields <EKEvlf> at about 2.5–2.8 3 1022
m2s22 for a 5 1 and Hrms 5 1.2 m (Table 2). Furthermore, <EKEvlf> broadly ranges between 1.89 3 1023
23.23 3 1022 m2s22 for Hrms 5 0.4–1.4 m (Figure 13). These in situ estimates agree with the present
model outcomes, ensuring that the intrinsic VLFs in the transient rip currents are substantially excited bybathymetric variability.

5. Decay of Unforced Eddies
5.1. Spin-Down Experiment
The intrinsically unstable rip currents comprise signiﬁcant three-dimensionality as shown in Figures 3 and 4,
whereas the barotropic analyses in the preceding sections do not fully capture it. To see the 3-D effects in
the rip currents more closely, we introduce a pseudo 2D (P2D) model based on the fully 3-D model by constraining the vertical eddy viscosity Kv not to be below 0.1 m2/s everywhere, while it may lead to inconsistent bottom stress through altered velocity vertical shear. Nevertheless, the P2D model enables us topractically eliminate the depth-dependency from the 3-D case, which mimics the 2-D barotropic model
quite well. Because the wave forcing keeps generating eddies through the vorticity source curl B (section
4.2) and hinders most of the three-dimensionality, we focus on initial development and subsequent decayof rip currents and associated coherent eddies. Hence, we carry out a spin-down experiment by imposing
the wave forcing until t 5 5 min to initiate offshore-directed rip jets and associated vortices, after which theforcing is abruptly turned off to let eddies decay spontaneously. This spin-down test is somewhat similar to
the 2-D numerical experiment for the self-excitation of the offshore shedding of a forced single vortex
dipole in Geiman and Kirby [2013], although we aim to retain three-dimensionality, the wave-induced rip
currents, and the realistic bathymetric information.
Figure 14 is sequential snapshots of the barotropic vorticity q in the 3-D and 2-D (i.e., P2D) cases with andwithout the wave forcing for t  5 min. At t 5 5 min, all the 4 cases show almost identical vorticity ﬁelds
with a rip current and an associated seaward elongated vortex pair around y 5 420 m. In the forced cases
the vortices are readily unstable, forming meanders and roll-up while the rip current develops. Apparently

UCHIYAMA ET AL.

3-D TRANSIENT RIP CURRENTS

5840

Journal of Geophysical Research: Oceans

10.1002/2017JC013005

Figure 14. Sequential snapshots of barotropic relative vorticity q5rH 3 ubt during the initial development of rip currents at t 5 5, 7, 10,
14 and 20 min for (a) the pseudo 2D (P2D) forced, (b) 3-D forced, (c) pseudo 2D spin-down, and (d) 3-D spin-down cases. The verticaldashed lines denote the approximate offshore extents of the inner surf zone (x < 340 m). The wave forcing vanishes after t  5 min for thetwo spin-down cases to let eddies decay spontaneously, whereas the steady forcing is imposed all the time for the two forced cases. The
P2D model is run with constrained eddy viscosity in the fully 3-D model to suppress depth-dependency, mimicking the 2-D barotropic
model. Notice that the forced 3-D vortices are signiﬁcantly elongated in the longitudinal direction, while the unforced 3-D vortices arepromptly collapsed into smaller eddies, as compared to more coherent, long-lived 2-D vortices.

the forced 2-D eddies are more coherent than the forced 3-D eddies, consistent with Figure 5. The 3-D
eddies are then stretched northeastward remarkably at t  14 min to evolve into a ﬁlament rather than circular vortices. In the unforced cases the 3-D eddies evidently break down into small eddies and decay
much more quickly than the 2-D eddies for t  5 min. The unforced 2-D eddies retain coherent structureseven at t 5 20 min, indicating longer lifetimes and larger wavelengths.
5.2. Kinetic Energy Wavenumber Spectra
Spatiotemporal evolution of coherent eddies in transient rip currents are diagnosed with the alongshorewavenumber barotropic kinetic energy (KEbt) spectra, Sbt, at x 5 140, 240, and 360 m, and at t 5 5, 20, and50 min (Figure 15). We focus on the unforced spin-down cases to observe the fate of eddies to avoid thecomplexities associated with the wave forcing that keeps recharging kinetic energy in the lower wavenumbers in the forced cases. At t 5 5 min, the spectral shapes of the P2D and 3-D models are almost identical as intended. The peak wavelength is evident, ranging between 90 and 160 m. In the P2D model, as
the KEbt gradually decreases with time, the peak wavelength increases due to the KEbt transferred to thelower wavenumbers. In contrast, the 3-D model more rapidly loses the KEbt particularly at the lower wavenumbers. The KEbt is transferred to the higher wavenumbers as most notably seen in (c), making mean
spectral slopes gentler than the P2D model. Hence, the 2-D eddies are dominated by the inverse kinetic
energy cascade, resulting from growth and merging of initially smaller nearshore vortices that travel offshore and evolve into coherent rolled-up circular eddies. In contrast, the 3-D eddies readily break down
into smaller eddies with time and further collapse to promote the transfer of the KEbt to the higher
wavenumbers.

UCHIYAMA ET AL.

3-D TRANSIENT RIP CURRENTS

5841

Journal of Geophysical Research: Oceans

10.1002/2017JC013005

Figure 15. Alongshore-wavenumber barotropic kinetic energy spectra Sbt for the P2D (thick solid curves) and 3-D (thin dotted curves)spin-down cases at t 5 5, 20, and 50 min denoted by the legend in Figure 15a. The selected cross-shore locations are at (a) x 5 140 m in
the feeder region, where the shore-face breaking occurs (Figure 2b), (b) x 5 240 m in the inner surf zone (Figure 2c), and (c) x 5 360 m inthe offshore rip-head region. The abscissae are the reciprocal wavelength in m21. Initially, the peak wavelength is evident, ranging
between 90 and 160 m. In the P2D model, the peak wavelength increases due to the kinetic energy transfer to the lower wavenumbers. In
contrast, the 3-D model rapidly loses the kinetic energy in the lower wavenumbers that is transferred to the higher wavenumbers as mostnotably seen in (c), resulting in gentler spectral slopes than the P2D model.

5.3. 3-D Vorticity Analysis
To further examine three-dimensionality of the rip eddies during the spin-down, we conduct a 3-D vorticity
analysis, in which the wave forcing terms (section 4.2) vanish to simplify the vorticity equation. For theincompressive (q5q0 : constant), nonrotating (fc 5 0), unforced 3-D ﬂuid, the vorticity equation for the vertical relative vorticity xz may be reduced to:@t xz 1ðu  rÞxz 5 xx @x w1xy @y w 2xz ð@x u1@y vÞ 1@z ðKv @z xz Þ
|ﬄ{zﬄ}|ﬄﬄﬄﬄﬄﬄﬄﬄﬄ{zﬄﬄﬄﬄﬄﬄﬄﬄﬄ} |ﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄ{zﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄ} |ﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄ{zﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄ} |ﬄﬄﬄﬄﬄﬄﬄﬄﬄ{zﬄﬄﬄﬄﬄﬄﬄﬄﬄ}

tendency

advectionvortex tilting

vortex stretching(8)

diffusion

where x5r 3 u5 xx ; xy ; xz : the vorticity. The ﬁrst RHS terms represent vortex tilting, which translates xx
and xy into xz due to inclination of the vertically oriented vortices. The second RHS term is the vortex stretchingterm that acts to alter xz due to the vertical compaction and expansion of the vortices associated with divergence of the lateral ﬂow in order to keep the circulation constant. The third RHS, the vertical diffusion term, isalso included to investigaete because we intentionally alter the vertical eddy viscosity Kv for the P2D model.
Figure 16 shows a snapshot of x and the tilting and stretching terms at z 5 20.8 m at t 5 10 min (i.e., 5 minafter turning off the wave forcing) when eddies start to diminish. The P2D model successfully reduces thedepth-dependency during the spin down, as xx and xy practically vanish (Figure 16). The horizontal vorticity isvery large in the 3-D model and a factor of 4–5 larger than the vertical vorticity, due primarily to the vertical
shear of the horizontal velocities. It is obvious that the vortex titling and stretching terms in the 3-D case dominate over the 2-D case, promoting the forward KE cascade (Figure 15), where the initially coherent eddies collapse into smaller eddies more promptly than the 2-D case. The depth-dependent rip currents are accompaniedby strong xx and xy as counter-rotating streamwise vortex pairs (Figures 4 and 16b) that easily yield the tilting
effects to modify xz . The stretching term is also signiﬁcant, while it is more intensiﬁed in the 3-D case than the2-D case, because the offshore eddy migration preconditions vortex stretching due to the increased depth.
The vertical cross-sectional plots of all the three RHS terms and the advection term in equation (8) at t 5 10
min are shown in Figure 17 along the transect at y 5 454 m (Figure 16). The P2D model generates thedepth-invariant xz , which is primarily maintained by the balance between the advection and the tendency(not shown). The tilting term vanishes, while the stretching term moderately contributes, for example, toenhance the negative xz centered at x 5 290 m traveling offshore. Evidently, the 3-D model yields considerable depth-dependency in xz and all the terms. Although the primary balance occurs between the advection and tendency terms, the depth-variant tilting and stretching terms concurrently contribute to the 3-D
vorticity balance. For instance, the negative xz on the seaward side of the vortex pair at x 5 200 m is

UCHIYAMA ET AL.

3-D TRANSIENT RIP CURRENTS

5842

Journal of Geophysical Research: Oceans

10.1002/2017JC013005

Figure 16. Snapshots of the subsurface vorticity xx, xy, xz, the vortex tilting and stretching terms in the xz conservation equation (equation (8)) at t 5 10 min (i.e., 5 min after vanishing
the wave forcing) at z 5 20.8 m for (a) the pseudo 2D and (b) 3-D spin-down cases. Horizontal velocity vectors uH at z 5 20.8 m subsampled every 4-grid points are overlaid (green).
The vertical dashed lines denote the approximate offshore extents of the inner surf zone (x < 340 m). The horizontal magenta lines deﬁne the cross-section at y 5 454 m used in
Figure 17. The constrained vertical eddy viscosity in the pseudo 2D model successfully diminishes xx and xy components of the 3D vortices arising from the depth-dependency, leading
to negligibly small xz tilting effects. Notice that the 3-D rip current is associated with counter-rotating longitudinal vorticity composed of xx and xy.

enhanced by the tilting and stretching terms near the surface, whereas reduced by the positive stretchingnear the bed, resulting in the seaward inclination of the vortex tube as shown by the advection. As a consequence, xz is converted into xx and xy to provoke the vertical circulations, and xz is diminished accordingly. In other words, the vortex tilting helps to make the 2-D eddies (xzxx ; xy ) transform the 3-D
eddies (xz  xx ; xy ) to readily collapse into smaller eddies due to the forward energy cascade. In contrast,
the diffusion term plays a minor role in the vorticity budget, even though the present model does notaccount for spatiotemporal evolution of Kv due to such as the offshore transport of EKE [e.g., Kumar and Feddersen, 2017a,2017b]. The nonwave, shear-induced Kv proﬁles evaluated by the standard KPP model provide
1–6 3 1023 m2/s (cf. Figure 3a in UMS10) with the incremental near-surface intensiﬁcation of an additional
7 3 1023 m2/s due to the wave breaking (cf. Figure 3c in UMS10). The parabolic vertical shape of the

Figure 17. Snapshots of xz and the terms in the unforced 3-D vorticity equation (8) at t 5 10 min along the cross section at y 5 454 m deﬁned in Figure 16 for (a) the pseudo 2D and (b)3-D spin-down cases. The P2D model yields the depth-invariant xz and the negligibly small tilting effects. The 3-D vorticity is generally depth-dependent, and so are the associatedvorticity terms. The diffusion terms for both models are relatively minor and do not play an essential role.

UCHIYAMA ET AL.

3-D TRANSIENT RIP CURRENTS

5843

Journal of Geophysical Research: Oceans

10.1002/2017JC013005

Figure 18. Time series plots of the volume-averaged (a) kinetic energy and (b) enstrophy of xz, and the lateral standard deviations of thedepth-averaged terms in the xz equation (8): (c) the tendency, (d) the advection, (e) the vortex tilting, and (f) the vortex stretching, for the
3-D (red) and pseudo 2D (blue) cases. Bold and thin lines represent the spin-down and forced cases. The vertical dotted lines are the timing when the forcing is shut off at t 5 5 min. The volume averages and lateral standard deviations are taken for x  600 m. The 3-D cases
indicate pronounced effects from the vortex titling and stretching on xz. The former is more substantially reduced than the latter in theP2D cases.

nonwave Kv is retained with a gradual decrease of its magnitude even after the forcing is shut off. The
imposed lower limit of 0.1 m2/s for the Kv in the P2D model acts to vertically redistribute the vorticity andconsequently, reduces the diffusion term.
To extend the above arguments to the temporal evolution of the vorticity balance, the 3-D vorticity terms
for both the forced and the unforced spin-down cases are depicted in Figure 18, along with the volumeaverages of the barotropic kinetic energy KEbt and enstrophy of xz . Note that the vorticity terms (Figures
18c–18f) are represented as the lateral standard deviations of the depth-averages of each term, as proxies
for the volume averaged magnitudes. The 2-D eddies are more energetic and long-lived as depicted in the
KEbt, the enstrophy, and the tendency term regardless of wave forcing. Although the differences of
the advection terms are subtle for each case, the tilting and stretching terms are concurrently intensiﬁed inthe 3-D cases and moderately contribute to the vorticity tendency. In particular, the tilting term is apparent
in the 3-D cases as intended, leading to the faster decay of the 3-D eddies. The tilting is initially magniﬁedby t 5 10 min, and subsequently, gradually decreases for about 20 min for the 3-D spin-down case, which is
sufﬁciently long and is comparable to the mean VLF period of about 18 min (section 4.3) to interact with
the VLFs. For the forced cases, each of the advection and stretching terms ﬂuctuates similarly in the 2-D
and 3-D models. The stretching term is ampliﬁed 25% in the 3-D forced case relative to the 2-D, due tothe vortex expansion and compaction, which eventually is compensated for by the tendency term.

6. Conclusions
The ROMS-WEC model [Uchiyama et al., 2010] based on an Eulerian wave-averaged vortex-force asymptotic
theory of MRL04 is applied to analyze the transient wave-driven 3-D rip currents that spontaneously occur
on a surveyed sandy-beach bathymetry at Duck, NC, under spatiotemporally uniform offshore incident
waves. Substantial depth-dependency is demonstrated in the vertical proﬁles of the rip currents that concentrate at 1–3 m below the surface due to the surf-zone recirculation associated with the surface-conﬁned
breaker acceleration and the near-bed undertow. The 3-D rip currents are composed of pairs of counterrotating longitudinal overturning roll cells that promote surface convergence along the rip center. A vorticity budget analysis exhibits that the vortex force torque plays the most dominant role other than the breaker acceleration torque that is the primary source of vorticity and, thus, the eddy kinetic energy (EKE)
production responsible for the rip eddy shedding. Overall, the vortex force torque preconditions vorticity
reduction from the rip eddies, while it conversely recharges vorticity to the rip eddies in the shoreward of
the wave breaking points where reformation of the broken waves occurs.

UCHIYAMA ET AL.

3-D TRANSIENT RIP CURRENTS

5844

Journal of Geophysical Research: Oceans

10.1002/2017JC013005

These rip currents are intrinsically unstable either in the 3-D or 2-D models. The ﬂuctuations of the transient rip currents are signiﬁcant, contributing about 70% of kinetic energy (KE) as EKE. The contribution to the generation of
the EKE and the Reynolds stress from the depth-variant velocity is evaluated for 5–6%. The dominant ﬂuctuation
period is estimated at about 18 min, but the energy is extensively distributed in broadband in the frequency space.
If we deﬁne our VLF (very low frequency) band for the periods between 5 and 100 min, the EKE contained in the
VLF explains about 70% of the total EKE. The modeled intrinsic EKE in the VLF band is quantitatively comparable tothat of the rotational VLFs observed on the same beach [e.g., MacMahan et al., 2010b; de Schipper et al., 2012], clarifying that bathymetric variability is the main source of the intrinsic VLFs. In addition, the VFLs due to group waves€[e.g., Long and Ozkan-Haller,
2009] and individual breaking waves [e.g., Spydell and Feddersen, 2009], both of which
are neglected in the present model, are anticipated to further excite the forced VLFs. The mean intrinsic VLF period
is consistent with those for the sinuous mode Bickely jet and the Karman vortex streets in wakes and shear turbulence, which suggests shear instability primarily controls the intrinsic VLF variability. Supplemental numericalexperiments illustrate that the feedback mechanism of the current effect on waves (CEW) [cf. Uchiyama et al., 2009]
alters not only the mean rip current structure but also the associated turbulence ﬁeld. The former has been wellknown as the attenuation effects of CEW on rip currents [cf. Yu and Slinn, 2003; Weir et al., 2011]. The latter modiﬁes
the cross-shore EKE proﬁle prominently with considerable accentuation of the eddying motions near the shore(feeder region) and in the inner surf zone centered at the inner bar (neck region), owing to shoreward contraction
effects of CEW that hinder the offshore EKE transport. Moreover, alongshore bathymetric variability is explicitly
demonstrated as an important factor that intensiﬁes KE and EKE with increased alongshore undulations.
The 2-D model generates more coherent, organized, long-lived, and circular eddies. The 3-D eddies arerather stretched laterally, resulting in elongated ﬁlaments and streaks. The underlying mechanism is examined through a spin-down experiment, where a pseudo 2D model is conﬁgured based on the 3-D modelwith an incremental background vertical eddy viscosity for quantitative comparisons with the fully 3-Dmodel. Because the horizontal vorticity components are inherent in the 3-D transient rip currents, the vortex
tilting effect enters vorticity budgets and plays an inevitable role, together with the vortex stretchingeffects, in collapsing 3-D eddies to enhance forward kinetic energy cascade. As a consequence, the 3-D
eddies are short-lived and decay more quickly than the 2-D eddies.

Appendix A: The ROMS-WEC Model
The prognostic variables in the ROMS-WEC are the wave-averaged 3-D Eulerian velocity u 5 (uH, w) 5 (u, v,w) and the wave-averaged sea surface elevation fc . The former is related to the 3-D Lagrangian velocity
ul 5 (uHl, wl) 5 (ul, vl, wl) and the nondivergent 3-D Stokes drift velocity uSt5 (uHSt, wSt) 5 (uSt, vSt, wSt) by
ul 5 u 1 uSt. All the governing equations are expressed in horizontal curvilinear coordinate and verticallystretched surface- and terrain-following s coordinate [Shchepetkin and McWilliams, 2005], which is deﬁned
by fc and the depth h with the transformed grid-cell thickness Hcz . The wave-averaged continuity equation
is expressed in the horizontal Cartesian and the vertical s coordinates for the brevity:@t Hcz 1rH  Hcz ulH 1@s xls 50;(A1)

where xls : the vertical motion passed s surface related to wl through equation (32) of UMS10; andrH 5 @x ; @y : the horizontal gradient operator. Similarly, the horizontal momentum equation is@t Hcz uH 1@x Hcz ul uH 1@y Hcz v l uH 1@s xls uH
h
i
w
~
;

r
1F1F
5Hcz 2fc ^z 3 ulH 2rH /c 1u~St
u
H
H
H

(A2)

where fc: the Coriolis parameter; ^z : the vertical unit vector pointing upward; and F: nonwave nonconservative bodyforces including parameterized momentum mixing terms. The tilde convention identiﬁes conjoined horizontal vectors in a dot product. Fw is the wave-induced nonconservative body forces from depth-induced breaking, surface
rollers, wave-induced boundary streaming, and the incremental breaking-enhanced vertical mixing in the KPP
model [Large et al., 1994; Durski et al., 2004]. In the present study, however, the ﬁrst three contributions are alternatively imposed as the equivalent surface and bottom stress boundary conditions of uH (Appendix C) instead of as the
body forces for satisfactory reproducibility with in situ data even though vertical grid resolution is not ﬁne enoughwith fewer tuning parameters (n.b., UMS10, sections 3.3, 4.4, and 4.6). The vertical momentum equation appears in
the wave-averaged hydrostatic balance, leading to the wave-enhanced geopotential function /c :UCHIYAMA ET AL.

3-D TRANSIENT RIP CURRENTS

5845

Journal of Geophysical Research: Oceans

10.1002/2017JC013005

ð0


qg
/c 5g fc 2^f 2ðP2KÞjfc 1
2K Hcz ds;
s q0

(A3)

where g: the gravity; q: the density of the sea water; q0 : the reference density (1027.5 kgm23); ^f: the quasistatic surface elevation known as wave set-down; P: the pressure correction from the interaction between
vertical current shear and depth-dependent Stokes drift; K: the Bernoulli head after subtracting ^f; and K:
the vertical vortex force (n.b., MRL04, section 9.6; UMS10, sections 2.1 and 2.2).

Appendix B: The Barotropic ROMS-WEC Model
The barotropic version of ROMS-WEC model has previously been developed and used in UMR09 for 2-D shear waves
with a vortex force formalism and its prototype was provided in Uchiyama and McWilliams [2008] in the context of
excitation and propagation of the oceanic infragravity long waves. The governing continuity and horizontal momentum equations of the wave-averaged barotropic model are expressed in the horizontal Cartesian coordinate as:@t fc 1rH  Dulbt 50;(A4)@t ubt 1ðubt  rH Þubtssur sbot
52fc ^z 3 ulbt 2grH fc 2^f 2^z 3 uSt
2;bt ðrH 3 ubt Þ1
q0 D q0 D

(A5)



St St
where ubt 5ðubt ; vbt Þ: depth-averaged, barotropic 2-D Eulerian current velocity; uSt
bt 5 ubt ; vbt : the depthlllSt
averaged Stokes transport; ubt 5 ubt ; vbt 5ubt 1ubt : barotropic 2-D Lagrangian current velocity; D5h1fc :the wave-averaged height of the water column; ssur : the surface stress; and sbot : the bed stress. The prognostic variables of the barotropic 2-D ROMS-WEC are ubt and fc . In the present study we discretize the barotropic advection term with a third-order upwind scheme equivalent to that used in the ROMS. Hence, themomentum equation slightly differs from that in UMR09 in which the centered advection scheme requires
an additional lateral harmonic mixing term to suppress dispersive artifacts (n.b., UMR09, equation (3)).

Appendix C: Boundary Conditions
Following UMS10 (section 3.3), the bed stress sbot used in the bottom boundary condition for uH is given as
cd
the sum of the bottom streaming stress swd
bot and the combined wave-current bottom drag sbot . In the precd
sent study we employ a formula for sbot , proposed originally by Feddersen et al. [2000], who optimized it forthe data set obtained during the SandyDuck experiment in Duck, NC.
 1=3 " #1=2kajubt j 2
1:351
juwrms jubt ;k10:015 q0
juwrms jrD
cd
sbot 5swd
bot 1sbot 5

wd

(A6)

where wd : the wave bottom drag dissipation [Thornton and Guza, 1983; UMS10, section 3.2]; r: the intrinsic
wave frequency; k5 kx ; ky : wavenumber vector of the primary wave whose magnitude is k; ka: the apparpﬃﬃﬃ
ent roughness (50.0125 m) [Ruessink et al., 2001]; uwrms 5uworb = 2: the RMS wave orbital velocity at the bed;
uworb 5rHrms =sinh kD  k=k: the wave orbital velocity at the bed; and Hrms : the RMS wave height. For better
reproducibility with the observation, we alternatively employ the bottom-most horizontal velocity insteadof ubt in the 3-D ROMS-WEC model to best ﬁt the ﬁeld data (not shown). Notice that Feddersen et al. [2000]
only parameterized the quadratic velocity moment for the alongshore momentum balance, while the€model applies to both horizontal balances as pointed out in Long and Ozkan-Haller[2009].
The surface boundary condition for uH is similarly given by the surface stress ssur as the sum of swind
sur : surface
wind stress, and sbsur : the stress-form of the nonconservative WEC body forces from depth-induced wave
breaking b and the surface roller r (see UMS10, section 3.2).
bwindssur 5swind
sur 1ssur 5ssur 1ð12ar Þb 1r
k;
r

(A7)

where ar : the fraction of b converted into surface rollers.
All the numerical simulations are carried out in rectangle domains with the shore to the west and the offshore open boundary to the east in the x direction. The alongshore-periodic condition is imposed in the y

UCHIYAMA ET AL.

3-D TRANSIENT RIP CURRENTS

5846

Journal of Geophysical Research: Oceans

10.1002/2017JC013005

direction. At the shoreline represented by the land-mask function, there is a zero cross-shore mass ﬂux: u1
St
StuSt 50 and ubt 1uSt
bt 50 with u  0 and ubt  0 as for strong dissipation by wave breaking; no tangential
ﬂow: v50 and vbt 50; and the Neumann condition: @x f50. Chapman-type radiation boundary conditions
are adapted for fc , u and ubt at the offshore open boundary with weak nudging toward ^f1f, 2uSt and
2uSt
bt , where f is a nonwave free surface elevation such as due to tides, which is assumed time-invariant in
the present study (Table 1). The Neumann condition is used for all the other variables at the offshore
boundary. The minimum water depth is set to 0.1 m, and movement of the shoreline, due to wetting anddrying, is not considered.

Appendix D: The WKB Spectrum Peak Wave Refraction Model
For a narrow-banded, spectrum peak wave ﬁeld, we exploit a set of ray equations based on the WKB
approximation [e.g., Mei, 1992] for refractive wave transformation by the conservation equations of wave
action density A5q0 gH2rms =8r and the wavenumber k as1
@t A1rH cg 1ubt A 52 b 1wd ;r~  rH u~bt 2@t k1 cg 1ubt  rH k52k

kr
rH D;
sihn 2kD

(A8)
(A9)

where cg : (non-Doppler shifted) group velocity, and the tilde convention is the same as that in equation(A2) in Appendix A. With the linear dispersion relation, the Doppler-shifted wave frequency rds is represented aspﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ
rds 5ubt  k1r5ubt  k1 gktanh kD:(A10)
To better estimate the surf zone wave ﬁeld, the surface roller model [Svendsen, 1984; UMS10, section 3.2]
can be exploited.
@t Ar 1rH ½ðc1ubt ÞAr 52

1 r
 2ar b ;
r

(A11)

where Ar : roller action density; c5rk 22 k: the phase velocity of the primary wave (notice not cg [Svendsen,
1984], but c cg anyway for kD ! 0); and r 5gkAr sinb: a parameterized roller dissipation rate. According
to Nairn et al. [1991] and Apotsos et al. [2007], we rely on ar 51 and sinb50:1. The depth-dependent Stokes
drift and its depth-average are then expressed as
uH St ðz Þ5

r2 cosh 2k ðz1hÞ
ðA1Ar Þ
r
St
ð
Þk;
u
5
k :
A1A
bt
q0 D
q0 gsinh2 kD(A12)The breaker dissipation b is determined by a quasi-empirical parameterization for Rayleigh distributed narrow banded waves described in Church and Thornton [1993]:
2 (
3
 


2 )25=2pﬃﬃﬃ
3 p B3b r 3
Hrms
H
rms
5;
H
21
 412 11
 5q0 g11tanh 8
cb D
cb D
16 2pD rms
b

(A13)

where Bb and cb are empirical constants depending on a beach proﬁle and the type of breaker.
The wave model is solved at every baroclinic time step for the 3-D cases, while it is solved at every barotropic time step for the 2-D cases as a subroutine in the ROMS-WEC with the same thread structure and,thus, is identically parallelized. The computational cost of the wave model is negligibly cheaper than the 3D ROMS-WEC, whereas it typically takes about 80% as much as the barotropic current model. We choose
the alongshore-periodic condition and the shoreline Neumann condition for all the wave variables. The offshore clamped conditions for A and k with prescribed incident wave parameters of Hrms , r, and the wave
direction hp are used. The feedback mechanisms of CEW participate in the wave model through the Doppler
shift due to ubt and wave-averaged surface elevation fc , hence in D. Note that the present wave model solutions do not generate caustics or other ray singularities with our particular choices of incident waves andbathymetries.

UCHIYAMA ET AL.

3-D TRANSIENT RIP CURRENTS

5847

Journal of Geophysical Research: Oceans
Acknowledgments
This research is supported by the JSPS
Grants-in-Aid for Scientiﬁc Research
15KK0207 and 15H04049 at Kobe
University, NSF OCE-1355970 and ONR
N00014–15-1–2645 at UCLA, and NSF
XSEDE computing resources.
M. Kawakami of Kobe University
carried out the preliminary simulationsin Secs 4.3 and 4.4. F. Feddersen of
Scripps Institution of Oceanography
provided knowledgeable commentson the manuscripts that improve theclarity of the study. We are grateful for
Ofﬁce of Naval Research, U.S. Army
Corps of Engineers, S. Elgar, R. T. Guza,
T. H. C. Herbers, P. A. Howd,
K. K. Hathaway, W. A. Birkemeier,C. E. Long and other scientists andstaff involved in the SandyDuck
experiment and the bathymetry
surveys around the Field Research
Facility for making their data set
available. The surveyed data used in
the present study are archived at
http://navigation.usace.army.mil/CHL_
Viewer/FRF/. The source code for the
ROMS-WEC model is available from theauthors upon request (uchiyama@
harbor.kobe-u.ac.jp).

UCHIYAMA ET AL.

10.1002/2017JC013005

References
Aagaard, T., B. Greenwood, and J. Nielsen (1997), Mean currents and sediment transport in a rip channel, Mar. Geol., 140, 25–45.
Allen, J. S., P. A. Newberger, and R. A. Holman (1996), Nonlinear shear instabilities of alongshore currents on plane beaches, J. Fluid Mech.,
310, 181 – 213.
Apotsos, A., B. Raubenheimer, S. Elgar, R. T. Guza, and J. A. Smith (2007), Effects of wave rollers and bottom stress on wave setup, J. Geophys. Res., 112, C02003, doi:10.1029/2006JC003549.
Ardhuin, F., N. Rascle, and K. A. Belibassakis (2008), Explicit wave-averaged primitive equations using a generalized Lagrangian mean,
Ocean Modell., 20, 35–60.
Bearman, P. W. (1967), On vortex street wakes, J. Fluid Mech., 28, 625–641.
Bonneton, P., N. Bruneau, B. Castelle, and F. Marche (2010), Large scale vorticity generation due to dissipating waves in the surf zone, Discrete Contin. Dyn. Syst., Ser. B, 13(4), 729–738, doi:10.3934/dcdsb.2010.13.729.
Bowen, A. J. (1969), Rip currents, 1, Theoretical investigations, J. Geophys. Res., 74, 5467–5478.
Bowen, A. J., and R. A. Holman (1989), Shear instabilities of the mean longshore current: 1. Theory, J. Geophys. Res., 94, 18,023–18,030.
Brander, R. W., and A. D. Short (2001), Flow kinematics of low-energy rip current systems, J. Coastal Res., 17(2), 468–481.
Bruneau, N., P. Bonneton, B. Castelle, and R. Pedreros (2011), Modeling rip current circulations and vorticity in a high-energy mesotidalmacrotidal environment, J. Geophys. Res., 116, C07026, doi:10.1029/2010JC006693.
Church, J.C. and E.B. Thornton (1993), Effects of breaking wave induced turbulence within a longshore current model, Coastal Eng., 20,
1–28.
de Schipper, M. A., A. Reniers, J. MacMahan, and R. Ranasinghe (2012), Vortical VLF motions under shore-normal incident waves, inProceedings of 33rd International Conference on Coastal Engineering, edited by P. Lynett and J. M. Smith, International Conference on
Coastal Engineering (ICCE), 13 pp., Santander, Spain, doi:10.9753/icce.v33.waves.58.
Dodd, N., J. Oltman-Shay, and E. B. Thornton (1992), Shear instabilities in the longshore current: A comparison of observation and theory,
J. Phys. Oceanogr., 22, 62–82.
Dodd, N., J. Oltman-Shay, and E. B. Thornton (1994), On the destabilization of a longshore current on a plane beach: Bottom shear stress,
critical condition and onset of instabilities, J. Geophys. Res., 99, 811–824.
Dodd, N., V. Iranzo, and A. Reniers (2000), Shear instabilities of wave-driven alongshore currents. Rev. Geophys., 38, 437–464.
Durski, S. M., S. M. Glenn, and D. Haidvogel (2004), Vertical mixing schemes in the coastal ocean: Comparison 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.
Feddersen, F. (2014), The generation of surfzone eddies in a strong alongshore current, J. Phys. Oceanogr., 44, 600–617, doi:10.1175/JPO-D13-051.1.
Feddersen, F., and R. T. Guza (2003), Observation of nearshore circulation: Alongshore uniformity, J. Geophys. Res., 108(C1), 3006, doi:
10.1029/2001JC001293.
Feddersen, F., R. T. Guza, S. Elgar, and T. H. C. Herbers (2000), Velocity moments in alongshore bottom stress parameterizations, J. Geophys.
Res., 105, 8673 – 8686.
Feddersen, F., D. B. Clark, and R. T. Guza (2011), Modeling surf zone tracer plumes: 1. Waves, mean currents, and low-frequency eddies,
J. Geophys. Res., 116, C11027, doi:10.1029/2011JC007210.
Fujimura, A. G., A. J. H. M. Reniers, C. B. Paris, A. L. Shanks, J. H. MacMahan, and S. G. Morgan (2014), Numerical simulations of larval transport into a rip-channeled surf zone, Limnol. Oceanogr. Methods, 59, 1434–1447, doi:10.4319/lo.2014.59.4.1434.
Geiman, J. D., and J. T. Kirby (2013), Unforced oscillation of rip-current vortex cells, J. Phys. Oceanogr., 43, 477–497, doi:10.1175/JPO-D-110164.1.
Geiman, J. D., J. T. Kirby, A. J. H. M. Reniers, and J. H. MacMahan (2011), Effects of wave averaging on estimates of ﬂuid mixing in the surfzone, J. Geophys. Res., 116, C04006, doi:10.1029/2010JC006678.
Haas, K. A., and I. A. Svendsen (2002), Laboratory measurements of the vertical structure of rip currents. J. Geophys. Res., 107(C5), 3047, doi:
10.1029/2001JC000911.
Haas, K. A., I. A. Svendsen, M. C. Haller, and G. Zhao (2003), Quasi-three dimensional modeling of rip current system, J. Geophys. Res.,
108(C7), 3217, doi:10.1029/2002JC001355.
Haller, M. C., and D. A. Dalrymple (2001), Rip current instabilities, J. Fluid Mech., 433, 161–192.
Hasselmann, K. (1971), On the mass and momentum transfer between short gravity waves and larger-scale motions, J. Fluid Mech., 50,
189–201.
Holman, R., and M. C. Haller (2013), Remote sensing of the nearshore, Annu. Rev. Mar. Sci., 5(1), 95–113.
Kennedy, A. B. (2003), A circulation description of a rip current neck, J. Fluid Mech., 497, 225–234.
Kumar, N., and F. Feddersen (2017a), The effect of Stokes drift and transient rip currents on the inner shelf. Part 1: No stratiﬁcation, J. Phys.
Oceanogr., 47, 227–241, doi:10.1175/JPO-D-16-0076.1.
Kumar, N., and F. Feddersen (2017b), The effect of Stokes drift and transient rip currents on the inner shelf. Part II: With Stratiﬁcation, J.
Phys. Oceanogr., 47, 243–260, doi:10.1175/JPO-D-16-0077.1.
Kumar, N., G. Voulgaris, J. C. Warner, and M. Olabarrieta (2012), Implementation of the vortex force formalism in the coupled oceanatmosphere-wave-sediment transport (COAWST) modeling system for inner shelf and surf zone applications, Ocean Modell., 47, 65–95.
Kumar, N., F. Feddersen, Y. Uchiyama, J. McWilliams, and W. O’Reilly (2015), Midshelf to surfzone coupled ROMS–SWAN model–data comparison of waves, currents, and temperature: Diagnosis of subtidal forcings and response, J. Phys. Oceanogr., 45, 1464–1490, doi:
10.1175/JPO-D-14-0151.1.
Large, W. G., J. C. McWilliams, and S. C. Doney (1994), Oceanic vertical mixing: A review and a model with nonlocal boundary layer parameterisation, Rev. Geophys., 32, 363–403.
Lesser, G. R., J. A. Roelvink, J. A. T. M. van Kester, and G. S. Stelling (2004), Validation and development of a three-dimensional morphological model, Coastal Eng., 51, 883–915.
Lippmann, T. C., and R. A. Holman (1990), The spatial and temporal variability of sand bar morphology, J. Geophys. Res., 95, 11,575–11,590.
€
Long, J. W., and H. T. Ozkan-Haller(2009), Low-frequency characteristics of wave group-forced vortices, J. Geophys. Res., 114, C08004, doi:
10.1029/2008JC004894.
Longuet-Higgins, M. S., and R. W. Stewart (1962), Radiation stress and mass transport in gravity waves, with application to ‘surf beats,’ J.
Fluid Mech., 13, 481–504.
MacMahan, J., A. J. H. M. Reniers, E. B. Thornton, and T. Stanton (2004a), Infragravity rip current pulsations, J. Geophys. Res., 109, C01033,
doi:10.1029/2003JC002068.

3-D TRANSIENT RIP CURRENTS

5848

Journal of Geophysical Research: Oceans

10.1002/2017JC013005

MacMahan, J., et al. (2010a), Mean Lagrangian ﬂow behavior on an open coast rip-channeled beach: A new perspective, Mar. Geol., 268,
1–15.
MacMahan, J. H., A. J. H. M. Reniers, E. B. Thornton, and T. P. Stanton (2004b), Surf zone eddies coupled with rip current morphology, J. Geophys. Res., 109, C07004, doi:10.1029/2003JC002083.
MacMahan, J. H., E. B. Thornton, T. P. Stanton, and A. J. H. M. Reniers (2005), RIPEX: Observations of a rip current system, Mar. Geol., 218,
113–134.
MacMahan, J. H., A. J. H. M. Reniers, and E. B. Thornton (2010b), Vortical surf zone velocity ﬂuctuations with 0(10) min period, J. Geophys.
Res., 115, C06007, doi:10.1029/2009JC005383.
Marchesiello, P., R. Benshila, R. Almar, Y. Uchiyama, J. McWilliams, and A. Shchepetkin (2015), On tridimensional rip current modeling,
Ocean Modell., 96(1), 36–48, doi:10.1016/j.ocemod.2015.07.003.
McWilliams, J. C., J. M. Restrepo, and E. M. Lane (2004), An asymptotic theory for the interaction of waves and currents in coastal waters, J.
Fluid Mech., 511, 135–178.
Mei, C. C. (1992), The Applied Dynamics of Ocean Surface Waves, 768 pp., World Sci., Singapore.
Munk, W. H. (1949), Surf beats, Eos Trans. AGU, 30, 849–854.
Nairn, R. B., J. A. Roelvink, and H. N. Southgate (1991), Transition zone width and implications for modelling surfzone hydrodynamics, inProceedings of 22nd Int. Conf. Coastal Eng. 1990, edited by B. L. Edge, pp. 68–82, Delft, Netherlands, Am. Soc. Civ. Eng., New York.
Newberger, P. A., and J. S. Allen (2007), 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.
Noyes, T. J., R. T. Guza, S. Elgar, and T. H. C. Herbers (2004), Field observations of shear waves in the surf zone, J. Geophys. Res., 109, C01031,
doi:10.1029/2002JC001761.
Oltman-Shay, J., P. A. Howd, and W. A. Birkemeier (1989), Shear instabilities of the mean longshore current: 2. Field observations, J. Geophys. Res., 94, 18,031–18,042.
€
Ozkan-Haller
H. T., and J. T. Kirby (1999), Nonlinear evolution of shear instabilities of the longshore current: A comparison of observations
and computations, J. Geophys. Res., 104, 25,953–25,984.
€
Ozkan-Haller,
H. T., and Y. Li (2003), Effects of wave-current interaction on shear instabilities of longshore currents, J. Geophys. Res., 108(C5),3139, doi:10.1029/2001JC001287.
Phillips, O. M. (1977), The Dynamics of the Upper Ocean. Cambridge Univ. Press, Cambridge, U. K.
Reniers, A. J. H. M., J. A. Battjes, A. Falques, and D. A. Huntley (1997), A laboratory study on the shear instability of logngshore currents,
J. Geophys. Res., 102, 8597–8609.
Reniers, A. J. H. M., J. H. MacMahan, E. B. Thornton, and T. P. Stanton (2007), Modeling of very low frequency motions during RIPEX, J. Geophys. Res., 112, C07013, doi:10.1029/2005JC003122.
Reniers, A. J. H. M., J. H. MacMahan, E. B. Thornton, T. P. Stanton, M. Henriquez, J. W. Brown, J. A. Brown, and E. Gallagher (2009), Surfzone
surface retention on a rip channeled beach, J. Geophys. Res. 114, C10010, doi:10.1029/2008JC005153.
Reniers, A. J. H. M., J. H. MacMahan, F. J. Beron-Vera, and M. J. Olascoaga (2010), Rip-current pulses tied to Lagrangian coherent structures,
Geophys. Res. Lett., 37, L05605, doi:10.1029/2009GL041443.
Ruessink, B., J. Miles, F. Feddersen, R. T. Guza, and S. Elgar (2001), Modeling the alongshore current on barred beach, J. Geophys. Res., 106,
22,451 – 22,463.
Shchepetkin, A. F., and J. C. McWilliams (2005), The Regional Oceanic Modeling System: A split-explicit, free-surface, topography-followingcoordinate oceanic model. Ocean Modell., 9, 347–404.
Shchepetkin, A. F., and J. C. McWilliams (2008), Computational kernel algorithms for ﬁne-scale, multiprocess, longtime oceanic simulations,in Handbook of Numerical Analysis: Computational Methods for the Ocean and the Atmosphere, edited by R. Temam and J. Tribbia,
pp. 119–181, Elsevier, Amsterdam.
Smith, J. A. (2006), Wave–current interactions in ﬁnite-depth, J. Phys. Oceanogr., 36, 1403–1419.
Sonu, C. J. (1972), Field observations of nearshore circulation and meandering currents, J. Geophys. Res., 77, 3232–3247.
Spydell, M., and F. Feddersen (2009), Lagrangian drifter dispersion in the surf zone: Directionally spread, normally incident waves, J. Phys.
Oceanogr., 39, 809–830.
Svendsen, I. A. (1984), Mass ﬂux and undertow in a surf zone, Coastal Eng., 8, 347– 365.
Thornton, E. B., and R. T. Guza (1983), Transformation of wave height distribution, J. Geophys. Res., 88, 5925–5938.
Uchiyama, Y., and J. C. McWilliams (2008), Infragravity waves in the deep ocean: Generation, propagation, and seismic hum excitation, J.
Geophys. Res., 113, C07029, doi:10.1029/2007JC004562.
Uchiyama, Y., J. C. McWilliams, and J. M. Restrepo (2009), Wave-current interaction in nearshore shear instability analyzed with a vortexforce formalism, J. Geophys. Res., 114, C06021, doi:10.1029/2008JC005135.
Uchiyama, Y., J. C. McWilliams, and A. F. Shchepetkin (2010), Wave-current interaction in an oceanic circulation model with a vortex force
formalism: Application to the surf zone, Ocean Modell., 34, 16–35.
Walstra, D. J. R., J. A. Roelvink, and J. Groeneweg (2000), Calculation of wave-driven currents in a mean ﬂow model, in Coastal Engineering,edited by B. L. Edge, pp. 1050–1063, Am. Soc. of Civ. Eng., Reston, Va.
Weir, B., Y. Uchiyama, E. M. Lane, J. M. Restrepo, and J. C. McWilliams (2011), A vortex-force analysis of the interaction of rip currents andgravity waves, J. Geophys. Res., 116, C05001, doi:10.1029/2010JC006232.
Yu, J. (2006), On the instability leading to rip currents due to wave-current interaction, J. Fluid Mech., 549, 403–428.
Yu, J., and D. N. Slinn (2003), Effects of wave-current interaction on rip currents, J. Geopys. Res., 108(C3), 3088, doi:10.1029/2001JC001105.
Zhao, Q., I. A. Svendsen, and K. Haas (2003), Three-dimensional effects in shear waves, J. Geophys. Res., 108(C8), 3270, doi:10.1029/
2002JC001306.

UCHIYAMA ET AL.

3-D TRANSIENT RIP CURRENTS

5849