Observations of Nearbed Turbulence over Mobile Bedforms in Combined, Collinear Wave-Current Flows

Collinear wave-current shear interactions are often assumed to be the same for currents following or opposing the direction of regular wave propagation; with momentum and mass exchanges restricted to the thin oscillating boundary layer (zero-flux condition) and enhanced but equal wave-averaged bed shear stresses. To examine these assumptions, a prototype-scale experiment investigated the nature of turbulent exchanges in flows with currents aligned to, and opposing, wave propagation over a mobile sandy bed. Estimated mean and maximum stresses from measurements above the bed exceeded predictions by models of bed shear stress subscribing to the assumptions above, suggesting the combined boundary layer is larger than predicted by theory. The core flow experiences upward turbulent fluxes in aligned flows, coupled with sediment entrainment by vortex shedding at flow reversal, whilst downward fluxes of eddies generated by the core flow, and strong adverse shear can enhance near-bed mass transport, in opposing currents. Current-aligned coherent structures contribute significantly to the stress and energy dissipation, and display characteristics of wall-attached eddies formed by the pairing of counter-rotating vortices. These preliminary findings suggest a notable difference in wave-following and wave-opposing wave-current interactions, and highlight the need to account for intermittent momentum-exchanges in predicting stress, boundary layer thickness and sediment transport.


