Mapping turbulent diffusivity associated with oceanic internal lee waves offshore Costa Rica

Breaking internal waves play a primary role in maintaining the meridional overturning circulation. Oceanic lee waves are known to be a significant contributor to diapycnal mixing associated with internal wave dissipation, but direct measurement is difficult with standard oceanographic sampling methods due to the limited spatial extent of standing lee waves. Here, we present an analysis of oceanic internal lee waves observed offshore eastern Costa Rica using seismic imaging and estimate the turbulent diffusivity via a new seismic slope spectrum method that extracts diffusivities directly from seismic images, using tracked reflections only to scale diffusivity values. The result provides estimates of turbulent diffusivities throughout the water column at scales of a few hundred meters laterally and 10 m vertically. Synthetic tests demonstrate the method’s ability to resolve turbulent structures and reproduce accurate diffusivities. A turbulence map of our seismic section in the western Caribbean shows elevated turbulent diffusivities near rough seafloor topography as well as in the mid-water column where observed lee wave propagation terminates. Mid-water column hotspots of turbulent diffusivity show levels 5 times higher than surrounding waters and 50 times greater than typical open-ocean diffusivities. This site has steady currents that make it an exceptionally accessible laboratory for the study of lee-wave generation, propagation, and decay.


Introduction
Understanding the spatial distribution of diapycnal mixing associated with turbulent ocean phenomena is a fundamental, yet incomplete, component of understanding the meridional overturning circulation (MOC).Measured open-ocean diffusivities, on average, fall an order of magnitude short of values required to account for the upwelling necessary to balance deep water formation in polar regions (Ledwell et al., 1993).This suggests that moving water masses and smaller "hotspot" regions of elevated mixing likely account for the difference between measured open-ocean diffusivities and those required to balance the MOC (Munk and Wunsch, 1998;Stewart, 2005).Work in this area has shown tides and deep currents interacting with topography to be an effective source of elevated turbulent diffusivity near seamounts (Lueck and Mudge, 1997;Toole et al., 1997), mid-ocean ridges (St. Laurent, 2002), island chains (Rudnick et al., 2003), and rough bathymetry (Polzin et al., 1997;Naveira Garabato et al., 2004).
While internal tides are thought to account for the majority of ocean mixing (Garrett and Kunze, 2007;Munk and Wunsch, 1998), recent work has shown that lee waves dissipate significant amounts of energy in the global ocean.Estimates by Nikurashin and Ferrari (2011) and Scott et al. (2011) predict 0.2-0.4TW of energy is dissipated due to lee waves, especially in the Southern Ocean.These figures represent approximately 10-20 % of the 2 TW required to maintain the observed MOC (Munk and Wunsch, 1998).Here we show evidence for increased turbulent diffusivity associated with Published by Copernicus Publications on behalf of the European Geosciences Union.lee waves offshore Costa Rica using methods from seismic oceanography (herein SO).In particular, we observe lee wave breakdown and associated higher levels of turbulence above the lee wave generation site and subsequent propagation.
SO is well suited for estimating regional heterogeneities in turbulence (Holbrook et al., 2003(Holbrook et al., , 2013;;Sheen et al., 2009) as seismic images yield information across large 2-D and 3-D swaths of the oceanic water column rapidly and in high horizontal resolution.Commonly, turbulent dissipation is measured in vertical casts that have great vertical resolution but are spaced many kilometers apart (Polzin et al., 1997;St. Laurent and Thurnherr, 2007).Sampling in vertical casts has difficulty characterizing extent and shape of turbulent hotspots as turbulence is highly spatially variable and intermittent in time (Ivey et al., 2008).Images produced by SO, in contrast, typically have a horizontal sampling of 6.25 m, so that the method developed by Klymak and Moum (2007b) described in Sect.2.1 can be used to estimate turbulent diffusivity as described in Sect. 3 below.In this paper we use slope spectral methods to produce a map of turbulent diffusivity along a seismic transect in the western Caribbean that shows elevated diffusivity associated with lee wave structures.

Oceanic lee waves
Oceanic internal lee waves are formed by the interaction of a steady flow over varying seafloor topography when the advective length scale of the flow is much larger than the scale of the topographic variation (Thorpe, 2005).The result of such flow-topography interaction is a standing, near-vertical internal wave that is phase locked to the seafloor topography (Baines, 1995).Low mode lee waves, such as those we observe here, appear as vertical disruptions to isopycnal surfaces as a result of the disturbed flow over the raised bathymetry (Eakin et al., 2011;Klymak et al., 2010).
Lee waves are less well studied than near-inertial internal waves and internal tides.Studies such as Ocean Storms and HOME used frequency-specific time series methods to link measured internal wave fine structure to turbulent mixing processes (Rudnick et al., 2003).Such methods cannot be easily used in the case of lee waves, as they are not band limited in frequency in the same manner.Moreover, time series analyses of lee waves are problematic due to their fixed nature relative to the seabed topography, rendering the absolute frequency zero (Klymak et al., 2010).
Despite the relative lack of study, lee waves are thought to be an important driver for ocean mixing (Nash et al., 2007;Scott et al., 2011).In particular, there may be significant mixing along the Antarctic circumpolar current (ACC) in the Drake Passage, where microstructure data suggest mixing associated with lee waves is not limited to hydraulic boundary layers as radiated energy drives mixing in the stratified oceanic interior (St. Laurent et al., 2012).However, in previous studies, this association between upward radiation of energy and lee waves is inferred and not measured di-rectly.Studies in the Camarinal Sill in the Strait of Gibraltar (Farmer and Armi, 1988;Wesson and Gregg, 1994), the Bosphorus (Gregg, 2002), and the Knight Inlet (Farmer and Smith, 1980;Klymak and Gregg, 2004) show elevated levels of dissipative turbulence associated with oceanic lee waves.In these studies, high-frequency acoustic backscatter data imaged dissipative turbulence, but the images are limited to the upper 250 m due to acoustic attenuation.For large, deep features, lower frequency seismic imaging is the only technique we are aware of that can image oceanic lee waves.

Seismic observation of internal waves
Seismic oceanography uses marine seismic reflection methods to acoustically image the ocean interior (Holbrook et al., 2003).Seismic reflections are produced when a highamplitude, low-frequency (20-200 Hz) acoustic signal is partially reflected by density and sound speed boundaries created by oceanic fine structure, producing a return that is essentially a convolution of the seismic source wavelet with the reflectivity of the water column.Such reflections are sensitive to temperature changes as small as 0.03 • C (Nandi et al., 2004;Sallares et al., 2009).SO is unique among oceanographic techniques in its capability to rapidly image on O(10m) horizontally and vertically over large areas.
The product of SO imaging is a cross sectional view of thermohaline fine structure in the oceanic interior.Data are collected with great spatial redundancy to increase signal to noise and is sorted into common mid-point (CMP) gathers that represent one vertical profile in the cross section.Traces in each CMP gather are then processed and summed together to create one seismic trace that essentially represents a convolution of the acoustic source wavelet with the normal-incidence reflection coefficient, which is a function of the acoustic impedance at that point along the profile.To produce the final image, we typically shade the positive seismic returns one color and the negative seismic returns another (black and white are a common choice) since drawing a wiggle-trace seismogram at each location would blacken the image due to the high horizontal sampling.Additional seismic processing method details can be found in Yilmaz (1987), andRuddick et al. (2009) provide a review of seismic images as they relate to the water column.
Beyond producing snapshots of oceanic fine structure, SO can provide quantitative oceanographic information.Holbrook and Fer (2005) showed that seismic images produce internal wave spectra consistent with the Garrett-Munk 76 tow spectrum (Katz and Briscoe, 1979) and noted enhanced internal wave energies near the continental slope.The technique has since been used to estimate internal wave energies off the Iberian Peninsula (Krahmann et al., 2008) and Costa Rica (Eakin et al., 2011).
Similar spectral techniques can estimate turbulence from SO data by utilizing the isopycnal slope spectrum method of Klymak and Moum (2007b).Turbulent diffusivities can be estimated by matching spectral levels with the predicted spectral slope across the horizontal wavenumber range of tracked seismic reflections (Holbrook et al., 2013;Klymak and Moum, 2007b;Sheen et al., 2009).Our work here expands upon these advances to estimate turbulence even in regions where seismic reflections are difficult to track due to high levels of turbulence that has disrupted continuous fine structure.Eakin et al. (2011) used SO to produce the first seismic image of an oceanic lee wave.Their data were collected along a transect coincident with the one used in this study but with different acquisition parameters.Their analysis concluded that lee waves are formed by the interaction of the Panama-Columbia Gyre and the rough bathymetry of the region.Additionally, Eakin et al. show the lee wave structure to have anomalously high internal wave energy over the horizontal wavenumber range corresponding to wavelengths measured in the lee wave structure.Here, we develop a new technique, described below in Sect.3, to estimate turbulent diffusivities associated with oceanic internal lee waves in the same region of the western Caribbean (Fig. 1).

Data set used in this study
The seismic data presented here (Fig. 2) come from a transect running parallel to the shore of the Caribbean coast of northern Costa Rica.The interaction of the Panama-Columbia Gyre (Richardson, 2005) with rough bathymetry creates lee waves that can be seen in the seismic data.These lee waves are stationary and persistent as evidenced by their repeated observation during the field expedition.Here, we present a second observation along the same transect as the image in Eakin at al. (2011)   The lee wave at 45 km has a smaller vertical disruption to fine structure, but exhibits a larger lateral extent displaying coherent vertical disruptions from ∼ 44 to 46 km.Nearer the sea surface, as the lee waves approach the well-stratified lower thermocline, between 200 and 350 m depth, the vertically continuous nature of the lee wave breaks down into higher horizontal wavenumber reflections that we show to be turbulent breakdown.
isopycnal displacements, ∼ 60 m, and a more striking breakdown above the generation site.
The most prominent lee wave in the image is at 40 km along the line, with a smaller lee wave at 45 km.These lee waves are seen in the data as vertically coherent displacements in otherwise laterally continuous seismic reflections above local seafloor highs.In contrast, where the seafloor is smooth (e.g. at locations 32-38 and 50-52 km along the line) no coherent vertical disruptions are observed.We show the section of line starting at ∼ 31 km as the seismic transect continues into much shallower water and out of the Panama-Columbia Gyre resulting in a lack of lee wave activity outside of our presented region.
Above the generation site at the seafloor and the imaged lee wave disruption at both 40 and 45 km are regions of more laterally discontinuous seismic reflections at 225-375 m depth.These regions coincide with the predicted termination of the vertical propagation of lee waves generated under the conditions observed in our data as described below in Sect.2.3.The method we develop here shows these regions to exhibit elevated levels of mid-water turbulence associated with the breakdown of the observed lee waves.

Turbulent isopycnal slope spectra
Oceanic turbulence is commonly measured directly by microstructure profilers in vertical casts.Such an approach yields highly detailed estimates vertically, but stations are often separated by kilometers.Given the spatial and tempowww.ocean-sci.net/12/601/2016/Ocean Sci., 12, 601-612, 2016 ral variability of turbulence in the ocean, comprehensive representations of turbulence in the ocean interior are difficult to produce using standard microstructure profiling.Horizontally towed instruments provide another approach to measuring turbulence and the results from such studies can be used to inform a seismic approach to measuring turbulence in the ocean.Klymak and Moum 2007b (herein KM07b) used horizontally towed thermistor data to calculate isopycnal slope spectra and examine levels of turbulence and internal wave energy present in the oceanic interior.Slope spectra are calculated from the horizontal gradient of the vertical isopycnal displacement (ζ ), which can be calculated from temperature measurements by ζ = z(T )−z o (T ), where z o (T ) is the mean depth-temperature relationship (Klymak and Moum, 2007a, b).Here, we are interested in the turbulence subrange of isopycnal slope spectra, which KM07b show to extend over four orders of magnitude in horizontal wavenumber (k x ).These wavenumbers include horizontal wavelengths of over 100 m, well within the resolution of SO techniques.The findings of KM07b closely match the predicted Batchelor model spectrum, yielding accurate estimates of turbulent dissipation (ε).Specifically, they find that the turbulence spectrum (ϕ ζ x ) can be written as a combination of inertial convective and inertial diffusive terms: Here the empirical mixing efficiency is 0.2 (Osborn and Cox, 1972), the constant C T is set to 0.4 (Sreenivasan, 1996), the empirical constant q is 2.3, ν is the viscosity of seawater, and N o is the mean buoyancy frequency.
At seismic resolution we can ignore the inertial diffusive term leaving: Equation ( 2) yields a predicted slope of one-third in logarithmic slope spectrum -wavenumber space with the intercept directly relating to the level of turbulent dissipation.Diffusivities (K ρ ) are calculated using K ρ = 0.2 ε / N 2 (Osborn, 1980).

Turbulent seismic slope spectra
Levels of turbulence from seismic data rely on the methods developed in KM07b.There are two techniques to examine slope spectra using seismic data that each have unique advantages and drawbacks: (a) displacement spectra from tracked seismic reflections and (b) seismic amplitude spectra taken along a depth slice of a seismic image (herein "data transforms").Displacement spectra from tracked reflections have been shown to accurately estimate both internal waves and turbulence (Eakin et al., 2011;Holbrook and Fer, 2005;Holbrook et al., 2013;Krahmann et al., 2008;Sheen et al., 2009Sheen et al., , 2011)).This method allows for estimation of turbulent diffusivity but relies on the presence of continuous seismic reflections, at least 0.8-1.2km long, which may not be found in highly turbulent areas where continuous fine structure is disrupted.Data transforms, in contrast, include all horizontal wavenumber information present in the data, can be applied equally throughout a data volume, and often show clear boundaries between internal wave, turbulence, and noise subranges (Holbrook et al., 2013).However, data transforms alone cannot provide estimates of diffusivity, as the spectral energy is affected by the amplitude of the seismic reflection, so a separate scaling step is required to link the spectral levels to isopycnal displacements.Since seismic amplitude is heavily dependent on the processing workflow, careful consideration must be given to building a sound speed model and data processing before seismic slope spectra are produced (Fortin and Holbrook, 2009;Holbrook et al., 2013).Previous studies using tracked reflections have mapped turbulent variability on a coarse scale.Sheen et al. (2009) mapped turbulent diffusivities across the Falkland Plateau using the tracked reflector method over areas 12 km long and 200 m deep to show the utility of SO to recover the uneven, patchy nature of turbulence over very large areas.Recent work by Holbrook et al. (2013) (herein H13) examines, in detail, the many issues in using seismic data to estimate turbulence.
The H13 analysis outlines many requirements and pitfalls in using SO data to estimate turbulence by analysis of both synthetic and real data.Our analysis requires that the considerations of H13 be met, so we briefly outline them here.First, careful processing of the seismic data is necessary in order to increase signal-to-noise ratios in the data, accurately capture the seismic amplitudes, and honor true subsurface reflector shapes.Second, seismic slope spectra must be calculated from Fourier data transforms over the data set in order to determine that the data are not noise dominated and that a suitable turbulent subrange exists.Estimates of diffusivity are later extracted from an interpreted quantity, tracked seismic reflections, so this step is necessary to demonstrate interpreted turbulence is not a result of spectral leakage associated with noise in the seismic data.In Figs. 3 and 4 we verify our turbulent subrange as 0.00625 to 0.0225 cpm as the data in this horizontal wavenumber subrange closely match the predicted spectra for turbulent energy in both the data transform and displacement spectra of tracked reflections.At large lateral scales in K x space, > 100 s of meters, where internal waves dominate, our image does not have continuity and does not show a clear internal wave subrange.Using only tracked reflectors can enforce continuity but our data are dominated by energy at higher wavenumbers thereby extending the turbulent subrange and masking internal wave signal in both the data transform and displacement spectra.Third, random noise is suppressed by bandpass filtering, and har- monic shot-generated noise is eliminated in the k x domain.Finally, seismic reflections are automatically tracked in order to produce displacement spectra.The reflector slope spectra are calculated by where ϕ R is the displacement power spectrum as in Holbrook and Fer (2005).Eddy diffusivity is estimated as Osborn, 1980).Additionally, H13 addresses necessary approximations and assumptions implicit in using a reflector slope spectra method for estimating turbulence from seismic images.A signal-to-noise ratio of about 4 or higher is required for spectral analysis.A larger issue is the assumption that seismic reflections, in fact, follow isopycnals.H13 finds that fine structure that persists for a sufficient timescale will record ambient displacement fields.In real data, they show that isotherms typically closely follow seismic reflections as far as 5 km and suggest that SO recovers a more highly detailed isotherm contour than a densely spaced XBT (eXpendable BathyThermograph) survey.Because our analysis windows are 400 m, we expect the seismic reflector -isopycnal assumption to be valid for our study.

