Surf Zone Turbulence and Suspended Sediment Dynamics—A Review

Abstract

The existence of sandy beaches relies on the onshore transport of sand by waves during post-storm conditions. Most operational sediment transport models employ wave-averaged terms, and/or the instantaneous cross-shore velocity signal, but the models often fail in predictions of the onshore-directed transport rates. An important reason is that they rarely consider the phase relationships between wave orbital velocity and the suspended sediment concentration. This relationship depends on the intra-wave structure of the bed shear stress and hence on the timing and magnitude of turbulence production in the water column. This paper provides an up-to-date review of recent experimental advances on intra-wave turbulence characteristics, sediment mobilization, and suspended sediment transport in laboratory and natural surf zones. Experimental results generally show that peaks in the suspended sediment concentration are shifted forward on the wave phase with increasing turbulence levels and instantaneous near-bed sediment concentration scales with instantaneous turbulent kinetic energy. The magnitude and intra-wave phase of turbulence production and sediment concentration are shown to depend on wave (breaker) type, seabed configuration, and relative wave height, which opens up the possibility of more robust predictions of transport rates for different wave and beach conditions.

1. Introduction

The coexistence and interaction of hydrodynamic motions over a wide range of timescales create a rather chaotic surf zone environment. Incoming wind waves and swell operate across a broad spectrum of frequencies, they shoal and break, and they generate time-averaged cross- and longshore currents, long-period infragravity waves, as well as turbulence that emanates from both the sea surface and from the seabed. While surf zone processes can often be described in a deterministic manner, many of the processes interact in complex and sometimes non-linear ways with each other and the seabed, making it challenging to describe aspects of the system other than stochastically. Consequently, detailed model simulations and accurate predictions are difficult. For example, the theoretical foundations for the process of wave breaking are still weak, and the role of surface-generated turbulence on vertical and horizontal mixing and sediment mobilization is not well understood for natural situations.

When surface gravity waves propagate into shallow water, wave energy is dissipated and momentum is transferred into the bottom boundary layer by friction between the fluid and the seabed. Further landward, the waves break in the outer surf zone; the free stream flow changes from (nearly) irrotational to rotational and vortices generated by wave breaking transfer momentum into the water column, sometimes reaching into the bottom boundary layer and affecting the seabed. The force per unit area that the fluid exerts on the sediment bed is the bed shear stress (τ_b). It mobilizes sediment grains and causes them to roll, slide, or saltate along the bed as bedload. With increasing stress, particles collide in the bed layer and some particles may gain sufficient vertical momentum to be suspended in the water column. Alternatively, sediment may be suspended through flow separation around seabed irregularities such as wave ripples, or when upward-directed pressure components within the bed induced by fluid turbulence become so large that the resulting upward force is sufficient for the particles to rise and detach from the bed [1]. This is the case, for example, when breaking-wave generated vortices induce large upward-directed pressure gradients in the seabed [2]. Because wave boundary layers are much thinner than boundary layers associated with time-averaged mean currents, wave-generated bed shear stresses are significantly larger than current-derived stresses. Understanding of the temporal variation of wave-generated τ_b is therefore essential for the prediction of sediment transport under wave motion [3].

Numerical models for nearshore sediment transport are typically underpinned by three assumptions: (i) Sediment is mobilized at the fluid/bed interface by small-scale turbulence generated by shear between the fluid and the seabed. (ii) Consequently, τ_b scales with the horizontal wave orbital velocity (u) and it is in phase with u, and (iii) the sediment load (bedload and suspended load) scales with τ_b. Therefore, net sediment transport, q, also scales with τ_b, $q\propto \left|u{\tau }_b\right|$. However, as pointed out by Foster et al. [4], because τ_b is unknown in most cases, the models often rely on simple quadratic assumptions (stress scales with velocity-squared) or eddy viscosity models to describe shear stress and sediment load.

Present-day operational sediment transport models typically predict sediment transport due to time-averaged flows quite well, also under breaking-wave conditions. Time-averaged sediment load indeed does tend to scale with shear stress, mean current speeds are, if not trivial, then relatively simple to simulate with some accuracy, and a mean (current-driven) sediment transport is not reliant on time-dependencies. However, the models encounter significant difficulties with predicting oscillatory transport due to wave motions [5,6,7], which is not ideal since this is the most important mechanism for onshore sediment transport involved in beach recovery. In contrast to mean transport, oscillatory transport is time-dependent and must be resolved on intra-wave scale, and intra-wave quantities such as turbulence and bed shear stress are not well defined in most models. More detailed and sophisticated but computationally demanding single-phase RANS-based approaches do exist [8] and these models can resolve flow, turbulence, and sediment transport over wave phase (e.g., the SedWaveFoam two-phase model for sheet flow; [9]) but they have only been applied and tested for nonbreaking waves outside the surf zone. In general, numerical models assume that turbulent shear stresses derive from friction between the (horizontal) oscillatory flow and the seabed, and sediment load and transport therefore scale with higher order moments of u(t). Other sources of bed shear stress, for example, turbulence generated by surface wave breaking are typically not included and presently available models are not, therefore, particularly well-suited for simulation of surf zone sediment transport.

Longo et al. [10] reviewed surf zone turbulence and Sumer and Fuhrman [11] (their chapter 7) recently published an exhaustive treatment of turbulence in coastal waters, while aspects of nearshore sediment transport were reviewed, for example, by Aagaard and Hughes [12]. The present paper aims to provide an up-to-date review of recent experimental advances on turbulence characteristics and generation in the surf zone and its effects on sediment mobilization, suspended sediment load, and wave-driven (i.e., time-dependent) suspended sediment transport. The focus is mainly on intra-wave properties due to their significance to (peak) bed shear stresses and onshore sediment transport, and the review is mainly concerned with experimental results rather than with numerical model simulations. Following this introduction (1), we discuss the following themes: (2) the definition and parameterization of shear stress and its relation to turbulent quantities; (3) current methods to measure turbulence under waves with a particular focus on field conditions; (4) a review of turbulence generation and its spatio-temporal structure with emphasis on surface-generated turbulence in the surf zone; (5) the effect of turbulence on sediment mobilization and sediment concentrations; (6) the intra-wave properties of turbulent kinetic energy production and sediment suspension under wind waves, again with an emphasis on conditions within the surf zone; (7) wave-driven sediment transport as affected by turbulence, its magnitude, and timing within the wave cycle; and finally (8), a summary and discussion of future perspectives in empirical studies and modeling of turbulence and sediment transport.