Introduction
Fluid motion in the coastal zone is often characterised by the simultaneous presence of both oscillatory (e.g., surface wind and swell waves) and unidirectional currents (e.g., tidal, fluvial, winddriven, rip currents, etc.). The interaction between currents and waves is extremely complicated [1], and has direct implications on sediment transport and bedform dynamics, as well as engineering problems such as the design of offshore structures and beach replenishment. Sediment entrainment and transport are often described as a function of bed shear stress, a measure of the fluid force per unit area. In oscillatory flows, bed shear stresses may be an order of magnitude larger than encountered in a current of comparable magnitude, due to sharper vertical velocity gradients associated with much smaller wave boundary layer thicknesses [2]. Non-linear interactions between the steady shear of the current and the oscillating flow often result in enhanced stresses beyond the simple summation of the current-only and wave-only components [1,3]. Thus, a current below the threshold for motion may move significant amounts of sediments in the presence of waves [4]. As a result, combined wave-current flows often result in highly dynamic bedforms and active sediment transport across the continental shelf and shallow coastal waters [5,6].
To model the transport of sediment under the combined action of waves and currents, the maximum bed shear stress is often used to determine the threshold for motion and entrainment of particles, while the flow velocity and suspension of sediment into the upper part of the flow are determined by the time-averaged, combined stresses within a wave cycle [3,7]. The oscillating component of stress is commonly presumed to be confined within a thin wave boundary layer in the direct vicinity of the bed, within which intense turbulence is generated (Figure 1). Its effect on the current profile is felt as enhanced bottom roughness [8,9]. Kemp and Simons observed, experimentally, an enhancement of turbulent intensity within the inner bed layer (which they defined to be within 1-2 roughness lengths) when currents are linearly superposed on waves for rough, rippled beds [10,11]. They showed that maximum turbulence intensities and Reynolds stresses depend upon the magnitude of current and wave-induced velocities rather than the relative direction; whereas mean stresses are enhanced by changes in vortex patterns formed between the bedforms. This implied an increase in sediment pick-up heights and diffusion beyond the combined layer. They concluded, therefore, that the relative direction only influences the outer layer of flow, but not the wave-current interaction near the bed. Following similar arguments, most numerical models used to predict bed shear stress in combined flow employ a linear, vector addition of the individual, enhanced steady (current) and unsteady (wave) components within the near bed region, accounting for the relative angle, , between the two [4,8,[12][13][14][15][16][17][18]. This means that the magnitudes of mean (and of maximum) bed shear stresses are equal for collinear flows, since | cos | = 1 for = 0 or 180°; irrespective of whether the current follows, or opposes, the direction of wave propagation. Strong non-linear interactions between the wave and current are at their peak in collinear flows in such models [19]. Nonetheless, comparisons between the various models used often reveal significant differences between estimated and measured stresses [1,20,21] and thus the estimate of sediment transport can differ significantly. Figure 1. An illustration of boundary layer structure under regular waves ( ) and steady currents ( ), and the associated shear stresses. The vectorial addition of mean ( ) and maximum ( , ) stresses in combined oscillatory (Uw) and steady (Uz) flows is shown, as is the zero-flux boundary condition, whereby : kinetic energy, : energy dissipation, : eddy viscosity, : sediment settling velocity, and : sediment concentration.
In addition, such models also impose a zero-flux condition limiting turbulent quantities; and by consequence sediment concentrations, at the interface of the wave boundary layer with the core flow [22], despite documented evidence of vertical momentum exchanges extending beyond the oscillating layer in the nearshore [23][24][25]. Holmedal et al. contend that such simplification is problematic given production and dissipation of turbulent kinetic energy are taking place within the current boundary layer at the interface with the wave boundary layer [22]. However, they argue that this simplification can still be used as the vertical gradients in both boundary layers are confined near the bed and are very small, exhibiting slow variations with time and decaying rapidly with distance from the bed. Under nonbreaking, regular waves following and opposing a current, combined Reynolds stresses (momentum exchange) decrease slowly with increasing height above the bed in the opposing flow case; while a more rapid decrease and change in sign is observed at mid-depth for the aligned flow case [26]. Such differences between the aligned and opposing cases have been attributed to turbulent mixing by secondary local vortices [27][28][29]. Lumley and Terray linked the spatial structure of such turbulent motions to temporal fluctuations due to advection by the wave motion [30]. Shaw and Trowbridge suggested that wave and turbulence induced fluctuations are incoherent with one another outside the thin wave boundary layer in the nearshore to obtain direct estimates of near-bed turbulent fluxes, which they found to be complicated by the intermittent nature of turbulent motions near the bed [25]. Trowbridge and Elgar measured imbalances in the momentum equation and deviations from the Prandtl-von Kármán logarithmic velocity law near the bed (below the reach of breaking induced turbulence) associated with inaccurate estimates of Reynolds stresses [31]. Such deviations have also been attributed to stratification by suspended sediments [32]. It follows that determining the vertical spatial scales of near-bed momentum-bearing motions is an important research question for modelling currents and sediment transport in the nearshore and on continental shelves.
Little is known about the spatial and temporal scales of momentum fluxes by coherent turbulence structures in combined flows, or about their interaction with bedforms and suspended sediment under such conditions [33]. Accurate estimates of combined wave-current bed shear stress and ensuing sediment transport require robust characterisation of the momentum exchanges by these motions in the near bed region in space and time. This paper examines the two fundamental assumptions used in describing/modelling fluid shear due to wave-current interactions. We assess whether: (i) the mean and maximum combined bed shear stresses-and the nature of intermittent momentum fluxes-are the same irrespective of whether a collinear current flows along or opposing the direction of regular wave propagation; and (ii) the highly-nonlinear interaction is restricted to the thin oscillating boundary layer thus limiting vertical exchanges of mass and moment to the direct vicinity of the bed (zero-flux boundary condition). Prototype-scale flume measurements of near bed flow, turbulence and bed morphology are presented under the combined effects of waves and currents. Two cases are considered: collinear currents aligned with; and opposed to, the direction of wave propagation. We highlight differences in the structural nature of turbulent exchanges and stress-bearing coherent structures as opposed to the deterministic view based on time-averaged quantities used to predict sediment transport.

Materials and Methods
Results are presented from an experiment conducted at the HR Wallingford Fast Flow Facility (FFF), UK, a large, dual-channel (racetrack-shaped) flume where combined wave-current-sediment interactions can be studied ( Figure 2). A flat, sandy bed section, 0.3 m thick, was constructed in the main (4 m wide) channel, starting at a distance of 27 m from the wave paddles and spanning a length of ~21 m, with the remaining bed lined with smooth concrete. A 10-paddle, HR Wallingford bottomhinged wave maker was used to generate waves. Currents were generated by two Bedford Pumps Ltd. SA-80.04.08 pumps, operating alternatively in forward (producing currents aligned with direction of wave propagation, counter-clockwise in flume) and reverse (opposing progressive waves, clockwise) modes. The sediment bed was composed of well-sorted, fine to medium fluvial silica sand, with a normal (Gaussian) grain size distribution and a median grain diameter, d50, of 0.247 ± 0.008 mm. A deeper section (the scour pit, 1.0 m thick, 4 × 4 m wide, Figure 2a) within the bed allowed the stable installation of an instrumented frame, equipped to measure flow hydrodynamics and bed morphology in combined wave-current flows at high frequency.
The three-dimensional velocity field was recorded by two frame-mounted, downward-looking, 10-MHz Nortek acoustic Doppler velocimeters (ADV) offset vertically by 10 cm, and horizontally by 8 cm. Each ADV samples a cylindrical volume (7 mm in length), situated 50 mm below the probe, at a rate of 25 Hz. The ADVs logged their sensor heights above the initial bed elevation at the start of each run, and the position of the sampling volume with respect to the dynamic bed morphology was assessed from the centreline of the sonar scans. The local bed morphology was surveyed by a Marine Electronics sector scanning sonar (SSS, 500 kHz) and a single-beam sediment image profiling sonar (SIS, 1.1 MHz), affording a three-dimensional view of the bed. The incident wave field was measured using twin-wire resistance wave probes recording immersed depth through conductivity, with channels located just before, parallel to, and just after the instrumented frame location within the sandy bed section of the flume.

The Hydrodynamic Forcing
As most models of wave-current interaction require a value for a depth-averaged current velocity, these were derived from pump calibration curves based on discharge, pump frequency and water depth. The non-directional wave parameters (peak wave periods, mean, significant, and maximum wave heights) were computed using zero-crossing and spectral analysis of the surface water elevation following Tucker and Pitt [34]. The wave heights were given as immersed water level variations around the still water level measured at the start of each test, derived from calibrated voltage across each of the twin-wire probes (channels). Two probe channels were available for A2 (Table 1), while three channels were available for the other runs, located just before, at, and just after the instrumented frame. It is suspected, however, that the probes in Runs W2, A1 and A2 were not fully submerged at the troughs of some of the waves, underestimating the wave heights. To this end, a correction was applied to compensate for the apparent truncation in the measurements in these runs, whereby an approximation for water surface elevations was inferred from the measured streamwise velocity, by assuming linear wave theory and therefore re-arranging the wave dispersion equation, whereby the wave celerity is computed (to within 0.1%) following Hunt's [35] approximation. This approach was applied to cases where no truncation is observed, yielding comparable results to within 4%.
The instantaneous, three-dimensional velocity field ( , , ) was measured by the ADVs at two discrete points above the bed. This comprises time series of the streamwise ( , defined as the mean direction of wave propagation along the flume); the cross-wise ( , perpendicular horizontally to , transverse across the flume); and the vertical ( , perpendicular to the ( ) plane, positive upwards) velocity components of the flow. Turbulence is extracted from these time series by Reynolds decomposition, following the approach outlined in Kassem et al. [24]. This encompasses a quality check such that no more than 20% of the record in a given ensemble falls below the correlation threshold, set at 70% following Elgar et al. [36] to account for noise arising from signal aliasing [37][38][39]. Values falling below the threshold were replaced by interpolation through a moving average algorithm [40]. Subsequently, the resulting velocity vectors were geometrically rotated to account for sensor misalignment [41], as even a small disparity may subject estimates of stress inferred from velocity fluctuations to significant errors [42,43]. The resulting signals for each run thus comply with the criteria set by Feddersen [44] to estimate turbulence properties from ADV measurements. The records were then zero-meaned and de-trended; with the mean velocity calculated from the three components. The mean direction was verified by calculating the relative angle between the two horizontal components to ensure the flow is two-dimensional. The wave signal, defined by an adaptive moving average, zero-phase low-pass digital filter, was subtracted from each velocity vector to obtain the fluctuations, , , ′ [40]. These filtered signals were de-spiked using the 3D phase-space method [37,38], resulting in time series of the three turbulent components (velocity fluctuations ', ', and ′). Less than 5.6% of any given data set was affected by the despiking process, and despiking was realised in less than 8 iterations in all runs.

Bed Morphology
The bedform geometry was inferred from the acoustic images of the sonars and underwater photography. The SSS backscatter intensity was used in raw acoustic form (in dB) to assess ripple crest shape and orientation, providing a top view of variations in planform. Mean bedform wavelengths, following slant range correction, and rotation of the image to align with the flume walls, were evaluated from the locations of the peaks in backscatter along the central line of the echograms at the start and end of each wave run. The SIS backscatter intensity was range-corrected along each beam of a given pencil-sweep (excluding the near-field), then normalised against the maximum intensity of that beam, with the background noise, defined as the mean backscatter within the water column, removed. The bed was detected as the range-corrected height at which the maximum normalised backscatter occurs, and bed relief thus inferred after de-trending the bed, and despiking to remove outliers resulting from high concentration suspension events. A zero-crossing algorithm was then used to determine ripple heights [45,46].

Near-Bed, Mean and Maximum Shear Stress
The resulting turbulence components from the ADV processing are used to estimate the turbulent properties of the flow to infer the local shear stresses. The turbulent kinetic energy, ( / ) per unit mass (for a fluid of density, ), is defined as: In dynamic terms, the turbulent kinetic energy (TKE) is thus defined as, = ; and the Reynolds stress ( ) in the principal plane of motion ( ) is defined by the covariance of the velocity fluctuations [47]: The shear stresses can be defined by the measured turbulent kinetic energy at each ADV [40,48]: The mean shear stress at the height of the ADV's sampling volume for each test run is then defined as the time-average, τ = τ ; and τ as the maximum stress within the time series. The TKE method is used to infer stresses as it provides a more robust measure compared to the Reynolds stress approach, which is particularly sensitive to sensor misalignment, and often underestimates the stresses [40,42,43]. The mean thickness of the combined, wave-current boundary layer, δ , , is estimated as δ , = 2κ u * /ω where κ is von Kármán's constant, and ω is the orbital frequency; with u * = τ /ρ defining the combined wave-current shear (friction) velocity [49]. The latter is derived from the bed shear stress using three semi-analytical/empirical models: modified Grant and Madsen (1986) [49,50] [20] modification of SC05 (hereafter MD12). MD12 present a modification of SC05, accounting for non-linearity in the maximum stress, whereby τ = τ + τ + 2τ τ |cos ϕ|, which assumes the strongest non-linear interaction is at a relative angle of ϕ = 0 , or ϕ = 180 . For the GM86 model, a modification to account for transport and ripple enhanced shear velocity is used [50], as implemented in the numerical model Sedtrans05 [51] whereby u * = u * /(1 − πη /λ ), with η and λ representing height and wavelength of ripples.

Turbulent Fluxes and Vertical Eddy Scales
The Reynolds stress describes the turbulent vertical advection of streamwise turbulent moment (i.e., the momentum exchange in the main flow plane). The intermittency of these exchanges is quantified using 'quadrant analysis', which distributes the momentum transport within a single plane of motion between four types of event structures in four quadrants: violent "Sweeps (S or Q4)" and "Ejections (E or Q2)", and weaker "Inward (II or Q3)" and "Outward Interactions (OI or Q1)", collectively termed as "bursting" [52,53]. This does not account for the entire/mean turbulent stress as crosswise fluctuations may contribute to the mean streamwise momentum (ibid.). The fractional contribution to stress of each of the four structural features of flow commonly focusses on the intense events restricted to values that lie above a critical threshold, H, which prescribes the hole between a set of hyperbolae defined by: The contributions of individual u'w' events in the (xz) plane are the ones which occur in each quadrant outside the central "hole" region bounded by these hyperbolae. The choice of H is often arbitrary, and ideally chosen so that it yields the same number of ejections as identified by visualization studies [53].
The number of events (concentration, C ) within each quadrant (i) in the primary plane (xz) is calculated as follows [54,55]: ϕ , = 1 if |u w | > H. u . w and belongs to n quadrant 0 otherwise (6) Meanwhile, the spatial scales of the vertical turbulent eddies can be derived from the spectrum of measured vertical fluctuations at a fixed point by assuming they advect faster than they develop in time, in line with Taylor's "frozen turbulence" theory [56,57]. Accordingly, the rate of strain of an eddy is a function of its wavenumber (k = 2πf/U), and the spectral energy of all eddies of a size 2π/k is proportional to the energy (E(k)) multiplied by the spectrum's width [58]. The normalised, premultiplied turbulence power spectra can thus be presented in non-dimensional wave-number space (k * = kz, where z is the height of measurement) in energy-preserving form such that an equal area under the curve represents equal spectral energy [59]. As the strain rate of a particular eddy is a function of its wavenumber, the peaks of these plots correspond to the dominant spatial scales of motions as they are advected past the fixed ADV sensor [60], and thus these peaks are used to infer the spatial scales of turbulent motion. The length scales, λ, are thereby inferred using λ = 2πz/k * with k * taken at the peak of the spectrum. The mean vertical turbulence flux is defined as k w (m /s ), and the vertical gradient can therefore be inferred given the vertical and horizontal offsets between the two ADVs were close enough to allow estimates of structure functions (statistical moments) of the flow field following Kolmogorov (1941) theory [23,61].

Temporal Scales of Turbulent Motions and Active Periods of Flow
To detect coherent turbulent motions in time, we apply a continuous wavelet transform (CWT), to analyse the local variations in power in time-frequency space [62][63][64]. We utilise the common Morlet wavelet function [65] as a "mother" wavelet (wavelet vectorψ = e . e | | ⁄ ; k = 6) which is then stretched, compressed and shifted to derive "daughter" wavelets that are then used to decompose the signal [66,67]. The continuous and complex nature of CWT allows detection of both amplitude and phase for different frequencies in the time series [68]. This continuous transform maps the time series of the instantaneous Reynolds stress into a 2D portrait of the evolving scales and frequencies in time. The transform can be sped up and the edge effects near large scales removed by padding the time series with zeroes at the start, and removing these afterwards, denoted as a cone of influence [63,64]. The limitation given by the cone of influence is an expression of Heisenberg's uncertainty principle, which implies that one cannot obtain good localisation in both time and frequency, and a trade-off thus exists with bad temporal resolution at the larger scales, and bad frequency resolution at the smaller scales [69,70]. The temporal scales of turbulence can therefore be inferred from the active periods of flow, derived from the high variance regions within the timefrequency domain. The global spectral power is also calculated as the time average of all local wavelet spectra, and is tested in a time-averaged sense using autocorrelation of background noise to determine the 95% confidence limits. This is a reasonable test of the null-hypothesis at the longer scales, whereby the distribution tends towards Gaussian as we convolute longer wavelets according to the central limit theorem [63].

Flow Properties and the Incident Wave Field
The hydrodynamic forcing, non-directional spectral properties of the incident wave field, namely the peak wave period, T , and the mean (H ); significant (H ); and maximum (H ) wave heights inferred from the immersed water depths, and the average flow velocities measured, at each ADV elevation, are presented in Table 1. Runs W1 and W2, represent wave only conditions, A1 -A3 indicate aligned conditions, and O1-O4 denote opposing conditions; with variations in current magnitude, and wave period. Runs W2, A3, and O4 were in 'deeper' water (1.55 m) with higher "design" waves of 0.3 m, whereas all other runs were forced by 0.2 m design waves in 1.05 m of water. The estimated wave field varied in height from these input (design) values in all the runs, due to paddle and channel geometry. The wave periods were consistently reproduced, with the surface wave spectra dominated by the driving wave frequency, and a secondary peak at half cycle, with small, random deviations at higher frequencies from a regular wave field. The mean wave heights, H , were largely consistent across all runs, while the significant wave heights, H , derived from the zero-order moment (integral of variance), m , of the wave spectrum; showed reduction in height in the aligned cases (up to ~9% in A3) by comparison to corresponding wave only flow conditions. Meanwhile, significant heights in the opposing conditions were 10-30 % higher than the wave-only conditions, particularly when stronger opposing currents were present. A similar, more pronounced trend was observed for maximum wave height, H . This deviation from the design value can be attributed to wave steepening due to the opposing current, whereby the opposing current causes the wavelengths to shorten and the wave heights to become larger in order to conserve energy, whereas aligned currents cause waves to flatten [71]. Table 1. Hydrodynamic forcing during the experimental runs. Symbols h: mean water dept U h;: depth-averaged, mean current speed (from pump calibration); ϕ relative angle between the mean current and the direction of wave propagation; H : mean wave height; H : zero-moment significant wave height; H : maximum wave height (* indicates values derived from measured velocities using linear wave theory); T : peak wave period; Z : height of ADV measurement (sampling volume) above bed at start of run; 〈U 〉: mean (ensemble-averaged) flow velocity at height of ADV sampling volume.

Bed Morphology and Sediment Motion
The observed bed morphology in each run was inferred from the acoustic sonar images and underwater photographs, revealing mobile bedforms and active sediment transport in all the runs. The dimensions of these bedforms, and modes of transport are presented in Table 2, with plan view sonar images of bed morphology at the end of each run, shown in Figure 3. The wave-only test runs were characterised by two-dimensional, symmetric wave ripples, with vortex shedding in shallower water in Run W1 (λ ~ 0.08 m), and much larger, rolling grain ripples in deeper water (W2) with larger waves, with wavelengths (λ ~ 0.16 m) scaling with the bed orbital diameter (d ~ 0.14 m). Run A1 was characterised by 2.5D symmetric [72,73], rolling-grain ripples with equilibrium wavelengths λ ~ 0.12 ± 0.001 m, scaling to the wave orbital diameter d ~ 0.1 m), and continuous but not straight crest lines with relatively rounded crests and ripple heights η ~0.056 m. This run was largely limited to bedload transport, active only when the orbital velocity and current were aligned, with vortex shedding entraining small suspensions (~6 cm) at flow reversal. Runs A2 and O2 were characterised by two distinct scales of bedforms: 2D, relatively symmetric combined flow ripples with intermediate equilibrium wavelengths (λ ~ 0.19-0.21 m), rounded crests with local deep scour at the lower end of the stoss side and discontinuous crest lines of irregular planform (η ~ 0.056 − 0.065), over which very small current ripples were superimposed. The wavelength of the dominant bedforms is clearly dictated by near-bed orbital diameter d , ~ 0.2 m. The dimensions of the smaller, superimposed ripples were inferred from a continuous wavelet transform applied at the centre-line of the SSS sonar image, revealing small wavelengths of ~λ ~ 0.03 ± 0.01 . Such superposition of bedform scales has been described for combined wave and current flow experimentally [72]; and in the field [74,75]. In both test runs, large, intense suspension clouds were observed at flow reversal, indicating entrainment by coherent eddies through vortex shedding. These were more persistent, and extended higher within the water column in the aligned case (A2), while O2 featured more defined, and shorter-lived suspension events. A3 is dominated by current-induced ripples, with pronounced three-dimensionality. This suggests that the strong current (U = 0.5 m) in A3 dominated the transport processes in deeper water, despite the slightly larger waves. In Run O1, the bed featured small, regular bifurcating asymmetric ripples, skewed on the current lee side, with a flatter upstream slope and steeper downstream face, with λ ~ 0.09 ± 0.008 m; and η ~ 0.035 ± 0.007 m. In Run O3, where the wave orbital velocities increase in comparison to O1 for similar opposing current, these transition towards more rounded, largely 2.5D with much longer wavelengths λ ~ 0.165 ± 0.008 m (scaling to d ~ 0.16 m), and much more pronounced scour on the lower end of the stoss, with ripple height η ~ 0.078 ± 0.009 m. The weaker currents in O2 and O4 have a less significant impact on the morphology, with 2D ripples as opposed to the currentmodified geometries in A3 and O3. Sediment transport appears as avalanching down the ripple face when the flows coincide (opposing waves), and reverts to intense vortex shedding at flow reversal entraining large suspension clouds extending several ripple heights into the water column. Despite deeper water and weaker currents in Run O4, enhanced waves resulted in more intense transport, with suspended sediments turning into intermittent sheet flow every half cycle, akin to plug flow conditions, where sediment starts to move instantly as a solid block (few mm deep) due to strong pressure gradients associated with accelerations at flow reversal [76,77].