Expected lee wave behavior
Lee wave behavior is governed by environmental variables related to the local bathymetry, flow velocities, and buoyancy frequency.The topographic Froude number is described by N h m /U o , where N is the buoyancy frequency, h m is the height of the bathymetry above the channel bottom, and U o is the flow speed far from the obstacle (Klymak et al., 2010).In our study location N = 5.17 × 10 −3 1 s −1 and U o avg.= 0.15 m s −1 (Eakin et al., 2011) and bathymetric heights, h m , are 100 m for the ridge at 40 km and 130 m for the ridge at 45 km.These values yield topographic Froude numbers of 3.5 and 4.9 for the lee waves seen at 40 and 45 km, respectively.The observed lee waves are thus low mode and appropriate to compare to the modeling of Klymak et al. (2010) and we can calculate the extent of lee wave propagation.
The vertical wavelength of the lee wave propagation, λ o =2π U m /N , is governed by the flow velocity, U m = U o H /(H − h m ), above the bathymetric rise over the extent of the disturbance where H is the water depth to the local low in bathymetry (Klymak et al., 2010).Using these equations, we find a flow velocity of 0.175 and 0.18 m s −1 above the lee wave generating bathymetry at 40 and 45 km, respectively.This results in a predicted vertical wavelength of the lee waves of 215 and 220 m upward into the water column from the generation site at the bathymetric highs at 40 and 45 km along the seismic transect.Using the full range of possible flow velocities reported in Eakin et al. (2011), we calculate a propagation range of 141-283 m for the lee wave at 40 km and a range of 145-292 m for the lee wave at 45 km.

