Click Here

JOURNAL OF GEOPHYSICAL RESEARCH, VOL. 113, C07022, doi:10.1029/2007JC004374, 2008

for

Full Article

Sandbar formation under surface waves: Theory and experiments M. J. Hancock,1,3 B. J. Landry,2,3 and C. C. Mei4 Received 5 June 2007; revised 7 September 2007; accepted 31 January 2008; published 18 July 2008.

[1] We report a combined theoretical and experimental study of sandbar formation under

simple-harmonic surface waves. For coarse grains and weak waves, an established empirical rule of bed load transport is used with an asymptotic theory for the fluid flow. The surface waves are governed by potential theory and a depth-linear eddy viscosity is employed in the turbulent boundary layer at the seabed. The derived bed stress is used to predict the sand-bed evolution. Laboratory experiments and corresponding numerical simulations for both high and low beach reflection are discussed. For weak reflection, the shear stress associated with the return current is found to be important. Partial simulation of a field record in Cape Cod Bay is also described. Citation: Hancock, M. J., B. J. Landry, and C. C. Mei (2008), Sandbar formation under surface waves: Theory and experiments, J. Geophys. Res., 113, C07022, doi:10.1029/2007JC004374.

1. Introduction [2] The presence of almost periodic sandbars along many beaches of gentle slope has been frequently recorded by aerial photography and acoustic soundings [Dolan and Dean, 1985; Elgar et al., 2003; and references therein]. The physical origin and oceanographic consequence of these sandbars are of considerable engineering and scientific interest. A host of sandbar generation mechanisms have been proposed, including vortex action by plunging breakers, steady currents induced by breaking waves, edge waves, harmonic decomposition of shoaling waves, the combined effect of waves and undertow, and partially standing waves [O’Hare and Davies, 1993; Yu and Mei, 2000a; and references therein]. In this article, we focus on longshore bars in intermediate depths, outside the surf zone, generated by non-breaking weakly nonlinear partially standing surface waves. [3] Based on the mechanics of induced streaming in the bottom wave boundary layer, Carter et al. [1973] gave theoretical and experimental evidence that the horizontally varying Reynolds stress in the bottom boundary layer under partially standing waves induces steady circulations which can initiate the accumulation of sand particles at a spacing equal to one half of the surface wavelength. Experiments for rigid bars by Heathershaw [1982] demonstrated that once periodic bars are present, waves twice as long as the bar spacing are strongly reflected by Bragg resonance. Bragg resonance over rigid bars has been explained theoretically 1 Department of Mathematics, Massachusetts Institute of Technology, Cambridge, Massachusetts, USA. 2 Department of Civil and Environmental Engineering, University of Illinois, Urbana-Champaign, Illinois, USA. 3 Formerly at Department of Civil and Environmental Engineering, Massachusetts Institute of Technology, Cambridge, Massachusetts, USA. 4 Department of Civil and Environmental Engineering, Massachusetts Institute of Technology, Cambridge, Massachusetts, USA.

Copyright 2008 by the American Geophysical Union. 0148-0227/08/2007JC004374$09.00

[Davies, 1982; Mei, 1985; Hara and Mei, 1987; Naciri and Mei, 1988]. For sandbars generated by surface gravity waves, Boczar-Karakiewicz et al. [1995] and Restrepo and Bona [1995] proposed theories for long waves in shallow water beneath which the sediment is transported almost entirely in suspension. In shallow water, higher harmonics are produced which interact over the recurrence ¨ nlu¨ata, 1972]. The resulting nonunidistance [Mei and U formities in the mass transport velocity generate sandbars whose lengths are close to the recurrence distance. The coupled evolution of partially standing surface waves and sandbars was first studied by O’Hare and Davies [1993] using the numerical technique of Devillard et al. [1988]. The wavefield was computed numerically by discretizing the bed into small horizontal steps. [4] Focusing on coarse sand, intermediate wavelengths, and constant depth, Yu and Mei [2000a] presented a theory for partially standing wave generated sandbars for cases affected strongly by Bragg resonance. By using a model of constant eddy viscosity, Yu and Mei showed that sandbar formation by surface waves is a process of forced diffusion, rather than instability as in the case for short ripples [Blondeaux, 1990; Vittori and Blondeaux, 1990; Mei and Yu, 1997]. The forcing is caused by nonuniformity in the wave envelope, and hence the bottom shear stress and bed load transport. The diffusivity is due to modifications in local bed stress caused by gravitational forces on sediment grains on a sloping bed. Because of the nonlinearity of the sediment dynamics, both the mean flow (Eulerian streaming) and the second time harmonic contribute to the forcing. This suggests that sandbars can be generated for any finite seaward reflection, as O’Hare and Davies [1993] found numerically and experimentally. Yu and Mei [2000a] have also shown that if shoreline reflection is not present, sandbars (on the scale of the wavelength) cannot be generated and any residual sandbars present initially are eventually washed out. Therefore, under the conditions of the study, finite reflection from the shore is both necessary and sufficient to generate and maintain sandbars. Yu and Mei

C07022

1 of 23

C07022

HANCOCK ET AL.: SANDBAR FORMATION UNDER SURFACE WAVES

C07022

Figure 1. Definitional sketch of sandbar formation under ocean surface waves in water of intermediate depth. The waves and seabed are assumed to be gently sloping. The water column is divided into three distinct regions: the inviscid core, the bottom boundary layer, and the sediment transport on the surface of the seabed. also explained the observation of O’Hare and Davies [1993] that the relative position of bar crests and envelope nodes affects the wave response over the bar patch. [5] Laboratory confirmation of Yu and Mei’s model was limited to existing experiments performed in small tanks which lack details of the slow evolution of both waves and sandbars. The earliest experiments were performed by Herbich et al. [1965] in a long wave flume with a seawall inclined at different angles. Half-wavelength sandbars with ripples superposed on each were noted to form. A detailed account of sandbar geometry and crest and trough locations was not given. Additional observations were made of partially standing waves over a monolayer of sand by Carter et al. [1973]. Several small scale laboratory experiments of sandbar formation under partially standing waves have been made by de Best et al. [1971], Xie [1981], O’Hare and Davies [1993], and Seaman and O’Donoghue [1996]. de Best et al. [1971] were the first to observe that for coarse sand, bar crests tended to form near wave envelope nodes, while for fine sand, crests formed near wave antinodes. Xie [1981] and O’Hare and Davies [1993] later repeated this finding. de Best et al. [1971] also demonstrated sediment sorting, by generating standing waves over a seabed initially consisting of a well mixed sand of two grain sizes. The finer sand was transported toward the antinodes as suspended load, and the coarser sand toward the nodes as bed load. Jan and Lin [1998] observed the formation of sandbars and ripples under oblique standing waves in front of a seawall in a laboratory tank with a horizontal bed. Dulou et al. [2000] studied sandbar formation under partially standing waves on a gently sloping bed in a small wave tank. Scaling effects limit the similarity of these small scale sandbars to those observed in the field [Yu and Mei, 2000a]. Also, fine grains were used so that sediment was transported largely in suspension. Recently, the phenomenon of sediment sorting has been reexamined experimentally in a relatively large wave flume with two sands of significantly different grain sizes [Landry et al., 2007].

[6] In this article, we report an improved theory of sandbar generation under partially standing waves, as sketched in Figure 1. We also report new laboratory tests in a wave flume with an initially horizontal sand bed. Long duration records of free-surface displacements and bed profiles due to waves of different intensity and degrees of reflection are compared with the theory. We focus on monochromatic partially standing waves and coarse grains, so that sediment is transported primarily as bed load. In the theoretical part, a major change over the work of Yu and Mei is the use of a depth-linear and time-invariant eddy viscosity model which makes possible agreement with laboratory data. Established empirical formulas of sediment transport are used to limit the number of fitting parameters. [7] We shall begin by citing relevant empirical information on sediment transport under oscillatory flows. After summarizing the theory of waves in the inviscid core, the turbulent boundary layer at the seabed is analyzed to the second-order in wave steepness. A depth-linear eddy viscosity model is used with a friction factor determined by treating the ripples present on the bed as a known roughness. Combining the calculated bed shear stress with the discharge rule of bed load transport yields a forced diffusion equation for the ripple-averaged bed profile. Three sets of laboratory experiments are then described and compared to numerical simulations of the forced diffusion equation. Finally, the transient evolution of sandbars in Cape Cod Bay is simulated and the predicted bed profile is compared with that recorded by Elgar et al. [2003] after a storm. This provides a partial test of our theoretical model for field-scale applications. [8] For convenience, a list of symbols is given in Appendix E.

2. Rate of Sediment Transport [9] In steady flow over a plane sediment bed, the Shields parameter is usually defined by

2 of 23

Qðx0 ; t 0 Þ ¼

t 0b ; rðs  1Þgd

ð1Þ

HANCOCK ET AL.: SANDBAR FORMATION UNDER SURFACE WAVES

C07022

C07022

where ym is the angle of repose (taken as 30°), h0 is the ripple-averaged seabed profile, and hats denote the maximum amplitude of a time periodic quantity. The factor

3=2

h pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffii b  QC b  QC ð4Þ Q0B ðx0 Þ ¼ 3 ðs  1Þgd 3 rd Q Hv rd Q

is the rate of sediment transport over a half-cycle, where QC denotes the threshold Shields parameter for incipient sand motion, found from the modified Shields diagram, and Hv is the Heaviside unit step function. Let e denote the characteristic wave-steepness, e ¼ k Ab 1; * *

Figure 2. Correction factor rd versus d/hb*. Numbers adjacent to curves indicate kAb*/d. where t b0(x0, t0) denotes the shear stress on the bed surface, r the fluid density, s the ratio of the sand to fluid density, g the acceleration of gravity, and d the sand diameter. Throughout this article superscript primes indicate dimensional variables; the corresponding dimensionless quantities are written without primes. Characteristic scales are distinguished by subscript asterisks. [10] In uniform but oscillatory flows of a single frequency w, with horizontal velocity U / cos wt0, the Shields parameter is monochromatic at leading order, Q / cos(wt0 + v), and Sleath [1978] found the bed load transport rate to be proportional to the fourth power of the time factor, i.e., jcos(wt0 + v)j3 cos(wt0 + v). The phase v is the lag between the shear stress and the flow. Thus, the bed load transport rate is proportional to jQj3Q. Taking the time average over a wave period 2p/w, denoted by overlines, it can be shown that j cosðwt 0 þ vÞj3 cosðwt0 þ vÞ ¼ 0:

ð2Þ

[11] Hence it follows that q0B vanishes at order unity and the net bed load transport rate after one period is at most of order O(e). A theory for sandbars generated under monochromatic waves must therefore consider higher order corrections to the velocity field and sediment transport rate. [12] Under sinusoidal flows, Sleath [1978] and Nielsen [1992] also found it necessary to introduce the skin-friction Shields parameter, which is the fraction of the bed friction due to grain roughness, by replacing t b0 by rd t b0 and Q by rdQ. The correction factor rd depends empirically on the ratio of the sand diameter d to the characteristic roughness height hb* as plotted in Figure 2 and explained in Appendix A. [13] On an inclined plane with finite slope, a further correction factor has been proposed by King [1991] to account for the effect of gravity. The modified rate of sediment discharge is " q0B

¼

1 1  tan1y

m

# @h0 Q @x0 jQj

  8 0 Q4 Q Q ; b  jQj 3 B Q

ð5Þ

where k* is the characteristic wave number and Ab* is the orbital amplitude just above the bed. With sufficient generality we assume that the ripple-averaged bed slope is comparable to the wave-steepness, i.e., @h0/@x0 = O(e) 1. Taking the time average of (3) over a wave period, denoted by overlines, gives the period-averaged discharge rate  8 jQj3 Q 8 Q0B jQj4 @h0 q0B ¼ Q0B þ þ O e2 : 0 4 4 b b 3 3 tan ym Q @x Q

ð6Þ

We have used the fact that the ripple-averaged bed profile h0 changes slowly in time relative to the wave period, hence depends on a slow time variable
t 0 , i.e., h0 = h0(x0,
t 0 ), where
t 0 =t 0 1. It is shown later that
t 0 = O(e6t0). [14] The bed load transport is parallel to the rippleaveraged bed. Since the slope of the bed is mild, @h0/@x0 = O(e) 1, the time period-averaged law of sediment conservation reads ð1  N Þ

 @h0 @q0B  0 ¼ O e2 ; 0 @t @x

ð7Þ

where N denotes the porosity. With (6), equation (7) leads to a diffusion-like equation for the ripple-averaged depth h0(x0,
t 0 ) below the still water level. [15] By reasoning that gravity reduces (increases) the threshold stress if the fluid moves upward (downward) along a slope, Fredsøe [1974] proposed a relation which can be reduced to the form of (6), with an extra factor of QC 0.05 in the second term, and has been used by Yu and Mei [2000a] and by Komarova and Newell [2000] in a theoretical model of tidal bars. [16] In order to minimize the subsequent algebra, it is advantageous to first estimate which parts of the velocity field in the boundary layer are needed to calculate the period-averaged sediment transport rate. In view of the assumptions of small wave steepness and bed slope, characterized by the same small parameter e, we anticipate that the longitudinal velocity in the seabed boundary layer can be expanded as

ð3Þ

3 of 23

 u0 ¼ u01 þ eu02 þ O e2 :

ð8Þ

HANCOCK ET AL.: SANDBAR FORMATION UNDER SURFACE WAVES

C07022

[17] Because of nonlinearity, the time dependence of the first two orders should be

0 ½10 u01 ¼ < u1 eiwt ;

ð9Þ



0 0 ½00 ½10 ½20 u02 ¼ < u2 þ u2 eiwt þ u2 e2iwt ;

ð10Þ

C07022

due to the spatial variation of bed shear stresses. A similar equation was derived earlier by Yu and Mei [2000a]. [19] What remains is to calculate the shear stress components present in (15) and (16). We first calculate the velocity field by separating the water column into an inviscid core and a thin boundary layer of thickness O(d*) just above the seabed.

3. Waves and Current in the Inviscid Core where the superscript [n] identifies the particular harmonic, and the integer subscript indicates the order of e. The corresponding bed stress must therefore be of the form



0 0 0 ½10 ½00 ½10 ½20 t 0b ¼ < t b1 eiwt þ e< t b2 þ t b2 eiwt þ t b2 e2iwt :

ð11Þ

From (1) and (11) we calculate the terms in the bedload transport (6), jQj3 Q 4

b Q

¼e

4 b rðs  1Þgd Q

n ½00 < t b2 j cosðwt 0 þ vÞj3

½10

½20

þ t b2 j cosðwt 0 þ vÞj3 eiwt0 þ t b2 j cosðwt0 þ vÞj3 e2iwt0  þ O e2    16 5 ½00 2iv ½20 < t b2 þ e t b2 þ O e2 ; ¼e b 3 5prðs  1Þgd Q

[20] From simple experiments [Carter et al., 1973], sandbars of length lB can be initiated by partially standing waves whose wavelength is twice the sandbar length, l = 2lB. If the barred region extends for many wavelengths, bars can augment reflection through Bragg resonance and alter the first-order waves [Mei, 1985], and thus their interaction is in general mutual [Yu and Mei, 2000a]. We shall verify that while unimportant for high reflection, Bragg resonance plays a crucial role in sandbar evolution for low to moderate reflection. [21] We expand the free surface displacement and velocity potential in the inviscid core as

o

ð12Þ

 z 0 ¼ z 01 þ ez 02 þ O e2 ;

 f0 ¼ f01 þ ef02 þ O e2 :

Assuming that the depth contours and shoreline are straight and parallel and the surface waves are normally incident in the x-direction, we write the first-order free surface displacement as

0 ½10 ½10 z 01 ¼ < z 1 eiwt ; with z 1 ¼ A0 eiS þ B0 eiS ;

and jQj4 3 ¼ j cosðwt 0 þ vÞj4 þ OðeÞ ¼ þ OðeÞ: 4 b 8 Q

ð13Þ

[18] Thus, for predicting the mean sediment transport up to O(e), the first harmonics of the second-order velocity and t [1]0 field and bed shear stress, u[1]0 2 2 , are not needed. In addition to the zeroth harmonic (the steady drift), the second harmonic is important and is derived later. Furthermore, since q0B = O(e), only the leading order (O(e0)) terms in b are needed to the accuracy of O(e). In QB0 (and hence Q) summary, substituting (6), (12), and (13) into (7) and dropping terms of order O(e2) and higher yields ð1  N Þ

  0 @h0 @ @q0 0 @h ¼ t0 ;  D n 0 0 0 @t @x @x @x

D0n ¼

Q0B ; tan ym

Z S¼

  5 ½00 ½20 < t b2 þ e2iv t b2 : b 3 15prðs  1Þgd Q 128eQ0B

k 0 ðex0 Þdx0 :

½00  0 0 x1 ;
t

f01 ¼ f1

ð19Þ



0 ½10 þ < f1 eiwt ;

ð20Þ

where
t 0 signifies the parametric effects of the slow sandbar growth. Note that the velocity associated with the long0 wave potential is e@f[0]0 1 /@x1 and is of the second order. The spatial amplitude of the first harmonic potential,

ð15Þ ½10

ð16Þ

x0

[22] The slow coordinate x10 = ex0 describes the gentle mean depth H0(x10 ) of the beach. The leading-order velocity potential consists of a mean current and a first-harmonic short wave,

f1 ¼  q0t ¼

ð18Þ

where A0 and B0 represent the slowly varying amplitudes of the incident (right-going or shoreward) and reflected (leftgoing or seaward) waves, and S denotes the wave phase

ð14Þ

which is a forced-diffusion equation for the seabed elevation h0, where

ð17Þ

ig cosh k 0 ðz0 þ H 0 Þ  0 iS A e þ B0 eiS w cosh k 0 H 0

ð21Þ

is the solution of a homogeneous boundary-value problem and w0 and k0(ex0) satisfy the dispersion relation

The effective diffusivity Dv0 , given by (15), results from the gravitational effect on bed load transport on a slope. The forcing term @qt0/@x0, where qt0 is given in equation (16), is 4 of 23

w2 ¼ gk 0 tanh k 0 H 0 :

ð22Þ

HANCOCK ET AL.: SANDBAR FORMATION UNDER SURFACE WAVES

C07022

Note that under our assumptions the sandbar height is at most O(e) of the total mean depth H0, and hence the sandbar height does not appear explicitly in (21). h0(x0,
t 0 ) denote the sandbar height above the mean [23] Let e 0 depth below seabed H , so that the total (ripple-averaged) 0 the still water level is h0 = H0  e h . We assume that the 0 sandbar height is small relative to the mean depth, e h /H0 = O(e). Anticipating the sandbars to be spatially periodic with half the local wavelength of the surface waves, we express 0 e h (x0,
t 0 ) as a Fourier series of spatial harmonics, 0 e h ¼<

1

X ½n0 e h e2inS :

ð23Þ

n¼1

[24] At the second order, the velocity potential and free surface elevation contain zeroth, first, and second harmonics in time, ½00  0 0 x1 ;
t

z 02 ¼ z 2 f02

¼

½00  f2 x01 ;
t 0



0 0 ½10 ½20 þ < z 2 eiwt þ z 2 e2iwt ;

0 0 ½10 ½20 þ < f2 eiwt þ f2 e2iwt :

Cg0

Cg0

½10 @A0 A0 @Cg0 þ ¼ iW0 e h B0 ; @x01 2 @x01 0

0

@B B  @x01 2

@Cg0 @x01

½10

*A0 ; ¼ iW0 e h

The last term in (31) is independent of all fast coordinates; it has negligible effect on the flow velocity, though it corresponds to an oscillating pressure affecting microseisms [Longuet-Higgins, 1950]. 0 t 0 ) corre[26] The zeroth-harmonic potential f[0]0 2 (x1,
sponds to an inviscid current above the bottom boundary layer. Once the sinusoidal incident waves have been maintained for a long time in a wave flume of finite length, a depth-averaged return current must exist to cancel the Stokes drift. While this current may be unimportant in an open sea, in a laboratory flume (of finite length) it can induce a shear stress in the bottom boundary layer, which can in turn affect the sediment transport. To find the return current, let us first calculate the period-average of the horizontal volume flux across a vertical cross-section. First, the instantaneous volume flux rate across a vertical crosssection is

ð24Þ Mþ0 ð25Þ

Mei [1985] has shown that only the first harmonic of the and f[1]0 second order flow, z [1]0 2 2 , are affected by the bar profile, and furthermore that the inhomogeneous boundary is solvable only if A0 and B0 value problem governing f[1]0 2 are coupled by ð26Þ

Cg0 ¼

0

0

w 2k H 1þ sinh 2k 0 H 0 2k 0

¼

Z

H 0

H 0 þ~ h0

Z

 0 f1x0 þ ef02x0 dz0  þ z 01 f01x0 jz0 ¼0 þ O e2 ;

½20

z2

½20

3iw cosh 2k 0 ðz0 þ H 0 Þ  02 2iS A e þ B02 e2iS 8k Ab sinh4 k 0 H 0 * * iwA0 B0 ð1  2 cosh 2k 0 H 0 Þ :  4k Ab sinh2 k 0 H 0 * *

ð32Þ

 ½00 ef02x0 ¼ e2 f2x0 ¼ O e2 ; 1

we obtain  ½00 Mþ0 ¼ eH 0 f1x0 þ z 01 f01x0 jz0 ¼0 þ O e2 :

ð33Þ

1

ð28Þ

0 Use is made of the fact that the sandbar elevation e h varies little during a few wave periods and is O(e) of the mean depth H0. The quadratic product above is the Stokes drift, which can be calculated from (18) and (21):

ð29Þ

has the dimensions of frequency. [25] The second harmonic is unaffected by the sandbar profile and is formally the same as that for a barless bed of mean depth H0,  k 0 1 þ 2 cosh2 k 0 H 0 cosh k 0 H 0 ¼ 4k Ab sinh3 k 0 H 0 * *   A02 e2iS þ B02 e2iS ;

0

where coordinates in subscripts indicate partial derivatives. Taking the time average over a wave period 2p/w, denoted by overlines, and noting that



wk 0 e h 2k Ab sinh 2k 0 H 0 * *

f0x0 dz0

H 0

z0 ¼H 0

½10

¼

0

)

