Mixing Driven by Breaking Nonlinear Internal Waves

Non‐linear internal waves (NLIW) are important to processes such as heat transfer, nutrient replenishment and sediment transport on continental shelves. Our unique field observations of shoaling NLIW of elevation revealed a variety of different wave shapes, varying from relatively symmetric waves, to waves with either steepened leading‐ or trailing‐faces; many had evidence of trapped cores. The wave shape was related to the position of maximum density overturns and diapycnal mixing. We observed both shear (where sheared currents overcome the stabilizing effects of stratification) and convective (where the local velocity exceeds the wave propagation speed) instabilities. The elevated diapycnal mixing (>10−3 m2s−1) and heat flux (>500 Wm−2) were predominantly local to the NLIW of elevation packets, and were transported onshore 10s kilometers with the wave packets. We demonstrate that wave steepness may be a useful bulk property for the parameterization of wave‐averaged diapycnal heat flux.


Introduction
The dynamics of nonlinear internal waves (NLIW) can be critically important to both the vertical and crossshore transport of constituents such as heat, sediment, nutrients, plankton and contaminants (e.g., Boegman & Stastna, 2019;Pineda, 1994). NLIW of elevation (that travel on a thermocline below the mid depth) are particularly important to transport on continental shelves. NLIW of elevation can be generated directly, evolve from waves of depression through the largely energy-conserving process of fission on gradual slopes or through more dramatic breaking processes on steeper slopes (Aghsaee et al., 2010). Despite their widespread importance, waves of elevation remain poorly understood, and there are very limited field observations of sufficient detail to reveal the wave shape and resultant magnitudes of diapycnal mixing (e.g., Klymak & Moum, 2003;Shroyer et al., 2009). This lack of understanding currently restricts our ability to assess, for example, seasonal heat or nutrient budgets on the shelf.
Idealized numerical modeling and laboratory experiments have studied the evolution of NLIW of elevation via the interaction of the internal waves with idealized continental slopes and shelves (e.g., Aghsaee et al., 2010;Venayagamoorthy & Fringer, 2007). These studies have revealed that both the wave and topographic steepness determine the fate of a shoaling wave of depression; for example, a collapsing internal wave is characterized by large overturning at the trailing-face due to convective instability (sometimes referred to as a static or kinematic instability) (Aghsaee et al., 2010). In contrast, a bolus has large overturns at the leading-face due to a lobe and cleft type of instability (Venayagamoorthy & Fringer, 2012). Field observations of NLIW of elevation have begun to capture the range of different wave shapes, from symmetric internal solitary waves with closed cores (Scotti & Pineda, 2004) to more asymmetric waves with enhanced dissipation on the leading- (Xie et al., 2017) or trailing-face (Moum, Farmer, et al., 2007) of the wave. However, in the ocean we have not yet observed the full range of wave types produced in numerical studies and very few field studies have resolved the wave structure alongside turbulence estimates. Here we use unique field observations to explore the relationship between the shape and character of NLIW of elevation, and the resultant diapycnal mixing and heat flux.

Site Description
Non-linear internal waves (NLIW) are frequently observed on the Australian Northwest Shelf (NWS), including near our study site (Figure 1) (e.g., Gong et al., 2019;Holloway et al., 2001;Van Gastel et al., 2009). The topographical slope at this site is aligned at a bearing of 152 o from true north (Figure 1 inset). The topographic slope is relatively mild (0.15%) from the shelf-break (at depth of 135 m of water) to a midshelf rise (of slope 0.55%), approximately 25 km inshore of the shelf-break. In shore of the mid-shelf rise, the slope is very small (0.08%). On the shelf, the M2 semi-diurnal constituent dominates tidal forcing and the major axis of the barotropic tidal ellipse is approximately oriented in the crossshore direction at a bearing of 125 o true north .