Overview
Estimates of turbulence from isopycnal slope spectra rely on large volumes of data with which to build statistically significant results.In the case of seismic data, this has previously required averaging spectra from seismic reflections tracked across large regions.However, mapping oceanic turbulence related to localized processes, such as lee waves, requires reswww.ocean-sci.net/12/601/2016/Ocean Sci., 12, 601-612, 2016 olution at the scale of the feature of interest.The method we develop here uses the unique advantages of seismic data to provide such spatial resolution.By calculating slope spectra along isobaths from the seismic data itself, we do not require laterally continuous fine structure and the resultant long, clear seismic reflections to reach statistically meaningful results.Using slope spectra obtained via a Fourier transform of the seismic data in small, normalized subsections of the 2-D image, we map relative turbulent diffusivity throughout a data volume by integrating spectral energy across a known turbulent subrange.We then scale the relative turbulent energies obtained in the data transforms to absolute measures of turbulent diffusivity determined from the displacement spectra of tracked reflections.
The workflow we present here follows five steps.First, careful processing of the seismic data must be done, including suppression of harmonic noise, as outlined above.Second, the turbulent wavenumber subrange in the data set must be identified in both the data transform and displacement spectra across the entire data set.Third, a map of relative diffusivity is made by integrating spectral energy in the turbulent subrange across many small subsets in the seismic transect.Fourth, measures of turbulent spectral energy and diffusivity are calculated from data transforms and displacement spectra to determine regional scaling factors for larger subsets of the data.Lastly, the regional scaling factors are applied to the map of relative levels of spectral energy to create a high-resolution map of turbulent diffusivity across the seismic transect.

