Tidally Forced Lee Waves Drive Turbulent Mixing Along the Arctic Ocean Margins

In the Arctic Ocean, limited measurements indicate that the strongest mixing below the atmospherically forced surface mixed layer occurs where tidal currents are strong. However, mechanisms of energy conversion from tides to turbulence and the overall contribution of tidally driven mixing to Arctic Ocean state are poorly understood. We present measurements from the shelf north of Svalbard that show abrupt isopycnal vertical displacements of 10–50 m and intense dissipation associated with cross‐isobath diurnal tidal currents of ∼0.15 m s−1. Energy from the barotropic tide accumulated in a trapped baroclinic lee wave during maximum downslope flow and was released around slack water. During a 6‐hr turbulent event, high‐frequency internal waves were present, the full 300‐m depth water column became turbulent, dissipation rates increased by a factor of 100, and turbulent heat flux averaged 15 W m−2 compared with the background rate of 1 W m−2.


Introduction
Much of the interior of the Arctic Ocean is quiescent (Fer, 2009;Lincoln et al., 2016), insulating the surface from heat imported into the Arctic Ocean by Atlantic Water intruding at intermediate depths (Carmack et al., 2015). In the Arctic basins, diapycnal heat fluxes out of this main oceanic heat source are weak (<1 W m −2 ) and typically driven by double diffusive convection (Polyakov et al., 2019;Shibley et al., 2017). However, in the critical region around Svalbard where warm water enters the Arctic Ocean, tides appear to play a key role in transporting this heat toward the surface. Mixing hotspots have been identified over regions of steep topography (Fer et al., 2010;Meyer et al., 2017;Padman & Dillon, 1991;Rippeth et al., 2015), where enhanced levels of turbulence cause diapycnal heat fluxes as large as 50 W m −2 . The regions of energetic turbulent dissipation correspond with areas of conversion of barotropic tidal energy to baroclinic waves and mixing Padman & Dillon, 1991;Rippeth et al., 2015).
Planetary rotation constrains the conversion of barotropic tidal energy to freely propagating linear internal tides. For almost all of the major tidal constituents, most of the Arctic Ocean is poleward of their critical latitude, ϕ c , the latitude at which the tidal frequency matches the local inertial period. In consequence, the linear radiation pathway that the energy would normally follow from tides to dissipation, equatorward of ϕ c , is not permitted (Falahat & Nycander, 2014;Vlasenko et al., 2003). The hypothesized energy pathway at high latitudes is, instead, that the tidal flow over sloping topography induces a lee wave, which is trapped and becomes unsteady as the flow becomes supercritical Rippeth et al., 2017;Vlasenko et al., 2003), that is, Froude number Fr = u/c > 1, where u is the flow speed and c is the internal wave speed. The energy is dissipated locally and by breaking nonlinear internal waves (NLIWs) evolving with the changing phase of the tide (Rippeth et al., 2017).
The generation of NLIWs by the lee wave mechanism (e.g., Jackson et al., 2012;Maxworthy, 1979) occurs when flow over sloping seabed vertically displaces the pycnocline. Over the upper slope and shelf edge, downslope (ebb) flow depresses the pycnocline, developing an initial lee wave (Sharples et al., 2007). The spatial scale of the generated wave is of the same order as the horizontal scale of the topographic feature, typically 5-30 km over the shelf edge. When Fr = 1, the disturbance remains stationary, accumulates energy, and attains phase speed equal to the flow speed. As the flow slackens and Fr falls below 1, the wave propagates toward the generation region. Waves propagating on-shelf steepen as they radiate into shallower water and further evolve and shorten (typically to about 1 km) as the tide turns, similar to internal hydraulic jumps that develop over steep slopes and deep-ocean ridges after the maximum downslope flow and propagate upslope as the flow relaxes (e.g., Legg & Klymak, 2008;van Haren, 2019). The large vertical displacements of the density surfaces alter the sea surface roughness (Jackson et al., 2012), allowing the surface signature of NLIWs to be detected in remote sensing imagery (e.g., Alford et al., 2015;Jackson et al., 2012;Kozlov et al., 2017).
At midlatitudes, intense turbulence close to topography has been linked to tidal lee waves in deep sea environments (Alford et al., 2015;Klymak et al., 2008;Li & Farmer, 2011;, continental shelf breaks (Sharples et al., 2007) and in fjords (Cummins et al., 2003;Klymak & Gregg, 2003), and to lee waves associated with geostrophic flows (Thorpe et al., 2018). In the Arctic, The inset shows the probability distribution function (PDF) for Froude number calculated using |u x − max | and the first two baroclinic (evanescent) wave speeds for the K 1 wave frequency and the MIMOC August stratification at the locations with |u x − max |>0.2 m s −1 (inside red contours) and total water depth greater than 50 m.
high-frequency NLIWs have been observed in the eastern Arctic (Kurkina & Talipova, 2011;Morozov et al., 2017), including the northern flanks of the Yermak Plateau (Czipott et al., 1991;Padman & Dillon, 1991) (YP in Figure 1). The recent decline in Arctic seasonal sea ice has facilitated the mapping of the surface signatures of NLIWs from space using high-resolution synthetic aperture radar (SAR), which shows the widespread occurrence of these features in the eastern Arctic Basin . Observations linking tides, NLIWs, and turbulence in the Arctic, however, are lacking despite their hypothesized importance in the magnitude and distribution of mixing.
Here we report detailed concurrent measurements of current structure, turbulence, and mixing in the Arctic Ocean that show intense dissipation of lee waves driven by the tidally induced stratified flow over sloping topography. We focus on the diurnal K 1 tidal current because it is directed approximately perpendicular to the shelf break at the study site, giving favorable forcing conditions for energy conversion. The aim of this paper is to identify the pathway of tidal energy to turbulence at a mixing hotspot in the Arctic, and to highlight the role of tidally-forced mixing along the Arctic margins.

Data and Methods
Data were collected from R.V. Kristine Bonnevie during 27 June to 10 July 2018 using a conductivity-temperature-depth (CTD) and lowered acoustic Doppler current profiler (LADCP) system, a vertical microstructure profiler (VMP), a shipboard acoustic Doppler current profiler (SADCP), and a scientific echosounder (Simrad EK80). Observations presented here are from a 100-km-long transect across the shelf break at 18°E and from a 24-hr process station near the 300-m isobath. The CTD and LADCP profiles were obtained concurrently. The VMP was deployed immediately after the CTD/LADCP profiling. During the process station, a set of 4-5 VMP profiles (10-15 min between profiles) were taken, while the ship drifted about 0.5-1 km, after which there was a gap when the ship was repositioned. In addition, CTD/LADCP profiles were collected every 4 hr. The cruise data are available from Fer et al. (2020).
The CTD profiles were acquired using a Sea-Bird Scientific, SBE 911plus system. Accuracies of the pressure, temperature, and salinity sensors are ±0.5 dbar, ±2 × 10 −3°C , and ±3 × 10 −3 , respectively. The SADCP (150-kHz Teledyne RDI Ocean Surveyor) velocity profiles (narrowband mode, 8-m vertical bins) were processed into 3-min averages using the University of Hawaii postprocessing software. The LADCP data (from two 300-kHz Teledyne RDI Workhorse ADCPs attached to the CTD frame) were processed using the LDEO software version IX-13 based on Visbeck (2002), using the constraints from ship navigation, bottom tracking, and SADCP. Error in horizontal velocity measurements is less than 0.03 m s −1 (Thurnherr, 2010).
Microstructure measurements were made using the loosely tethered VMP2000 (Rockland Scientific International). The fall speed was 0.6 m s −1 . The CTD data, averaged to 0.5 dbar, were taken from the pumped SBE-CT sensor. Microscale velocity shear from shear probes was processed into 2-dbar averaged dissipation rate estimates assuming isotropic turbulence, ϵ ¼ 7:5ν ∂u ′ =∂z ð Þ 2 , where ν is the viscosity (∼ 1.6 × 10 −6 m 2 s −1 for T = 3°C). The small-scale shear variance (the term with the overline) was obtained by integrating the wavenumber spectrum of shear from half-overlapping 2-s segments in a wavenumber range that is relatively unaffected by noise and corrected for the unresolved variance using the empirical model from Nasmyth (1970).
Acoustic data were collected using the EK80, operated in the narrowband mode at discrete frequencies. The volume backscatter data were reprocessed to remove noise and correct for spatial offsets and time delays between the transducers using the Large Scale Survey System (Korneliussen et al., 2016). Each of the EK80 frequencies displayed discrete targets across the water column and a continuous volume backscatter layer located midwater that followed the pycnocline variability. We used the acoustic image at 38 kHz, which best matched the isopycnal displacement pattern observed by the VMP profiles.
The squared buoyancy frequency was computed using N 2 = −(g/ρ 0 )∂σ θ /∂z, where g = 9.8 m s −2 , ρ 0 = 1,027 kg m −3 , and σ θ is the potential density anomaly. Turbulent heat flux was obtained from where T is the background temperature profile, C p is the specific heat capacity of seawater, and the diapycnal eddy diffusivity was estimated from K ρ = 0.2ϵ/N 2 (Gregg et al., 2018;Osborn, 1980). Richardson number, Ri = N 2 /S 2 , at the time of a dissipation profile, was calculated using the 3-min averaged velocity 10.1029/2020GL088083

Geophysical Research Letters
profile from the SADCP closest in time to the VMP cast. The vertical gradients of density in N 2 and horizontal currents (u,v) in shear-squared, S 2 = (∂u/∂z) 2 + (∂v/∂z) 2 , were calculated at an 8-m vertical scale.
Surface manifestations of NLIWs were derived from optical Landsat-8 and Sentinel-1 A/B SAR images (Kozlov et al., 2015). We show the identified wave crests but not the actual images. We used Level 1 near-infrared Band 5 and panchromatic Band 8 images from Landsat-8, with spatial resolutions of 30 and 15 m, respectively, to detect NLIW signatures on 3 July 2018. We used Sentinel-1A/B Level-1 mediumresolution Ground Range Detected images taken in Extra Wide swath mode with a spatial resolution of 90 m to identify NLIWs during 4-5 July 2018.
We obtained tidal currents from Arc5km2018 (Erofeeva & Egbert, 2020), a barotropic inverse tidal model on a 5-km grid. At each grid point, we predicted the total tidal current over a 30-day window using all constituents. We rotated the tidal currents to local along-and cross-isobath components and calculated statistics of the total tidal current speed and cross-isobath component.
The Froude number was computed as Fr = u/c, where u is the flow speed and c is the internal wave speed. We obtained the phase speeds of the subinertial evanescent response for the first (c 1 ) and second (c 2 ) baroclinic modes by solving the flat-bottom eigenvalue equation using the observed N(z) (e.g., Chapter 5, Phillips, 1977) and a relation above the critical latitude (Equation 2.10 in , with local f and K 1 wave frequency. For the time-averaged N 2 profile at the process station, c 1 = 0.14 m s −1 and c 2 = 0.06 m s −1 . Fits of the first two baroclinic modal shapes to the cross-isobath baroclinic (depth-average removed) currents and the pressure anomaly explained most of the variance (section 3). The second mode dominated the currents during the turbulent period.
To calculate Fr from tidal currents over the Arctic Ocean, we used the 30-day maximum cross-isobath tidal current, |u x − max |, and the evanescent wave phase speed for the K 1 frequency. We obtained the phase speeds of the first two baroclinic modes using the August stratification from the 0.5°resolution monthly isopycnal mixed-layer ocean climatology (MIMOC; Schmidtko et al., 2013), linearly interpolated to the Arc5km2018 grid.

Results
In July 2018, sea ice melt resulted in strong ocean stratification north of Svalbard, with a buoyancy period of 20 min in the pycnocline at 50 m. Profiles collected along the transect showed a depression of isopycnals of about 50 m near the 300-m isobath (Figure 2a). The transect started from the shelf on 3 July at 12:30 (all times are given in coordinated universal time UTC) and progressed north to deeper water, ending 4 July 02:30. The water column was highly turbulent near the 300-m isobath (57 km), where ϵ reached 10 −6 W kg −1 .
Current measurements were converted to along-isobath, u a , and cross-isobath, u x , components using the orientation of the 500-m isobath (60°counterclockwise from east; see Figure 2b). The large values of depth-averaged u a between 40 and 60 km mark the core of the warm boundary current (Figure 2c), whereas the upslope (negative) values of u x are primarily from the diurnal tide, as will be shown later. Fr attains near-critical (Mode 1) and supercritical (Mode 2) values between 40 and 60 km along the transect. We detected signatures of NLIWs from the SAR image from 3 July 16:45, approximately when the station at 50 km with Fr > 1 was taken (Figure 2b). The ship returned to 57 km after the transect was completed to carry out detailed time-series measurements (the process station).
Tidal ellipses from Arc5km2018 are approximately normal to the isobaths, particularly for the diurnal constituent K 1 , giving favorable forcing conditions for energy conversion ( Figure 3a). As we found for the dayof the transect, SAR images near the start and end of the process station show the presence of NLIWs (Figure 3b). The wave crests were densely spaced in the portion of the transect where the flow was observed to be critical or supercritical (see also Figure 2b). The location, however, is not unusually energetic, and such NLIWs are frequently detected in SAR images from the Arctic shelves .

Geophysical Research Letters
The process station was occupied for 24 hr starting from 07:00 on 4 July, about 4 hr after the slope transect was completed. Repeat profiling revealed the strong variability of isopycnal depths and the highly intermittent nature of turbulence (Figure 4a). Isopycnal depth variations were consistent with a dominant diurnal cycle but with substantial changes on shorter time scales. An abrupt depression of isopycnals was observed near 17:00, similar to the depression of the isopycnal surfaces observed at 17:00 the previous day at 300-m depth (Figure 2a). The event was highly turbulent, with ϵ ∼ 10 −6 W kg −1 across the entire water column, similar to the profile taken 24 hr earlier (Figure 2a).
The vertical structure of the temporal variability in u x was complex but showed the overall diurnal cycle (Figure 4b). The turbulence event coincided with the slackening and reversal of u x . The downslope current reached a subsurface peak of 0.3 m s −1 approximately 1 hr before the isopycnal dip, and the current direction reversed in the short duration of the abrupt isopycnal perturbation (Figure 4b). The acoustic trace of the internal wave observed in the midwater column with the EK80 (Figure 4c) was better temporally resolved than the other profiling measurements. In particular, the isopycnal dip of about 50 m near 17:00 at the time of the midwater flow reversal was more abrupt than could be resolved by the two consecutive microstructure profiles at that time. However, the EK80 data quality is compromised by the ship's drift and frequent repositioning by using propellers, and a detailed quantitative (e.g., spectral) or qualitative description is not possible.
The depth-averaged u x at the process station was dominated by the diurnal tide with amplitude 0.15 m s −1 , while u a was predominantly semidiurnal (Figure 5a). The observations are in reasonable agreement with the predicted currents using Arc5km2018, with amplitudes within 5% for diurnals (1.5-hr phase offset) and 20% for semidiurnals (not shown). During maximum ebb, the second baroclinic mode explained 90% of the cross-isobath current structure and remained comparable to the first baroclinic mode at slack water ( Figure 5b). The corresponding vertical structure can be seen in Figure 4b between 12:00 and 17:00. Coincident with the Mode 2 wave, Fr 2 was supercritical (Figure 5c), peaking during the maximum offslope flow accompanied by energetic turbulence (Figure 5d).
Dissipation rate averaged between 50 and 250 m increased by 2 orders of magnitude starting after the ebb maximum (after 12:00) until it decayed after the peak of the flood at 21:00 (Figure 5d). We chose this averaging depth range to minimize contributions from turbulence in the surface and bottom boundary layers; however, the pattern for the full-depth averaged dissipation was similar. The onset of the turbulent event followed the supercritical conditions for Mode 2. At this time, the shear associated with the subsurface maximum in u x (Figure 4b) reduced Ri. During the turbulent event, Ri < 1 in 40-90% of the water column, and Ri < 0.25 in 10-30%, continuously for a period of 6 hr, despite the relatively coarse vertical resolution of the Ri estimates.

Discussion
Our observations are from a site north of Svalbard, where diurnal tides dominate cross-isobath barotropic currents. The region experiences strong, tide-modulated production of NLIWs and turbulence which increases the time-averaged heat flux four times above the background values of 1 W m −2 .
Conversion of barotropic tides to baroclinic waves will be strongest when the cross-isobath component of barotropic tidal currents is largest. Energetic cross-isobath tidal currents are common along the Arctic Ocean margins (Figure 1). Near-critical or supercritical Froude numbers in these regions (Figure 1, inset) will allow trapping and unsteady evolution of lee waves and NLIWs. These calculations from the climatology likely underestimate the critical conditions (e.g., Fr inferred from the August climatology at the measurement location is 30% less than the value from observations). SAR images indeed show a widespread occurrence of NLIWs in the eastern Arctic Basin .
We estimate the relative contribution of double-diffusive and the tidally driven heat flux using the approximate surface area in the deep Canada (3 × 10 12 m 2 ) and Eurasian (2 × 10 12 m 2 ) basins and the areas with |u x − max |>0.2 m s −1 (inside red contours above 74°N in Figure 1, 0.5 × 10 12 m 2 ). The fraction of area covered by staircases and average double-diffusive heat fluxes are summarized in Shibley et al. (2017). Using 0.3 W m −2 in 80% of the Canada Basin and 1 W m −2 in 25% of the Eurasian Basin, we obtain 1.2 × 10 12 W. If 15 W m −2 (average over the turbulent event) occurs at 50% of the |u x − max | area (approximate probability of Fr 2 > 1 using the PDF from Figure 1), 25% of the time (based on our observations), or a continuous 4 W m −2 (our station average) over the same area, the Arctic-wide contribution of the tides to the diapcynal heat flux is comparable to that due to double diffusion. This coarse estimate suggests that, despite the limited spatial extent and temporal occurrence of turbulent events, the contribution of tidal-driven mixing is substantial. Bulk estimates of heat loss from the Atlantic Water during its transformation from the Yermak Plateau to the Lomonosov Ridge (e.g., see Ivanov & Timokhov, 2019) suggest that of about 16 W m −2 total heat loss, roughly 40% is released to the atmosphere, 40% warms the deeper water, and 20% is lost laterally. This implies the tidally driven mixing is an important component of this budget.
Along the Arctic Ocean margins, the pathway for the energy from the tide to turbulence is nonlinear and dominated by breaking unsteady lee waves and critical flow. Here we have shown the enhanced turbulence associated with supercritical flow of a trapped tidal lee wave. The vertical scale of sharp vertical undulations (10-50 m), duration (hours) and amplitude of bursts of dissipation (10 −7 to 10 −6 W kg −1 ), their vertical extent (100-400 m), and timing relative to the tidal currents (at the end of the ebb tide) are similar to 10.1029/2020GL088083

Geophysical Research Letters
observations from lower latitudes. In contrast to the lower latitudes where linear internal tidal waves disintegrate into packets of short NLIWs, linear response above critical latitudes is evanescent but nonlinear response has properties inherent to lee waves (at wavelengths comparable to the length scales of the bottom topography), independent of the rotation (Vlasenko et al., 2003).
The Arctic sea ice cover variability is sensitive to small changes in ocean-ice heat flux (Carmack et al., 2015). Sea ice and stratification changes, e.g., reduced sea ice, weakening of the halocline, and shoaling of the intermediate-depth Atlantic Water layer observed in the eastern Eurasian Basin (Polyakov et al., 2017(Polyakov et al., , 2020 will affect the generation, radiation, and dissipation of the internal tide. The combination of these changes together with climate-driven variability of the warm boundary current will alter both the mean flow and the internal wave speed, c. The polar regions are unique with relatively small values of c; hence, even modest tidal currents can lead to near-critical states. For constant stratification, c scales as NH, where H is the total water depth. In the changing Arctic with weaker stratification, we expect smaller c, hence larger values of Fr. In turn, the transition toward near-critical and supercritical states will be changed, altering the rate of tidal energy conversion, particularly through the nonlinear processes identified in this research. Understanding the pathway for the energy from tides to turbulence in the Arctic Ocean, the magnitude and distribution of the associated ocean mixing rates, and the role of feedbacks between mixing rates, stratification, sea ice, and the tide is key to predicting the fate of the Atlantic water in the Arctic and the evolution of the Arctic Ocean in a warming world.