Field Experiment
We deployed the Benthic Underwater Boundary Surveyor (BUBS) mooring from 4-14 April 2012 (day 95-105) over the bottom 35 m of the water column at the 102 m depth contour (Bluteau et al., 2017). We measured the mean velocity profile from 9.25-40.25 m above bed (mAB) with an upward looking acoustic Doppler current profiler (ADCP, 300 kHz Teledyne RDI) sampling with 60 s averages and 1 m vertical bins. An acoustic Doppler velocimeter (ADV; Vector, Nortek AS) was located at 7.5 mAB. The near-bottom vertical density structure was measured with 22 thermistors (SBE39 and SBE56, Seabird Electronics) spaced at approximately 1 m intervals from 3.65 to 31.4 mAB. Of these thermistors, the bottom-most 15 sampled at 2 Hz, while the remainder sampled at 0.1 Hz. We used data from the Australian Integrated Marine Observation System through water column mooring at PIL100 (collocated with BUBS), which during the period of interest measured mean currents (ADCP, Teledyne RDI) with 600 s averages and 1 m vertical bins. Density was measured using a combination of CTDs (SBE37, Seabird) and thermistors (SBE39, Seabird) at a vertical spacing of roughly 10 m . Finally, we used data from the Naval Research Laboratory N4I ADCP mooring, which measured currents at 150 s intervals and 1 m vertical bins.

Data Analysis
A proxy for the nonlinear wave amplitude Amp was extracted from the PIL100 data by least-squares fitting the 1 minute density data to modal structure functions . We estimated the linear phase speed c 1 in a similar way, but using the low-pass filtered density data. The NLIW speed c nl and direction θ was found from the BUBS ADCP data by using the backscatter phase-lag method of Scotti et al. (2005). The ADCP velocities were then corrected for beam separation using a phase-lag approach (Scotti et al., 2005). For a selection of 16 waves, we used the vertical displacement of the temperature contours at PIL100 to estimate the amplitude a. The wavenumber was estimated as k ¼ 2πTc nl, where T is the wave period corrected for any Doppler-shift effect due to mean currents advecting the waves. The Richardson number, Ri ¼ N 2 /S 2 , where N is the buoyancy frequency and S is the velocity shear, was estimated from 60 s averages of the resorted vertical temperature data and velocity data collected at ∼1 m vertical spacing. We combined the ADV data from 7.5 mAB with the ADCP data to increase the extent of the near-bed data. We interpolated the temperature data onto the depths of the velocity data and both data sets were smoothed before centered-differencing was used to find N and S. We estimated the Thorpe length scale L T , as a visual aid for the vertical extent and position of mixing events, from 10 s averaged temperature data. There was a linear density-temperature relationship at the site and salinity varied by only 0.14. L T is the root mean square (RMS) of the displacement d needed to adiabatically reorder a profile containing density inversions due to overturns, , where the ⟨⟩ indicates the mean displacement for each time step.
We estimated the diapycnal diffusivity as K ρ ¼ 0:09L 2S E , following Ivey et al. (2018), in which evaluation of the method was performed with data from the BUBS site. We estimated the Ellison length where θ′ is the temperature fluctuation around the (time) mean value θ, from the 0.5 Hz SBE56 data. To separate the internal waves from the turbulence, whilst also considering the full range of turbulent wavelengths, θ ′2 1=2 was determined from the time-frequency decomposition of the temperature variance (i.e., wavelet transforms) (Cimatoribus et al., 2014). The total temperature variance was determined by integrating from the local minimum buoyancy period T buoy (determined using 1 minute averaged data in the 14 minutes surrounding the wavelet estimate) to the Nyquist frequency (1 Hz). T buoy ranged from 3-16 minutes, with a median value of 4 minutes. The gradients dθ=dz and S were estimated as the time-mean over T buoy . Due to limited instrument resolution we removed estimates of L E when dθ=dz 0.002 o Cm −1 and S 0.004 s −1 .
Last, we estimated the vertical heat flux J q ¼ ρC p K ρ dθ=dz using the measured, average density ρ, and heat capacity C p , where downward J q is positive.