Seismic acquisition and processing
Data were acquired aboard the R/V Marcus Langseth in February 2008 as part of the Caribbean leg of the TICO-CAVA project (Van Avendonk et al., 2011, 2010).Acoustic signal was produced from a 36-element, 108 L airgun array and recorded on an 8 km, 636-channel hydrophone array, yielding subsurface sampling of 6.25 m at a sample rate of 2 ms.For our seismic line, Limon1, the acoustic source was fired every 150 m.Buoyancy frequencies were calculated from a nearby XCTD (eXpendable Conductivity, Temperature, and Depth) sample taken ∼ 120 km to the northeast 13 days after the seismic acquisition.
Data processing was completed in the manner recommended by Holbrook et al. 2013 (see Sect. 2.2) to produce reliable turbulent diffusivities from the seismic data.Standard processing steps were applied, including geometry analysis, trace editing, CMP sorting, detailed velocity analysis (Fortin and Holbrook, 2009), stacking, filtering, and migration (for further information see Yilmaz, 1987).Specialized post-processing to remove harmonic noise caused by geometric shot and receiver move-up along the transect were applied (Holbrook et al., 2013).

Seismic slope spectra
Using a data transform to supplement displacement spectra from tracked seismic reflections produced three benefits.First, we can be confident that our data are not noise dominated, as all horizontal wavenumbers are examined and noise has an easily identifiable spectral character.Second, seismic slope spectra often closely match the predicted k 1/3 x slope for turbulence and have distinct boundaries between internal wave, turbulent, and noise subranges (Holbrook et al., 2013), providing confidence in our analysis as well as clear wavenumber boundaries for the turbulent subrange.Lastly, by not relying on a limited number of tracked reflections, every sample recorded is used for spectral analysis, thus providing a comprehensive turbulence map, improved resolution, and ample statistics.
The drawback of a data transform analysis is that spectral energies are scaled by the seismic amplitude and are therefore not solely measures of absolute dissipation or diffusivity.However, within a seismic data set, relative spectral energies are preserved, assuming the data are carefully processed and the data subsets are normalized such that the maximum amplitude is one.Taking slope spectra of increasingly small subsections of a seismic data set ensures that neighboring regions were nearly identically processed, as the major control, sound speed, was quite similar.However, smaller subsections come at a loss of statistics and limit on lower wavenumber energies, so a compromise between resolution and accuracy must be made.Here, we use an analysis window of 400 m laterally by 10 m vertically to encompass the entire turbulent subrange (44-160 m), which provided enough samples to yield meaningful statistics.Figure 5 demonstrates the ability to shrink the analysis window down to small sizes yet retain the appropriate level of turbulent energy represented by the y intercept of the characteristic slope.In theory, our analysis window could shrink to the longest horizontal wavelength in the turbulent subrange, here 160 m, but we found that an analysis window ∼ 2 times larger than the maximum horizontal wavelength allows a more complete estimate of energy at lower wavenumbers.
To produce the map of relative turbulent energy, we took the data transform of normalized seismic data over regions corresponding to our chosen resolution, 400 m laterally by 10 m vertically, at non-overlapping points within the seismic data.We then integrated the spectral energy over the horizontal wavenumber range determined from the analysis of both seismic and displacement slope spectra.The resultant map of spectral energy is smoothed with a 3 × 3 sample boxcar function.Figure 6 shows the resultant map of relative seismic spectral energies in the turbulent subrange.
Where seismic reflections are difficult to track, this hybrid approach still allows for estimations of diffusivity.When the source wave field is propagated throughout the ocean, it interacts with temperature and salinity gradients of every magnitude.When we plot seismic data, the brightest reflec-  2009), who presented the first turbulence estimates from seismic data using tracked reflectors.Such a method would be inadequate to examine relatively small features, such as lee waves.The middle, green box is shown to demonstrate the match of seismic slope spectra to the characteristic slope at the scale of our tracked reflector analysis.tions dominate the image and represent the steepest T and S gradients.If the image were to be resized to include only "dim" reflections and scaled to show the brightest amplitudes in the new section, previously unseen reflections would appear.Tracking these reflections is often problematic, as they are not long and continuous.However, the data still contains information about the T and S gradients in the dim region.To make estimates in dim regions, it is important that analysis windows are normalized in the data transform method as the peak energies are in the noise subrange at high horizontal wavenumbers above our utilized turbulent subrange.This re- sults in lowered relative, and thus final, measures of turbulent diffusivity where appropriate.