þ

1

is the group velocity, e h is the first 0spatial harmonic ½1 amplitude of the sandbar profile, and (e h )* is its complex conjugate. The coupling coefficient ½10

z0

þ H 0

½00

½10

h W0e

Z

0

þ

  0 0  e ¼ h f1x0 

where 

(Z

f01x0 ¼ ef1x0 ; ð27Þ

C07022

z 01 f01x0 jz0 ¼0 ¼

ð34Þ

In order for the net volume flux Mþ0 to vanish, the return current must be, from (33) and (34), ½00

ef1x0 ¼  1

ð30Þ

gk 0 0 2 jA j  jB0 j2 : 2w

gk 0 0 2 jA j  jB0 j2 : 0 2wH

ð35Þ

In the special case of perfect reflection, jA0j = jB0j, there is no return current. [27] For later computations, all variables are made dimensionless according to:

f2 ¼ 

ð x; h; H Þ ¼ k ðx0 ; h0 ; H 0 Þ; *

ðt;
t Þ ¼ wðt 0 ;
t 0 Þ;

ð36Þ

ð31Þ ðA0 ; B0 Þ ¼ A ð A; BÞ; * 5 of 23

Cg0 ¼

w Cg ; k *

W¼

Ab * W0 ; w

ð37Þ

HANCOCK ET AL.: SANDBAR FORMATION UNDER SURFACE WAVES

C07022

where all normalizing scales denoted by asterisks are constants. The slow time variable
t describing the sandbar growth is found later. We also choose to normalize the sandbar height by the orbital amplitude Ab* just above the bed, i.e., 0 e hAb ; h ¼e *

A ð1 þ jRL ð0ÞjÞ where Ab ¼ * ; * sinh k H * *

ð38Þ

and RL(
t = 0) is the initial value of the reflection coefficient at the shoreward end (x = L) of the bar patch. In this work, RL is assumed to be a given constant or function of the long timescale
t. [28] For convenience, we rewrite the evolution equations (26) and (27) in terms of the dimensionless incident amplitude A(x1,
t) and the local reflection coefficient R(x1,
t) = B/A = jRjeiqR,   ½1 @A A 1 @Cg ¼ þ iWe h R ; @x1 Cg 2 @x1   @R iW e½1* e½1 2 þh R ; ¼ h @x1 Cg ½ 1

½10

where e h =e h /Ab* is the dimensionless amplitude of the ½1 first harmonic sandbar height and e h * is its complex conjugate. [29] We now cite a result important later for understanding our experiments. For a flat mean seabed (dH/dx1 = 0), (26) and (27) reduce to, in dimensionless variables, Cg

½1 @A ¼ iWe h B; @x1

Cg

½1 @B ¼ iWe h *A: @x1

ð41Þ

[30] Multiplying the first equation in (41) by A* and the second by B* and adding the resulting equations to their complex conjugates, we obtain Cg

Cg

½ 1 ½ 1 h eiqR g = 0, there is no reflection, then e h = eiqR, =fe coupling, and Bragg resonance is not effective.

4. Turbulent Boundary Layer [31] In what follows, we use the flow in the inviscid core to calculate the flow in the seabed boundary layer, and then calculate the bed shear stress components that appear in the forcing terms in the forced-diffusion equation (14) for the seabed elevation. [32] Inside the oscillatory turbulent boundary layer, we adopt the simple eddy viscosity model relating the shear stress and strain rate,

n ½1 o @jAj2 ¼ 2W= e h eiqR jAjjBj; @x1

ð42Þ

n ½1 o @jBj2 ¼ 2W= e h eiqR jAjjBj: @x1

ð43Þ

From these equations, Yu and Mei [2000a] observed that if ½1 =fe h eiqR g > 0, jAj and jBj increase shoreward as energy is transferred from the reflected wave to the incident wave, ½1 i.e., from B to A. Conversely, if =fe h eiqR g < 0, jAj and jBj decrease shoreward as energy is transferred from A to B. Note that the energy jBj2/2 of the reflected wave propagates seaward. Furthermore, Yu and Mei showed that the relative position of the bar crest and the wave node is related to the ½ 1 sign of =fe h eiqR g. In particular, if the bar crest is shoreward of the antinode and seaward of the node, ½1 =fe h eiqR g < 0, and vice versa. If bar crests appear directly under the wave nodes, as will be shown for perfect

½10

½10

t 1 ¼ rn 0e

@u1 ; @h0

½20

½20

t 2 ¼ rn 0e

@u2 ; @h0

ð44Þ

where

ð39Þ

ð40Þ

C07022

h0 ðx0 ;
t 0 Þ ¼ z0 þ h0 ðx0 ;
t 0 Þ

ð45Þ

denotes the height above the mean surface of the rippleaveraged seabed. The same eddy viscosity n e0 is used for the oscillatory parts of the first and second order shear stresses, since momentum diffusion is due to eddies created by the first-order orbital velocity. The eddy viscosity increases linearly with height h0, n 0e ¼ ku0f h0 ;

ð46Þ

where the local friction velocity uf0(x0,
t 0 ) varies slowly in time as the sandbars grow, and k = 0.41 is the Ka´rma´n constant. [33] The eddy viscosity model, (44) and (46), is one of many that have been proposed for simple oscillatory flows [Sleath, 1990]. Although models of transient eddy viscosity exist [Trowbridge and Madsen, 1984a, 1984b; Davies and Villaret, 1999], the simple model (44) and (46) is known to give reliable prediction of the flow velocity and makes analytical calculation possible, and hence is adopted here. [34] We shall assume that u10 vanishes at the roughness height hb0, which is related empirically to the local ripple height hr0 by Grant and Madsen [1982] h0b ¼ 4h0r =30:

ð47Þ

Because ripples are generated by waves, their local height and hence h0b depend on the wave envelope, which varies in x0 and
t 0 . The bed shear stress components t [n]0 bm (n = 0, 1, 2 and m = 1, 2) used above are simply the corresponding shear 0 0 stress components t [n]0 m evaluated at h = h b. 0 0
0 [35] The friction velocity u f (x , t ) is defined by the magnitude of the leading-order bed stress,

6 of 23

u02 f ¼

   ½10  t   b1  r

0    @u½1   ¼ ku0f h0 10  : @h h0 b

ð48Þ

HANCOCK ET AL.: SANDBAR FORMATION UNDER SURFACE WAVES

C07022

[36] The characteristic thickness d* of the oscillatory boundary layer is chosen to be sffiffiffiffiffiffiffiffiffiffi 2n e *: d ¼ * w

ð49Þ

is the normalized inviscid velocity just above the boundary layer, i.e.,

½1 ½0 ½1 ½2 U1 ; eU2 ; eU2 ; eU2 1 ¼ U *

The velocity scale is taken as the characteristic orbital velocity above the bed, U* ¼ Ab* w:

ð50Þ

*

w

ðZ; Zb Þ ¼

 0 0 2 h ; hb ; uf ð x;
t Þd *

ðu; wÞ ¼

u0f

0

ðu ; w Þ ; U *

uf ¼

uf

;

*

 0 0 t ; tb ðt; t b Þ ¼ ; rn e U =d * * *

! ð59Þ z0 ¼H 0

ð52Þ

ð53Þ

ð54Þ

where u, w, t, t b are functions of x, Z, t,
t and uf = uf (x,
t ). Note that the normalized coordinates Z, Zb vary slowly in time with the growth of the sandbars and vary in space with the local wave amplitude. In the boundary layer coordinates (x, Z), the component of velocity orthogonal to the rippleaveraged seabed is w0N ¼ w0 þ

½20

A similar expansion holds for the dimensionless bed shear stress t b. Let us now examine the boundary layer flow by a perturbation analysis. 4.1. First-Order Flow 4.1.1. Velocity and Bed Shear [38] Expanding (u, wN) in ascending powers of e, one derives from (57), after using (8), (9), (55), and (58), the horizontal momentum equation at the leading order, n o @u1 1 @t 1 ½1 ¼ þ < iU1 eit ;
@t uf ð x; t Þ @Z

and define the normalized velocities and shear stresses by 0

½10

ð60Þ

[37] For further analysis we introduce the following normalized coordinates in the oscillatory boundary layer: x ¼ k x0 ; *

½10

ð51Þ

It follows from (49) and (51) that 2kuf

½00

 t ¼ t 1nþ et 2 þoO e2 n o  ½1 ½0 ½1 ½2 ¼ < t 1 eit þ e< t 2 þ t 2 eit þ t 2 e2it þ O e2 :

Let the scale of the friction velocity be uf *, which is deduced later. The scale of the eddy viscosity is inferred from (46),

d ¼ *

½10

0 @2f @f1 @f1 @f2 @f ; ; þe h 0 1 0 ; 20 0 0 0 @x @x1 @x @x @z @x

and the dimensionless shear stress is

u f0

n e ¼ kuf d* : * *

C07022

ð61Þ

where t 1 ¼ uf ð x;
tÞZ

@u1 ; @Z

ð62Þ

in accordance with (44), (53), and (54). The boundary conditions are u1 ¼ 0

at Z ¼ Zb ;

ð63Þ

0

@h 0 u  k d U wN ; * * * @x0

ð55Þ

where wN is dimensionless. The normalized continuity equation becomes @u Z @uf @u 2 @wN ¼ 0:  þ @x uf @x @Z uf @Z

½1

U1 ¼

ð57Þ

A A  iS * e  ReiS : Ab sinh kH *

ð65Þ

[39] The solution for this boundary-value problem has been found by Kajiura [1968],     n o K0 ðZ Þ it ½1 ½1 e ; u1 ¼ < u1 eit ¼ < U1 1  K0 ðZb Þ

ð66Þ

where pffiffiffi

pffiffiffi

pffiffiffi

K0 ðZ Þ ¼ ker0 2 Z  i kei0 2 Z ¼ K0 2eip=4 Z ; ð67Þ

where n o n o ½1 ½0 ½1 ½2 UI ¼ < U1 eit þ e< U2 þ U2 eit þ U2 e2it

ð64Þ

where, from (21) and (58),

ð56Þ

We now assume that the boundary layer thickness is no more than O(e) of the amplitude of the bars or waves, d */ Ab* = O(e). Ignoring terms of order O(k*d *) = O(e2), the normalized horizontal momentum equation in the boundary layer reads   @u @u Z @uf @u 2wN @u þe þ eu  uf @Z @t @x uf @x @Z 1 @ @UI @UI ¼ þ eUI ; þ @t @x uf @Z

n o ½1 u1 ! < U1 eit as Z ! 1;

ð58Þ

kern(z) and kein(z) are the Kelvin functions of order n, and Kn(z) is the modified Bessel function of the second kind of

7 of 23

HANCOCK ET AL.: SANDBAR FORMATION UNDER SURFACE WAVES

C07022

C07022

and   K0 ðZ 0 Þ 1 dZ 0 K0 ðZb Þ Zb pffiffiffi pffiffiffiffiffi 1  i Z K1 ðZ Þ  Zb K1 ðZb Þ ¼ Z  Zb þ pffiffiffi : K0 ðZb Þ 2 Z

Z

F1 ðZ; Zb Þ ¼

ð73Þ

The solution can be checked directly by substituting (66) and (72) into the continuity equation (56), noting that Zb = Zb(x,
t ) and using (70) and (73) for the derivatives of K0 and F1, respectively. 4.1.2. Friction Velocity, Friction Factor, and Shields Parameter [ 42] Following Kajiura [1968], Grant and Madsen [1979], and others, the local friction velocity is calculated via the friction factor f in terms of the flow above the boundary layer: Figure 3. Friction factor f(y) (solid) from equation (75) compared with the fitted formulas equation (77) (dash) and that of Madsen [1994] (dash-dot). In Madsen’s formula, ubr/ (kNwr) is replaced by y/(30k).

 2 2  1  ½10  1 2  ½1  U u02 ¼ ¼ : f U fU f 2  1  2 * 1 

Combining (48), (53), (54), (68), and (74), one obtains a transcendental relation for the friction factor f in terms of a single dimensionless variable y,  h  pffiffiffiffiffiffiffiffii1     K y f =2 0 pffiffiffiffiffiffiffiffi 3=2    h  ky ¼ y f =2  pffiffiffiffiffiffiffiffii1  ; K1 y f =2   

order n [Abramowitz and Stegun, 1972]. The corresponding shear stress amplitude on the bed surface is evaluated from (62) and (66), ½1

t b1 ¼ uf Z

pffiffiffiffiffi ½1 ½1  @u1  1 þ i uf Zb K1 ðZb ÞU1 p ffiffi ffi ¼  ;  @Z Zb K0 ðZb Þ 2

ð68Þ

where

   ½1  U U1  * 1 : yð x;
t Þ ¼ pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ¼ ¼ Zb u0f 2uf h0b =d Zb f ð yð x;
t ÞÞ=2 * *

where pffiffiffi

pffiffiffi

pffiffiffi

K1 ðZ Þ ¼ ker1 2 Z  i kei1 2 Z ¼ iK1 2eip=4 Z : ð69Þ

To obtain (68), the derivative of K0 (Z) was evaluated from Abramowitz and Stegun [1972] (equation (9.6.27)), pffiffiffi

dK0 eip=4 K1 ðZ Þ ¼  pffiffiffi K1 2eip=4 Z ¼ ð1 þ iÞ pffiffiffiffiffiffi : dZ Z 2Z

ð70Þ

[40] The phase lag v between the shear stress and the flow, introduced earlier as Q / cos (t + v) and present in the forced-diffusion equation (14) for the seabed elevation, can now be deduced from (1), (54), and (62):  ½1  ½1    @u1  @u1   iv  e : uf Z ¼ uf Z @Z Zb  @Z Zb 

ð71Þ

[41] Next, by integrating the dimensionless continuity equation (56) and using (66) and the no-slip condition u1, w1 = 0 on Z = Zb(x,
t ), we obtain the normal velocity wN at ½1 the leading order, w1 ¼
w1 ¼ 

ð72Þ

   ½10  U   1 

ð75Þ

ð76Þ

The middle equality in (76) is a simple rearrangement of (74). The last equality follows from (53) and (54). Equation (75) for the friction factor was found in similar forms by Kajiura [1968] and Grant and Madsen [1979], and is plotted in Figure 3. For the convenience of numerical computations, this curve is fitted by f ð yÞ ¼ 100:0433ð log10 yÞ

2

0:637 log10 y0:131

ð1 þ errð yÞÞ:

ð77Þ

The error is jerr(y)j < 0.02 for 2 < y < 105. [43] We now deduce the scale of the friction velocity uf * from that of y. From (38) and (65), the dimensionless orbital velocity amplitude jU[1] 1 j = O(1), and from (47), the roughness height hb0 = O(4hr*/30). Thus from (52) and (76), we choose the characteristic value of y to be U kU kA * * ¼ b* : ¼ y ¼ * 2uf hb =d whb hb * * * * *

ð78Þ

From (74) and (78), we choose 1

u2f ¼ f y U 2 ; * * * 2

½1 F1 ðZ; Zb Þuf @U1

2 @x     ½1 U1 @uf K0 ðZ Þ þ  F1 ðZ; Zb Þ Z 1 2 @x K0 ðZb Þ rffiffiffiffiffi  ½1 U1 uf K1 ðZb Þ Z @Zb þ K1 ðZ Þ  K1 ðZb Þ @x Zb 2K20 ðZb Þ

ð74Þ

ð79Þ

where f (y*) can be found from (75) or (77). [44] The local dimensionless friction velocity uf (x,
t), used to calculate the local Shields parameter and the bed load transport rate, is found as follows. Once uf * is found from (79), the local friction factor f (y(x,
t)) may be found for any x from (75) or (77) and the last equality in (76). The

8 of 23

HANCOCK ET AL.: SANDBAR FORMATION UNDER SURFACE WAVES

C07022

local dimensionless friction velocity uf (x,
t ) = uf0(x0,
t 0 )/uf * can then be calculated from (74). [45] The timescale of bed load transport and the sandbar growth rate depend on the scale of the Shields parameter. Using the normalization so far introduced, the natural scale of Shields parameter is

u2f f y U2 * * * Q ¼ ¼ : * ðs  1Þgd 2ðs  1Þgd

   ½10  t   b1 

b ð x;
t Þ Q ¼ 2 ¼ u2f ð x;
t Þ: Q ruf * *

Substituting (66) and (72) gives 8    ½1 2 > >   > U ½1 1 < ½1 @U1 @c1  1  @Zb @c2 RHS½0 ð x; Z;
t Þ ¼ < U1 * þ @x @Z Zb @x @Z 2 > > > :

ð80Þ þ

From (11), (48), and (54), the leading-order time-amplitude b ð x;
t Þ of the Shields parameter Q b ð x; t;
t Þ is related to the Q dimensionless friction velocity uf via

þ

@u2 1 @t 2 @u1 u1 Z @uf @u1 2w1 @u1  ¼  u1 þ  @t @x uf @Z uf @Z u @x @Z ( ½1 f ) ½ 1  ½1 ½1 * U1 @U1 U1 @U1 2it ½1 it þ<  iU2 e þ e 2 @x 2 @x n o ½2  2< iU2 e2it ;



½2

U2

is the dimensionless amplitude of the second harmonic of the second-order potential flow at z = H. As pointed out at the end of section 2, the first harmonic u[1] 2 of the secondorder flow does not contribute to the first-order bed load [1] are not pursued here. transport; the details of u[1] 2 and U2 The forcing terms on the right hand side of (82) can be formally written as n o ½1 ½2 RHS ¼ RHS½0 þ < iU2 eit  2iU2 e2it n o þ < RHS½1 eit þ RHS½2 e2it ;

ð87Þ

ð89Þ

The functions cn(Z, Zb) depend on x through both Z and Zb. As shown shortly, only their values along Z = Zb are needed here, and found in Appendix B by integrating equations (87) – (89). [47] Omitting the details of the first harmonic, which is not needed for calculating the bed load transport, the depthdependent part of the second harmonic forcing is ½1

RHS½2 ð x; Z;
t Þ ¼

½1

½1

½1

U1 @U1 u @u1  1 2 @x 2 @x! ½1 ½1 ½1 Zu1 @uf w1 @u1 þ  : 2uf @x uf @Z

ð90Þ

Substituting (66) and (72) gives ð84Þ

where the linear part corresponding to @U2/@t from the pressure gradient of the inviscid flow above the boundary layer is separated from the rest. The zeroth harmonic forcing corresponds to the Reynolds stress and can be evaluated to give ½1 @U1

ð86Þ

 ! iZb dK0  1 dK*0 Z  Zb þ K0 ðZb Þ dZ Zb K*0 ðZb Þ dZ  !  Zb dK0  F1 ðZ; Zb Þ dK*0 Zb dK0  1þ  K0 ðZb Þ dZ Zb K*0 ðZb Þ dZ K0 ðZb Þ dZ Zb !   Z  Zb dK*0 K0 ðZ Þ  K0 ðZ Þ 2 ; ð88Þ þ  K*0 ðZb Þ dZ K0 ðZb Þ K0 ðZb Þ 

RHS½2 ð x; Z;
t Þ ¼

(

@uf @c3 ; @x @Z > > > ;

 ! @c3 iZb dK0  1 dK*0 : ¼  Z  Zb þ @Z K0 ðZb Þ dZ Zb K*0 ðZb Þ dZ

ð82Þ



2 ½1 ½1 3k A =Ab @U1 * A2 e2iS  R2 e2iS ¼  3iU1 * ¼ 4 sinh4 kH 4 sinh2 kH @x ð83Þ

uf

9 > > > =

   @c1 K0 ðZ Þ 2 F1 ðZ; Zb Þ dK*0  ; ¼ 1  1  @Z K0 ðZb Þ  K*0 ðZb Þ dZ @c2 ¼ @Z

where, from (31) and (58),

   ½1 2 U   1 

where

ð81Þ

4.2. Second-Order Flow 4.2.1. Boundary Layer Equation and Forcing Terms [46] At O(e) the horizontal momentum equation reads

C07022

½1 c4 U1

2



½1 2   U1 c5 @Zb c6 @uf þ þ ; @x 2 Zb @x uf @x

½1 @U1

ð91Þ

where

½1 @u1

1 ½1* ½1* RHS½0 ð x; Z;
t Þ ¼ < U1  u1 @x @x 2

) ½1* ½1 ½1 ½1* u1 Z @uf @u1 2w1 @u1 þ  : ð85Þ @x @Z uf @Z uf

9 of 23

  K0 ðZ Þ 2 F1 ðZ; Zb Þ dK0 c4 ðZ; Zb Þ ¼ 1  1   ; K0 ðZb Þ K0 ðZb Þ dZ

ð92Þ

HANCOCK ET AL.: SANDBAR FORMATION UNDER SURFACE WAVES

C07022

c5 ðZ; Zb Þ ¼

   Zb dK0  K0 ðZ Þ 2 K0 ðZ Þ  K0 ðZb Þ dZ Zb K0 ðZb Þ K0 ðZb Þ ! F1 ðZ; Zb Þ  ðZ  Zb Þ dK0 þ ; dZ K0 ðZb Þ

where L is a coefficient of order unity and hb* = 4hr*/30 from (47). Solving for t [0]0 2C , nondimensionalizing, and using (51) gives 

ð93Þ

 ½0  t 2C 

¼

Zb

c6 ðZ; Zb Þ ¼ 

F1 ðZ; Zb Þ dK0 ; K0 ðZb Þ dZ

ð94Þ

and F1 is given in (73). We display explicitly the dependence on Z and Zb to keep track of the dependence on x and
t. [48] The solution for u2 is expected to be of the form n o ½0 ½1 ½2 u2 ¼ < u2 þ u2 eit þ u2 e2it :

ð95Þ

ð96Þ

in accordance with (82). Integrating (96) once and imposing the shear-free condition at the outer edge of the boundary layer yields the shear stress at Zb  ½0  t 2W 

Z

1

¼ uf

RHS½0 ð x; Z;
t ÞdZ:

ð97Þ

Zb

Zb

Substituting (86) gives  ½0  t 2W 

Zb

(

¼



½0 U2



RHS½2 ð x; Z;
t Þ @2 1 @ 2i ½2 ½2 U ¼ þ  u þ ; 2 2 @Z 2 Z @Z Z Z


ð100Þ

½2

u2 ¼ 0

on Z ¼ Zb ;

ð101Þ

[50] The dimensional bed stress t [0]0 2C due to the second-order return flow U[0]0 2 is found by assuming a logarithmic current velocity very close to the bed that takes the value U[0]0 2 at an elevation O(d *) from the bed [Grant and Madsen, 1979]:

¼L

h0b rku0f

d ln * ; hb *

½2

as

Z ! 1:

ð107Þ

ð102Þ

½2  @u2  @Z Zb

½2 pffiffiffiffiffi ð1 þ iÞuf U2 Zb K1 ð2Zb Þ K0 ð2Zb Þ ½1 Z 1 uf ½1 @U1 K0 ð2Z Þ c4 ðZ; Zb Þ þ U1 dZ 2 @x Zb K0 ð2Zb Þ Z uf ½1 2 1 @Zb 1 K0 ð2Z Þ c5 ðZ; Zb Þ þ U1 dZ 2 Zb @x Zb K0 ð2Zb Þ Z

1 2 @uf 1 K0 ð2Z Þ ½1 c6 ðZ; Zb Þ dZ; þ U1 @x Zb 2 K0 ð2Zb Þ

0

½00

½2

u2 ! U2

¼

Substituting (99) and (100) into (98) gives

U2

ð106Þ

subject to the boundary conditions

½2

ð99Þ

 ½00  t 2C 

ð105Þ



ð98Þ

1 þ i pffiffiffiffiffi K* ðZb Þ ; c1 ðZb ; Zb Þ ¼ pffiffiffi Zb 1 2 K0* ðZb Þ

( ) ½1 uf 1 þ i pffiffiffiffiffi K1* ðZb Þ ½1* @U1 U1 ¼ < pffiffiffi Zb : 2 @x 2 K* ðZb Þ

2 A =Ab  2 * * jAj  jBj2 : ¼ 2H tanh kH

4.2.3. Second Harmonic Bed Stress [52] From (44), (54), (82), and (84), the momentum equation for the second harmonic at the second order can be written as

In Appendix B, equations (87) – (89) are integrated to find



ð103Þ

where the dimensionless return flow velocity can be derived from (35),

t b2 ¼ uf Z

½0  t 2W  Zb

ern e U =d * * *

uf U2

: L ln d =hb * *

The straightforward, but lengthy procedure to solve (106) and (107) is outlined in Appendix C; the solution u[2] 2 and are given by equations (C5) and (C6), shear stress t [2] 2 respectively. Substituting (54) and (C6) into (44) gives the corresponding bed stress

½1 @U1

uf ½1 < U1 * c1 ðZb ; Zb Þ 2 @x    )  ½1 2 c2 ðZb ; Zb Þ @Zb c3 ðZb ; Zb Þ @uf : þ U1  þ @x @x Zb uf

½0

¼

( ) ½0 ½1 uf uf U2 1 þ i pffiffiffiffiffi K* ðZb Þ ½1* @U1

; U1 < pffiffiffi Zb 1 þ 2 @x 2 K0* ðZb Þ L ln d =hb * * ð104Þ

½0

½0

½00  t 2C  h0b

We shall later choose L to fit the experiments. [51] In summary, the total mean bed stress follows from (101) and (103), t b2 ¼

[2] Only the shearing rates of u[0] 2 and u2 are needed and calculated below. 4.2.2. Mean Bed Shear Stress [49] There are two sources for the mean shear on the sandbar surface. One is the Reynolds stress inside the boundary layer and the other is due to the inviscid return current U[0] 2 above the boundary layer. Let the mean shear [0] [0] [0] stress t [0] 2 be decomposed as t 2 = t 2W + t 2C. The first part must satisfy

@t 2W ¼ uf RHS½0 ð x; Z;
t Þ @Z

C07022

ð108Þ

where cn(Z, Zb), n = 4, 5, 6, are listed in equations (92) – (94). The integrals in (108) must be computed numerically; the details are found in Appendix B.

5. Equation for Sandbar Evolution [53] Under our nondimensionalization in section 3, including equation (38), the dimensionless ripple-averaged depth h can be expressed as the barless mean H minus the

10 of 23

HANCOCK ET AL.: SANDBAR FORMATION UNDER SURFACE WAVES

C07022

bar elevation ee h,

[56] Substituting (71), (81), (83), (104), and (108) into (116), we finally obtain hð x;
t Þ ¼ H ðx1 Þ  ee hð x;
t Þ;

ð109Þ qt ¼

where H(x1) depends on the slow coordinate x1 =ex. Note that, from (36) and (109), ! @h @h0 dH @e h ¼e  ¼ : @x @x0 dx1 @x

Q0B QB ¼ pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi

3=2 : 3 ðs  1Þgd 3 rd Q *

ð110Þ

ð111Þ

Z iK1* ðZb ÞK0 ðZb Þ 1 c4 ðZ; Zb ÞK0 ð2Z ÞdZ K0 ð2Zb Þ 2K1 ðZb ÞK0* ðZb Þ Zb pffiffiffiffiffiffiffi * 5ð1 þ iÞ 2Zb K1 ðZb Þ þ ; 12K* ðZ Þ

M1 ðZb Þ ¼ 

ð114Þ

which is a forced diffusion equation for the sandbar elevation e h(x,
t), where 3QB ; 8ð1  N Þ tan ym

  16kU QB 5 ½0 2iv ½2 * < þ e t t b2 ; 3 b2 5pð1  N Þuf u2f *

QB ¼

b  QC rd Q rd Q *

!3=2



b  QC : Hv rd Q

b

1

b

iK1* ðZb ÞK0 ðZb Þ 2K ðZ ÞK* ðZ Þ b

0

0

Z

1 Zb

b

ð120Þ

b

 c5 ðZ; Zb ÞK0 ð2Z ÞdZ ; K0 ð2Zb Þ ð121Þ

[55] Substituting the dimensionless variables in (110), (111), and (113) into the dimensional seabed evolution equation (14) gives

Dn ¼

M3 ðZb Þ ¼ < 

1

ð113Þ

!   @e h @ @e h @ dH ;  Dn ¼ qt þ Dn @
t @x @x @x dx1

pffiffiffiffiffi ð1  iÞ Zb K1 ð2Zb ÞK1* ðZb ÞK0 ðZb Þ ; K ð2Z ÞK ðZ ÞK* ðZ Þ (

ð119Þ

b

0

0

ð112Þ

ð118Þ

where

M2 ðZb Þ ¼ 

The second equality follows from (52) and (80). In the parameter regime of interest, the depth-to-wavelength ratio is k*H* = O(1) and the boundary layer thickness-towavelength ratio is k*d *  O(e2). Hence a = O((k*d*)3)  O(e6) and @h/@t  O(e6), which implies
t = O(e6t) and defines the slow sandbar growth time
t ¼ t=a:

5U2

3L ln d =hb * *

) ½1 ½1* ½2 @U1 U1 U2 þ< þ M2 ðZb Þ ½1 @x U1 ! ½1 ½1 jU1 j2 @Zb jU1 j2 @uf þ M3 ðZb Þ þ M4 ðZb Þ ; Zb @x uf @x

Substituting (110) and (111) along with the dimensionless variables defined in (36) into the dimensional seabed evolution equation (14) results in @h/@t  O(1/a), where Ab w * a¼

3=2 pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 8e rd Q k ðs  1Þgd 3 * * ðs  1Þk3 1 ¼ 3=2

3 : rd tanh k H k d * * * *

½0

16kU QB * 5pð1  N Þuf uf * (

½1* M1 ðZb ÞU1

[54] The time rate of sandbar growth is found in terms of the maximum rate of bed load transport, from (4),

qt ¼

C07022

ð115Þ

( M4 ðZb Þ ¼ < 

iK1* ðZb ÞK0 ðZb Þ 2K ðZ ÞK* ðZ Þ 1

b

0

b

Z

1 Zb

 c6 ðZ; Zb ÞK0 ð2ZÞdZ : K0 ð2Zb Þ ð122Þ

[57] The functions Mn(Zb) have been plotted by Hancock [2005]. Numerically the functions M1, M2 are much larger than M3, M4. From (53), the derivative of Zb can be written in terms of the derivatives of the friction velocity uf and the dimensionless roughness height hb0/d *. [58] Yu and Mei [2000a] derived a forced-diffusion equation similar to (114) for the limiting case of constant eddy viscosity and constant H. The improvements made here include the depth-linear eddy viscosity in the boundary layer, the skin-friction Shields parameter in the bed load transport, and our inclusion of the shear stress induced by the return current. Furthermore, by using the depth-modification of King [1991] to the bed load transport rate, we avoid an additional fitting parameter (unrelated to L) present in Yu and Mei’s diffusivity Dv that was directly proportional to their sandbar heights for all degrees of reflection. The consequence of these changes will be assessed by comparison with experiments.

ð116Þ

6. Local Properties of Diffusivity and Forcing

ð117Þ

[59] To elucidate the physical picture, let us consider the local properties of the sandbar evolution equation for one sandbar wavelength at some stage of the slow evolution. Without loss of generality, we take the wave parameters to be locally constant: H = k*H*, k = 1, A = 1, R = RL. Computed

11 of 23

C07022

HANCOCK ET AL.: SANDBAR FORMATION UNDER SURFACE WAVES

C07022

greatest under the wave node and smallest under the antinodes. Figure 4 shows that the effective diffusivity Dn, which limits the steepness of the sand bed, is greatest under the wave node and smallest under the wave antinode. On the other hand, the bed load forcing, @qt/@x, is greatest and positive near the node implying that deposition occurs to form sandbar crests. The bed load forcing is negative somewhere between the node and antinode; scouring occurs there to form sandbar troughs. [62] In the limit of perfect reflection under a standing wave, the forcing is zero over a significant portion of the sand surface where the horizontal motion of the flow is small and the Shields parameter is below the threshold for sand transport. This implies the presence of a flat plateau with no sand motion. We point out that in the earlier theory of Yu and Mei [2000a], the full shear stress was used instead of the skin-friction shear stress (see Appendix A) to estimate the bed load transport. Consequently, the shear stress under the antinodes was over-estimated, leading to much narrower plateaus. Laboratory evidence supporting the present prediction will be discussed shortly. [63] Also, under perfect reflection, the forcing and diffusivity are symmetric about the wave node, and are positive under the wave node. Therefore, sandbar crests will form directly under the wave nodes which, as discussed at the end of section 3, renders Bragg scattering ineffective so that the incident and reflected waves are decoupled as if on a flat bottom.

7. Numerical Procedure for Multibar Evolution

Figure 4. Dimensionless leading-order wave amplitude jz 1[1]j, diffusivity Dv, and forcing @qt/@x across one bar length for various reflection coefficients, RL = 0.1 (solid), RL = 0.3 (dash), RL = 0.5 (dash-dot), RL = 1.0 (solid). Numbers adjacent to curves indicate the reflection coefficient. results are shown for inputs typical of field cases: T = 8 s, H* = 6 m, A*(1 +jRLj) = 50 cm, d = 0.5 mm, and hr* =5 cm. In dimensionless terms, results for typical laboratory cases are quite similar and are thus omitted. The value chosen for the fitting coefficient L is 0.3, close to the values later chosen to give best agreement with our laboratory data. [60] In Figure 4, the dimensionless leading-order wave amplitude   qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi  ½1  z  ¼ jAj 1 þ jRj2 þ 2jRj cosð2cÞ;  1 

c ¼ x  qR =2;

ð123Þ

diffusivity Dn , and bed load forcing @qt/@x are plotted across a single bar for various reflection coefficients RL. In this comparison, the maximum wave height A*(1 + jRLj) is fixed while the reflection coefficient RL is varied. Hence a higher reflection coefficient RL corresponds to a smaller incident wave amplitude A*. [61] From (4), (48), (65), and (81) the horizontal orbital velocity, bed shear stress, and bed load transport rate are

[64] In all computations, we consider the evolution of sandbars in a finite patch 0 < x < L parallel to a straight shore. Seaward of the patch (x < 0), we assume the seabed is either too deep or nonerodible. On the shoreward side of the patch, we assume the existence of either a vertical wall at x = L or a partially reflective shoreline somewhere beyond x > L. The dimensionless incident amplitude A(x1 = 0,
t ) and the reflection coefficient R(x1 = eL,
t ) = RL(
t ) are prescribed. The sandbar elevation e h is assumed to vanish at the ends of the sandbar patch x = 0, L for all
t > 0. For strong reflection, the ends x = 0, L are chosen to coincide with sub-critical regions where the forcing and e h are zero. For weak reflection, artificial damping is imposed within a few grid points of x = L to maintain numerical stability and a continuous profile. The boundary condition at x = L could be improved using information about the sediment transport at x = L, if available. Initially the bed is flat, e h(x,
t = 0) = 0=e h[1](x1,
t = 0). Near the transition region where the local b= Shields parameter is near the threshold of erosion, i.e., rdQ QC, the bed slope may be large. We have added artificial diffusion to achieve local smoothness and eliminate numerical instability. [65] For comparison with our laboratory experiments, where the mean depth H0 was constant, the measured dimensionless incident wave amplitude A(x1 = 0,
t ) = A0(x10 = 0,
t 0 )/ A* and reflection coefficient RL(
t) = R(x1 = eL,
t) provided boundary values for the dimensionless local amplitude A(x1,
t) of the incident wave and the local reflection coefficient R(x1,
t) governed by the steady-state Bragg resonance equations (39) and (40). For each time
t of the slow evolution of the sandbars, a fixed-step fourth-order Runge-Kutta

12 of 23

C07022

HANCOCK ET AL.: SANDBAR FORMATION UNDER SURFACE WAVES

scheme integrated (40) from x1 =eL to 0 to find R(x1,
t) and then integrated (39) from x1 = 0 to eL to find A(x1,
t). If the mean depth varies, H = H(x1), the local wave number k(x1), group velocity Cg(x1), and phase S(x) are found from the dimensionless ray theory equations (19), (22), and (28). [66] The wave parameters A and R are then interpolated over the short scale x and used to compute the local forcing qt and diffusivity Dn in the sandbar equation (114), where the spatial derivatives are discretized by finite differences. The sandbar elevation e h is then advanced in time
t by an adaptive-step 4th –5th order Runge-Kutta scheme subject to the boundary conditions e h = 0 at x1 = 0, eL. The sandbar ½1 e amplitude h is then calculated from the updated sandbar elevation e h using the Fast Fourier Transform (in Matlab). The wave parameters A, R are then updated from the Bragg resonance equations, followed by the sandbar elevation e h. This process is repeated until the desired time is reached. [67] Two types of simulations are discussed next. The first is for the transient bar evolution in three laboratory experiments in a wave flume of constant depth. The second is to use measured wave data as input to reconstruct sandbars observed in Cape Cod Bay. Because of the absence of data on the time-history of the Cape Cod beach, the second comparison may only suggest partial relevance of the theoretical model.

8. Laboratory Experiments and Simulation [68] Experiments were conducted in a horizontal wave flume of length 28 m, width 76 cm, and depth 90 cm with glass sidewalls and large holding tanks on either end [Landry, 2004]. Waves were generated by a programmable piston wavemaker at one end of the flume. Weak reflection was achieved by installing a 1/10 sloping absorber beach at the opposite end. The distance between the mean position of the wavemaker and the toe of the beach was 26.5 m. For strong (perfect) reflection, a vertical seawall was installed instead at 18.95 m from the mean position of the wavemaker. Beginning at 3.75 m from the wavemaker, a 10 cm thick layer of Ottawa silica sand was spread at the bottom for a length of 15.2 m. The sand layer extended to the end wall for the strong reflection test, and to 0.55 m before the 1/10 beach slope for the weak reflection tests. The sand density was rs = 2650 kg/m3 with specific gravity s = rs/r = 2.65. The mean grain diameter was d = 0.20 mm and the range of variation was not large: d16 = 0.14 mm and d84 = 0.28 mm, where dn stands for the diameter of the grain that is larger than n%, by mass, of all the sampled grains. Before each test, the sand bed was leveled. The still water depth over the sand layer was 60 cm. During each test, which lasted several days, time series of the free surface elevation were recorded, three stations at a time, by conductivity probes mounted on a movable cart. At each station, the free surface was sampled 2048 times at 24.1 Hz and 24.3 Hz for wave periods of 2.50 s and 2.63 s, respectively. Sixty stations, separated by 25 cm, covered the entire sand bed. One complete set of wave measurements required 49 min, which is considerably shorter than the timescale of sandbar evolution. [69] Sandbar profiles were recorded by a digital camera mounted on a side arm of the same cart and positioned at the level of the initial surface of the sand bed. Images of

C07022

successive portions of the seabed were taken and later stitched together to obtain the profile of the entire sand bed. Each set of seabed images required 10 min to collect. [70] In each test the wave intensity was chosen to be strong enough to move sediment without appreciable suspension, yet weak enough to avoid significant nonlinear effects which would lead to slow variations of the wave envelope with time and to seiching. The wave intensities were well above the threshold for creating a fully roughturbulent boundary layer, according to the following empirical criterion of Sleath [1990], 0 1 ! 2 lr hr 1:16 wA * G@ b *  108:2A n Ab l r * 0 1 !1:16 h r A b * @  0:042A  0:58; lr lr * *

ð124Þ

where hr*, lr* 5hr* are the typical ripple height and length, respectively. As shown in Table 3, Sleath’s criterion for turbulence is satisfied in our experiments. [71] As purely progressive waves proportional to exp (ikx  iwt) propagate, a bound second harmonic proportional to exp (2ikx  2iwt) is generated by nonlinearity (see equations (30) and (31)). Because of the presence and the motion of the wavemaker piston, a free second harmonic proportional to exp(ik2x  2iwt) with (2w)2 = gk2tanh k2H is also generated [Madsen, 1971]. The bound second harmonic of wave number 2k and the free second harmonic of wave number k2 periodically modulate the free surface. The free second harmonic can also cause seiching in the tank. By suitable programming of the piston motion (see Appendix D), the free second harmonic was essentially eliminated in the wave flume. The success of this elimination can be seen in the top panels of Figures 6 through 9 where the amplitude of the second harmonic is not modulated over long spatial scales. [72] Initially the sand bed was smooth and horizontal. Shortly after the wavemaker started, ripples began to form on the initially flat bed beneath the wave nodes and spread outward toward the antinodes. After approximately 30 to 60 min, the ripples reached their maximum heights and ceased to grow or expand. The observed characteristic ripple height confirms the empirical formula due to Nielsen [1981], hr

U * * ¼ 0:275  0:022 pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi : Ab ð s  1Þgd *

ð125Þ

[73] By removing the spatial average, the measured horizontal variations of the ripples over the sandbars are plotted in Figure 5. For test R10 with perfect reflection (RL = 1.0), Figure 5 shows that the highest ripples were found directly under the wave nodes where the ripple wavelength was almost constant, while flat plateaus representing the uneroded regions lay under the antinodes. In test R02a reflection was weak with jRLj = 0.240. The ripple amplitude was essentially uniform, hr* = 1.8 cm, which agrees with Nielsen’s empirical formula (125). In test R02b reflection was also weak with jRLj = 0.236. The ripple envelope undulated slightly. All measured ripple envelopes

13 of 23

HANCOCK ET AL.: SANDBAR FORMATION UNDER SURFACE WAVES

C07022

C07022

Table 2. Calculated Quantities in Physical Dimensions for the Laboratory Tests R10, R02a,b and for the Field Observations of Elgar et al. [2003]a Test

R10

R02a

R02b

Field

l*, m d*, cm a/w, days U*, cm/s uf *, cm/s wS, cm/s

5.67 2.23 1.55 25.9 6.82 2.46

6.01 2.26 1.92 25.1 6.59 2.46

5.67 2.52 1.06 30.1 7.74 2.46

124.3 31.4 162 62.2 10.1 4.48

a

The fall velocity wS was predicted by the formula of Ahrens [2000].

those for spherical grains. For the wave and sediment parameters in our laboratory, uf /wS = 2.7 to 3.1 (Table 2).

Figure 5. Ripple elevations in tests R10, R02a,b. Numbers adjacent to ripple profiles indicate the corresponding elapsed time in days. Solid curves are fitted envelopes of ripple height. For test R02a, hr0 = hr* = 1.8 cm. The origin is 4.78 m (R10), 2.80 m (R02a), and 3.50 m (R02b) from the mean wavemaker position. Each increment on the vertical axis corresponds to 4 cm. are fitted by smooth curves shown in Figure 5 for later numerical modeling of sandbars. [74] After two to three hours, bars began to emerge. The initial inputs over the horizontal bed are shown in Table 1. The calculated parameters in physical dimensions are given in Table 2. The salient dimensionless parameters are listed in Table 3. [75] In all laboratory tests presented here, the sediment was observed to be transported exclusively as bed load. The maximum excursion of the sediment grains from the bed occurred in vortices shedding off ripple crests; the grains were lifted at most a few millimeters before quickly falling back to the bed. The fall velocity wS was estimated using the formula of Ahrens [2000], based on the work of Hallermeier [1981], for natural sediment grains. The estimates are close to

Table 1. Physical Input Parameters for the Laboratory Testsa R10, R02a,b and for the Field Observations of Elgar et al. [2003]: Input Wave Period 2p/w; Water Depth H*; Incident Wave Amplitude A*; Beach Reflection Coefficient RL; Sand Diameter d; Characteristic Ripple Height hr* Predicted by Nielsen’s Formula (125) for R10 and R02a,b; Fitting Coefficient L for the Return Current; and Duration of Experiment Tf Test

R10

R02a

R02b

Field

2 p/w, s H*, cm A*, cm RL d, mm hr*, cm L Tf, days

2.50 60 3.68 1.0 0.20 1.8 n/a 4.04

2.63 60 5.68 0.24e1.05i 0.20 1.8 0.28 3.04

2.50 60 6.92 0.236e0.071i 0.20 1.9 0.31 3.98

23.8 280 23.1 pffiffiffiffiffiffi ffi 0:2e2.51i 0.33 5.8 0 n/a

a For test R02b, the phase of RL varied in time (see equation (128)); the value listed herein is that at
t = 0.

8.1. Test R10 Under a Standing Wave [76] The evolution of the sand bed under complete reflection jRLj =1 (test R10) over four days is shown in Figures 6 and 7. Wave records at discrete stations are marked by circles, while the computed values are shown by continuous curves. The normalized incident wave amplitude was constant, A(0,
t) = 1. The sandbar crests appeared directly under the wave nodes. Throughout the entire test the wave envelope remained essentially unchanged as seen from the initial and final wave envelopes of the first and second harmonics of the free surface shown in Figure 6. All these observations are consistent with our theory which predicts symmetric sandbar crests directly under wave nodes and therefore zero Bragg resonance. To fully test our numerical model, we still used the Bragg resonance equations (41) to calculate A and B for our simulations. As the sandbars grew in time, the maximum ½ 1 of the coupling coefficient =fe h eiqR g in the energy exchange equations (42) and (43) increased from 0 to 0.081 due to numerical error. Hence A and B were unaffected by the sandbars and could be calculated from the simple theory for a horizontal bottom. [77] With the parameters given in Tables 1 – 3 and the fitted ripple amplitude profile, predictions of the mean bed profiles were made and are shown in Figure 6 by the heavy and continuous curves. After averaging over the ripples, the agreement between the sandbar predictions and measurements is good at all times and is evident not only in the

Table 3. Calculated Dimensionless Parameters for the Laboratory Tests R10, R02a,b and for the Field Observations of Elgar et al. [2003]a Test

R10

R02a

R02b

Field

k*H*

0.66 0.11 0.025 0.0002 0.053 1.44 0.11 0.14 0.22 471

0.63 0.11 0.024 0.0002 0.053 1.34 0.12 0.14 0.21 470

0.66 0.13 0.028 0.0002 0.053 1.85 0.11 0.13 0.20 660

0.14 0.12 0.016 0.000017 0.041 2.02 0.14 0.053 0.05 40,808

e k*d* k*d QC Q* rd f(y*) Zb G

a The values of the l.h.s. of equation (124), G, indicate the boundary layers in all tests were fully rough turbulent.

14 of 23

C07022

HANCOCK ET AL.: SANDBAR FORMATION UNDER SURFACE WAVES

C07022

Figure 6. Measured and predicted wave amplitudes and sandbar elevations for test R10 at various times. The measured first and second harmonic wave amplitudes jz 1[1]0j (open circle) and jz 2[2]0j (solid circle) are plotted in the top and bottom rows. Between them are the recorded seabed elevations h0 (jagged 0 curves) and corresponding sandbar predictions eh (smooth curves without ripples). The origin is 4.78 m from the mean wavemaker position. Numbers adjacent to profiles indicate the elapsed time in days. Each increment on the vertical axis corresponds to 5 cm. positions and magnitudes of the maxima and minima, but also in the horizontal extents of the uneroded plateaus. Since there is no return current (U[0] 2 = 0), the value of the fitting parameter L is immaterial. Figure 7 compares the predicted and measured histories of sandbar heights. The sandbar growth rate is large initially and diminishes as the sandbars approach their steady state. [78] In Hancock [2005], the earlier theory of Yu and Mei [2000a] has also been used to simulate test R10. The measured equilibrium sandbar height was used in Yu and Mei’s model to set the fitting coefficient in their diffusivity. The width of the plateau is grossly underpredicted, being roughly one-tenth of that observed. 8.2. R02 Tests Under Weakly Reflected Waves [79] Two tests were performed for sandbar evolution under weakly reflected waves. For nearly spatially uniform envelopes, it was necessary to use larger wave amplitudes in order to cause significant granular movement. In tests R02a and R02b, the piston displacement was found to slowly increase and decrease with time for unknown reasons. From measured data, the normalized incident wave amplitudes A(0,
t) = A0(0,
t)/A* near the piston were fitted to the following empirical formulas: Test R02a : Að0;
t Þ ¼ 1 þ 0:0263
t; Test R02b : Að0;
t Þ ¼ 1  0:0870
t þ 0:0113
t 2 :

ð126Þ ð127Þ

For test R02a, the shoreline reflection coefficient was constant, R(L,
t ) = 0.240e1.05i. However, for test R02b, the phase of the shoreline reflection coefficient RL was observed to change with time according to the empirical formula qR ð L;
t Þ ¼ 0:0706 þ 0:2214
t  0:0284
t 2 ;

ð128Þ

while the magnitude of RL remained constant (jRLj = 0.236) in time. [80] In earlier numerical simulations without the return current, the extents of crests and troughs were reasonably well predicted by the present model, but not the position of the sandbars relative to the wave envelope. For best fit with measurements, we took the coefficient L to be 0.28 for test R02a and 0.31 for test R02b. With the fitting coefficient L, the parameters in Tables 1 –3, the measured ripple amplitude profile, and the specified forms of A(0,
t ) and RL (
t ), predictions were made for each laboratory test. The predicted and measured sandbar profiles for both tests R02a and R02b are in satisfactory agreement as shown in Figures 8 and 9, respectively. [81] The Bragg resonance mechanism is important in tests R02a and R02b.( The predicted maximum of the ) coupling coefficient = e h½1 eiqR in the energy exchange equations (42) and (43) increased from 0 (initially with no bars) to approximately 0.29 for test R02a and 0.15 for test R02b. The physical effect of the bed stress due to the return flow is to move the sandbar crests seaward of the wave nodes, rendering the Bragg resonance coupling

15 of 23

C07022

HANCOCK ET AL.: SANDBAR FORMATION UNDER SURFACE WAVES

C07022

ancies in Figure 10 between the predicted and observed incident wave amplitudes. Another indicator of Bragg scattering is that it causes the reflection coefficient R to increase seaward, as is seen in Figure 10, whereas dissipation causes R to decrease seaward. [83] For comparison, the earlier theory of Yu and Mei [2000a] has also been used by Hancock [2005] to simulate tests R02a and R02b. The measured final sandbar height was used in Yu and Mei’s model to find the fitted diffusivity. The return current was not accounted for and the predicted sandbar crests appeared shoreward of the nodes of the surface wave envelope, rather than seaward. [84] With the inclusion of the mean bed stress due to the return current, the predicted sandbar crests are indeed seaward of the wave nodes, as observed. Further improvement to the model for the return current and the associated mean shear stress is desirable in order to improve quantitative agreement, which is poorer near the beach where the presence of an undertow due to wave breaking is likely.

9. Reconstruction of Sandbars in Cape Cod Bay Figure 7. Measured and predicted sandbar height evolution in time for tests R10, R02a, and R02b. The measured sandbar heights correspond to the first (open circle), second (plus), third (asterisk), fourth (open square), and fifth (open triangle) bars, counted from the left in Figures 6, 8, and 9, with the first crest near x0 = 2 m for test R10 and near x0 = 4 m for tests R02a,b. Solid lines are the predicted sandbar 0 0 h ). For test R10, the solid line is the height max(e h )– min(e predicted maximum sandbar height across the sandbar patch. For tests R02a and R02b, the solid lines are the height predictions for the first and fourth bars. The variation in the experimental data is of the order of the ripple height hr* listed in Table 1. ( ½1 iq ) h e R negative and causing energy to be coefficient = e transferred from the incident wave to the reflected wave. Thus, once sandbars begin to form, the incident wave amplitude A and reflection coefficient R increase in magnitude seaward, as is evident in both the observations and predictions shown in Figure 10. Initially, the reflection coefficient R = RL and incident amplitude A = 1 were constant across the entire sand bed 0 < x < L. At the end of tests R02a and R02b, the reflection coefficients at the seaward end of the sandbar patch (x = 0) had increased by 36% and 29%, respectively, above jRLj near the beach (x = L). To summarize, Bragg scattering causes the seaward increase of A and R, and hence the total wave envelope of the standing waves, which in turn causes the sandbar amplitudes to increase seaward. [82] Shoreward attenuation of waves is also due to dissipation in the bottom boundary layer. Kajiura [1968] estimated that waves attenuate on the spatial scale l*/(k*d*), where k*d* = O(e2). Since we have neglected terms of O(e2) and higher, we have not included this effect in our model. In our experiments (see Table 3), k*d * ranged from 0.02 to 0.03 and thus the attenuation distance due to bottom dissipation was 33 to 50 wavelengths. Fewer than three wavelengths spanned the sand bed in our wave tank. Weak attenuation due to dissipation likely accounts for the slight discrep-

[85] Thus far we have been unable to find comprehensive field records of concurrent measurements of both waves and sandbar growth. Dolan and Dean [1985] made an extensive survey of longshore sandbars in Chesapeake Bay, Maryland, but without records of either waves or sandbar evolution. Recently, Elgar et al. [2003] reported wave measurements over existing sandbars in Cape Cod Bay, near Truro, Massachusetts, and compared the wave records with Yu and Mei [2000b]’s theory of Bragg resonance over fixed bars. The wave records were taken with bottommounted pressure gages and acoustic Doppler current meters at five cross-shore locations during August to October of 2001. The sandbars are known to have been present for a long time and remained unchanged during the wave measurements. [86] During a storm, Elgar et al. obtained a 1 hour record of very large waves of frequency 0.042 Hz. The waves broke over the first sandbar and lost 30% of p their ffiffiffiffiffiffiffi energy (and hence their height decreased by a factor of 0:7). They then proceeded unbroken until the shoreline where 80% of the energy was dissipated (implying jRLj2 0.2). The significant wave height, defined as four times the standard deviation of the surface fluctuations, was 0.8 m at the seaward sensor in water of depth H* = 2.80 m. [87] In order to see the possible relevance of the present model, we used the measured wave parameters of Elgar et al. to simulate sandbar growth on the mean topography of the Truro beach. We fit a smooth profile to the observed beach bathymetry to obtain the initial barless depth H0(x0). There were no sandbars in water deeper than 3.2 m. Since our model is only valid for non-breaking waves, we chose the initial incident wave height for our predictions to be that of the waves after breaking over the first sandbar, i.e., pffiffiffiffiffiffi ffi 2A*(1 + jRLj) = 80 0:7 cm at the depth H* = 2.80 m. Sandbar formation is simulated over a horizontal distance L0 = 364 m from an offshore depth of H0(0) = 4.02 m to a shallow depth of H0(L0) = 2.25 m. Invoking the conservation of wave energy flux we estimate the dimensionless incident wave amplitude at the offshore depth H0(0) to be qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 0 0 A(0,
t ) = A (0,
t )/A* = Cg =Cg0 = 0.92, where Cg* and Cg0 *

16 of 23

C07022

HANCOCK ET AL.: SANDBAR FORMATION UNDER SURFACE WAVES

Figure 8. Measured and predicted wave amplitudes and sandbar elevation for test R02a at various times. Caption details are given in Figure 6, except the origin is 2.80 m from the mean wavemaker position.

Figure 9. Measured and predicted wave amplitudes and sandbar elevation for test R02b at various times. Caption details are given in Figure 6, except the origin is 3.50 m from the mean wavemaker position. 17 of 23

C07022

C07022

HANCOCK ET AL.: SANDBAR FORMATION UNDER SURFACE WAVES

Figure 10. Measured (open circle, cross) and predicted (dash, solid) dimensionless incident wave amplitude jAj = jA0j/A* and reflection coefficient jRj at the beginning (open circle, dash), t0 = 0, and end (cross, solid), t0 = 3.02 and t0 = 3.96 days, of tests R02a and R02b, respectively. The origin is 2.80 m (R02a) and 3.50 m (R02b) from the mean wavemaker position.

Figure 11. Comparison of our predictions with the observations of Elgar et al. [2003] (open circle) of sandbars in Cape Cod Bay. The predictions are based on the parameters listed in Tables 1 – 3. The 0 h (bottom) are amplitude of the first wave harmonic jz 1[1]0j (top) and seabed profiles z0 = h0 = H0 + e given after 8 hours (dash) and 16 hours (solid) of wave action. The initial beach profile z0 = H0 is indicated by the dotted line. 18 of 23

C07022

C07022

HANCOCK ET AL.: SANDBAR FORMATION UNDER SURFACE WAVES

are the group velocities at the characteristic depth H* and at the offshore depth H0(0), respectively. The phase of the shoreline reflection coefficient RL was not measured; we adjusted its value to align the predicted and measured sandbar crest positions. The dimensionless offshore ampli
t )ffi = 0.92 and the shoreline reflection coefficient tude A(0, pffiffiffiffiffiffi RL = 0:2e2.51i provided boundary values for the wave amplitude A and reflection coefficient R in our sandbar simulation. [88] Elgar et al. reported that small sand ripples comparable in length to wave orbits were observed by SCUBA divers on the sandbar crests, and again when the sandbar crests were exposed during the spring low tide. From a photograph of a sandbar crest (Figure 1c in Elgar et al.), the typical ripple height was approximately 4 – 6 cm and length 10– 20 cm. This is consistent with the empirical formulas of Nielsen [1981] and Wikramanayake and Madsen [1990] for ripples in the field under irregular waves. The larger of these predictions gives hr0 = hr* = 5.8 cm at the characteristic depth H* [Hancock, 2005], which we use here. [89] The following values representative of seawater are used in all the predictions: sediment specific gravity s = 2.57, fluid viscosity n = 0.0115 cm2/s, and water density r =1.03 g/cm3. [90] The water depth is relatively shallow (k*H* = 0.14), and hence the term proportional to 1/sinh2 kH dominates the bed load forcing in equation (114). In particular, the fitting parameter L in the mean stress becomes immaterial since the return flow is overshadowed by the factor 1/sinh2 kH. [91] On the basis of the parameters listed in Tables 1 –3 and the fitted mean depth H(x), the predictions made of the beach profiles after 8 and 16 hours of wave action are shown in Figure 11. No computations have been made for longer times since the assumption of constant wave amplitude is unlikely valid during a storm typical in Cape Cod. Otherwise, the agreement between the predicted and observed beach profiles in the length, height, and location of the sandbars is remarkable, despite the crude estimates of the wave parameters, the uncertainty of the phase of the shoreline reflection, and the duration of wave action. This agreement suggests that waves of similar characteristics, which may occur frequently in this area, may be responsible for the sustained presence of these sandbars.

10. Concluding Remarks [92] In this work we have developed a theory and performed new experiments on the wave-induced formation of sandbars. The work makes significant improvements to the wave-induced forced-diffusion sandbar model of Yu and Mei [2000a] to achieve better quantitative agreement between theory and laboratory experiments. We neglect suspended sediment and limit our study to cases where the sediment transport is dominated by bed load. While still relying on empirical relations between bed shear and sediment discharge, we have replaced the constant eddy viscosity model used by Yu and Mei with a depth-linear model to analyze the boundary layer flow to the second order in wave slope. Ripples, which form quickly, are treated as the known roughness parameter in the boundary layer.

C07022

[93] In order to match predictions with observations, two further improvements were necessary. First, use of the skinfriction Shields parameter in the empirical formula for bed load transport enabled accurate prediction of the flat subcritical regions present under high shoreline reflection. Second, making the slope adjustment to the bed load transport rate instead of the bottom shear stress allowed us to accurately predict sandbar heights and avoid the use of an additional fitting coefficient proportional to sandbar height present by Yu and Mei [2000a]’s model. [94] For high shoreline reflection, our model needs no fitting parameters and predictions agree well with laboratory experiments. For weakly reflected waves, the return current alters the bed shear stress and affects the position of the sandbars relative to the profile of the wave envelope. A better turbulence model for modeling this current and its shear stress is desirable. [95] Yu and Mei [2000a] showed that Bragg resonance between the surface waves and sandbars was tied to the size and sign of the quantity =fe h½1 eiqR g. We have shown that ½1 iqR e while =fh e g is small under strongly reflected standing waves (e.g., test R10) where Bragg resonance is weak, it is larger for weakly reflected waves (tests R02a and R02b) where Bragg resonance causes the shoreward attenuation of both incident and reflected waves and sandbars. The effect of Bragg resonance is expected to be more dramatic, even for high reflection, over longer barred regions. [96] In addition to testing our model against laboratory experiments, numerical simulation of a field record from Cape Cod Bay provides partial evidence that the present model can be relevant to nature. [97] Further laboratory data on a sloping seabed are desirable for testing our theory. So far only Dulou et al. [2000] have measured sandbars on a gently sloping bed. Their short tank of length 4.7 m, width 39 cm, and depth 15 cm housed a bed of sand grains of diameter 0.08 mm. The mean water depth over the bed varied between 4 to 8 cm and the total wave height ranged from 1 to 2 cm. Under monochromatic waves, half wavelength sandbars formed along the bed, with depth-dependent length ranging from 20 to 23 cm. Under the wave antinodes, sandbar crests of approximate height 1.2 cm appeared. Under the wave nodes, sandbar troughs covered with 0.12 to 0.14 cm ripples were found. These observed features are characteristic of the transport of fine grains in suspension [de Best et al., 1971; Xie, 1981; O’Hare and Davies, 1993] and can be anticipated by the theory of mass transport [Carter et al., 1973]. Furthermore, based on the wave and sediment data, the ratio of the friction velocity uf* to the sediment fall velocity wS was between 5.4 and 5.7, almost twice that (2.7 to 3.1) in our laboratory experiments for coarse grains. Since our theory neglects suspended sediment transport, valid comparison must await new experiments on sloping beds of coarse grains in a large flume. [98] Further effort on the modeling of sandbar formation in nature should of course include the transport of suspended sediment, the randomness of sea waves, and the additional effects of currents. A greater challenge is to construct a fundamental theory describing the two-phase mechanics of particles and water with less reliance on

19 of 23

HANCOCK ET AL.: SANDBAR FORMATION UNDER SURFACE WAVES

C07022

purely empirical relations between the sediment transport rate and the bed shear.

Z

Zb

Appendix A: Correction Factor in the Skin-Friction Shields Parameter [99] The empirically fitted bed load transport formula devised by Sleath [1978] and Nielsen [1992] for use in oscillatory flow uses a Shields parameter based on the grain roughness, where the Nikuradse roughness kN0 = 30hb* is replaced by 2.5d and consequently y* in (79) is redefined as kU 12kU * *: ¼ y ¼ * ð2:5d=30Þw dw

ðA1Þ

From (80), the ratio of the grain roughness Shields parameter to Q* is 

12kU

f

dw



rd ¼ f

*

kU

* *

 d 12hb

 ’

hb w

!0:6370:0433 log10

 12d hb

*

Z

Z

1

1 Zb

Z

Z

1

Zb

1  i pffiffiffi K0 ðZ
Þd Z
¼ pffiffiffi Z K1 ðZ Þ; 2

ðB7Þ

b

c3 ðZb ; Zb Þ ¼

pffiffiffiffiffi pffiffiffiffiffi ð1 þ iÞ Zb K*1 ðZb Þ ð1  iÞ Zb K1 ðZb Þ pffiffiffi  : pffiffiffi * 2K0 ðZb Þ 2K0 ðZb Þ

b

dK0 ðZ  Zb Þ dZ ¼  dZ

0

Z

ðB8Þ

Equations (B7) and (B8) are clearly imaginary, as written in equation (100) in section 4.2.2. [102] The integrals in the second harmonic of the second order bed shear stress, equation (108), are now simplified and written in terms of three integrals An (Zb) that must be computed numerically: ! pffiffiffiffiffi K1 ð2Zb Þ K0 ð2Z Þ K1 ðZb Þ c4 ðZ; Zb Þ dZ ¼ 2ð1  iÞ Zb  pffiffiffi K0 ð2Zb Þ K0 ð2Zb Þ 2K0 ðZb Þ

1 Zb

 A1 ðZb Þ  A2 ðZb Þ; Z

1 Zb

ðB9Þ

pffiffiffi 2Zb K1 ðZb Þ K0 ð2Z Þ c5 ðZ; Zb Þ dZ ¼  K0 ð2Zb Þ K0 ðZb Þ

! K1 ð2Zb Þ K1 ðZb Þ   pffiffiffi K0 ð2Zb Þ 2K0 ðZb Þ pffiffiffiffiffi ð1 þ iÞ Zb K1 ðZb Þ pffiffiffi þ 2K0 ðZb Þ  ðA1 ðZb Þ þ A2 ðZb Þ  A3 ðZb ÞÞ; ðB10Þ

ðB1Þ

!  K0 ðZ Þ 2 K21 ðZb Þ ; dZ ¼ Zb 1 þ 2 K0 ðZb Þ K0 ðZb Þ

0

Zb

0

Z

    pffiffiffiffiffi  K0 ðZ Þ 2  dZ ¼ Zb < eip=4 K1 ðZb Þ ;   K ðZ Þ  K ðZ Þ 1

pffiffiffiffiffi pffiffiffiffiffi ð1 þ iÞ Zb K*1 ðZb Þ ð1  iÞ Zb K1 ðZb Þ p ffiffi ffi  pffiffiffi 2 2K0 ðZb Þ 2 2K0* ðZb Þ    K1 ðZb Þ 2  ; þ iZb  K ðZ Þ 



Zb

Z

c2 ðZb ; Zb Þ ¼

*

[101] In this appendix, we calculate the functions cn(Zb, Zb), n = 1, 2, 3, present in equation (98) for the Reynolds stress inside the boundary layer, and integrals of the functions cn(Z, Zb), n = 4, 5, 6, found in equation (108) for the second order second harmonic shear stress. First, equations (87) – (89) are integrated from Z = Zb to 1 to obtain the functions cn(Zb, Zb), n = 1, 2, 3. The integration involves the following six integrals: 1

The integrals (B4) and (B6) are calculated by parts. The definition of F(Z, Zb), equation (73), was used in (B6). Integrating equations (87) – (89) from Z = Zb to 1 and substituting (B1) – (B6) gives equation (99) for c1(Zb, Zb) in section 4.2.2, and

ðA2Þ

The Functions cn(Z, Zb)

Z

1 F1 ðZ; Zb Þ d*K0* F1 ðZ; Zb ÞK*0 ðZ Þ  dZ ¼  Zb K* ðZb Þ K0* ðZb Þ dZ Z 10 @F1 ðZ; Zb Þ K0 ðZ Þ dZ  @Z Z K0* ðZb Þ !  Z 1   K0 ðZ Þ 2 K0* ðZ Þ   ¼ dZ:ðB6Þ  K ðZ Þ   * 0 b Zb K0 ðZb Þ

kU 2  * dw

[100] The approximation (77) to f(y) has been used to b obtain the approximate formula on the right. The term rdQ is therefore used as the amplitude of the grain roughness (or skin-friction) Shields parameter in the empirical formula (4) for the bed load transport rate. The ratio rd is plotted vs. d/hb* in Figure 2 for typical values of kU* /(dw) ranging from 50 to 500. Since rd can be as low as 0.05, this correction can affect the quantitative predictions significantly.

Appendix B:

1

C07022

ðB2Þ Z

1

c6 ðZ; Zb Þ

ðB3Þ

Zb

K0 ð2Z Þ dZ ¼ A2 ðZb Þ; K0 ð2Zb Þ

ðB11Þ

b

where 1

K0 ðZ ÞdZ;

ðB4Þ

Zb

Z

1

K20 ðZ ÞK0 ð2Z Þ dZ; K20 ðZb ÞK0 ð2Zb Þ

ðB12Þ

F1 ðZ; Zb ÞK0 ð2Z Þ dK0 dZ; K0 ðZb ÞK0 ð2Zb Þ dZ

ðB13Þ

A1 ðZb Þ ¼ Zb

!

pffiffiffiffiffi K1 ð2Zb Þ K0 ðZ ÞK0 ð2Z Þ K1 ðZb Þ ; dZ ¼ ð1  iÞ Zb  pffiffiffi K0 ðZb ÞK0 ð2Zb Þ K0 ð2Zb Þ 2K0 ðZb Þ ðB5Þ 20 of 23

Z

1

A2 ðZb Þ ¼ Zb

HANCOCK ET AL.: SANDBAR FORMATION UNDER SURFACE WAVES

C07022

Z

1

A3 ðZb Þ ¼ Zb

ðZ  Zb ÞK0 ð2Z Þ dK0 dZ: K0 ðZb ÞK0 ð2Zb Þ dZ

ðB14Þ

Equation (B5) was also used to derive equations (B9) – (B11). The range 0.00019 < Zb < 1.02 corresponds to the range of validity 2 < y < 105 of the approximate formula (77) for the friction factor f(y). For this range of Zb, the integrals in An (Zb) need only be computed to an upper limit of 25 to obtain relative errors less than 5  105, and 100 for relative errors less than 1010.

Appendix C: Stress

Second Harmonic Velocity and

[103] The homogeneous solutions of (106) are K0 (2Z) and I 0 (2Z), pffiffiffiwhere K0 is defined in (67) and I 0 (Z) = I0(2 eip/4 Z ), where I0 is the modified Bessel function of the first kind of order zero. By the method of variation of parameters [e.g., Hildebrand, 1964, section 1.9], the full solution of (106) is ½2 U2



½2 u2

Z



t ÞI 0 ð2Z
Þ RHS½2 ð x; Z;

½K0 ð2Z
Þ; I 0 ð2Z
Þ d Z ZW

Z

¼ K0 ð2Z Þ Zb

Z

Z

þ I 0 ð2Z Þ Zb

ðC1Þ

where the Wronskian is W ½K0 ð2Z Þ; I 0 ð2Z Þ ¼ 2K0 ð2Z ÞI 00 ð2Z Þ  2K00 ð2Z ÞI 0 ð2Z Þ R Z 1 1 1 ¼ e Z dZ ¼ : ðC2Þ 2 2Z

The integration constants c7 and c8 are found by imposing the boundary conditions (107). As Z!1, K0 (2Z)!0 and [2] as jI 0 (2Z)j!1 and the matching condition u[2] 2 !U2 Z!1 implies, with (C2), Z

1

c8 ¼ 2

RHS½2 ð x; Z;
t ÞK0 ð2Z ÞdZ:

ðC3Þ

Zb

As Z ! Zb, the integrals in (C1) vanish and the no-slip condition u[2] 2 = 0 at Z = Zb implies ½2

c7 ¼

U2  c8 I 0 ð2Zb Þ : K0 ð2Zb Þ

ðC4Þ

Substituting c7, c8 from (C3), (C4), and the Wronskian from (C2) into (C1) gives ½2

u2 ¼

 K0 ð2Z Þ ½2 U K0 ð2Zb Þ 2 Z 2K0 ð2Z ÞI 0 ð2Zb Þ 1  RHS½2 ð x; Z;
t ÞK0 ð2Z ÞdZ K0 ð2Zb Þ Zb Z Z

t ÞI 0 ð2Z
Þd Z
þ 2K0 ð2Z Þ RHS½2 ð x; Z; Zb Z 1

t ÞK0 ð2Z
Þd Z:
þ 2I 0 ð2Z Þ RHS½2 ð x; Z;

Note that the lower limits of the first and third integrals are different, Z and Zb. The integrals converge as Z!1 since RHS[2](x, Z,
t ) contains only Kn (2Z) and not I n (2Z). [104] To find the bottom shear stress, we differentiate (C5) in Z, multiply by Z, and set Z = Zb to obtain ½2  ½2 pffiffiffiffiffi @u2  ð1 þ iÞU2 Zb K1 ð2Zb Þ Z ¼ @Z Zb K0 ð2Zb Þ   þ 2Zb W ½K0 ð2Z Þ; I 0 ð2Z Þ Zb Z 1 K ð 2Z Þ 0  RHS½2 ð x; Z;
t Þ dZ ð 2Z K 0 bÞ Zb ½2 pffiffiffiffiffi ð1 þ iÞU2 Zb K1 ð2Zb Þ ¼ K0 ð2Zb Þ Z 1 K0 ð2Z Þ þ RHS½2 ð x; Z;
t Þ dZ: K0 ð2Zb Þ Zb

ðC6Þ

The second equality in (C6) follows from (C2). Substituting RHS[2](x, Z,
t ) from (91) into (C6) gives (108), where the integrals are calculated in (B9) – (B11).

Appendix D: Cancellation of Free Second Wave Harmonic



t ÞK0 ð2Z
Þ RHS½2 ð x; Z;

½K0 ð2Z
Þ; I 0 ð2Z
Þ d Z ZW

þ c7 K0 ð2Z Þ þ c8 I 0 ð2Z Þ;

C07022



1

[105] For a flume of uniform depth, the free second harmonic can be eliminated by adding a second harmonic component to the displacement of the wavemaker piston [Madsen, 1971; Mei et al., 2005, Part II, p 567ff]. [106] In our wave flume, a sloping ramp was present near the wavemaker where the depth was larger than the constant depth of the main section. To avoid a complicated theory for the uneven bottom, we empirically eliminated the free second harmonic in the main section by adjusting, in sequence, the amplitude and phase of the second harmonic piston motion. During the elimination, the depth in the main section was kept constant so that the wave parameters there were also constant. the free second [2]0 harmonic [107] The complex amplitude of [2]0 ik[2]0x0 [2] ik x0 + R e ). The wave has the form z [2]0 2F = A (e 0 complex incident amplitude A[2] and reflection coefficient R[2] were calculated by combining wave gage records at several stations along the tank with equation (30) for the bound second harmonic wave. three steps. First, the free [108] The elimination followed 0 second harmonic wave z [2] 2(ORIG) generated from monochromatic piston motion was measured. Second, a second harmonic component was added to the piston motion and 0 the free second harmonic wave z [2] 2(NEW) was remeasured. The amplitude of the second harmonic piston motion was adjusted to match the magnitudes of the incident amplitudes 0 [2]0 [2]0 second harmonic z  z (jA[2] j) of the added 2(NEW) 2(ORIG) and [2]0 the original z 2(ORIG). Finally, the phase of the second harmonic piston motion was adjusted to cancel the free second harmonic wave. The success of this process can be seen in the top panels of Figures 6 through 9 where the amplitude of the second harmonic is not modulated along the tank. Further details are given by Hancock [2005].

ðC5Þ

Z

21 of 23

C07022

Appendix E:

HANCOCK ET AL.: SANDBAR FORMATION UNDER SURFACE WAVES

Hv k0(x10) qB0(x0, t0) QB0(x0,
t 0 )

Symbols

C07022

Heaviside unit step function wave number rate of bed load transport rate of bed load transport over a half cycle dimensionless forcing term in the sandbar equation (114) local reflection coefficient reflection coefficient at the shoreward end of the sandbar patch ratio of skin friction Shields parameter to the full Shields parameter (see Appendix A) sediment specific gravity wave phase time coordinate slow sandbar growth time amplitude of horizontal component of velocity just above the boundary layer friction velocity horizontal and vertical velocity components in the boundary layer component of velocity orthogonal to ripple-averaged seabed horizontal coordinate slow horizontal coordinate normalized boundary layer coordinate normalized roughness height

X0, X*, X = X0/X* respectively the dimensional form, characteristic value, and normalized form of X b amplitude of a time periodic X qt (x,
t ) quantity X X[n] the n’th time harmonic amplitude R(x1,
t ) = B/A m of the m’th order perturbation of X RL(
t ) X* the complex conjugate of X < and = the real and imaginary parts of a complex quantity rd RHS right-hand side d* characteristic boundary layer thickness s e = k*Ab* dimensionless measure of wave S(x0,
t 0 ) slope t0
t 0 h0(x0,
t 0 ) = z0  h0(x0,
t 0 ) height from ripple-averaged bed surface U0(x0,
t 0 ) 0 0
0 hr (x , t ) ripple height hb0(x0,
t 0 ) = 4hr0/30 roughness height k = 0.41 Ka´rma´n constant uf0(x0,
t 0 ) 0 0 0 0 0 0 0 0 L fitting coefficient of order unity u (x , z , t ), w (x , z , t ) in the second order mean bottom shear stress due to the return wN0(x0, z0, t0) current N seabed porosity, taken as 0.3 x0 0 0
0 n e (x , t ) eddy viscosity x1 = ex w wave angular frequency Z W0(x0) Bragg resonance coupling coefficient Zb(x,
t ) f0(x0, z0, t0) velocity potential v(x,
t ) phase lag of bed stress behind flow [109] Acknowledgments. The authors wish to thank the US Office of Naval Research for sustained financial support through several grants ym angle of repose, taken as 30° (N00014-89-J-3128, N00014-04-1-0077, and N00014-95-1-0840 (ASr fluid density SERT) directed by Thomas Swean). Partial support from the US National [0]
t 2W(x, t ) wave-induced mean shear stress Science Foundation (CTS 0075713 directed by M. Plesniak, Roger Arndt, and John Foss) is also acknowledged. We thank Ole Madsen for a number at second order t ) return-current induced mean shear of valuable discussions, including the estimate of L, and Steve Elgar and t [0] 2C(x,
Britt Raubenheimer for providing additional field data. We also thank stress at second order Vladimir Barzov and Yile Li for their help with our laboratory experiments. 0 0 0 0 0 0 0 t (x , z , t ), t b (x , t ) shear stress and bottom shear References stress qR(x1,
t ) phase of local reflection coefficient Abramowitz, M., and I. A. Stegun (1972), Handbook of Mathematical Functions, Dover, Mineola, N. Y. Q(x0, t0) Shields parameter Ahrens, J. P. (2000), A fall-velocity equation, J. Waterw. Port Coastal QC critical Shields parameter for Ocean Eng., 126, 99 – 102. incipient sediment motion Blondeaux, P. (1990), Sand ripples under sea waves. Part I: Ripple forma0 0 0 tion, J. Fluid Mech., 218, 1 – 17. z (x , t ) free surface displacement Boczar-Karakiewicz, B., D. L. Forbes, and G. Drapeau (1995), Nearshore Ab* characteristic bottom orbital bar development in southern Gulf of St. Laurence, J. Waterw. Port amplitude Coastal Ocean Eng., 121(1), 49 – 60. 0 0
0 0 0
0 Carter, T. G., L.-F. Liu, and C. C. Mei (1973), Mass transport by waves A (x1, t ), B (x1, t ) local incident and reflected and offshore sand bedforms, J. Waterw. Port Coastal Ocean Eng., 99, complex wave amplitudes 165 – 184. 0 0 Cg(x1) group velocity Davies, A. G. (1982), The reflection of wave energy by undulations on the d mean sediment diameter seabed, Dyn. Atmos. Oceans, 6, 207 – 232. Davies, A. G., and C. Villaret (1999), Eulerian drift induced by progressive Dn(x,
t ) dimensionless diffusivity in the waves above rippled and very rough beds, J. Geophys. Res., 104, 1465 – sandbar equation (114) 1488. f friction factor de Best, A., E. W. Bijker, and J. E. W. Wichers (1971), Scouring of a sand bed in front of a vertical breakwater, in Proc. Int. Conf. Port Ocean g acceleration of gravity Engineering Under Arctic Conditions, vol. II, pp. 1077 – 1086, Tech. h00(x0,
t 0 ) ripple-averaged seabed elevation Univ. of Norway, Trondheim, Norway. 0
0 e h0 (x , t ) sandbar elevation Devillard, P., F. Dunlop, and B. Souillard (1988), Localization of gravity e waves on a channel with a random bottom, J. Fluid Mech., 186, 521 – h[1] (x10,
t 0 ) first spatial0 harmonic of sandbar 538. elevation e h Dolan, T. J., and R. G. Dean (1985), Multiple longshore sand bars in the H0(x10) mean (bar-averaged) depth upper Chesapeake Bay, Estuarine Coastal Shelf Sci., 21, 727 – 743. 22 of 23

C07022

HANCOCK ET AL.: SANDBAR FORMATION UNDER SURFACE WAVES

Dulou, C., M. Belzons, and V. Rey (2000), Laboratory study of wave bottom interaction in the bar formation on an erodible sloping bed, J. Geophys. Res., 105, 19,745 – 19,762. Elgar, S., B. Raubenheimer, and T. H. C. Herbers (2003), Bragg reflection of ocean waves from sandbars, Geophys. Res. Lett., 30(1), 1016, doi:10.1029/2002GL016351. Fredsøe, J. (1974), On the development of dunes in erodible channels, J. Fluid Mech., 64, 1 – 16. Grant, W. D., and O. S. Madsen (1979), Combined wave and current interaction with a rough bottom, J. Geophys. Res., 84, 1797 – 1808. Grant, W. D., and O. S. Madsen (1982), Movable bed roughness in unsteady oscillatory flow, J. Geophys. Res., 87, 469 – 481. Hallermeier, R. J. (1981), Terminal settling velocity of commonly occurring sand grains, Sedimentology, 28, 859 – 865. Hancock, M. J. (2005), Generation of sand bars under surface waves, Ph.D. thesis, Dept. of Civil and Environmental Engineering, Mass. Inst. Tech., Cambridge, Mass. Hara, T., and C. C. Mei (1987), Bragg scattering of surface waves by periodic bars: Theory and experiment, J. Fluid Mech., 178, 221 – 241. Heathershaw, A. D. (1982), Seabed-wave resonance and sand bar growth, Nature, 296, 343 – 345. Herbich, J. B., H. D. Murphy, and B. Van Weele (1965), Scour of flat sand beaches due to wave action in front of sea walls, in Int. Conf. Coastal Eng. Santa Barbara Specialty Conference, pp. 705 – 726, ASCE, Reston, Va. Hildebrand, F. B. (1964), Advanced Calculus for Engineers, Prentice-Hall, Englewood Cliffs, N. J. Jan, C.-D., and M.-C. Lin (1998), Bed forms generated on sandy bottom by oblique standing waves, J. Waterw. Port Coastal Ocean Eng., 124, 295 – 302. Kajiura, K. (1968), A model of the bottom boundary layer in water waves, Bull. Earthquake Res. Inst. Univ. Tokyo, 46, 75 – 123. King, D. B. (1991), The effect of beach slope on oscillatory flow bedload transport, in Coastal Sediments ’91, pp. 734 – 744, ASCE, Reston, Va. Komarova, N. L., and A. C. Newell (2000), Nonlinear dynamics of sand banks and sand waves, J. Fluid Mech., 415, 285 – 321. Landry, B. J. (2004), Bathymetric evolution of sand bed forms under partially standing waves, S.M. thesis, Dept. of Civil and Environmental Engineering, Mass. Inst. Tech., Cambridge, Mass. Landry, B. J., M. J. Hancock, and C. C. Mei (2007), Note on sediment sorting in a sandy bed under standing water waves, Coastal Eng., 54, 694 – 699. Longuet-Higgins, M. S. (1950), Theory of origin of microseisms, Proc. R. Soc., Ser. A, 243, 1 – 35. Madsen, O. S. (1971), On the generation of long waves, J. Geophys. Res., 76, 8672 – 8683. Madsen, O. S. (1994), Spectral wave-current bottom boundary layer flows, in Int. Conf. Coastal Eng., pp. 384 – 398, ASCE, Reston, Va. Mei, C. C. (1985), Resonant reflection of surface water waves by periodic sandbars, J. Fluid Mech., 152, 315 – 335. ¨ nlu¨ata (1972), Harmonic generation in shallow water Mei, C.C., and U. U waves, in Waves on Beaches, edited by R. E. Meyer, pp. 181 – 202, Elsevier, New York. Mei, C. C., and J. Yu (1997), The instability of sand ripples under partially standing surface waves, Phys. Fluids, 9, 1606 – 1620.

C07022

Mei, C. C., M. Stiassnie, and D. K.-P. Yue (2005), Theory and Applications of Ocean Surface Waves, World Sci., Hackensack, N. J. Naciri, M., and C. C. Mei (1988), Bragg scattering of water waves by a doubly periodic seabed, J. Fluid Mech., 192, 51 – 74. Nielsen, P. (1981), Dynamics and geometry of wave generated ripples, J. Geophys. Res., 86, 6467 – 6472. Nielsen, P. (1992), Coastal Bottom Boundary Layers and Sediment Transport, World Sci., Hackensack, N. J. O’Hare, T. J., and A. G. Davies (1993), Sand bar evolution beneath partially-standing waves: Laboratory experiments and model simulations, Cont. Shelf Res., 13, 1149 – 1181. Restrepo, J. M., and J. L. Bona (1995), Three-dimensional model for the formation of longshore sand structures on the continental shelf, Nonlinearity, 8, 781 – 820. Seaman, R. C., and T. O’Donoghue (1996), Beach response in front of wave-reflecting structures, in Int. Conf. Coastal Eng., pp. 2284 – 2297, ASCE, Reston, Va. Sleath, J. F. A. (1978), Measurements of bed load in oscillatory flow, J. Waterw. Port Coastal Ocean Eng., 104, 291 – 307. Sleath, J. F. A. (1990), Seabed boundary layers, in The Sea, vol. 9, edited by B. Le Me´haute´ and D. M. Hanes, pp. 693 – 727, John Wiley, N. Y. Trowbridge, J., and O. S. Madsen (1984a), Turbulent wave boundary layers: 1. Model formulation and first-order solution, J. Geophys. Res., 89, 7989 – 7997. Trowbridge, J., and O. S. Madsen (1984b), Turbulent wave boundary layers: 2. Second-order theory and mass transport, J. Geophys. Res., 89, 7999 – 8007. Vittori, G., and P. Blondeaux (1990), Sand ripples under sea waves. Part II: Finite amplitude development, J. Fluid Mech., 218, 19 – 39. Wikramanayake, P. N., and O. S. Madsen (1990), Calculation of movable bed friction factors, in Technical Progress Report 2, Dredging Research Program, U. S. Army Corps of Eng., Coastal Eng. Res. Center, Vicksburg, Miss. Xie, S.-L. (1981), Scouring patterns in front of vertical breakwaters and their influences on the stability of the foundations of the breakwaters, in Technical Report, Coastal Engineering Group, Dept. of Civil Engineering, Delft Univ. of Technol., Delft, Netherlands. Yu, J., and C. C. Mei (2000a), Formation of sand bars under surface waves, J. Fluid Mech., 416, 315 – 348. Yu, J., and C. C. Mei (2000b), Do longshore bars protect the shore?, J. Fluid Mech., 404, 251 – 268. 

M. J. Hancock, Department of Mathematics, Massachusetts Institute of Technology, Room 2-236, Cambridge, MA 02139, USA. (hancock@ alum.mit.edu) B. J. Landry, Department of Civil and Environmental Engineering, University of Illinois, Urbana-Champaign, Hydrosystems Laboratory, Room 1534, 205 North Mathews Ave., Urbana, IL 61801, USA. ([email protected]) C. C. Mei, Department of Civil and Environmental Engineering, Massachusetts Institute of Technology, Room 48-413, Cambridge, MA 02139, USA. ([email protected])

23 of 23