Site Conditions
The through-water column observations at PIL100 demonstrated the highly variable and energetic conditions throughout the 10-day experiment (day 95 to 105 of 2012). The surface tide amplitude increased to a maximum in the middle of the observations at day 100 ( Figure 2(a)). The background alongshore current (low-pass filtered at 30 h) was always towards the southwest (Figure 2(b)), but the cross-shore current varied in direction (Figure 2(c)). At the beginning of the record, the alongshore current was moderate, vertically sheared and accompanied by a predominantly onshore velocity and relatively strong and depth-uniform stratification (low-pass filtered at 30 h). The alongshore current significantly weakened (to 0.1 ms −1 ) from day 97 to 99 and the cross-shore current transitioned to predominantly offshore, except for the near-surface. After day 101, the alongshore current strengthened but had less shear, and the crossshore current became onshore again, coinciding with a significant decrease in stratification strength (Figure 2(d)). The mode-1 linear phase speed c 1 decreased from 0.4 to 0.3 ms −1 as the background stratification weakened (Figure 2(e)).
Typically, a cold bore with onshore, near-bed baroclinic currents arrived every 12 h, with the leading front accompanied by a packet of 5-10 large amplitude NLIW of elevation (Figures 2(e) and 2(g)). Similar features have been termed solibores, having both bore-like and wave-like features (Henyey & Hoering, 1997). The cold-water bore typically persisted for approximately 6 h before the water warmed, as the pycnocline was depressed and the near-bed baroclinic current reversed to offshore (Figure 2(g)). This "relaxation" phase was either gradual (Figure 2(g) day 103.1) or characterized by a sudden warmfront accompanied by waves of depression (Figure 2(g), day 102.575), which have previously been termed canonical and non-canonical bores, respectively (Walter et al., 2012). The bores were not present from day 97.5 to 100 (except for one on day 99) when the alongshore current significantly weakened (from 0.2 to 0.1 ms −1 ) and the cross-shore current transitioned to predominantly offshore, except for the near-surface. We observed the largest amplitude NLIW of elevation after day 101.5 when the background stratification was weakest (Figure 2(e)).
The periods between the NLIW ranged from 6 to 23 minutes. Herein we will focus on the elevation waves associated with cold water bore arrival as they have the largest amplitude and influence the near-bed dynamics most dramatically.

Characterizing the Instabilities of NLIW of Elevation
During the 10-day observation period, we captured a range of shapes of NLIW of elevation with variation in position and type of instability. We selected four examples for detailed analysis that demonstrate the range of variation (Figure 3).

Example1: Solitary-like wave with trapped core
The first example was the leading wave in a packet and was relatively symmetric (Figures 3(a), 3(b), 3(c)). It was characterized by a clockwise-recirculating core (u > c nl close to the seabed) with the center 18 mAB and biased towards the leading-face of the wave. The center of the recirculation corresponded with a maximum backscatter. Both Moum, Klymak, et al. (2007) and Scotti and Pineda (2004) observed similar rotating trapped cores in NLIW of elevation. We observed small to medium overturns in the core of the wave and ; positive onshore); d) background temperature (°C); e) Mode-1 buoyancy amplitude from modal fitting to observed temperature data (black) and mode-1 linear phase speed c 1 (red). f) Zoom in on 24 h: Surface tide onshore (black) and alongshore (red) velocity and mode-1 linear phase speed c 1 (dashed cyan lines); g) 1 minute averaged contoured temperature (°C).

Geophysical Research Letters
on the trailing-face of the wave coincident with strained isotherms. NLIW of depression with elevated turbulence in the core have also been observed; driven by a convective instability where u > c nl (Shroyer et al., 2010a;Zhang & Alford, 2015). Example 1 has similar properties to these depression waves, although it is likely that shear at the sea-bed also plays a role (Moum, Klymak, et al., 2007). The temperature structure in the wave core indicated that the instabilities were three-dimensional as regions of warmer water suddenly appeared within the core.

Example 2: Bolus with lobe and cleft instability
This asymmetric wave, had a steep leading-face and more gradual sloping trailing-face and had many features of an internal bolus, similar to those produced by numerical models of collapsing and surging NLIW of depression (Figures 3(d), 3(e), 3(f ),) (Aghsaee et al., 2010). It was the first wave in the train. Less dense water was entrained under the leading-face of the wave, leading to very large overturns similar to the lobe and cleft instabilities observed at the nose of gravity currents (Simpson & Britter, 1979). Wallace and Wilkinson (1988) observed a similar feature in their laboratory experiments of shoaling waves of elevation, while Xie et al. (2017) observed similar NLIW generated by an estuarine gravity current. The peak wave velocity was located towards the front-face of the wave at 16-20 mAB and exceeded c nl by 0.2 ms −1 , suggesting a convective instability may have caused overturning towards the bed where the velocities were slower. A recirculating region was centered ∼30 mAB in the apex of the wave. Sub-buoyancy frequency "oscillations" were observed at the trailing-face of the wave where the isopycnals were strained, similar to the "finger-like" motions observed by van Haren and Gostiaux (2010) on the trailing-face of the internal tide, where S ∼ N. We observed the largest backscatter at the bottom extent of the observations and towards the trailing-face of the wave.