Scaling to diffusivity
To scale the relative turbulence map to absolute diffusivities, it is necessary to use displacement slope spectra from tracked reflections.Since turbulent diffusivity can vary significantly across distances and depths spanned in seismic images, our scaling method uses regional "windows" that span the lateral and vertical extent of the seismic data.These windows are larger than the resolution of the data transforms described above, but significantly smaller than the entire seismic transect.Here we use windows that span 3.2 km laterally and 50 m vertically.Figure 7 shows the data transform of a representative window and Fig. 8 the corresponding displacement spectra.The consonance between the data transform and displacement spectra in our turbulent subrange allow for scaling between the two methods.
Applying the diffusivity level determined from the displacement spectra to the higher resolution data transform map is problematic due to the large size of the regional windows.The window location plays a major role in the final result since localized regions of elevated turbulence, such as those observed above the lee waves, can dominate or be overwhelmed by diffusivities estimated in the large regional window.Additionally, few tracked reflectors exist in each regional window (average = 14.9), meaning slight variances in window location, and thus included tracked reflectors, can affect the calculated diffusivities.To eliminate this effect, we calculated overlapping regional windows in a 5 × 5 grid spanning a half window length in each principal direction and assigned the average value to the regional window, the result of which is shown in Fig. 9.
Diffusivities were calculated by scaling the median value of the data transform map over each regional window to the diffusivity derived via the displacement spectra of that regional window by K ρ = (ϕ DT / DT )× K ρ (ϕ R x ), where ϕ DT  is the relative seismic spectral energy, DT is the median of the regional window of the data transform, and K ρ (ϕ R x ) is the diffusivity calculated from the regional window of tracked reflections.This procedure translates the regional window of spectral energy to diffusivity values while preserving relative window-internal differences present in the data transform map.A 3 × 3 boxcar smoothing function was then applied producing the final map of diffusivities at the resolution of the data transform map, shown in Fig. 10.