Near-Bed Shear Stresses
The time-averaged (test run mean) and maximum, local shear stresses at both ADV elevations are presented together with stresses estimated using the models of GM86, SC05, MD12 and the corresponding maximum oscillating boundary layer thicknesses in Table 3. The maximum, theoretical thickness of the modified oscillating boundary layer, hereafter referred to as the combined wave-current boundary layer, δ , was of the order of O(δ ) = 1-2 mm for all runs using SC05 and MD12; smaller than estimates using GM86, where O(δ ) = 10-30 mm. Note that SC05 and MD12 predict skin-friction only stresses (and hence u*), while GM86 accounts for bedforms through a transport-enhancement term in Nielsen's modification. Nonetheless, these predicted thicknesses suggest that both ADVs fall outside the theoretical oscillating boundary layer in all runs. Therefore, the measured stresses at the corresponding ADV elevations are smaller than the predicted bed shear stresses, as expected, and increase towards the bed (at ADV2) in all the runs. Predictably, runs with higher mean waves registered larger stresses than those of comparable geometry in the aligned runs. Run O3, characterised by strongest currents, resulted in higher mean and maximum stresses, compared to O1, which, otherwise, has identical forcing, and flow geometry. The measured stresses agree with the observed transport modes, with higher stresses resulting in enhanced transport. In Runs A1, A3, O1 and O3, with similar mean and maximum shear stresses recorded by ADV2, bedload and vortex shedding were observed. Run O4, which registered the highest stresses, featured the most intense transport, with intermittent, plug flow at half cycles. Table 3. Estimated and predicted mean and maximum shear stresses, and corresponding maximum oscillating (/combined wave-current) boundary layer thickness, δ , derived from shear velocities inferred through the turbulent kinetic energy method for ADV1 and ADV2; and the models of GM86 (effective: τ , , and transportenhanced, τ , , ), SC05, and MD12.