2. Bed Shear Stress

Bed shear stress is defined as the tangential force that the fluid exerts on the seabed per unit area. While the bed shear stress can rarely be measured directly, various methods exist for its estimation and parameterization [13,14]. The classic expression for τ_b is based on the log-law derived from unidirectional flows:

$$\frac{u}{u_{*}}=\frac{1}{\kappa }\ln \frac{z}{z_0}$$

(1)

where $u_{*}$ is the friction velocity, κ is the von Karman constant, z is the elevation (within, or at the top of the boundary layer) at which the two-dimensional horizontal fluid velocity vector, u, is specified and z₀ is the hydraulic bed roughness. Measuring the horizontal velocity profile within the boundary layer provides an estimate of $u_{*}$, and $u_{*}$ is related to shear stress through

$${\tau }_b=\rho u_{*}\left|u_{*}\right|$$

(2)

where ρ is fluid density. Due to the unsteady nature of oscillatory flow under waves, the wave boundary layer is typically thin, in the order of a few centimeters, and estimating the friction velocity from a measured velocity profile using (1) and (2) is generally difficult, especially in the field where highly energetic conditions prevent the use of fragile high-resolution instruments. Alternatively, shear stress (τ) can be estimated from a point measurement of the Reynolds stress associated with fluid shear (sometimes termed the velocity covariance method):

$$\tau =-\rho \langle u^{\prime }{w}^{\prime }\rangle$$

(3)

where w is vertical velocity, primes denote turbulent quantities, and the brackets indicate time-averages. To obtain the bed shear stress from Equation (3), the turbulent velocities must be measured at a point within a constant stress layer near the bed. Reynolds stress estimates include only the coherent part of the turbulence and are not always simple to obtain because of potential errors associated with sensor misalignment relative to the axes of the flow. Further, the assumption of constant stress may be difficult to justify in accelerating/decelerating wave-driven flows.

Shear stress can also be estimated from measurements of turbulent kinetic energy (TKE; [15,16]):

$$\tau =C\rho k$$

(4)

where C is a constant in the order of 0.2, and k is turbulent kinetic energy per unit mass, k = 0.5(<u′²> + <v′²> + <w′²>), where v is the transverse velocity component. Kim et al. [13] and Pieterse et al. [14] evaluated Equations (2)–(4) and both found that the TKE method yielded the most consistent results of τ, although more work is required to determine the magnitude of the constant of proportionality, especially under wave-driven flows.

A fourth way of estimating turbulence and (bed) shear stress is through estimation of the turbulent dissipation (ε) using, for example, the velocity frequency spectrum. This approach requires the assumption that turbulent production and dissipation are balanced everywhere, as well as application of Taylor’s hypothesis of frozen turbulence. The latter relates the turbulent length scale to turbulent frequency within the inertial subrange (e.g., [17]):

$$E\left(f\right)=\frac{\alpha }{2{\pi }^{2/3}}{u}^{2/3}{\epsilon }^{2/3}{f}^{5/3}$$

(5)

where α is the Kolmogorov constant (≈0.65), ε is the turbulent dissipation rate (defined as ε = −<u′w′> du/dz; [18]), and E(f) is the spectral density at frequency f. The assumption is made that turbulent vortices are advected by currents (or ‘long’ waves) past the point of measurement, but the criteria for its use, namely that u >> u′ and u >> Δu (where Δu is the variation of u over the turbulent length scale) are often grossly violated in the field [19].

3. Turbulence Measurement

Turbulent velocities under waves can be measured with a range of different instruments. Hot-film sensors, Laser Doppler Anemometers (LDA; [20,21]), and Particle Image Velocimetry (PIV; [22,23,24]) have all found favor in laboratory settings where they can provide detailed and precise measurements. However, these instruments are all quite fragile and measure velocities in only two dimensions. That may not be a critical issue in laboratory wave flumes, but it certainly is in surf zone field experiments where waves are directionally spread and strong currents may be oriented perpendicularly to the incoming waves. In field settings, more robust—but often perhaps less accurate—instruments are used routinely for turbulent measurement such as 3D acoustic Doppler velocimeters (ADV), which may be mounted in vertical arrays (e.g., [25]), acoustic Doppler profilers (ADP; [26,27]), or electromagnetic current meters ([28]).

Both the Reynolds-stress (Equation (3)) and TKE (Equation (4)) methods for estimating stress rely on isolation of turbulent quantities from the total flow vector(s), which may be less than straightforward in wave-driven oscillatory flows. The instantaneous velocity field can be resolved into three components:

$$\begin{gathered} u=\overline{u}+\tilde{u}+u^{\prime } \\ \mathrm{and} \\ w=\overline{w}+\tilde{w}+w^{\prime } \end{gathered}$$

(6)

where overbars indicate mean components, tildes represent the oscillatory flow, and primes are turbulent quantities. While it is trivial to subtract the mean from a time-series, separation of oscillatory and turbulent components is less so, partly because they often overlap in frequency space (e.g., [28,29]).

For regular waves, a common way to remove the wave component is to use ensemble averaging (e.g., [21,30]) by calculating the quantity:

$$\overline{X}(t)=\frac{1}{N}{\displaystyle \sum _{n=0}^{N-1}X(t+nT)}$$

(7)

where X represents the variable in question, in this case u or w, t represents time, T is wave period, and N is the number of ensembles. The turbulence component is calculated as the residual between the velocity time-series and the ensemble-averaged product. However, since waves are rarely regular in nature, other techniques must be used here, for example frequency cut-off methods, involving high-pass filtering to isolate turbulence. Thornton [31] calculated the co-spectrum between u and surface elevation (η) to identify the frequency limit between wave motion and turbulence. This method turned out to be less useful as large turbulent vortices do impose a surface signal [32].

Alternatively, horizontal/vertical velocity spectral cut-offs have been defined on the basis of breaks in spectral slope within the high-frequency part of the spectrum (e.g., [4,33,34]). Surface wave energy should decay as f⁻³ and the inertial subrange is expected to exhibit a spectral decay of f^−5/3 under the assumption that Taylor’s hypothesis of frozen turbulence (Equation (5)) applies.