Synthetic seismic analysis
To test the accuracy of the data transform method, we created synthetic seismic data with known structures and levels of internal wave energy and turbulent diffusivity.Our synthetic data were developed from a sound speed cross section calculated from XCTD 51 from the Norwegian sea (as used in Nandi et al., 2004) over the depth range of 202-400 m and lateral extent of 24 km.This uniform background model was then subjected to internal wave and turbulence displacements, as specified in KM07, applied in a sheared checkerboard pattern.The range of turbulent diffusivities was set from 1-10 × 10 −4 m 2 s −1 to capture a range of diffusivity levels similar to those observed in our real data.The sheared checkerboard pattern was chosen to test the method's ability to reliably capture both levels of turbulent diffusivity and structural changes.Synthetic seismic data were created by convolving a 45 Hz Ricker wavelet with reflection coefficients calculated from the sound speed model described above, assuming constant density.The result is an artificial seismic data set with known levels and structures of turbulent diffusivity designed to resemble observations in our real data.
We then applied our turbulence mapping procedure to this synthetic seismic data to test the method's accuracy.Figures 11 and 12 show the data transform and displacement spectra of the data and were used to determine the turbulent subrange.Figure 13 shows the seismic data, turbulence perturbation applied, and recovered result of our turbulent mapping method.While the general checkerboard pattern and levels of diffusivity are generally recovered, the shear direction is not.By not resolving the direction of hotspot shear, coupled with the low amplitude of one hotspot, these results may indicate that features of this size and diffusivity are near the resolution of our method, likely due to our choice of smoothing parameters.While we recover the checkerboard and accurately estimate diffusivity levels, further method testing would be required to attempt to image turbulent features significantly smaller than ∼ 2 km laterally by ∼ 30 m vertically.

Uncertainties for real data
Uncertainties in our method mainly stem from fitting the characteristic slope of turbulent spectra to the displacement spectra of tracked reflectors and finding the intercept, or diffusivity, associated with the characteristic.We estimate uncertainty by calculating the 95 % bootstrap confidence interval for our displacement spectra.For the final turbulence map, we calculate confidence intervals for each regional window.

Seismic observations
Our observations of lee waves in the seismic data are wellmatched to the predicted behavior of lee waves calculated in Sect.2.3.For the lee wave at 40 km along the transect we see a clear end to the vertical propagation of the lee wave and a vertical wavelength of 225 m, close to the predicted 215 m.The lee wave at 45 km is less pronounced in amplitude and has a less abrupt end to its vertical propagation but appears to change behavior 290 m above its generation site, slightly shallower than the predicted 220 m wavelength.This discrepancy is likely due to local complications concerning the downstream, secondary lee wave.Alternatively, a current velocity of 0.2 m s −1 predicts an extent of 292 m and is within the range of velocities reported in Eakin et al. (2011), which may suggest a faster current velocity at 45 km.In both cases, the seismic data show the termination of coherent vertical lee wave disruption near the calculated wavelength.

Synthetics
Our synthetic seismic analysis shows that turbulence estimates on features as small as 2-5 km wide with vertical extents of only tens of meters can be extracted from seismic data.Our recovered result resembles the known perturbation; the checkerboard pattern is recovered, but the direction of shear is less well resolved (Fig. 13).Edge effects and the application of a smoothing window decimate the upper left hotspot and have some amplitude-decreasing effect on the other two upper hotspots.Levels of diffusivity are generally well recovered, with the exception of the lower right hotspot, demonstrating the method's efficacy to provide accurate estimates of turbulent diffusivity at this scale.However, for the underestimated diffusivities in the lower right hotspot, the known maximum diffusivity of 10 × 10 −4 m 2 s −1 is within the uncertainty calculated over the regional window, as the diffusivity recovered is ∼ 6.5 × 10 −4 m 2 s −1 and the 95 % confidence interval is 3.5-12 × 10 −4 m 2 s −1 .While the percent error seems high, this method is best suited to investigate relative turbulent differences across regions tens to hundreds of kilometers long and hundreds of meters deep.Additionally, a factor of less than 2 (or < 100 % uncertainty) is relatively small considering turbulent diffusivities can span many orders of magnitude across similar regions.