Intermittent Coherent Motions
Quadrant analysis is used to define the intermittent nature of momentum exchange, manifest through the Reynolds stress. Hereafter, 'intermittency' refers to the alternate passage of energetic coherent structures, which denote active flow periods identified through quadrant and wavelet analysis [78,79]; rather than the classical phenomenological definition [80,81]. The passage of such motions may result in observed deviations from the classical (Kolmogorov, 1941) description of turbulence in the inertial range [82], but we focus specifically on the implications to shear stresses and momentum exchanges in the vertical. The distribution of the streamwise-vertical Reynolds stresses (u w in (x-z)) plane, across the four different types of quadrant event structures for each of the test runs is presented in Figure 4, for the exemplary runs A1 and O1. For both cases, the joint distribution of streamwise and vertical velocity fluctuations is presented at the two ADVs, without (blue) and with imposed quadrant thresholds with hyperbolic hole sizes H = 1 (black) and H = 4 (orange), as defined by Equation 5. H = 0 defines all detected events, H = 1 delimits values exceeding the root-mean-square stress, and H = 4 is used to identify extremes (which constitute 1-3% of recorded events, in all runs). From the outset, there are clear differences between the aligned and the opposing cases, and between the ADVs at both elevations in a given run. In Run A1, the streamwise component visibly dominates the flow structure at the higher ADV1, and the most extreme events, highlighted to be outside H = 4, mostly occur in the 1 st and 4 th quadrants (Outward Interactions and Sweeps), characterised by positive velocity fluctuations (i.e., along the direction of wave propagation and current). This pattern is reproduced in the other aligned cases, and the most extreme events often feature in a chronological sequence manifest within Q2 (Ejection) and Q4 (Sweep) pairs [counterrotating motions], suggesting intense turbulent events are associated with the passage of spatially defined energetic structure past the sensor. Closer to the bed, the instantaneous Reynolds stresses appear more evenly distributed amongst the 4 quadrant types (isotropic). In the opposing cases (O1 to O4), Inward and Outward Interactions are prevalent in all cases at the higher elevation, while Sweeps and Ejections feature more prominently closer to the bed. The opposing cases exhibit more variance between one another, due to the different current conditions, with runs O1 and O3 (having the strongest currents), featuring more frequent extreme Inner and Outer Interactions, whist O2 and O4 (lowest currents) are more evenly spread across all quadrants. Collectively, the quadrant analysis suggests anisotropic distribution of turbulent Reynolds stresses at the higher elevations, and more isotropic turbulent fluctuations at the lower elevation. This implies that the two ADVs fall into distinct "layers" within the flow, with a stronger current influence at the higher elevation.
The contribution of different types of quadrant events to Reynolds stress (in the streamwise vertical plane) is shown in Figure 5. Without applying a threshold, Outward Interactions (OI) and Sweeps contribute more (29-37%, and 28-43.7% respectively) to the time averaged Reynolds shear stress than the Inward Interactions and Ejections (9-24.1 and 13-19%) in all runs and both elevations above bed. Given that all events appear in nearly equal proportions, then Sweeps and Inward Interactions (both characterised by positive fluctuations in the direction of wave propagation) are therefore the more violent, stress-bearing structures. With a threshold H = 1 applied, the aligned cases are dominated by flow-aligned motions (Q1 and Q4) which contribute 58-93% of stress at both elevations in these runs. Sweeps are also the major contributors to stress in the opposing runs at the lower locations, but collectively, Outward and Inward Interactions comprise 62-77% of the time averaged stress at the higher elevations, and 40-53% in the opposing wave-current cases. Currentaligned motions are more prevalent and dominant at ADV1 elevation, confirming the different nature of vortices between the two heights.