Example 3: Backward breaking wave
In this asymmetric wave, the trailing-face of the wave had steepened and was overturning similar to the laboratory experiments of Wallace and Wilkinson (1988), where a solitary wave of elevation steepened and broke (Figures 3(g), 3(h), 3(i)). It was the first wave in the train. Plunging warmer water was entrained up into the core of the wave in a similar fashion to the front-face of plunging surface waves. In the reference frame of the wave, we observed the largest onshore wave velocity close to the seabed (exceeding c nl by 0.15 ms −1 ) and then diminishing to 0 ms −1 at 20 mAB and then becoming offshore above this. This led to the rear-face of the wave steepening until it broke due to a convective instability. Higher temporal resolution velocity observations would allow more detailed analysis. The wave form was similar to the result of a collapsing wave of depression as modeled by Nakayama et al. (2019) and Arthur and Fringer (2014). The largest overturns were found at the trailing-face of this wave. Ri > 1 in the leading-edge of the wave, but Ri 0.25 in the back-face of the wave (not shown). An isolated region of low backscatter water was located at a height of 18-28 mAB in the center of the wave.

Example 4: Shear instability on trailing-face
This wave immediately followed the backward breaking wave of Example 3 (Figures 3( j), 3(k), 3(l)). The wave was asymmetrical, with an extended trailing-face where the isotherms remained elevated well beyond their initial position. The minimum wave velocity was located at 32.5 mAB. On the trailing-face of the wave, braids or fingers were present with the first two forming overturns in the direction of wave propagation, as expected for a Kelvin-Helmholtz billow forming from a shear instability. Geyer et al. (2010) and van Haren and Gostiaux (2012) observed similar features in the deep ocean and an estuary pycnocline, respectively. Numerical and laboratory experiments of NLIW of depression have demonstrated that enhanced shear at the apex of the wave can result in very small Ri, that if sustained long enough can lead to the growth of a shear instability resulting in Kelvin-Helmholtz billows on the trailing-face of the wave (e.g., Barad & Fringer, 2010). We did not observe N at the wave crest as it was beyond the height of our observations, but typically N would be smaller at the crest due to vertical straining, and S was enhanced at the crest, thus likely leading to small Ri. We observed a series of medium overturns in the billow region. In our observations, the peak wave velocity occurred at 5 minutes and a height of 25-30 mAB.
We classified the dominant form (sometimes more than one feature was present) of a total of 28 NLIW of elevation (the first two in each packet) during the ten-day observation period. The solitary-like wave was the most common (13), followed by backward breaking waves (7), bolus with lobe and cleft instability (5) and shear-instability on trailing-face (3).

Diapycnal Diffusivty
The diapycnal diffusivity K ρ varied by six orders of magnitude (10 −7 to 10 −1 m 2 s −1 ) throughout the 10-day observation period, ranging from near-molecular diffusion up to values at ocean mixing "hot spots" (e.g., Waterhouse et al., 2014). To evaluate the influence of NLIW on the near-bottom mixing, we averaged log 10 (K ρ ) over the measurement depth (8-20 mAB) (K ρ ) (Figure 4(b)). For the majority of the NLIW of elevation events the first two waves in each packet had the highest K ρ and overall K ρ followed a log-normal distribution with a mean of 0.5x10 −4 m 2 s −1 (Figure 4(a)). We identified two other periods where K ρ >10 −3 , which were associated with large waves of depression (Figure 4(b)). K ρ was on average five times larger during the 13 NLIW of elevation events compared to the 2 h period just prior (Figure 4(d)). During the NLIW events, 11% of the individual K ρ estimates were greater than 1x10 −3 m 2 s −1 , compared with just 4% in the 2 h prior.
We observed the highest K ρ in the unstable region/s of the wave coinciding with the largest Thorpe overturns, agreeing with the modeling results of Arthur and Fringer (2014). For example, the "backward breaking wave" (Example 3) recorded the largest K ρ in the trailing-face of the wave (Figure 4(c)). K ρ decreased to smaller, but not background levels, at the front-face of the following wave. The second wave also produced the largest K ρ in the trailing-face of the wave, coinciding with the location of the shear-instabilities on the trailing edge (Figures 3(d) and 4(c)). Moum, Klymak, et al. (2007) observed elevated dissipation at 1 mAB in the trailing-face of what they identified as a solitary-like wave, which they attributed to the contribution from the near-bed stress, but their density measurements were restricted to one height and therefore wave-breaking may not have been detected. For the "solitary-like wave" (Example 1) the largest K ρ occurred in the wave core, whilst for the "Bolus lobe and cleft instability wave" (Example 2) the largest K ρ was located on the leading-face of the wave (not shown). For 16 clearly characterized waves, we did not observe a relationship between the type of wave breaking and the magnitude of the median K ρ for the wave.
The elevated K ρ generally remained for the initial phase of the bore, coincident with the largest decrease in temperature and the largest amplitude waves (2-3 h) before returning to background levels as the water temperature became approximately constant (Figure 4(c)). Within the wave packets K ρ oscillated, dropping by an order of magnitude between each wave, but not to background levels. The turbulence must have been continuously produced throughout the passage of the bore as the timescale of elevated K ρ was orders of magnitude larger (O(2 h)) than the timescale for the decay of unforced, stratified turbulence 0.6(2π/N) ∼ 5 minutes (Lienhard & Van Atta, 1990).
Our observations suggest that the elevated mixing was largely moving onshore with the bore instead of decaying locally. It is likely the bore continuously produced turbulence as it traveled over some distance and the breaking waves evolved, before eventually fully dissipating. We observed the wave packets 8.1 km inshore in ∼75 m water depth (N4I mooring). However, the observed peak vertical wave velocity was diminished to 1/3 rd of the value at PIL100/BUBS, indicating dissipation was exceeding shoaling (Sandstrom & Oakey, 1995), and the waves had on-average approximately halved their propagation speed (estimated from the time of wave propagation between the two sites). This indicated substantial energy was dissipated over this distance. Further studies are required to determine the evolution and fate of breaking waves.