Near-seafloor turbulence
We observe elevated levels of turbulent diffusivity just above the seafloor over rough bathymetry (Fig. 10), a common feature where flow interacts with varying topography (Polzin et al., 1997;St. Laurent and Thurnherr, 2007).The highest levels of seafloor turbulence are in the troughs between seafloor highs; turbulence there is likely created as the flow of the gyre spills over the ridge.Here, we observe turbulence similar to that of St. Laurent and Thurnherr (2007) and confirm the speculation of Eakin et al. (2011) that these troughs are regions of elevated turbulence.
The levels we observe in these troughs indicate strong, but not extreme, ocean mixing of 4-5 × 10 −4 m 2 s −1 .However, combined over roughly 50 m depth and 10 km extent, this mixing layer likely represents a significant contributor to regional diapycnal mixing.While caused primarily by the relatively steady current of the Panama-Columbia Gyre, there could also be some contribution of 3-D effects related to bottom flow moving near orthogonal to our seismic transect.

Mid-water turbulence
The most significant finding we present is the observation of mid-water turbulence non-local to the generation site of large amplitude lee waves associated with lee wave breakdown near the lower thermocline, a result consistent with Nikurashin and Ferrari (2010) as our diffusivity levels are not extreme.In the seismic image, it is visually clear that the horizontal wavenumber content of the seismic reflections is enhanced between 300 and 400 m depths directly above the lee wave generation sites.Our slope spectra analysis quantifies this, showing significantly more energy in the turbulent subrange in this region, revealing 4-5 times more turbulent diffusivity than adjacent waters both vertically and laterally.We note that this is a ∼ 3 times greater turbulent diffusivity than estimated from tracked seismic reflections.This is expected, as turbulent energy disrupts laterally continuous oceanic fine structure, thereby eliminating continuous seismic reflections necessary for the tracking algorithm to follow.
At 44-46 km at mid-water depth of 225-350 m, turbulent diffusivities are even higher than those observed near the seafloor.Both lee waves exhibit elevated levels of turbulent diffusivity above their generation sites at the predicted wavelengths of ∼ 215 and ∼ 220 m above the seafloor.Peak diffusivity for both lee waves on this transect lies around 5 × 10 −4 m 2 s −1 , tapering down to background levels 100-125 m below the maximum.Given levels of uncertainty, our elevated regions of turbulence could be as low as 3.5 × 10 −4 m 2 s −1 , or as high as 1.3 × 10 −3 m 2 s −1 .
In similar regions where ocean ridges, islands, and seamounts exist, mid-water turbulence associated with lee waves could play an important role in ocean mixing.Here, we observe lee waves in intermediate waters cascading into turbulence at the bottom of the thermocline.Global non-local lee wave mixing likely contributes significantly to meridional overturning (Scott et al., 2011) and could play a significant role in the downward diffusion of heat in the ocean, thereby impacting global climate.

Conclusions
Using seismic techniques, we demonstrate that turbulent diffusivity can be mapped across ocean basins and full water depth at a scale of 400 m laterally by 10 m vertically.Applied to synthetic seismic data, our technique reliably reproduced both locations and diffusivity levels of synthetic "hotspots".The value of the method's improved resolution is demonstrated in lee waves observed in the western Caribbean, a feature often difficult to study with traditional oceanographic tools due to their limited spatial extent and stationary nature.We observe two regions of elevated turbulence associated with the generation and breakdown of each lee wave.Offshore Costa Rica we see elevated turbulence around local bathymetric highs that generate lee waves.Turbulence near rough bathymetry is elevated to ∼ 3 times that of the adjacent water column.We see high seafloor turbulence in troughs between seafloor highs, particularly when the seafloor is strongly sloped.
In addition to elevated turbulence along the seabed, we observe strong mid-water turbulent diffusivity where the lee wave propagation ends at the lower thermocline for both lee waves observed in our data.In the seismic image, this is seen as a sudden increase in the horizontal wavenumber content of seismic reflections at the upper limit of the lee wave's coherent vertical structure between 250 and 400 m depth.Both data transforms and displacement spectra show the characteristic behavior for turbulence and elevated spectral energies in the turbulent subrange above the lee wave.We estimate turbulent diffusivities associated with the lee waves Ocean Sci., 12, 601-612, 2016 www.ocean-sci.net/12/601/2016/ of ∼ 5 × 10 −4 m 2 s −1 , which is 4-5 times higher than background levels along the transect.

Figure 1 .
Figure 1.Location map of seismic line Limon1.The flow interaction of the Panama-Columbia Gyre (Richardson, 2005) with the rough bathymetry creates lee waves along the transect in two locations.