Vertical Turbulent Fluxes and Eddy Scales
To describe the vertical exchanges within the flow, the vertical length scales of turbulent motions are derived from the normalised, wave-number weighted, turbulence power spectra of the vertical turbulent component, w ; constructed using Welch's method with 50% overlapping Hann windows [83]. These are presented in non-dimensional wave-number space (k * = kz, where z is the height of measurement) in Figure 6; in an energy preserving form such that an equal area under the curve represents equal energy [59]. The peaks are then used to derive spatial scales as outlined in §2.4. These spectra show that the bulk of energy in the vertical motions is contained at non-dimensional wavenumbers in the range: 10 < k * < 10 ; in agreement with previous work in shallow settings [59,84]. For the heights of measurements at ADV2 (between 0.07 and 0.13 m above bed), this corresponds to length scales between 0.01 < λ < 0.3 m, exceeding the inferred and theoretical thicknesses of the combined wave-current boundary layer predicted by all models. The peak energies correspond to similar length scales in all the runs (presented in Table 4), in the order of λ = 0.03 − 0.07 m; with the exception of Run O1. These integral scales of vertical fluxes exceed the theoretical thicknesses of the combined wave-current boundary layer predicted by all models; scaling to ~1.2 -3 times the thicknesses δ , inferred from measurements and the GM86 model, and an order of magnitude larger than predicted by SC05 and MD12. This implies that vertical fluxes are not restricted to within the oscillating region, as suggested by the zero-flux condition imposed in combined flows. Similar scaling is observed in the vertical at the higher elevation, with larger peaks ( λ = 0.027 − 0.12 m ). To assess the nature of these exchanges, the rate of vertical turbulent transport (turbulent kinetic energy flux in the vertical), T = ∂ k′w′ (where ∂ denotes the vertical gradient, [23]) is also presented (Figure 6, c). T represents an estimate of turbulent advection, given molecular diffusion is negligible at high Reynolds numbers (and we cannot calculate the pressuredriven transport). The vertical turbulent flux in the deeper runs (W2, A3 and O4) is largely negligible at both elevations. A clear distinction between the aligned and opposing cases is evident at the higher ADV location, whereby [shallower] aligned cases are characterised by upward diffusion (negative T ), whereas the opposing cases show downward net transport of momentum (positive T ). Closer to the bed, all runs display upward diffusion, except those in deeper water suggesting bed generated turbulence dominates. W1 shows significant negative flux at lower ADV2 (upwards transfer of momentum from bed), i.e., eddies from the wave boundary layer transferring momentum into upper flow. This is expected given W1 is characterised by vortex shedding; and these vortices do not reach ADV1 (higher); whereas W2 is characterised by rolling grain ripples. Similarly, runs A2, O1, O3 and O4 exhibit significant negative turbulent fluxes at the lower elevation, indicating upward transport from the bed through vortex shedding. This is not evident at O3, which also features vortex shedding.