Christensen et al. [35] introduced a modified frequency cut-off method using the co-spectral phase between u and w to define the frequency that separates wave orbital motion from turbulent motion. For wave orbital motions outside the wave boundary layer, the phase between u and w is expected to be ±π/2 and with statistically significant coherence, whereas for turbulence, the u-w phase is expected to exhibit continuous and erratic phase jumps with low coherence. Figure 1 shows an example of the u-w co-spectral technique, where the transition from dominance of wave orbital motion to turbulence (the cut-off frequency) was estimated at f = 0.62 Hz, based on the u-w phase jumps and the drop to statistically insignificant coherence. It is clear, however, that such imposed cut-off frequencies do not uniquely define, for example, the gravity wave and turbulence frequency bands since gravity wave motion may leak into frequencies above the cut-off, and large-scale turbulent motions may exist at frequencies below the cut-off such that an overlap region must be expected. However, the cut-off should identify the spectral bandwidths where the respective motions (wave orbital vs. turbulent) are dominant [36].

Alternatively, Bricker and Monismith [37] devised a phase-method operating in the spectral domain that produced satisfactory results in cases when the wave-induced strain field was weaker than the turbulence strain field. However, recently a standard method has emerged where two or more velocity sensors are deployed, separated in either the vertical or the horizontal dimension. This velocity-differencing method is based on work by Trowbridge [38], which was later extended by [39,40]. The method relies on the assumption that a turbulence decorrelation length scale exists and that correlations between two sensors separated in space will include wave, but not turbulence components; hence, wave bias in the time series can be removed. As a check on Reynolds-stress estimates derived from this method, Feddersen and Williams [40] introduced the non-dimensional integrated co-spectrum Og_u′w′(f) (the ogive curve) for u′w′, which was defined as

$$O{g}_{u^{\prime }{w}^{\prime }}(f)={\displaystyle \underset{f}{\int }\frac{C{o}_{u^{\prime }{w}^{\prime }}(f)df}{\overline{u^{\prime }{w}^{\prime }}}}$$

(8)

where Co_u′w′ is the co-spectrum of u′w′. When co-spectral wave bias is minimal, ogive curves are expected to increase gently from 0 to 1 over a wide frequency range [41]. On the other hand, if the curves increase sharply over a narrow frequency range, or fluctuate erratically with frequency, the stress estimates are likely wave-bias contaminated. However, a robust test to reject bad Reynolds stress estimates does not yet exist [25].

Feddersen and Williams [40] concluded that the velocity-differencing technique appears to work well in the presence of strong currents. That would imply that it may be most suitable for estimating turbulence in the longshore dimension where the oscillatory velocity component is relatively weak, while, on the other hand, turbulence estimates may contain wave bias in cases when currents are relatively weak. A further question mark relating to this technique is the uncertainty with respect to the optimal sensor separation distance; Brinkkemper et al. [42] found that turbulence magnitudes increased with sensor separation. In addition, it is unclear if, and how, a potential optimum separation distance might change with wave conditions.

Christensen et al. [35] evaluated the u-w co-spectral technique relative to the velocity-differencing method and found that turbulent kinetic energy (k) estimated by the u-w method was about a factor 2 smaller than estimates from velocity-differencing. The simplest explanation would seem to be that the former may tend to exclude large turbulent eddies with low frequencies [34,43], and/or that the latter may contain wave bias related to uncertainties with respect to optimum sensor separation distance. Importantly, however, the k-profiles, and the temporal variation of k within the time series were qualitatively similar using the two methods.

4. Turbulence Generation, Spatial Structure, and Magnitude

When waves are breaking in a surf zone, turbulence may be generated in several ways: through fluid shear, through flow separation around roughness elements, and by injection of turbulent kinetic energy from breaking waves. Over a smooth seabed, turbulence starts to appear in the boundary layer flow when the Reynolds number Re > 1.5 × 10⁵ (Re = Au/ν, where A is orbital amplitude and ν is kinematic viscosity) but this number is expected to be smaller for rough beds [11] (p. 260). Shear between the flow and the seabed creates micro-turbulent vortices such as vortex tubes and turbulent spots (e.g., [44,45]) that originate at the bottom boundary and spread upward through diffusion. Horizontal flow velocity and turbulent production are more or less in-phase and production of k scales with u² [46]. When incoming waves are skewed, maximum k-production therefore occurs under the wave crest [23]. In the field, the situation is rarely that simple because strong mean currents and/or infragravity waves may distort the in-phase relationship between u and k [47]. For example, a strong seaward directed undertow may enhance trough velocities to an extent where maximum k (and thus bed shear stress) occurs under the wave trough.