Figure 2 .
Figure 2. Seismic image of line Limon1 shows two lee wave structures as continuous vertical displacements at 40 and 45 km.The lee wave at 40 km is significantly more pronounced and has a vertical displacement of 35 m.The lee wave at 45 km has a smaller vertical disruption to fine structure, but exhibits a larger lateral extent displaying coherent vertical disruptions from ∼ 44 to 46 km.Nearer the sea surface, as the lee waves approach the well-stratified lower thermocline, between 200 and 350 m depth, the vertically continuous nature of the lee wave breaks down into higher horizontal wavenumber reflections that we show to be turbulent breakdown.

Figure 3 .
Figure 3. (a) Region of analysis of seismic slope spectra (data transform) of Limon1.(b) Spectra of Limon1 matches the characteristic turbulence slope in horizontal wavenumber -slope spectra space over the same wavenumber range as the displacement spectra of Fig. 4 as marked by gray vertical lines.The signal-to-noise ratio here is 5.54.The drop-off at high wavenumbers is due to frequency filtering in the seismic data from 30 to 80 Hz.

Figure 4 .
Figure 4. (a) Tracks of 553 continuous seismic reflections longer than 800 meters are shown in black over seismic line Limon1.(b) Displacement slope spectra with 95 % bootstrap confidence intervals of the tracked seismic reflections show line Limon1 to be dominated by the characteristic slope of turbulence, shown in cyan, over the horizontal wavenumber range (0.00625-0.0225)cpm, or (44-160) m horizontal wavelengths.Average turbulent diffusivity calculated for Limon1 over this image is 1.54 × 10 −4 m 2 s −1 .The turbulence subrange of this data set is emphasized herein by gray vertical lines marking the upper and lower wavenumber limits.Characteristic slopes of internal waves (GM77) and noise (KM07b) are shown in yellow and green.

Figure 5 .
Figure 5. (a) Line Limon1 with three analysis boxes of varying size that correspond to normalized seismic slope spectra of the same color in (b).The three analysis regions match the characteristic turbulence slope of the Batchelor spectrum (cyan line) over our prescribed wavenumber range marked by gray vertical lines.The smallest, red box is our analysis window and resolution for turbulence mapping via data transform and sized 400 m by 10 m.The largest, blue box has the best statistics and fit to the characteristic slope and is 12 000 m by 200 m in size.This large, blue box was chosen to mimic the analysis of Sheen et al. (2009), who presented the first turbulence estimates from seismic data using tracked reflectors.Such a method would be inadequate to examine relatively small features, such as lee waves.The middle, green box is shown to demonstrate the match of seismic slope spectra to the characteristic slope at the scale of our tracked reflector analysis.

Figure 6 .
Figure 6.Map of seismic spectral energy in the turbulent subrange calculated via data transform.Here, the color-map does not show diffusivities.Instead, this map gives a qualitative view of relative turbulent energy across the section.White lines track the location of the lee wave features, traced from the seismic data.

Figure 7 .
Figure 7. Seismic data and data transform of a representative window from which relative turbulent energies are calculated and scaled.

Figure 8 .
Figure 8. Tracks of seismic data and displacement spectra with 95 % bootstrap confidence intervals (vertical bars) of the representative window over which scaling from relative turbulent energy to estimates of diffusivity are calculated.At this scale, the same wavenumber subrange exhibits the characteristic slope of turbulence (cyan line).This segment yields an average diffusivity of 9.0 × 10 −5 m 2 s −1 with a 95 % confidence range of 5.6 × 10 −5 -1.5 × 10 −4 m 2 s −1 .Characteristic slopes of internal waves (GM77) and noise (KM07b) are shown in yellow and green.

Figure 9 .
Figure 9. Turbulent diffusivity map created by displacement spectra calculated over windows 3.2 km wide by 50 m deep.At this rough scale, elevated turbulent diffusivity is observed around the rough bathymetry as well as above the lee wave generation sites.

Figure 10 .
Figure 10.Final map of turbulent diffusivity across line Limon1 at resolution 400 m by 10 m. White tracks of seismic reflections along the large amplitude lee waves are shown for reference.We observe elevated turbulence near the rough bathymetry and above the lee waves at 40 and 45 km.Levels of turbulent diffusivity are significantly higher above the lee waves structures and ∼ 50 times higher than the abyssal open ocean.

Figure 11 .
Figure 11.Synthetic seismic data and data transform.

Figure 12 .
Figure 12.Synthetic seismic section overlaid with tracked reflections and displacement spectra with 95 % bootstrap confidence intervals (vertical bars).Displacement spectra match the characteristic slope of turbulence (cyan line) across horizontal wavenumbers 0.017-0.043cpm.Areas where we introduced turbulence are clearly seen as they contain few tracked reflections highlighting the difficulty in using displacement spectra to estimate turbulent energy in seismic data.Characteristic slopes of internal waves (GM77) and noise (KM07b) are shown in yellow and green.

Figure 13 .
Figure 13.(a) Synthetic seismic data where reflections represent fine structure calculated from an XCTD density profile disrupted by an internal wave field with energy in the turbulent subrange applied fitting the sheared checkerboard pattern shown in (b).The recovered turbulent diffusivity from the seismic slope spectra method developed here is shown in (c).
but collected approximately 24 h earlier.The lee waves captured in this data set show larger vertical