Heat Flux
The heat flux J q was strongly enhanced during the wave trains associated with the bores, similar to the enhancement by waves of depression observed by Shroyer et al. (2010b). On average J q was five times larger during the NLIW events (2 h periods) compared with the 2 h prior and J q > 500 Wm −2 for 18% of the time during the events, compared with just 1.5% in the 2 h period prior (Figure 4(e)). The enhanced K ρ during NLIW events was largely responsible for the increase in J q . We determined the average heat flux within a wave J q,wave (from log 10 (J q )) for 16 waves of elevation with well-defined characteristics. These waves were all very energetic with the wave Reynolds number Re w ¼ ac 1 /ν ¼ O(10 7 ), where ν is the molecular viscosity. A strong positive relationship was found between J q,wave and the wave steepness ka, which is commonly used to characterize interfacial waves in deep water (Figure 4(f)). For an interfacial wave, ka is approximately the maximum wave Froude number and (ka) 2 is inversely proportional to the minimum Richardson number of a wave (e.g., Fringer & Street, 2003). The trend in our observations supports the relationship developed by Thorpe (2018), demonstrating an increased transfer of energy to turbulence with larger ka.

Conclusions
The observations demonstrated the varied characteristics of the NLIW of elevation at one site and even within one wave packet over a 10-day period. Many of the waves had evidence of a closed core that could 10.1029/2020GL089591

Geophysical Research Letters
lead to the horizontal transport of trapped material over large distances. The NLIW of elevation were all associated with elevated diapycnal mixing and heat flux. 20% of our observations had K ρ > 1x10 −3 m 2 s −1 , largely attributable to NLIW. We found a positive relationship between the average heat flux within a wave and the wave steepness ka. Although there are too few data points to confidently model the process parametrically, the relationship suggests that for energetic waves we can use bulk properties of NLIW of elevation (from observations or numerical models) to predict vertical fluxes of heat and other scalars such as dissolved nutrients.
As NLIW propagate 10s kilometers they can both contribute to net onshore transport and drive vertical mixing of important constituents such as heat, nutrients and sediment. It is important for future studies to capture the propagation and dissipation of these waves so that the cross-shore impact is better defined. In many regions, NLIW of elevation dynamics are likely to dominate transport mechanisms, thereby capturing these extremely variable waves and their resultant turbulent mixing in regional estimates of heat budgets, sediment transport and primary production will be crucial. Our study has highlighted the relatively small length and time scales that need to be resolved (both in observations and numerical models) in order to capture the NLIW dynamics, with large variations in vertical mixing occurring over timescales much less than the wave period. Alternatively, wave steepness may be a useful bulk property for future parameterization of waveaveraged heat flux.