Observations of Enhanced Sediment Transport by Nonlinear Internal Waves

The mechanisms responsible for sediment resuspension and transport by nonlinear internal waves (NLIWs) remain poorly understood largely due to a dearth of detailed field measurements. We present novel observations of the turbulent benthic boundary‐layer (BBL) beneath trains of NLIWs of depression in the ocean. At the 250 m deep, low‐gradient (<0.2%) continental shelf site the BBL was near well mixed to an average height of about 10 m above the bottom. Above this bottom mixing‐layer, stratification constrained the extent of vertical sediment transport. NLIWs drove sediment transport by a combination of bed‐stress intensification, turbulent transport, and a vertical pumping mechanism associated with the compression and subsequent expansion of the mixing‐layer. There was no evidence that the observed dynamics were associated with a global instability, as proposed in previous studies. The results have implications for cross‐shelf mass transport and highlight future challenges for measuring and modeling boundary‐layer processes within shelf seas.

The instability is global as it depends on shear in the near-bed vertical profiles of horizontal velocity, which, driven by pressure gradients beneath the wave, vary along the direction of wave propagation. That is, it is a 2-D or three-dimensional (3-D) mechanism driven by properties of the global flow field (for a more detailed description, see Diamessis & Redekopp, 2006;Huerre & Monkewitz, 1990;Sakai et al., 2020). This is distinct from the more canonical local instability, which is dependent solely on the shape of the local vertical profile of velocity. Such concepts are predominantly applied to the study of laminar to turbulent transition and arise less commonly in the study of fully turbulent flows. Large NLIWs are typically found, however, on energetic continental shelves with highly turbulent BBLs, often with near-bed stratification. Modeling NLIW influence on a fully turbulent, stratified BBL is highly challenging in both laboratory and numerical simulations, and no such studies have been performed to date. The influence of NLIWs and NLIW-induced pressure gradients on a fully turbulent BBL, including flow separation and the role of a GI or other mechanisms for transition to turbulence, is thus largely unknown. The fate of a GI as it ultimately decays into turbulence is also poorly understood as numerical simulations are, with the notable exception of Sakai et al. (2020), typically halted before the wave-induced stirring cascades energy to fine-scale motions. Further, only the studies of Harnanan et al. (2015Harnanan et al. ( , 2017, Sakai et al. (2020), andSoontiens et al. (2015) are 3-D, as required to resolve a turbulent cascade.
As highlighted in the recent review by Boegman and Stastna (2019), our limited understanding of sediment transport mechanisms in the turbulent ocean is due in part to a lack of appropriate field observations from which to either (1) extract process understanding directly or (2) guide process-oriented simulations. Here we present new, highly resolved field observations of large-amplitude Mode-1 NLIWs of depression propagating over near-flat topography. The observations include simultaneous time-series measurements of the NLIWs in the overlying water column, BBL velocity profiles, comprehensive near-bed turbulence statistics, and observations of suspended sediments.

Experiment
This study was part of the Kimberley Internal Solitons Sediment and Mixing Experiment (KISSME), conducted in April-May 2017. The study site was a broad region of low-gradient (<0.2%) topography in the Browse Basin on Australia's Northwest Shelf (NWS) (Figure 1). We deployed three through-water-column moorings in a triangular array with side length of 500 m, each measuring temperature, currents, and conductivity to characterize the NLIWs (see Rayson et al., 2019, for more details). In addition, a bottom landing frame (herein lander) was deployed approximately 100 m from the southernmost mooring, to characterize the boundary-layer's response to the NLIW forcing. The instrumentation on the lander included a 1 MHz ADCP sampling at 8 Hz with 0.2 m vertical resolution from 0.55-20 m ASB and two 6 MHz ADVs (at 0.49 and 1.4 m ASB) with coupled thermistors sampling at 64 Hz. Two additional fast thermistors (at 0.1 and 1.1 m ASB) were mounted to the lander, and further 10 (3.0, 3.6, 4.5, 5.5, 6.5, 8.1, 10.1, 12.1, 14.2, and 16.2 m ASB) were mounted to a buoyant thermistor string (herein T string), each sampling at 2 Hz. One additional logger attached to the top of the T string (at 20.3 m ASB) sampled both temperature and pressure every 15 s, and optical backscatter (NTU) was recorded at 1.1 m ASB every 90 s. Mooring displacement (herein knockdown) was evaluated using the upper pressure sensor. In addition to the moorings, ship-based CTD profiling took place on 3 April.

NLIW Overview
We observed numerous Mode-1 and Mode-2 NLIWs throughout the deployment, though this paper focuses on the large Mode-1 waves of depression observed in early April (Figure 2a). Rayson et al. (2019) demonstrated that these waves likely formed at or beyond the outer shelf break and propagated toward the southeast ( Figure 1). The NLIWs had Mode-1 amplitudes (A 1 ), estimated from the observed density perturbation using weakly nonlinear (WNL) theory (Rayson et al., 2019), of up to 70 m and were accompanied by significant sediment resuspension events (Figures 2a and 2b).

Boundary-Layer Overview
A bottom mixing-layer with weak density stratification and small density fluctuations was a persistent feature of our observations (Figure 2c). This mixing-layer was separated from the stratified interior by a thin layer of locally elevated buoyancy frequency (N), herein referred to as the mixing-layer pycnocline (MLP). Similar layers form in numerical studies of turbulent stratified flows over flat beds (e.g., Gayen et al., 2010). We defined the thickness of the mixing-layer (δ ML ) as the height above the bottom where N first reached a value of 0.009 s −1 . The record of δ ML is not continuous as the mixing-layer thickness occasionally exceeded the height of the T string, which varied due to current-induced T-string motion. δ ML and A 1 were well correlated, with δ ML ranging from 4 m under forcing of NLIWs of depression to over 20 m when A 1 was positive (Figures 2b and 2c).
Backscatter intensity and the thickness of the near-bed layer of elevated backscatter (δ BS ) were greatest immediately following the passage of the NLIWs of depression ( Figure 2b). We defined the scale δ BS | 45 (δ BS | 60 ) as the bottom sediment-layer thickness corresponding to an echo intensity threshold of 45 dB (60 dB). When δ ML could be estimated it was highly correlated with δ BS , with δ ML typically falling between δ BS | 45 and δ BS | 60 (Figures 2b and 2d). Assuming that sediment concentration remained detectable across the entire mixing-layer, δ BS | 45 and δ BS | 60 provide continuous estimates of the mixing-layer thickness and indicate an average value of~10 m over the period (Figure 2d).

NLIWs of Depression of 3 April
We now focus on the train of Mode-1 waves of depression observed around 8 a.m. on 3 April, which had leading wave amplitude A 1 of~70 m and a (Doppler-shifted) period of~20 min ( Figure 3). The nondimensional amplitude of A 1 /H~0.3 exceeds the formal limits of WNL theory, but we use WNL theory only to calculate a characteristic wave height. At the moment the wave train passed the mooring, the barotropic current speed was 0.17 m s −1 with a bearing of 324°. The nonlinear wave velocity (c), estimated using a lagged correlation of the A 1 time series from each mooring in the triangular array, and accounting for advection of the waves by barotropic currents, was 0.95 m s −1 at 134°. The associated wave Reynolds number (Re w = cA 1 /ν) was O(10 8 ).
A CTD profile prior to the arrival of the wave confirmed the approximate location and temperature jump of the MLP derived from the T string. A second CTD profile taken after the arrival of the wave supports the relation between mixing-layer thickness and the sediment-layer thickness described in the previous section ( Figure 3b). The sediment boundary-layer thickness δ BS | 60 did not increase appreciably until after the passage of the trough, at which time it increased from about 4 m to almost 10 m in a time frame O(10) min. The layer thickness remained near 10 m during the subsequent passage of the trailing waves, with the small fluctuations in δ BS | 60 correlated with the amplitude of the wave train ( Figure 3c). Appreciable velocity gradients occurred only over the region with appreciable backscatter (Figure 3c), suggesting that sediments were transported to the very edge of the turbulent boundary-layer.
We calculated turbulence quantities over 5 min subsections of each ADV record. We calculated Reynolds stress ( u ′ 1 u ′ 3 ) and turbulent kinetic energy (TKE) by integrating under the ADV velocity spectra, and TKE dissipation (ε) using both an inertial subrange fitting method (ε IDM ; Bluteau et al., 2011) and by ε RS ¼ u ′ 1 u ′

3=2
=κz; where κ is von Karman's constant and z is the measurement height. At each ADV location the production (P) of TKE was estimated as the product of the local Reynolds stress and the vertical shear of the (5 min)

Geophysical Research Letters
ZULBERTI ET AL.
mean velocity from the nearest ADCP bins. The magnitude of all turbulence quantities grew as the wave approached, reached a peak near the trough, and then remained elevated during the passage of the subsequent wave train (Figures 3e-3h). The near-bed flow was highly turbulent, with TKE levels at the bed exceeding 4 × 10 −2 m 2 s −2 and dissipation rates ε exceeding 2 × 10 −4 m 2 s −3 (Figures 3e-3h). At 0.49 m ASB, dissipation ε and shear production (P) were in approximate balance (as expected in a local log layer) both prior to and after the passage of the wave. At 1.4 m ASB, however, dissipation exceeded P prior to the arrival of the wave, though an approximate balance was obtained once dissipation levels exceeded ε > 10 −5 m 2 s −3 (Figures 3e and 3g). We estimated the turbulent Reynolds number, defined as Re T = (δ BS | 45 /η K ) 4/3 , where η K is the Kolmogorov microscale, to be~O(10 5 ).
We estimated bed stress by the Reynolds stress method τ b ¼ ρ; u ′ 1 u ′ 3 and the quadratic drag method (τ bCd ¼ ρC d u 2 ) using measurements from the ADV at 0.49 m. Both methods assume that a log layer exists and extends to this measurement height. We used a C d,100 cm of 0.0016 as typical for a silt bed (confirmed by sediment grab samples) adjusted to a height of 0.49 m using log scaling (Soulsby, 1983). The two independent bed-stress estimates were in strong agreement, increasing from initially very small values to a peak under the leading trough and remaining elevated during the passage of the trailing waves ( Figure 3i). The bed stress exceeded the Shields critical stress for incipient sediment motion (τ crit ) more than 20 min prior to the arrival of the leading NLIW trough and remained greatly in excess of τ crit over the bulk of the wave train. The near-bed (1.1 m) turbidity signal responded almost immediately and continued to increase after the leading trough had passed (Figure 3i). The active resuspension prior to the arrival of the trough was not associated with growth of δ BS | 60 , and the sediment-layer growth observed after the passage of the trough occurred at a time when the bed stress was rapidly decreasing (Figures 3b and 3i).
Coherent structures within the sediment boundary-layer became apparent just prior to the arrival of the leading wave (Figure 4c). These structures had a quasiregular horizontal-vertical aspect ratio of approximately 3-2. As the thickness of the sediment layer grew after the passage of the trough the structures appeared to grow in size while maintaining their aspect ratio (Figures 4d-4f). Observable coherence was lost toward the end of the wave train (Figure 4g).

Discussion
The NLIWs reached wave Reynolds numbers 3 orders of magnitude greater than typical laboratory or numerical experiments (Re w ≈ 10 4 -10 5 ; Boegman & Stastna, 2019, and the references therein). The waves generated intense bottom stress resulting in vigorous sediment resuspension and subsequent transport high into the water column, in qualitative agreement with Bogucki et al. (1997). In this section, we discuss the nature of the turbulence responsible for sediment resuspension and the initial stages of vertical transport (section 4.1), and the role of the MLP in controlling vertical transport (section 4.2). We then introduce a two-phase pumping mechanism, unique to NLIWs of depression in the presence of a mixing-layer (section 4.3).

Sediment Resuspension Mechanisms
BBL turbulence prior to the arrival of the wave is a key feature that distinguishes these observations from past laboratory and numerical experiments of NLIW-BBL interaction. The approximate balance of P and ε (Figures 3e and 3f) and the agreement between τ b and τ bCd (Figure 3i) indicate the presence of a log layer of at least 0.5 m thickness for much of the wave duration. Conventional log dynamics suggest vigorous sediment resuspension occurred over the bulk of the wave train after τ b exceeded τ crit . If a GI mechanism was contributing to the resuspension, we might expect a peak in the near-bed velocity spectra, associated with the most unstable mode, to develop in the adverse pressure gradient (APG) region behind the trough and for this peak to continue to gain energy as the instability developed (e.g., Diamessis & Redekopp, 2006). We did not observe such a peak. Rather the velocity spectra were consistent with a conventional turbulent cascade (e.g., Tennekes & Lumley, 1972) at all times under the wave (Figures 4b-4g, Column 2). This suggests either that (1) no GI mode developed or (2) the GI energy levels were dwarfed by turbulence sustained directly by near-bed shear production.  Diamessis and Redekopp (2006) and others describe the GI mechanism forming vortical structures of consistent wavelength behind the wave and these vortices being subsequently shed from the bed. The large coherent structures evident in our backscatter observations, however, were present before the arrival of the trough, seemed to grow in wavelength with time, and appeared to remain attached to the bed (Figure 4). These characteristics are more representative of classical boundary-layer turbulence than GI. They also maintained an approximate horizontal to vertical aspect ratio of 3-2, typical of boundary-layer turbulence (Adrian, 2007;Adrian & Marusic, 2012).

Geophysical Research Letters
It is possible that the waves we observed simply did not satisfy the criteria for GI. It is generally thought that GI requires a threshold strength of APG (Diamessis & Redekopp, 2006), and this is influenced by wave shape and amplitude (Stastna & Lamb, 2008). The waves in previous GI studies have largely focused on large-amplitude steady soliton solutions to either WNL (e.g., Diamessis & Redekopp, 2006) or fully nonlinear equations (Stastna & Lamb, 2002) (a recent notable exception is Soontiens et al., 2015, who demonstrated instability under relatively low-amplitude developing waves near to their generation source). While our NLIWs were large amplitude, they were complex, time-evolving wave trains and not well described by such steady equations. These factors preclude direct evaluation of numerically predicted stability criteria for our field observations. Another feature of our observations is the lack of flow separation in the lee of the wave (velocity time series Figure 3c). Separation "bubbles" (where the boundary-layer separates, forming a region of flow reversal) are a common feature of GI simulations (though Stastna & Lamb, 2002, argued that the separation region was not a necessary precursor for GI provided particular conditions for upstream boundary vorticity were met). The breakdown of separation bubbles in strongly turbulent boundary-layers is well documented in other flow regimes (see, e.g., Kalter & Fernholz, 2001). This has not been documented for the NLIW boundarylayer, though Diamessis and Redekopp (2006) noted intermittent breakdown of the separation bubble at high Re w even in nonturbulent simulations. We do not yet understand the criteria for the existence of separation bubbles or coherent (shed) vortices in highly complex and energetic ocean environments, as no studies of GI to date have incorporated preexisting turbulence.

The Role of the MLP
Past laboratory and numerical studies of NLIW-BBL interaction (Boegman & Stastna, 2019, and references therein) have almost exclusively been initialized with zero or weak near-bed stratification (Harnanan et al., 2017, is a notable exception, though their boundary-layer rollup dynamic differs significantly from our observations). The vortex shedding, which occurs in these simulations, has thus been effectively unconstrained until the main thermocline. Another key feature of the present observations is the MLP, which limits the vertical transport of both momentum and suspended sediments (Figures 2 and 3).
Following the arguments of section 4.1, we assume attached boundary turbulence to dominate transport across the mixing-layer. Neglecting particle settling or inertia, the timescale for sediment transport from the bed to the top of the mixing-layer is dependent on the turbulence intensity and the mixing-layer thickness: Here u * ¼ ffiffiffiffiffiffiffiffiffi τ b =ρ p was O(1 cm s −1 ) under the leading trough. T T was thus O(10 min) around the time of the leading trough, similar to the wave period (T NLIW ), indicating that appreciable turbulent transport across the entire mixing-layer occurred under a single wave.
The correlation of δ ML and A 1 suggests that interior motion displaces the MLP at internal wave timescales, shifting this transport barrier vertically, and modulating this transport timescale. To confirm this we estimated the wave-induced MLP displacement behind the leading trough (Δδ ML ) using the wave-induced vertical velocities (u 3 ): 10.1029/2020GL088499

Geophysical Research Letters
ZULBERTI ET AL.
Taking the ADCP vertical velocity at 10 m ASB (Figure 3c) as an approximation for u 3 z ¼ δ ML ð Þ , we estimated an NLIW-induced MLP displacement of~5 m upward over half an NLIW wave period, which was approximately equal to the growth of the sediment layer over this time. The timescale for the mixing-layer to grow by an equivalent amount due to turbulent entrainment alone is Here δ p ML is the mixing-layer thickness under the trough and N 0 describes stratification beyond the mixing-layer (Valipour et al., 2015). We estimate that T e~O (10 hr), much longer than the period over which the apparent mixing-layer expansion occurred, indicating that the expansion was dominated by NLIW pumping rather than true mixing-layer growth (entrainment). Modulation of mixing-layer thickness by the NLIW-induced vertical velocity was thus critical to the observed sediment dynamics.

Conceptual Pumping Mechanism for Fine Sediments
The influence of the MLP, the rate of turbulent transport across the mixing-layer, and the NLIW-induced deflection of the MLP offer a conceptual model whereby an NLIW of depression enhances the suspension of sediment in a turbulent mixing-layer: 1. Just prior to wave arrival. Flow accelerates as the wave approaches, leading to increased velocity shear, and consequently enhanced turbulence intensity and sediment resuspension. The efficacy of turbulent transport across the mixing-layer also increases due to increased turbulence intensity and a reduction in transport distance (Equation 1) owing to the NLIW compressing the mixing-layer. 2. Under the trough. Concentration of fine sediments across the mixing-layer is the greatest owing to enhanced transport prior to the wave arrival and reduced volume of the now thin mixing-layer. 3. After the trough. Turbulence intensities drop, yet the mixing-layer rapidly expands under rising isotherms in the interior, rapidly dispersing fine sediment.
The contraction and subsequent expansion of a turbulent mixing-layer thus leads to an effective sediment pumping mechanism. Sediments must be sufficiently fine that particle settling can be ignored over NLIW timescales. To the best of our knowledge this pumping process, which is unique to NLIWs of depression, has not been described previously and explains our field observations of enhanced vertical sediment transport in a high-energy ocean environment without recourse to any particular instability mechanisms.

Concluding Remarks/Implications
We observed greatly enhanced sediment transport associated with the passage of large-amplitude NLIWs of depression. Sediment resuspension was dominated by intense boundary-layer turbulence typical of marine boundary-layers, and the subsequent transport was controlled by turbulent diffusion and wave-induced dispersion. Recourse to novel instability mechanisms was not required to explain the observations. Owing to the onshore propagation direction of these NLIWs of depression, this enhanced transport was associated with off-shelf flow near bed, providing a mechanism for net off-shelf transport of bed load. The very high ambient turbulence levels prior to the arrival of the NLIW, the range of scales of turbulence in the BBL, the existence of a bottom mixing-layer capped by a thin pycnocline (MLP), and the low boundary-layer thickness (relative to wave amplitude) were important features that distinguish the present observations from previous laboratory and numerical modeling studies.
The interactions between the boundary-layer turbulence, the MLP, and the NLIW-induced flow field are complex (see, e.g., Lee, 2017) and warrant further process-oriented simulation. Key unknowns are how NLIWs modulate turbulent structures and mean momentum profiles and whether vortex-shedding mechanisms would be capable of overcoming an MLP formed by balance between realistic ambient stratification and ambient turbulence levels. Realistic mixing-layer thickness will be key to such simulations, as this constrains the BBL thickness, and consequently affects the hydrodynamic stability of the boundary-layer. Challenges for future observational studies include reducing uncertainty in estimation of turbulence statistics to resolve subtle NLIW modulations, achieving greater vertical resolution across the rapidly spatially varying boundary-layer, and measuring sufficiently close to the bed to observe NLIW-induced pressure gradient effects.
Basin-scale ocean simulations based on Reynolds averaged Navier-Stokes equations do not resolve the BBL processes described here. A typically lower grid cell thickness in such models may be a similar thickness to the mixing-layer we observed, which is an order of magnitude thicker than the MLP that controlled near-bottom momentum and sediment transport (M. Rayson, personal communication, March 2020).
Careful parameterization of such processes may be necessary in many applications.

Data Availability Statement
Through-the-water-column mooring data from this experiment are archived on the UWA library research data repository (https://doi.org/10.4225/23/5afbf8fc55ed1). Bottom lander data used in this manuscript are also archived on the UWA library research data repository (https://doi.org/10.26182/5f62f6e361872).