Bedforms, such as wave/current ripples, or megaripples may appear both outside and inside the surf zone, although ripples begin to decay when the mobility number $\psi =u^2/((s-1)gD$ > ≈200 (e.g., [48,49]), where (s − 1) denotes the relative sediment density, g is the acceleration of gravity, and D is the mean sediment grain size. When ripples are of moderate or large steepness, they often create flow separation and vortex shedding, as schematized in Figure 2. Turbulent vortices are generated on the lee slope of ripples during onshore and offshore wave phases and are ejected at flow reversals when the horizontal velocity is zero. In these coherent vortices, turbulence spreads upward by convection, rather than through diffusion [50]. Experimental evidence [51] and RANS modeling [52] have shown that vortex ejection enhances near-bed k and amplifies TKE to a level of about 4–8 ripple heights, in the model example corresponding to 25 cm above the bed.

Turbulent production caused by flow separation around bedforms scales with u and with bed roughness, with the latter depending on both sediment grain size and bedform steepness [50] and/or bedform height [54]. For skewed, shoaling waves, the production of k is maximum at the onshore-to-offshore flow reversal [55] such that u and k are in quadrature, as illustrated in Figure 2.

These two mechanisms for turbulent kinetic energy production may be relevant anywhere under wave motion, and an important added source appears within the surf zone where surface-generated turbulence due to wave breaking may become the dominant source of TKE, even close to the seabed [30,56]. Wave breaking is characterized by a sudden transition from irrotational to rotational flow, accompanied by a violent transformation of wave energy into turbulence and eventually heat [57,58]. Large-scale coherent vortices create strong vertical mixing, and surface-generated breaker turbulence spreads downward by convection [59]; however, in general, a relatively small part of the wave energy is dissipated below trough level with the majority being dissipated between wave crest and wave trough [60,61]. The relative wave height γ (γ = H/h, where H is wave height and h is water depth) is a useful parameter for characterizing, or scaling, a range of surf zone processes; for example, γ is used to predict the onset of wave breaking and breaking intensity, and wave energy dissipation scales with γ. If the relative wave height is sufficiently large (γ > ≈0.4; [62]), breaker turbulence may invade the wave boundary layer and hit the seabed [63].

Since surface- and bottom-generated turbulence coexist within the surf zone and generate turbulence at different wave phases (as will be discussed in Section 6), it is often useful to distinguish and separate the different sources of turbulence. Ruessink [25] argued that when surface-generated turbulence is dominant, the time-averaged turbulent momentum flux (ρ<u′w′>) should be predominantly negative (downward-directed) since vortices carry high-speed fluid towards the bed [64]. On the other hand, the net momentum flux may be expected to be positive when bottom-generated turbulence is dominant [21].

While wave breaking is a major turbulent source in the surf zone, the process is best described as stochastic. In an irregular wave field, some waves break while some do not. Turbulence production (and dissipation) is therefore highly intermittent, and instantaneous k-levels may be several orders of magnitude larger than wave-phase-averaged <k> [30]. Moreover, different types of wave breaking exist and the scale and intensity of the turbulent vortices depend on breaker type. In the outer surf zone, waves break through plunging or spilling (Figure 3). Both breaker types gradually transition into depth-saturated surf bores and this transition constitutes the limit between the outer and inner surf zones [65]. Sometimes it may be useful to distinguish these different wave types, and Grasso et al. [66] suggested that breakers and surf bores may be separated on account of the ratio between wave skewness and wave asymmetry.

When breakers are spilling, turbulence spreads downward towards the seabed behind the wave crest. The process was probably initially described by Nadaoka et al. [32] who coined the turbulent structures ‘obliquely-descending eddies’ (ODEs) that trail behind the wave crest. They were also observed in laboratory experiments with surf bores [67], and Zhou et al. [68] visualized them numerically using a 3D Large Eddy Simulation Model (LES), see Figure 4. This figure shows that turbulence is created at the front of the wave and it spreads slowly downward towards the bed via descending eddies. The coherent turbulent structures impinge onto the bed as spatially isolated events some distance behind the wave crest.

When waves break through plunging, turbulence penetrates deeper into the bottom region than for spilling breakers, and the turbulent mixing rate below trough level is much larger [57]. Plunging breakers create large vortices, or downbursts, that rotate about a horizontal axis parallel to the wave crest [10,30] and generate oscillations in the vertical velocity [69] which may sometimes produce scour holes on the seabed [70]. In the field, Aagaard and Hughes [28] observed instantaneous velocities up to w ≈ 1 ms⁻¹ at an elevation 15 cm above the bed, and w_rms was about a factor 2 larger for plunging compared to spilling breakers of the same height, while deSerio and Mossa [71] reported that u′ was a factor 3 larger for plunging compared to spilling breakers.

Turbulent kinetic energy may be non-dimensionalized using Froude-scaling, $k^{\prime }=\left(\sqrt{k/gh}\right)$ which makes intuitive sense because the rate of turbulence production should scale with wave celerity [59]. For irregular waves in field and large-scale laboratory settings, <k′> ≈ 0.02–0.06 [42,72,73] which is typically smaller than for regular laboratory waves of similar root-mean-square height [25,56]. The difference arises because each regular wave injects turbulence while only some waves are breaking in an irregular wave field. <k′> has been observed to scale with γ [25,42] which is also consistent with expectations since turbulent production from wave breaking must scale with the intensity of wave breaking, and γ is a proxy for that intensity.

Christensen et al. [74] measured the turbulent kinetic energy at two beaches dominated by plunging breakers (Durras Beach, NSW, Australia, which is exposed to swell waves) and spilling breakers (Vejers Beach, Denmark, exposed to mainly wind waves). <k′> scaled with γ at both beaches (Figure 5) but with a significant degree of scatter, which is partly due to the different measurement elevations in the water column and uncontrollable bed level changes within and between periods of recording.

A different perspective is offered in Figure 6a. Here, data from roughly 60,000 waves were extracted from about 500 irregular-wave time series using a zero-downcrossing algorithm; the crest-to-trough height (H_z), mean water depth (h_z), and k were calculated for each wave and turbulence data were then binned within γ_z-ranges. Wave-averaged <k′> was similar at the two beaches (for similar γ) and scaled with the relative wave height, but maximum k′ during individual wave cycles (k′_max) were larger for plunging waves than for spilling (Figure 6b), which indicates that instantaneous levels of k are larger, but of shorter duration for plunging waves. This is consistent with Hsu et al. [73] who noted that for plunging breakers, instantaneous k(t) is typically many times larger than mean and phase-averaged <k>, while k may be up to a factor 5 larger for plunging breakers compared to shoaling waves [63].

Turbulent kinetic energy is related to eddy viscosity through

$$\sqrt{k}={\nu }_t/l$$

(9)

where ν_t is the eddy viscosity and l is a turbulent length scale. Eddy viscosity is often used in model simulations to describe sediment pick-up and the shape of the suspended sediment concentration profile, and the turbulent length scale expresses a distance over which the turbulent velocity becomes decorrelated; in other words, the length scale of the turbulent vortices. Measurements suggest that l = 0.1–0.3 h under breaking waves in the field [19,22,75,76] but l would be expected to depend on both breaker type as well as relative wave height. For plunging breakers, Aagaard and Jensen [77] calculated significantly larger turbulent length scales, l up to ≈ 0.65 h.

Time-averaged turbulence (<k>) tends to exhibit vertical structure within the water column, and the structure depends on location within (or outside) the surf zone. In general, observations have shown that for non-breaking waves, <k> is maximum at the bed. For breaking waves, Figure 7 illustrates some generally observed tendencies. The data were obtained from field measurements of the turbulent dissipation rate at Truc Vert beach [66] and show that for breaking waves, ε was almost uniformly distributed in the vertical, with a slight increase towards the upper parts of the water column. This is consistent with other field and laboratory experiments that have documented that plunging breakers, in particular, produce a very well-mixed water column with an almost uniform vertical distribution of <k> [25,78,79]. For surf bores in the inner surf zone, the field data in Figure 7 show a large decrease of ε towards the seabed.

To illustrate near-bed tendencies in more detail, Figure 8 plots <k′>-profiles from two field experiments [74]. For strongly breaking wave conditions (γ_s = H_s/h > 0.5), <k′> was maximum at the highest measurement position (z = 50 cm) with a decrease of about 30% close to the bed (z ≈ 5 cm). For nonbreaking wave conditions (γ_s = H_s/h < 0.3), <k′> increased towards the bed due to the roughness created by small (mainly anorbital) wave ripples in fine-grained sand and <k′> was similar at the top and the bottom of the water column. Weakly breaking waves (0.3 < H_s/h < 0.5) exhibited a mixture of the two shapes with <k′> being largest at the highest measurement elevation and likely even larger near the sea surface.

The vertical turbulence structure depends not only on the relative wave height but also on bed roughness, breaker type, and mean current shear. Brinkkemper et al [42] reported results from laboratory tests with regular and irregular wave series on a medium-coarse grained bed (median grain diameter, D₅₀ = 0.42 mm) with large ripples (η = 2–16 cm). Even for strongly breaking waves (γ_s > 0.67), <k> was maximum at the seabed, thus overwhelming surface-generated turbulence. When wave-driven longshore currents are strong, wave/current interaction in the bottom boundary layer increases bed shear stress relative to situations without longshore currents [36] and may conceivably also lead to <k>-maxima near the seabed. Considering effects of wave breaker type, observations have shown that for spilling breakers and surf bores, the vertical mixing is slower compared to plunging breakers and <k> tends to decrease more rapidly towards the bed [42,78,79], see also Figure 7.

In the horizontal dimension, <k> tends to be maximum at the transition from the outer to the inner surf zone [73], i.e., at the transition from breaking waves to surf bores. This was illustrated by measurements with regular plunging breakers in a wave flume over a sandy bed [63]. Vertical profiles of <k′> were measured from trough level to the seabed with very high resolution in the boundary layer, using an Acoustic Concentration and Velocity Profiler (ACVP; Figure 9). The upper panel in the figure shows the vertical <k′> profiles which are largely consistent with the data shown in Figure 7 and Figure 8 except for the secondary maximum observed in some of the profiles very near the bed that could be resolved with the high-resolution sensor. The lower panel shows that near-bed k′_max peaked slightly landward of the main breakpoint (at x = 55–56 m) and exceeded k′_max in the shoaling and inner surf zones by a factor 3–4.

5. Sediment Mobilization and Time-Averaged Sediment Concentration

Turbulent vortices lift sediment particles off the seabed and into suspension. For non-breaking waves on a flat bed, the lift force on the particles is mainly caused by an upward-directed pressure gradient associated with the ejection of low-speed velocity streaks [11] (p. 132). When the lift force exceeds the gravity forces on the sediment particles, the grains are picked up from the bed and particle suspension is an essentially diffusive process. However, over steep vortex ripples, the suspension is fundamentally different from that associated with gradient diffusion since in this case, it is caused mainly by convective vortex shedding from the bedforms [80]. The formation of turbulent vortices over ripples was discussed briefly in Section 4 (Figure 2). Vortices form in the lee of steep ripples at orbital velocity maxima under wave crests and wave troughs, they pick up sediment from the bed and are subsequently ejected from the bed at velocity reversals [55,81,82,83,84].

Wave breaking is a particularly efficient mechanism of sediment suspension since the breaker vortices may penetrate into the wave boundary layer where sediment pick-up and hence sediment transport are significantly increased [69,79,85,86]. Overall, turbulence generated by breaking waves plays a key role for suspension in the surf zone [87]. For example, Aagaard and Hughes [28] observed that the impact of breaking waves on the seabed generated large sediment clouds, and under plunging breakers up to 85% of the suspended sediment load was associated with large breaker vortices. Clouds of sediment are lifted upward in a convective process [88], and the clouds may be advected several meters away from the breakpoint before sediment settles out of the water column [89,90]. Consequently, the suspended sediment concentration (SSC) in the water column may be several times larger under broken compared to nonbreaking waves [69].

Sumer and Fuhrman [11] explained the physical processes associated with sediment suspension under plunging breakers. Laboratory measurements of pore water pressure in a sandy bed [2] demonstrated large lift forces under plunging breakers; the lift forces created suction in the sediment bed and subsequent lifting of sediment particles. The lift was ascribed to a pressure minimum due to the breaker vortex moving in the onshore direction: At the bottom of the vortex, flow velocity is relatively small, while it is large at the top of the vortex. The pressure gradient lifts a sediment plume, such as illustrated in Figure 3 where isolated suspended sediment clouds are visible on the front face of small plunging breakers. Field observations [77,90] have shown that sediment clouds lifted by plunging breaker vortices can extend a significant distance above the bed (up to z > 50 cm) and Figure 10 shows that when plunging is violent, the clouds may occupy the entire water column.

Laboratory experiments in wave flumes have shown that inside the surf zone, mean sediment concentrations indeed scale with k rather than with u (or u²; [91,92]). Since k scales with γ (Figure 6), it might be expected that there is also a relationship between SSC and γ. This relationship is explored in Figure 11 using the same dataset and techniques as in Figure 6 and the figure shows that mean SSC was significantly larger at Durras Beach which was dominated by long-period swell waves and plunging breakers. While <k′> was largely similar for the two beaches (Figure 6a), k′_max was significantly larger at Durras Beach (Figure 6b), which again suggests that sediment suspension events are predominantly related to isolated turbulent events rather than to wave-averaged turbulence levels [93]. Unexpectedly, the effect of (mean) grain size did not appear to be critical for this dataset; at Vejers Beach, D = 0.21 mm while it was almost a factor two larger at Durras (D = 0.36 mm).

When averaged over several wave periods, the vertical distribution of suspended sediment may be described by classic concepts of gradient diffusion in the form of the time-averaged sediment convection-diffusion equation (e.g., [94]):

$$w_{c}C(z)+\overline{{\epsilon }_s\frac{dC}{dz}}=0$$

(10)

where C is the (mean) sediment concentration, z is vertical distance above the bed, and the overbar indicates wave-averaged terms. The equation states that the upward mixing of sediment, expressed through a ‘sediment diffusivity’, ε_s, is balanced by sediment settling, expressed by the sediment fall velocity, w_s. The sediment diffusivity is proportional to the vertical mixing length and a mixing velocity (${\epsilon }_s\approx l_m{w}_m$) and the property is related to eddy viscosity (ν_t; Equation (9)) which is used to model the transfer of momentum due to turbulent eddies [50]. Some doubt exists as to whether ε_s is in fact equal to ν_t; experimentally derived proportionality constants range between 0.43 [95] and 1 [96]. Assuming that the mixing velocity is expressed by $\sqrt{k}$, ε_s may be written [59]:

$${\epsilon }_s\approx l_m\sqrt{k}$$

(11)

where l_m can be estimated from measurements of the SSC-profile [97]. Using this approach, Aagaard and Jensen [77] found that for plunging breakers, ε_s was vertically invariant, while for spilling breakers and surf bores, ε_s was very small at the seabed and increased linearly upward, similar to findings by Ogston and Sternberg [95]. For nonbreaking waves over ripples, ε_s also increases upward as vortices ejected from the ripples expand and eventually dissipate [96].

As an alternative to Equation (10), vertical profiles of SSC can be predicted by modeling sediment pick-up, cw’ [50,72]. In numerical models, pick-up functions are often parameterized using proxies for c and w’, such as u and quadratic stress formulations. However, Amoudry et al. [52] found that these proxies did not provide reasonable results over rippled beds since pick-up maxima are predicted at times of maximum velocity. Instead, an approach involving shear production of k resulted in considerable improvement of their numerical model predictions. It is also doubtful whether the velocity/quadratic stress approach is appropriate for breaking wave conditions, as will be discussed in the next section.

6. Intra-Wave Properties of Turbulence and Sediment Suspension

The time-averaged bed shear stress magnitude—whether it is expressed by k, Reynolds stress, or any other relevant parameter—is crucially important to the time-averaged sediment load and therefore also to both bedload and suspended sediment transport (SST). However, in oscillatory wave-driven flows, the phase relationships between shear stress, sediment load, and u are equally important to the rates of transport. Net SST at a point (i.e., the sediment flux) under waves can be defined as:

$$\langle q_s\rangle =\langle uc\rangle =\langle u\rangle \langle c\rangle +\langle \tilde{u}\tilde{c}\rangle +\langle u^{\prime }{c}^{\prime }\rangle$$

(12)

where brackets denote time-averaged terms and tildes are non-steady (oscillatory) terms. The first term quantifies the transport associated with mean currents, the second term comprises the transport due to waves (both high-frequency sea/swell waves and infragravity waves), and the third term is turbulence transport which is usually negligible. In the following, we will consider only the high-frequency sea/swell wave component, which is critical to onshore transport of suspended sediment and depends on the co-variation (the phase coupling) between wave orbital velocity and the oscillating SSC. Numerical sediment transport models have notorious difficulties with predicting this second term. It is less than trivial to predict the wave phase where maximum shear stress (or k) production occurs—particularly since this appears to be different for different wave and seabed conditions (see earlier sections)—and even more problematic to predict the phase relationships of the sediment response. The latter involve sediment properties such as grain size and settling velocity that may result in significant lags between hydrodynamic forcing and sediment response. This section is concerned with examination of the phase-dependent intra-wave properties of turbulence and suspended sediment load.

In the relatively simple situation of non-breaking waves over flat beds, turbulent production scales with orbital velocity, and k-maxima therefore occur under wave crests and troughs, while wave nonlinearities cause asymmetries in k on wave crest and trough phases. These relationships were demonstrated, for example, in recent laboratory flume experiments [23,46] and they have also been observed in the field over low anorbital ripples [86]. Due to lags in response and lift of sediment grains, SSC is expected to lag u and k, and the lag is expected to increase vertically upward in the water column. The lag also depends on sediment grain size and wave period. Finer grains may respond more readily to forcing than coarser grains but will take longer to settle out of the water column, such that lags on increasing and decreasing SSC-phases may be different. On the other hand, long wave periods facilitate sediment mobilization and settling within the same half-wave period.

When steep wave ripples form on the seabed, bed-generated turbulence and vortex formation affect the phase relationships between SSC and wave orbital motion. For both regular [52,55] and irregular [84] laboratory waves, vortex ripples tend to produce two suspension events within the wave cycle: one at each flow reversal (Figure 2 and Figure 12). Wave skewness affects vortex strength such that positively skewed waves produce larger suspension clouds at velocity zero-downcrossing compared to clouds ejected after the trough phase (Figure 12). The regular shedding process at each flow reversal has, however, been difficult to confirm in the field, where it is in fact only rarely observed; the number of SSC-peaks within a wave cycle varies in response to both vortex strength and wave nonlinearity (e.g., [82,98]). This was exemplified by field observations under shoaling waves and low anorbital ripples [86], Figure 13. Only one SSC-peak was observed, and at z = 4 cm, the timing of the peak lagged the zero-downcrossing by 50–60° such that maximum SSC occurred on the following wave trough. For this (ensemble-averaged) irregular wave data set, there was no strong evidence for a secondary peak following the velocity zero-upcrossing.

O’Hara Murray et al. [84] suggested that vortex shedding occurs only when ripples are sufficiently steep (η/λ > 0.1), and when the ratio between orbital diameter and bedform wavelength, d_o/λ > 1.2, i.e., when the wave orbits cannot be accommodated within one ripple wavelength. A further complication in predictions of suspension events over ripples is that strong vortices favor longer life spans of the turbulent sediment-laden eddies, such that advection of suspension clouds from neighboring bed forms may occur. In this case, several suspension events may be observed within one wave cycle [98]. On the other hand, strongly skewed waves favor a single large suspension event at velocity down-crossing when orbital velocities are significantly larger beneath the wave crest compared to the wave trough [51,82]; Figure 13. All in all, prediction of turbulent production, sediment suspension and net (suspended) sediment transport rates, and directions over rippled beds is complicated due to the strong interactions between hydrodynamic processes and the still rather unpredictable geometry of the seabed.

Under breaking waves, the timing of turbulence production and the travel time from the surface to the seabed become additional critical factors in the timing of k-arrival at the seabed, SSC, and consequently SST. Recent evidence points to the fact that turbulence production and propagation associated with plunging and spilling breakers, and surf bores have different characters (Section 4) and therefore affect sediment mobilization and transport differently. Van der Zanden et al. [99] presented an instructive example of the effect of different wave types on intra-wave near-bed k using LDA-measurements on a fixed bed profile without bedforms and subjected to regular waves (Figure 14). The examples show wave-cycle time series for shoaling waves, plunging breakers, and surf bores. In the shoaling wave zone, k-production originated from the planar bed on the front face of the wave and propagated upward with increasing time lag. The origin of turbulence is less clear in the plunging wave example due to the limited vertical extent of the measurement section, but k is clearly much larger and was again produced on the front face of the wave. For the sawtooth bores in the inner surf zone (x = 60.2 m), k was not well resolved over wave phase with small k-maxima on both front face and back slope of the wave but of a magnitude much smaller than for the breaking wave.

A 45-s time series of k and SSC at the plunge point of three breakers over a mobile sandy bed is shown in Figure 15. At the zero-upcrossing of each wave, w became upward-directed, followed by a succession of down/upward directed vertical velocities with a velocity amplitude of w = 0.2–0.5 m/s, associated with the rotating breaker vortex. Closer inspection shows that a phase-lead occurred between k and u; k ramped up prior to velocity maxima and was followed almost instantly by sharply defined puffs of suspended sediment. The general signatures shown in Figure 15 echo several reports [71,73,100] that point to the fact that for plunging waves, k-arrival at the seabed typically occurs under the front face, or at the crest of the waves (the relative timing probably depending on water depth (h) and the turbulent length scale (l)). For spilling breakers, however, TKE moves more slowly downward and typically arrives either on the back of the wave crest [86] (see also Figure 4), or as late as the on/offshore velocity reversal [100].

Figure 16 shows intra-wave properties of SSC using ensemble-averaged data from field measurements with both plunging and spilling breakers, for similar wave heights/water depths and mean sediment grain size (D = 0.22 mm; [72]). In these examples, ensemble-averaged k arrived at the seabed at a relative wave phase t/T ≈ 0.2 earlier for plungers than for spillers; for the former, k peaked near the zero-upcrossing, while the maximum was delayed until the wave crest phase for the latter. Sediment suspension (at z = 0.05 m) slightly lagged maximum k and was again delayed by about t/T = 0.2 for spilling waves. Based on these examples, an in-phase relationship between u and c near the bed could be a reasonable assumption for spilling waves, while for plunging waves, c clearly leads u. Closer examination of Figure 16 also reveals a different crest/trough asymmetry of both k and SSC for plunging and spilling waves; k tends to be more rapidly dissipated under plungers than under spillers. For the latter, turbulence may be carried over from one wave cycle to the next [100]. Consequently, SSC-peaks were more well-defined under plunging waves. These observations carry implications for SST under different wave types, and this will be discussed in Section 7.

Extracting the ensemble-averaged data at each relative wave-phase bin for the six time series used in plotting Figure 16, we may examine the functional relationship between k′(t) and SSC(t) on instantaneous time scales and at different elevations above the bed (Figure 17). The plots suggest that k (or k′) may be used to predict SSC even at such instantaneous time scale; there was a statistically significant (α = 0.05) linear relationship of SSC on k′ at all elevations, but the slope of the functional relationship decreased with elevation, as would have been expected given that concentrations decrease with elevation above the bed.

Field measurements of k and SSC on intra-wave scale under asymmetric surf bores in the inner surf zone are rare. Christensen et al. [86] observed that k-peaks occurred near the velocity zero-upcrossing, similar to plunging breakers (Figure 18), but with poorer peak definition (compare with Figure 16) and slower mixing (and/or smaller turbulent length scales), which resulted in significant lags between u and SSC. The signatures are consistent with the detailed laboratory measurements of k shown in Figure 14, which also indicated that there may be little temporal structure in the turbulent velocity signal from surf bores.

Intra-wave properties of k and SSC in the inner surf zone may be influenced by strong mean cross-and longshore currents. Wave-current interaction becomes important since strong currents that oppose the wave motion affect the wave boundary layer and increase near-bed turbulence [50], but will also distort the phase relationships between u and SSC. For example, Ruessink et al. [101] observed an increasing boundary layer thickness for situations with an opposing mean current, however, the timing of k-maxima depended on whether currents were strong or not. For purely oscillatory-flow experiments, k_max emerged under the wave crest, while in cases with a counter-current, k_max shifted backwards to the trough phase, because the offshore-directed current added to the wave-generated (bed) shear stress on the trough phase and decreased the shear stress on the crest phase. The implication is that current shear may affect the phasing of sediment suspension and further delay SSC relative to u such that maximum sediment concentrations sometimes may appear at the trough phase under surf bores [74].

The general picture that emerges from examination of available field and laboratory data on the intra-wave characteristics of SSC is that SSC-peaks are shifted forward on the wave phase with increasing turbulence levels [86], corresponding to increasing relative wave height (Figure 5). For shoaling waves over mobile (rippled) sand beds, c lags u; the ensemble-averaged examples show that c_max(t/T) occurs near the velocity zero-downcrossing (Figure 12) or even on the trough phase (Figure 13) and sometimes (but not always) with a smaller secondary peak near, or after the ensuing upcrossing, depending on wave skewness. For spilling breakers, c_max(t/T) shifts forward and is typically located on the crest/back slope of the wave, while concentration tends to lead velocity for plunging breakers; c_max(t/T) is often located on the crest/front face of the wave, as exemplified in Figure 16. For surf bores, dissipation rates are often large and even though k-maxima may occur on the front face, the slower mixing and the adverse effects of mean current shear in the inner surf zone may shift SSC-peaks backward on the wave phase. These differences carry ramifications for wave-driven transport of suspended sediment, as will be discussed next.

7. Wave-Driven Suspended Sediment Transport

Turbulence affects the rate of wave-driven suspended sediment transport because of the way it affects the magnitude of the suspended sediment load and especially the phase relationship between u and SSC. The time-averaged net flux of SST through a two-dimensional vertical section of the water column under waves and currents is defined as:

$$\langle Q_s\rangle =\frac{1}{T}{\displaystyle \underset{0}{\overset{T}{\int }}{\displaystyle \underset{0}{\overset{h}{\int }}u(z,t)c(z,t)dzdt}}$$

(13)

and <Q_s> can be further separated into mean and oscillatory terms associated with mean current and wave motions, respectively (Equation (12)). Prior to a discussion of (high-frequency sea/swell) wave-driven transport, we note that these formulae rest on the assumption that the advection speed of the sediment grains is identical to that of the surrounding fluid, u_grain = u_fluid. While that is likely to be true for very fine sediment carried as wash-load, it is more doubtful whether it also holds for sand-sized sediment. The problem was examined by Frank-Gilchrist et al. [24], who used PIV and particle tracking velocimetry (PTV) techniques to track both coarse sand grains (0.7 mm diameter) and fluid tracers. They concluded that the measured advected sediment velocities—even for sediment this coarse—were of similar magnitude and phase to the fluid velocity.

As discussed earlier, intra-wave characteristics of SSC and its magnitude is a first-order control on wave-driven SST. SSC magnitude and the phase relationship with u depend significantly on wave (breaker) type but also on sediment grain size and wave period. When suspended sediment is mobilized on the wave crest phase, coarser grain sizes and longer wave periods favour larger onshore transport rates since sediment grains may be able to settle prior to the wave trough phase. On the other hand, reduction and/or reversal of wave-driven transport may occur when strong opposing mean currents and/or smaller mixing lengths in the inner surf zone tend to shift k_max(t/T) and c_max(t/T) backwards on the wave phase. Onshore wave-driven SST is therefore often smaller under (asymmetric) surf bores in the inner surf zone than under (plunging) breakers [93]. However, despite recent progress on the qualitative understanding of these complexities, accurate scaling/parameterization of SST in the surf zone still seems some way off, partly because of measurement difficulties. Field data, in particular, notoriously contain measurement inaccuracies, for example, with respect to the elevation above the bed where variables are measured, and the risks of bubble contamination of signals from optical and acoustic sensors. Moreover, robust prediction of wave (breaker) type is still not achievable, although recent efforts [102] represent a step forward in this respect.

Some general characteristics have emerged on suspended sediment transport driven by incoming sea/swell waves. Brinkkemper et al. [93] measured turbulence and suspended sediment transport rates in a large-scale laboratory experiment with irregular waves, over a mobile coarse-sediment bed. They reported a clear dependency of suspended sediment transport rate and direction on the relative wave height (Figure 19). Even though there is significant scatter in the diagram, <Q_s> increased rapidly when H_s/h increased above 0.65. The diagram also illustrates the fraction of breaking waves (Q_b) which increased from the outer to the inner surf zone and when Q_b > ≈ 0.6–0.7, the high-frequency wave transport decreased dramatically.

Similar results were obtained in field experiments reported by Christensen et al. [74]. In an attempt to obtain common scaling of sediment transport rates at two beaches where suspended sediment loads were very different (Figure 11), near-bed wave-driven sediment flux (q_s; Equation (12)) was normalized by the mean sediment concentration, and the ratio was termed a “flux efficiency”, F:

$$F=\frac{\langle q_s\rangle }{\langle c\rangle }$$

(14)

The dimension of F is ms⁻¹; it is a measure of transport velocity and expresses the efficiency by which the waves move a given sediment load onshore (or offshore). F was calculated for the ≈60,000 waves in the data set and averaged within γ-bins (Figure 20). Data from the two beaches follow the same trend; F was negative, signifying offshore-directed wave-driven transport under shoaling (nonbreaking) waves and became positive at the onset of wave breaking at γ ≈ 0.3. In agreement with Figure 19, the dataset with mainly plunging breakers exhibited a break-off point near γ = 0.7 where F decreased due to the strong undertow and infragravity wave activity [74] under strongly asymmetric surf bores in the inner surf zone. The dataset for mainly spilling breakers did not extend to sufficiently large γ to ascertain whether that was also the case here. Hence, given robust parameterization of <c>, it may eventually become possible to predict <qs> for different wave conditions on different beaches with different grain sizes, although much more work is required to test and refine the relationships in Figure 19 and Figure 20.

8. Summary and Future Perspectives

Beaches exist because waves restore sediment that has been lost offshore to the shoreface during storms. Wave-driven transport is of first-order importance for an onshore transport of sediment and thus beach recovery during post-storm conditions when offshore-directed mean currents weaken.

The prevailing opinion in coastal geomorphology is that beach recovery takes place during calm conditions when shoaling, non-breaking waves push sand onshore. Nevertheless, field and laboratory data overwhelmingly show that wave-driven transport in suspension is actually often directed offshore during such conditions (e.g., Figure 20), even when ripples are low and of the anorbital type (Figure 14). If significant onshore sediment transfer does occur under non-breaking calm weather conditions, it is more likely through bedload transport [103].

For sediment in suspension, onshore transport is clearly more efficient for breaking waves in the outer surf zone, particularly when waves are plunging and of long periods (corresponding to a swell wave climate). With increasing relative wave height, SSC shifts forward on the wave phase [86,93] and for plunging breakers, there is a phase lead of k and SSC relative to u; they both occur on the front face of the waves such that sediment is available for onshore transport under the entire wave crest phase.

The evidence reviewed here points to the fact that cross-shore velocity, u, or higher-order moments thereof, is a poor predictor of both suspended sediment load and wave-driven suspended sediment transport within the surf zone. Model calibration might be able to accommodate and predict the former, but not the latter, since transport depends on the phase relationship and intra-wave properties of u and SSC. Time-dependent SSC critically depends on the timing of TKE-production—which clearly depends strongly upon wave (breaker) type, and SSC is further modulated by seabed configuration and co-existing mean currents.

A robust model for wave-driven suspended sediment transport in the surf zone therefore does not yet exist. A parameterized model might be based on the observed proportionalities between normalized wave-driven SST and relative wave height (Figure 19 and Figure 20). However, these relationships clearly need further examination and testing using field (and laboratory) data from different types of beaches and with a wider range of sediment grain sizes. The critical quantity that needs to be predicted in any transport model is the SSC-magnitude, and while γ does appear to also hold some predictive capability in this respect (since it is a measure of wave energy dissipation and hence TKE-production), the scaling is different for different beaches (Figure 11). Apart from potential grain size effects, the reason is that SSC depends on the intra-wave properties of k (Figure 17) and scales with k_max rather than with $\overline{k}$.

Alternatively, more detailed process-based models of surf zone sediment transport would require resolution and prediction of k over wave phase. Encouraging progress has been made recently on simulation of time-varying turbulent properties using RANS-models [104,105] or Large Eddy Simulation [87,106] but computational demands are significant. Models for application to the surf zone would need to consider the type of wave breaking and the k/SSC lags that depend on wave (breaker) type, grain size as well as further lags associated with nonlinear interactions between waves and currents. RANS-models could be used as ‘numerical laboratories’ [105] and provided that they are sufficiently accurate, the output could be parameterized for inclusion in less complex operational models.