Temporal Scales and Active Periods of the Flow
The temporal scales of turbulent motions and active periods of flow are derived from the high variance regions within the time-frequency domain defined by a continuous wavelet transform (CWT) decomposing the instantaneous streamwise-vertical Reynolds stress ( u′w′ signal. This is presented for both ADV elevations, for the exemplary runs W2, A2 and O2 in Figure 7 as timefrequency plots, but we opt to represent period (inverse frequency) on the y-axes for clarity. In the CWT time-frequency plots, brighter colours indicate high power (variance), and all the plots were scaled to the same limits to allow comparisons of the magnitude of wavelet variance/power at local and global scales between the different runs and elevations. The global power, for all runs, is presented in Figure 8, and pertains to the time average of all local wavelet spectra (i.e., all the vertical slices through the wavelet plot); with the red lines in these subplots denoting the 95% confidence level for the global wavelet spectra, tested in a time-averaged sense using autocorrelation of background noise. Short-lived, intermittent events are evident at high frequencies (>1 Hz), in all runs, often spanning multiple frequencies and lasting up to a maximum of ~10 seconds, and featuring more prominently at the lower elevation (ADV2), with weaker power higher up (ADV1). While some of these exhibit frequencies that are significant globally within the series, most show significant peaks in the local wavelet power but are less significant globally. Closer to the bed, intermittent, short-term high frequency events occur more commonly throughout the time series and are of both global and local significance. In the wave only runs, lower frequencies (higher periods) pertaining to the wave forcing, occupy significant proportions of the time series; particularly at ADV2. In the aligned flow cases, these are even more prominent, persisting continuously for up to 300 seconds; and thus dominate the global normalised wavelet power. Many of the lower frequency (longer period) events stem from these scales, suggesting higher frequency motions result from breaking of larger coherent clusters, as they dissipate energy. For the opposing runs, the CWT of streamwise-vertical Reynolds stresses generally exhibit active periods with weaker power (lower variance) than in the aligned cases, and often rapidly diminishing. These are characterised by intermittent, very short-lived events across small ranges of high frequencies (1)(2)(3)(4)(5) with active durations of ~3-9 seconds. Low frequencies (periodicity ~1-3 s) feature for shorter durations than in the aligned cases (up to 50 seconds) at half wave cycles (vortex shedding), but last for longer (~200 s) in narrow frequency bands at these scales. The observed events are also weaker in the deeper runs (W2, A3 and O4) than the shallower runs, and are shorter lived in runs where the waves have shorter periods (A1, O1 and O3). Runs O2 and O4, with weaker currents, have the weakest motions. ; for exemplary runs W1 (waves-only); A2 (current aligned with waves); and O2 (current opposing wave). In each panel, the warmer colours (yellows) indicate higher spectral power (variance); the white shaded region represents the cone of influence, and the thicker, black contour lines delimit the 95% confidence limits (5% significance against red noise given by a Monte Carlo simulation). The y-axes defines the frequency scales of turbulent momentum exchanges in the streamwise-vertical plane, represented here as periods (inverse frequency) for clarity. for all test runs. The global power pertains to the time average of all local wavelet spectra (i.e., all the vertical slices through the wavelet plot); and the red lines in these subplots denote the 95% confidence level for the global wavelet spectra, tested in a time-averaged sense using autocorrelation of background noise.

Discussion
In combined wave-current flows, shear stresses near the boundary are often assumed to be independent of whether the current flows with or against the waves, yet some variations in turbulent intensities have been reported in the literature. Experimental results have reported currents aligned with waves result in reduced stresses in the outer flow, suggesting the generation of turbulence is periodic with the waves at the bed and thus "true turbulence" may be lost in the Reynolds decomposition needed to infer stresses [10,11,26]. In opposing flow, turbulence intensities depend on the outer flow velocities generating vortices which peak at the deceleration phase of the wave [11]. Flow reversal against a weak current could be attained each wave cycle, resulting in the generation and ejection of an upstream vortex in addition to the usual downstream vortex. However, stronger currents weaken or inhibit the formation of upstream vortices, ejecting one vortex only each cycle (ibid.). Our results show that an equal magnitude current opposing waves (Run O1) appears to have a more significant impact on time-averaged shear stresses than in the aligned case (Run A1) and high, mean and maximum shear stresses are outside the theoretical oscillatory boundary thickness, which the GM86, SC05 and MD12 models suggest to be confined to a few millimeters in the immediate vicinity of the bed. Yet, we report shear stresses exceeding those predicted by the models at the bed, with measured maximum stresses exceeding those estimated by SC05 and MD12 and the effective stress of GM86, albeit lower than the transport-enhanced bed shear stresses predicted by the latter. The three models use estimates of near-bed orbital velocity to define a wave-only shear stress, and a current velocity (depth-averaged or at a height within a presumed logarithmic velocity distribution) to derive a current only shear stress, then apply a vector summation of the two components to account for the relative angle. For the individual contributions of current and wave, the quadratic stress law is used to calculate mean and maximum bed shear stresses, using the mean current drag coefficient or currentmodified wave friction factor. The main difference among these models arises from the treatment of the current-modified wave (oscillatory) shear stress. The near-bed horizontal orbital velocities measured by the lower ADV2 ( U = (var(U ) + var(V ) ) . compared remarkably well (within 10%) to current-modified orbital velocities estimated from independently-measured depth-averaged current, wave height and period, using an efficient Newton-Raphson numerical approximation to solve the wave-current dispersion equation [71,85]. GM86 uses a simplified scaling through a generalised friction factor applied to wave-only shear (friction) velocity, to account for shear stress enhancement, such that u * = C u * ; where the coefficient C is slightly larger than unity (enhancement). The model of SC05 uses a combined shear stress defined by τ = τ + τ , with modified drag coefficients used in the quadratic stress laws, and defines the maximum shear stress τ , as the period zero-mean stress corresponding to the wave-induced oscillation about the mean shear stress τ , . A comparison on the period zero-mean stress, on a wave-by-wave basis, shows little difference to the ensemble-averaged stresses reported in Table 3. A comparison of wave-only (τ ) and current-only (τ ) stress components to a linear summation of both, derived following SC05 and MD12 suggest relatively equal contributions (τ /(τ + τ )~ 0.4 − 0.6); while using GM86 suggests wave dominance in all runs τ /(τ + τ ) > 0.5. Predictably, runs with higher mean waves registered larger stresses than those of comparable geometry in the aligned runs. Run O3, characterised by strongest currents, resulted in higher mean and maximum stresses, compared to O1, which, otherwise, has identical forcing, and flow geometry. In all runs, the lower ADV shows non-linear enhancement, suggesting interaction between the wave and current, and thus suggesting the combined oscillatory boundary layer must extend to that elevation. The observed stresses also tally with the observed modes of transport, with higher stresses manifest in enhanced transport. In Runs A1, A3, O1 and O3, with similar mean and maximum shear stresses recorded by ADV2, bedload and vortex shedding were observed. Run O4, which registered highest stresses pertaining to the largest recorded waves, despite a weaker current, featured the most intense transport, with intermittent mass transport (plug flow) at half cycles, and associated with strong pressure gradients at flow reversal.
All three models, used above to predict stresses, assume a two-layer, discontinuous timeinvariant eddy viscosity profile, accounting for near-bed modification of the wave orbital velocity by the current, but assuming velocity outside the combined boundary layer must equate to the current only velocity; hence subscribing to the zero-flux condition imposed on vertical momentum exchanges. Our results show momentum exchanges at both ADV elevations. The production of Reynolds stress (P ; calculated following [23]) was more prominent at the higher ADV1 (P ~0.4 − 1.8 × 10 m s in aligned cases; and lower (P ~ 6× 10 m s ) in the opposing cases; with opposite patterns closer to the bed at ADV2 (P ~2.8 m s in aligned runs, and negligible in the opposing runs). The rate of production is significantly higher for lower wave periods (A1, O1 and O3). This is reflected in a downward flux of turbulent energy observed in opposing flow in contrast with an upward transfer of momentum in the aligned case (from the oscillatory layer to the outer flow). Quadrant analysis also reveals that the aligned cases were dominated by intense, currentaligned stress bearing structures (Sweeps and Outward Interactions), with counter-rotating vortices (pairs of Ejections and Sweep) being more prominent in longer wave periods, with more time available for the formation of these vortices before flow reversal.
In the opposing case, current aligned motions were also prevalent, but the near bed flow was characterised by counter-rotating pairs of Inward and Outward Interactions, with a downward flux of Reynolds stresses produced within the outer flow and shorter-lived, intermittent active periods characterised by high frequency fluctuations. This could imply higher dissipation of energy imparted by the adverse current upon the eddies formed at flow reversal. An increase in current magnitude resulted in enhanced vertical turbulent fluxes from the boundary layer into the near bed region in the aligned case, and vice-versa in the opposing case. A stronger opposing current also resulted in shorter-lived energetic, turbulent eddies spread across more length and frequency scales; hinting at an enhanced dissipation of the wave energy. To assess this, the turbulent energy dissipation rate, ε, is calculated, at the lower ADV2, following the method of [44], and presented in Figure 9 for frequencies within the inertial and dissipation ranges. The figure shows stronger turbulent energy dissipation associated with the opposing cases compared to aligned flows, and enhanced dissipation linked to stronger currents. Dissipation is stronger in the vertical, and peak dissipation occurs at frequency bands linked to those identified by the CWT in Figure 7. This confirms that the intense motions are correlated with strong dissipation events. Spatial scales of turbulence indicate well-formed eddies in aligned and opposing flow scaling to the wave orbital motions with vertical excursions beyond the theoretical oscillating layer thickness, but an increase in current magnitude in reverse reduces the extent implying more dissipation of eddies formed at the bed. Wall-attached structures tend to be formed by paired Sweep-Ejection events, often with an associated vortex cluster (detached and isotropically-oriented) embedded within the ejection [86,87] measured the streamwise advection velocity of coherent structures and implied that these are deformed by the mean shear which also controls their lifetime (duration through which they persist before they are completely consumed by viscosity). [88] also suggest that the mean strain rate, which depends on the type of structure considered, is skewed along the positive flow direction, amplifying vorticity. We therefore suggest that in aligned flow of comparable wave and current strengths, the steady current dictates the mean direction of strain, favouring positive streamwise motions (Q1 and Q4 events). This agrees with the observed percentages of occurrence; and, with Q2 and Q4 motions being the most violent, these vortex pairs provide the most significant contributions to stress in the streamwise-vertical plane. This could also explain why, higher up, the maximum normalised shear stresses appear to be a linear superposition in the aligned cases, of nearly equal strengths of the current and waves, in comparison to the nearbed stresses.
Cross-examination of the active periods of flow with the time series of quadrant-event Reynolds stresses (delimited by threshold H = 1) allows us to define sequences of intense bursting motions contributing significantly to stress and momentum exchanges during these periods, and therefore transport of sand. In Run A1, the active, lower-frequency period between t = 230 -340 seconds coincides with a cluster of large-magnitude Outward Interactions and Sweeps, at the higher elevation. These are interspersed with several successions of very intense pairs of Sweeps and Ejections ( τ (S, E) ≅ 4.9 Pa ) followed by weaker interaction pairs (OI & II) between 238 -275 seconds, and then again between t = 330 -340 sec, and at t~350 s. Closer to the bed (at ADV2), Outward and Inward Interactions are featured throughout the time series (τ (OI, II) ≅ 2.8 Pa) yet all the significant regions within the CWT plot only coincide with the presence of Sweeps, and the lower frequency active regions are instigated by clusters of strong Ejections and Sweeps. All the Ejections are concurrent with active flows with a wavelet period of ≥1 s. Moreover, comparisons between A1 and A2 highlight the effect of wave period in aligned flows whereby the longer period waves (A2) resulted in a more discernible impact on observed trends, such that active periods of flow better mapped to the wave harmonics, with stronger upwards fluxes from the oscillatory boundary layer into the near bed region, and more symmetric, 2D bedforms associated with oscillatory flows. This is reflected in the observed sediment entrainment by vortex shedding at flow reversal in Runs A1, A2, O1 and O2, which were longer lived, and extended higher up, in the aligned cases; and were shorter, and quickly advected in the opposing cases, particularly under stronger advected currents in Run O1. The observed clouds of suspension corresponded to the inferred vertical turbulence length scales in §3.4.1. This possibly highlights the role of the superimposed current, whose mean direction governs the prevalence of a positive or negative streamwise fluctuation (with orbital motion switching signs every half cycle) in defining the predominant types of detected turbulence events and ensuing suspension, while the stronger waves dictated intense bursts of transport, as evident in plug flow transport at the intense stages of stress production where wave and current flows coincided in Run O4. [89] have shown Q4 (sweep) events to dominate near the wall (with Q2 away from the wall) as a result of a "spatial phase difference" between streamwise and wall-normal turbulent fluctuations caused by the advection of a large coherent structure which causes counter-gradient stresses (Q1 and Q3) events above and below the structure. This causes more negative peaks in fluctuations at higher wall units than positive ones, and thus more frequent Q2 motions higher up. Such intense sequences of sweeping motion therefore enhance the mobilisation of sediment, resulting in observed, intermittent sheet flow conditions for much of the wave cycle in Run O4. Collectively, these observations attest to the differences in behaviour between aligned and opposing, collinear wave-current flows. This implies that our present approach to dealing with the problem from a timeaveraged perspective is insufficient, because there is a need to capture the time-resolved momentum exchanges throughout the water column.
This work presents a study focusing on stress and momentum exchanges in combined wavecurrent flows. In the absence of three-dimensional, time-resolved information at field scales, advanced computational techniques are relied upon to extract temporal and spatial information on fluid motion (Reynolds stresses and quadrant events) and suspended particles, measured at discrete locations with the flow. This is warranted given that quadrant events, often considered manifestations of the passage of coherent structures, reside at the top of the turbulent energy cascade, and are important for the transfer of momentum and generation of turbulent energy, and hence sediment entrainment. An assessment of shear stress and its distribution across different quadrant events as well as its variability in time often provides a better account of sediment entrainment and transport [84,90]. What remains to be explained is the feedback of the suspended particles, once entrained, on the vortex structure itself. Unfortunately, the limitations of this experiment meant that no joint analysis of the temporal and spatial attributes of the turbulence and suspension signals could be undertaken.

Conclusions
In this study, turbulence measurements at two elevations, outside of the theoretical thickness of the oscillatory boundary layer, are used to examine fundamental assumptions on the nature of shear stresses and momentum fluxes in collinear flows aligned with and opposing the direction of wave propagation over mobile bedforms. Our observations suggest that the core flow (outside of the combined wave-current boundary layer) is not uniform (purely current dominated) with its lower portion comprising eddies both advected downwards/upwards depending on prevailing interaction between the steady and unsteady shear components; as well as locally generated by the bed morphology. This highlights the limitation on existing models which impose a zero-flux condition at the interface. Waves, contribute to pick-up by vortex entrainment, making material readily available for transport by net currents in the nearbed and outer flow regions. The entrainment of sediment and maintaining it in suspension is associated with the passage of convected or locally formed packets of eddies linked to large and intense coherent turbulence structures within the flow [24]. The intermittent nature of both mass and momentum transfers thus warrants a re-assessment of the existing models for wave-current boundary layer thickness (to account for intermittent exchanges of momentum) and the invalidity of the zero-flux condition commonly imposed at its the interface with the outer layer.
The key findings from our experiments can be summarised as follows: (1) Combined (2.5D) wave-current ripples were observed in all runs with comparable wave and current strengths, with intense suspensions at flow reversal that persist-longer in aligned flows.
(2) Near-bed interactions outside of the theoretical combined flow boundary layer were not uniform throughout the outer flow region, with discernible difference in flow structure and momentum exchange at the two elevations.
(3) The aligned flows are characterised by upward turbulent fluxes with Reynolds stresses produced throughout the outer layer, whilst opposing flows featured downward turbulent flux of eddies generated in the core flow.
(4) Spatial scales of turbulence indicate similar vertical scaling of eddies with fluxes extending beyond the theoretical thickness of the combined layer, often imposed as a zero-flux boundary in numerical models. It follows that the oscillatory (combined wave-current) boundary must be thicker than predicted by commonly used combined flow models.
(5) In aligned flow, active lower-frequency (wave-related) periods of flow are characterised by pairs of counter-rotating vortices implying longer-lived wall-aligned structures prevail. Opposing flows are characterised by shorter-lived intense stress-bearing structures, hinting at higher dissipation rates. Intense momentum exchanges (with vortex shedding inducing suspension clouds) and plug flow sediment transport conditions are observed when strong adverse currents correspond with the reversing wave motion within each cycle. Funding: This work was undertaken during the lead author's PhD candidacy at the University of Southampton, through a scholarship funded jointly by the Graduate School at the National Oceanography Centre Southampton, the Faculty of Engineering and the Environment and the Southampton Marine and Maritime Institute at the University of Southampton, and HR Wallingford Ltd. Access to the Fast Flow Facility was granted during its inaugural testing by HR Wallingford, as part of their contribution to the PhD.

Conflicts of Interest:
The authors declare no conflict of interest.