## Abstract

Broadband measurements of the internal wavefield will help to unlock an understanding of the energy cascade within the oceanic realm. However, there are challenges in acquiring observations with sufficient spatial resolution, especially in horizontal dimensions. Seismic reflection profiling can achieve a horizontal and vertical resolution of order meters. It is suitable for imaging thermohaline fine structure on scales that range from tens of meters to hundreds of kilometers. This range straddles the transition from internal wave to turbulent regimes. Here, the authors analyze an 80-km-long seismic image from the Falkland Plateau and calculate vertical displacement spectra of tracked reflections. First, they show that these spectra are consistent with the Garrett–Munk model at small horizontal wavenumbers (i.e., *k*_{x} ≲ 3 × 10^{−3} cpm). There is a transition to stratified turbulence at larger wavenumbers (i.e., *k*_{x} ≳ 2 × 10^{−1} cpm). This transition occurs at length scales that are significantly larger than the Ozmidov length scale above which stratification is expected to modify isotropic Kolmogorov turbulence. Second, the authors observe a rapid onset of this stratified turbulence over a narrow range of length scales. This onset is consistent with a characteristic energy injection scale of stratified turbulence with a forward cascade toward smaller scales through isotropic turbulence below the Ozmidov length scale culminating in microscale dissipation. Finally, they estimate the spatial pattern of diapycnal diffusivity and show that the existence of an injection scale can increase these estimates by a factor of 2.

## 1. Introduction

The oceanic internal wavefield probably arises from a forward cascade of energy from large-scale to small-scale processes (Thorpe 2005). Spectral analysis of this wavefield has played a useful role in developing quantitative models. For example, the power spectrum of vertical density displacements as a function of horizontal wavenumber, *ϕ*_{ζ}(*k*_{x}) (see Table 1 for a description of all variables), shows that distinctive regimes exist with different spectral slopes. At *k*_{x} < 5 × 10^{−3} cpm, corresponding to length scales of > *O*(10^{2}–10^{3}) m, the Garrett–Munk model provides an accurate empirical description of the behavior of internal waves (Garrett and Munk 1975). At higher values of *k*_{x}, a transition into what is conventionally assumed to be a turbulent regime is observed (Fig. 1). This transition is generally attributed to breaking of internal waves and to different kinds of convective and/or shear instabilities that can occur within a stratified fluid. In this turbulent regime, *ϕ*_{ζ}(*k*_{x}) varies as a function of , which distinguishes it from the internal wave regime. At sufficiently small length scales, an exponent of −5/3 is consistent with an inertial convective subrange that is based on isotropic turbulent models (Kolmogorov 1941; Obukhov 1949; Corrsin 1951; Batchelor et al. 1959).

It is increasingly evident that flow at horizontal length scales of *O*(10^{2}) m within a sufficiently stratified fluid does not always satisfy the underlying assumptions of these canonical models (Lindborg 2006; Riley and Lindborg 2008). For example, at horizontal scales greater than the Ozmidov length scale *l*_{O}, isotropic overturning can be strongly suppressed and the fundamental properties of turbulence are moderated by stratification. The length scale *l*_{O} is given by

where *ε* is the dissipation rate of turbulent kinetic energy per unit mass and *N* is the buoyancy frequency given by

where *ρ* is potential density.

Since *l*_{O} is typically *O*(10^{−2}–10^{0}) m, it is reasonable to infer that larger scales are associated with anisotropic flow, which fundamentally differs from that postulated by the Obukhov–Corrsin model (Gargett and Hendricks 1981). Horizontal flow is unconstrained by a stabilizing buoyancy force and so vertical fluctuations are expected to be smaller than horizontal fluctuations. Lindborg (2006) suggested that a horizontal energy spectrum with a power-law exponent of −5/3 is energetically consistent with a strongly anisotropic inertial flow regime, which is perhaps confusingly referred to as “stratified turbulence” (Riley and Lindborg 2010). To discriminate between turbulence within a stratified regime and stratified turbulence, we use the term layered anisotropic stratified turbulence (LAST) to define the regime referred to by Lindborg (2006). The existence of this LAST regime is supported by reinterpretation of published observations and by numerical simulations (Riley and Lindborg 2008; Brethouwer et al. 2007).

Here, we describe and analyze a seismic reflection experiment from the Falkland Plateau in the South Atlantic Ocean. Records from this experiment are used to construct a vertical image of the water column that reveals the detailed thermohaline structure at equal horizontal and vertical resolutions. We have four principal aims. First, we wish to demonstrate that meaningful information about the internal wave and turbulent regimes can be extracted by careful processing of seismic reflection datasets. In this regard, our approach builds on and complements the analysis and recommendations of Holbrook et al. (2013). Second, spectral analysis of vertical displacements of undulating reflections is carried out in order to investigate internal wave and turbulent regimes as a function of horizontal wavenumber (Holbrook and Fer 2005). Third, we use averaging and normalization methods to investigate the nature of the transition between internal waves and turbulence that has significant fluid dynamical implications. Fourthly, we estimate the spatial distribution of mixing and dissipation along a seismic image (Sheen et al. 2009; Holbrook et al. 2013). This approach complements global calculations made using one-dimensional microstructure profiling (e.g., Waterhouse et al. 2014).

## 2. Seismic imaging of thermohaline structure

Seismic reflection experiments use a controlled source to make well-resolved images of the earth’s subsurface. Acoustic energy is generated by priming tuned arrays of air guns with compressed air. These arrays are repeatedly fired to expel regular pulses of compressed air into the water column. Such arrays have total volumes of >150 L, and the vertically directed acoustic energy has a typical frequency bandwidth of 10–200 Hz. Energy from each pulse is transmitted through the subsurface and reflected at impedance contrasts. In the oceans, these contrasts are produced by temperature contrasts as small as 0.03°C over a few meters (Nandi et al. 2004). Salinity generally makes a minor contribution (Sallarès et al. 2009). Reflected acoustic energy is recorded by a towed streamer of hydrophones that is typically 2–12 km long. Since the reflected energy has a low signal-to-noise ratio, each point in the subsurface is recorded multiple times over a period of tens of minutes. This sampling redundancy enables signal stacking that is used to improve the signal-to-noise ratio. Following Holbrook et al. (2013), we estimate the signal-to-noise ratio for two adjacent seismograms to be

where *c* is the maximum value of the cross-correlation of both traces and *a* is the zero-lag autocorrelation of the first trace.

Although seismic reflection technology was developed to image the solid earth, Holbrook et al. (2003) demonstrated that this technology is eminently suitable for mapping thermohaline fine structure. In a typical two-dimensional seismic experiment, vertical slices extending from the sea surface down to the sea bed are acquired. The 80-km-long seismic image analyzed here is located ~100 km east of the Falkland Islands in the South Atlantic Ocean (Fig. 2). The original experiment was carried out by WesternGeco Ltd. in February 1993. Its geometric configuration is shown in Fig. 3. During this experiment, a tuned array of 36 guns with a total volume of 119 L was towed behind the vessel at an average depth of 7.5 m. Vessel speed was 2 m s^{−1} and the gun array was fired every 40 m (i.e., every 20 s). Further astern, a 4.8-km-long streamer consisting of 240 hydrophones spaced every 20 m was towed at a depth of 10 m. Horizontal offset between the air gun array and the start of the active streamer was 97 m. The common midpoint interval is 10 m, which yields a 60-fold redundancy of coverage. Note that the first second of two-way travel time was not recorded during acquisition.

This dataset was previously processed and analyzed by Sheen et al. (2009). Subsequently, Holbrook et al. (2013) have shown that significantly improved seismic images can be produced by paying particular attention to elements of the processing sequence (e.g., suppression of random and harmonic noise, poststack migration). Following Ruddick et al. (2009), Fortin and Holbrook (2009), and Holbrook et al. (2013), our refined processing methodology exploits standard techniques that are adapted from those used to construct seismic images of the solid earth (Yilmaz 2001). There are three particularly important steps. First, bandpass and wavenumber filtering is applied to ameliorate the influence of ambient and harmonic noise, respectively. Randomly generated ambient noise is suppressed using a zero-phase, bandpass (i.e., 12–100 Hz) Butterworth filter. As Holbrook et al. (2013) remark, harmonic noise can be especially significant when seismic images are spectrally analyzed in the horizontal wavenumber domain. This form of noise is shot-generated and occurs at integer multiples of the shot spacing (i.e., every 40 m or 0.025 cpm). These noise spikes are suppressed by applying a band-stop notch filter centered over each spike in the wavenumber domain.

Second, individual shot records are sorted into common midpoint (cmp) records that are stacked to generate a coherent seismic image with an optimal signal-to-noise ratio. Stacking is carried out by correcting for offset between each shot–receiver pair. This correction relies on carefully choosing the root-mean-square (rms) sound speed of seawater as a function of two-way travel time for shot–receiver pairs that share a common point of reflection at depth. Although sound speed generally varies only between 1470 and 1530 m s^{−1}, these rms functions must be chosen and applied with considerable care. It is also important that velocity picking is sufficiently dense (e.g., every 1–3 km) to allow for horizontal changes in sound speed (Fortin and Holbrook 2009).

Finally, seismic data are recorded as a function of the time elapsed between generation and detection of acoustic energy (i.e., two-way travel time). To correctly locate reflected signals within the spatial domain, seismic images are migrated from elapsed time into correct depth. This migration process is carried out either before or after a two-dimensional seismic image is constructed by stacking. It requires knowledge of sound speed as a function of two-way travel time. Sheen et al. (2009) carried out an iterative prestack depth migration. However, this form of prestack algorithm can degrade slope spectra at higher wavenumbers (Holbrook et al. 2013). Here, we have followed the recommendations of Holbrook et al. (2013) and carried out poststack time migration using a standard frequency–wavenumber algorithm (Stolt 1978). They also suggested that conversion to depth be carried out using a sound speed of 1500 m s^{−1}. We note that changing the sound speed used for depth conversion by ±30 m s^{−1} does not significantly affect the conclusions we draw from spectral analyses.

Coeval hydrographic measurements of temperature and salinity were not acquired during this seismic experiment. Here, we have chosen a legacy hydrographic database of meter-scale-resolution CTD casts acquired during December–April of 1972–2011 (www.nodc.noaa.gov). These casts are located less than 200 km from our seismic experiment (Fig. 2a). We chose to display calculated buoyancy frequency profiles as a function of mensal range (Figs. 2b–e). The average profile does not change significantly over ±4 months (note that a subset of CTD casts, shown in Fig. 2e and acquired in a single cruise, is offset to higher-than-expected values and is not used in our analysis). In this study, we use an average profile of *N* as a function of depth based on CTD casts acquired between December and March (i.e., ±2 months on either side of the seismic experiment). During this period, the standard deviation of the average *N* profile is ± 0.3 cph between 0.5 and 1.5 km.

## 3. Spectral analysis of fine structure

### a. Reflective event tracking

Seismic images of thermohaline fine structure reveal patterns of coherent undulating reflections. A substantial number of these reflections can be traced over distances of several kilometers (Fig. 4). Although these reflections occasionally occur as transgressive filaments, they often track isopycnal surfaces (Holbrook and Fer 2005; Krahmann et al. 2008, 2009; Sheen et al. 2009; Biescas et al. 2014). This observation is sufficient, but not strictly necessary, to make inferences about the internal wavefield. A more important requirement is that, over length scales of interest, these undulations are governed by the internal wavefield. This requirement is thought to be the case when 5 × 10^{−4} < *k*_{x} < 10^{−1} cpm (Krahmann et al. 2009). Most practitioners deem that it is reasonable to infer that seismic images are approximate snapshots of vertical isopycnal displacements.

To analyze stacked seismic images spectrally, it is necessary to track reflections (Holbrook and Fer 2005; Sheen et al. 2009). Accurate and automated tracking of discontinuous events with variable signal-to-noise ratios that variously grow, climb, descend, bifurcate, merge, and die is not straightforward. Here, automated tracking was carried out using the method described by Holbrook et al. (2013). First, the amplitude of each reflection is normalized to ±1 by calculating the cosine of the instantaneous phase angle. This angle is determined from the Hilbert transform of each individual vertical seismic trace. Second, the normalized reflections are contoured in order to identify and enclose individual continuous reflections. Third, individual tracks are identified using the average vertical position of each contour along its length. Holbrook et al. (2013) recommend using a contour value of ±0.6. We tested a range of values and found that a value of ±0.8 maximizes the number of tracks, while still yielding faithful tracking. To remove long wavelength features that may not be generated by the internal wavefield, tracked features were linearly detrended.

A total of 856 reflections were individually tracked across the seismic image (Fig. 4b). The total length of tracked reflections on this image is 1200 km, which is broadly comparable to 880 km of tracked internal waves from a typical hydrographic experiment using a towed instrument in the vicinity of Hawaii (Klymak and Moum 2007b). Subsequently, we have chosen to analyze a subset of the total tracked length consisting of tracks, each of which is longer than 2 km and has a signal-to-noise ratio of greater than 3.5. These chosen values fulfill the requirement for a large range of wavenumbers and are based on the recommendations of Holbrook et al. (2013). This subset has 88 tracks and a total track length of 270 km.

### b. Spectra of tracked reflections

Power spectra of the vertical displacement of detrended horizontal tracks were calculated using multitaper spectral analysis. This technique produces significantly less variability and bias than a standard periodogram (Thomson 1982). Vertical displacement power spectra are converted into horizontal slope spectra using (Klymak and Moum 2007a). This conversion emphasizes the transition from the internal wave to the turbulent regime, which now takes the form of a switch from negative to positive exponents.

We note in passing that there is little consensus on the exact value of the exponent for internal wave slope spectra, which is unlikely to be constant throughout the oceanic realm. For example, the GM75 model of Garrett and Munk (1975) has an exponent of −0.5 for the internal wave slope spectrum. In contrast, the GM76 model of Cairns and Williams (1976) has an exponent of zero. Other studies suggest that a roll-off occurs at an exponent of −1 toward higher wavenumbers (Gargett and Hendricks 1981). It is reasonable to infer that a range of values from 0 to −1 are consistent with slope spectra of the internal wave field. This range is qualitatively distinct from the turbulent spectrum that is expected to have an exponent of −5/3 + 2 = 1/3, where +2 comes from the multiplication by (2*πk*_{x})^{2} when converting vertical displacement spectra to slope spectra.

The suitability of a seismic image for spectral analysis is gauged by calculating its power–wavenumber spectrum (Holbrook et al. 2013). Figure 5 shows slope spectra that have been calculated for two panels of tracks shown in Fig. 4. These spectra demonstrate that internal wave and turbulent regimes are present with power-law exponents of −1 and ⅓, respectively. At wavenumbers >0.04 cpm, white noise starts to dominate and these higher wavenumbers were discarded. These spectral tests show that the turbulent regime is clearly identifiable at high wavenumbers.

Holbrook et al. (2013) emphasize the importance of identifying and removing harmonic noise, which can badly contaminate slope spectra, especially at higher wavenumbers. On the dataset presented here, a single harmonic noise spike occurs at *k*_{x} = 2.5 × 10^{−2} cpm, which has been excised using the method described by Holbrook et al. (2013). In Fig. 6, spectral analysis of a panel from Fig. 4 demonstrates that harmonic noise has been successfully removed.

### c. Temporal blurring

Finally, we tackle an issue that afflicts all hydrographic sampling technologies, namely, how to adequately sample moving fluid structure. Seismic images are constructed by stacking together shot–receiver pairs that are recorded over a finite period of time. Therefore, the resultant images are susceptible to blurring. This susceptibility might compromise our ability to adequately image internal wave and turbulent regimes. During stacking, multiple shot–receiver pairs (i.e., a cmp gather) that image the same portion of the subsurface are added together (Fig. 3). The time taken for a common midpoint gather to be acquired, *τ*, depends on the ship’s speed *V* and on the length of the streamer *L*, where

A finite duration of imaging will tend to blur structures that translate either vertically or horizontally by distances that are comparable to the spatial resolution of the seismic experiment. Inevitably, *V* is constrained by the technical requirements of towing a long streamer. However, *L* can effectively be changed by changing the length of streamer used during processing (i.e., discarding records from more distal portions of the streamer). A shorter streamer has a smaller imaging duration that will have the effect of sharpening the image of a moving structure at the expense of a lower signal-to-noise ratio. Conversely, a longer streamer yields an improved signal-to-noise ratio but has a greater susceptibility to blurring. In this seismic experiment, *L* = 4800 m and *V* = 2 m s^{−1}, which yields *τ* ≲ 17 min. If the geostrophic velocity is 0.1 m s^{−1}, structures could move horizontally by up to 100 m during this interval. Similarly, if *N* = 1 cph, 17 min represents more than one-quarter of the buoyancy period. In both cases, the stacked image may suffer from blurring. Thus, at the horizontal length scales of interest in this study, the vertical and horizontal motion of internal waves might, or might not, be significant compared with *τ*.

To estimate how spatial blurring could alter our spectral analyses, we have analyzed a series of partially stacked images that were constructed using different values of *L*. As *L* is progressively reduced from 4.8 to 1 km, *τ* correspondingly reduces from 17 to 3.5 min. The effect that decreasing values of *τ* have on calculated slope spectra is illustrated in Fig. 7. As *τ* is reduced, the transition between the internal wave and turbulent regimes sharpens (cf. Figs. 7a,c). For *τ* ≲ 3.5 minutes, spectral deterioration is caused by a decrease in the signal-to-noise ratio.

This result suggests that spatial blurring is not significant at the considered time scales. An alternative, but less plausible, possibility is that blurring is always significant. We support the first possibility for two reasons. First, synthetic seismic experiments, in which *τ* is varied, do not significantly distort spectra. Second, we do not think that the clear consistency between our observed spectral power-law exponents and those measured by other hydrographic techniques is fortuitous (Klymak and Moum 2007a,b). Here, we have used *L* = 4.8 km because the signal-to-noise ratio is marginally better than for *L* = 3 km.

### d. Grouped and averaged spectra

An important goal is identification of spectral subranges from their characteristic slopes. Unfortunately, individual slope spectra have low signal-to-noise ratios and some form of preliminary averaging is desirable. First, spectra of tracked reflections >2 km in length are sorted according to their estimated energy level, which is given by the median value of each spectrum for 0.004 ≤ *k*_{x} ≤ 0.024 cpm. Sorted spectra are then averaged into groups of four, yielding a total of 22 groups (Fig. 8).

At low wavenumbers (i.e., *k*_{x} < 0.002 cpm), observed exponents are consistently negative with a pronounced rollover at the lowest wavenumbers (i.e., a shallowing of the gradient of the reflection slope spectra). With increasing wavenumber, the steepest gradients of the slope spectra occur just before a crossover into positive exponents. These observations are consistent with slope spectral predictions of the GM76 model, which has a rollover of up to −1 (Cairns and Williams 1976; Gargett and Hendricks 1981; Gregg 1993).

At higher wavenumbers (i.e., *k*_{x} > 0.005 cpm), a positive exponent of ⅓ is observed. Klymak and Moum (2007b) demonstrated that the slope spectrum of the inertial convective turbulent regime, , is given by

where Γ = 0.2 is the classical turbulent flux coefficient that relates the kinetic energy dissipation rate *ε* to an appropriately averaged buoyancy flux (Osborn 1980), *C*_{T} = 0.4 is the Kolmogorov constant (Sreenivasan 1995), and *N* is the buoyancy frequency [Eq. (2)].

Here, we are less concerned with the inertial diffusive subrange where isotropic turbulence occurs at higher wavenumbers. At horizontal wavelengths that exceed 100 m, isotropic turbulence is unlikely to be the dominant process. Instead, it has been suggested that an inherently anisotropic and stratified (i.e., LAST) turbulent model applies. In this case, the horizontal kinetic energy spectrum is given by

where *k*_{x} is horizontal wavenumber and *C*_{k} ≃ 0.5 is an empirical constant estimated from numerical simulations of strongly stratified turbulent fluid flow (Lindborg 2006). This model also has a power-law exponent of −5/3 that is equivalent to a slope spectral gradient of ⅓.

The grouped slope spectra shown in Fig. 8 suggest that internal wave and turbulent regimes are identifiable and that spectra are displaced vertically and horizontally according to energy level. However, these grouped spectra are still quite noisy, and it is difficult to determine with confidence the nature of the crossover between the two regimes. Crossover from negative to positive gradients for slope spectra marks the transition from an internal wave regime to an appropriately defined turbulent regime. D’Asaro and Lien (2000) pointed out that the shape of this crossover ought to contain important information about the dynamics of the transition from one regime to another (e.g., Fig. 1). An additive model assumes that internal waves and layered anisotropic stratified turbulence coexist across a range of scales, whereas an onset model assumes that a significant change of behavior occurs at a crossover scale that triggers turbulence. This turbulence is still strongly affected by stratification since this crossover scale is assumed to be substantially larger than the Ozmidov scale *l*_{O} (cf. D’Asaro and Lien 2000). Thus, from a fluid dynamical perspective, an important goal is to determine the spectral shape of this crossover. For slope spectra, the crossover for an additive model is expected to be smooth and U-shaped without a sharply defined minimum, whereas the crossover for an onset model is expected to be sharp and V-shaped with no transitional subrange.

Here, we address the crossover imaging problem by calculating average normalized (i.e., stacked) spectra with a view to further improving the signal-to-noise ratios in the vicinity of the crossover locus. Simple averaging does not faithfully preserve crossover shape since the wavenumber at which crossover occurs varies as a function of both energy level and stratification (Figs. 10a,d,g). To bring the crossover region into better focus, we have developed and tested two different forms of normalization (Fig. 9). Both forms of normalization shift spectra with respect to each other. Although scaling along the *x* and *y* axes is preserved, absolute values are not. These values have been omitted from figure panels where appropriate.

Preliminary averaging into 22 groups helps to improve the signal-to-noise ratio and also allows the approximate crossover loci to be identified on grouped spectra. For each grouped spectrum, this approximate locus is determined by fitting a three-component model with subranges that have power-law exponents corresponding to internal waves, turbulence, and white noise. Intersections between internal wave and turbulent subranges yield a set of approximate crossover loci. To avoid bias, those parts of the spectra within ±0.2 logarithmic units of the predicted crossover wavenumber are ignored when fitting the three-component model.

In linear normalization, approximate crossover loci are fitted with a straight line using linear regression (e.g., Fig. 11c). Each crossover locus is projected orthogonally onto this line to give a projected crossover point. Grouped spectra are then averaged in a direction that is parallel to this best-fit line (i.e., all projected crossover points are collapsed in the direction of this line onto a single average value; Figs. 10b,e,h and 11d). Thus, linear normalization is equivalent to averaging parallel to a rotated *y* axis where the angle of rotation is that between the axis and the best-fit line. Note that linear normalization is not the same as point normalization where spectra are shifted so that the approximate crossover loci become coincident in *k*_{x}– space.

In nonlinear normalization, a value of *ε* is estimated from the turbulent subrange of each grouped spectrum using Eq. (5). Internal wave energy levels were then determined from values of *ε* using the Gregg–Henyey parameterization. Each energy level is used to calculate an internal wave spectrum for a GM76 model with a high wavenumber roll-off where *N* = 1.4 cph and is the bandwidth parameter (J. Klymak 2014, personal communication; Cairns and Williams 1976; Gargett and Hendricks 1981).

Intersections between predicted internal wave and turbulent slope spectra constrain a set of crossover points that lie along a curve in *k*_{x}– space. Normalization is achieved by sliding grouped spectra along this curve before summing and averaging (Figs. 10c,f,i). In other words, averaging is carried out along a curved rather than a straight line.

Figure 10 shows the resultant spectra for simple, linear, and nonlinear normalization of all 22 groups of slope spectra. Note that usage of the term “normalization” does not mean that there is a single normalization factor that relates these spectra and the original spectra. Quality of fit for all three forms of averaging with reference to the two competing models is quantitatively assessed in Figs. 10d–i. When simple averaging is carried out, it is difficult to discriminate between additive and onset models. With either linear or nonlinear normalization, a sharply defined crossover location is visible, which suggests that an onset model is more appropriate. It is important to emphasize that this result is not dependent on the use of multitaper spectral analysis. Thus, a method based on constructing periodograms also yields a sharp crossover, but the resultant spectra are noisier. Linear normalization is preferred since it does not require an internal wave model, apart from the choice of a representative power-law exponent. We note in passing that a sharp crossover between internal wave and turbulent regimes has also been observed on direct data transforms of seismic images (Holbrook et al. 2013).

An important consideration is that normalization is underpinned by fitting spectra using a fixed set of subranges. To address the possibility of bias, we carried out 4941 individual calculations for which power-law exponents of the internal wave and turbulent regimes were varied from −0.4 to 0.2 in 81 steps, and from −0.1 to 1.8 in 61 steps, respectively (Fig. 11). As before, linear normalization was carried out to determine an average spectrum in each case. All 4941 average spectra were used to produce a density plot that shows the resulting final averaged spectrum is robust with respect to model choice (Fig. 11e). This plot reinforces the observation that the transition between the internal wave and turbulent regimes is rapid and that the internal wave slope spectrum is consistent with a power-law exponent of −1 (Gargett and Hendricks 1981).

### e. Monte Carlo analysis

To further test the robustness of the normalization method, Monte Carlo analysis of synthetic spectra was performed. The purpose of this analysis is to address the following questions. First, can an underlying onset model be reliably recovered? Second, could an underlying additive model with a smooth crossover transition be artificially sharpened to mimic an onset model? By analyzing different synthetic datasets, we can assess the robustness and reliability of both linear and nonlinear normalization of spectra.

The normalization method uses a simple spectral model to identify the approximate position of the crossover between internal wave and turbulent regimes. This procedure is necessary because normalization requires observed spectra to be translated in a direction that is compatible with all crossover loci. It is important to ascertain whether or not this model-based translation biases the calculated average spectra in any way.

Two measures were employed to avoid artificially sharpening the crossover region. First, when fitting the model spectrum, regions within ±0.2 logarithmic units of the model’s crossover point were omitted. This omission prevents any single deviation from biasing crossover location or geometry. Second, once the crossover location is found, observed spectra are always normalized by translation in one direction that is either a straight line (i.e., linear normalization) or a curve (i.e., nonlinear normalization). Point normalization where all crossover locations are averaged to give a single point should be avoided.

Monte Carlo analysis was tested on a database of 88 individual synthetic spectra. These spectra were generated by adding normally distributed (1*σ* = 0.3) random noise to either additive or onset spectral models (Fig. 12a). Crossover loci of these synthetic spectra shift to lower wavenumbers with increasing power as expected. Consequently, a simple average of all 88 spectra will always yield an average spectrum with a smooth transition between the internal wave and turbulent regimes. As before, individual spectra were grouped according to median amplitude into 22 spectra that are shown in Fig. 12b. For each group spectrum, the approximate crossover location was found by fitting a model spectrum (Fig. 12c). Group spectra were then normalized to yield an average spectrum (Fig. 12d).

This procedure was repeated 500 times for different populations of random noise. The 500 calculated average spectra are summarized in the form of a density plot (Fig. 12). When either an onset or an additive model is used to generate synthetic spectra, the resultant density plots show that the correct spectral shape is reliably recovered, provided that a suitable averaging procedure is applied (Fig. 12e,i). The two most important features of this procedure are linear (or nonlinear) normalization and omission of the central portion of grouped spectra. These features strongly mitigate against “self-sharpening” of crossover loci.

If central portions of spectra in the vicinity of crossover loci are included, the expected spectral shapes are usually preserved (Figs. 12f,j). If point normalization is used instead of linear normalization, spectral shapes are also largely unchanged, although a small kink is visible on the additive model (Figs. 12g,h). However, if both of these features (i.e., retention of central portions and point normalization) are used, more noticeable spectral distortion can occur (Figs. 12h,l). It is clear that both onset and additive spectra are artificially sharpened. The greater the value of 1*σ*, the more pronounced this distortion becomes.

We conclude that appropriate normalization of spectra does not cause artificial sharpening of the crossover region. We have shown that a combination of linear normalization and omission of the central portion of spectra ensures that sharpening does not occur. It is particularly important not to use point normalization, which can result in self-sharpening of spectra.

## 4. Fluid dynamical implications

Careful analysis of slope spectra from seismic images demonstrates that the turbulent regime exists to horizontal wavenumbers as low as 10^{−2} cpm. The transition from the internal wave to the turbulent regime is sharp. We wish to outline the fluid dynamical implications of these observations. Lindborg (2006) argued that a turbulent regime, exhibiting horizontal spectra with characteristic power-law dependence at length scales that greatly exceed the Ozmidov scale, is energetically consistent with a strongly anisotropic, yet still inertial, flow regime. The existence of such a regime is supported by atmospheric and oceanographic observations with some underpinning provided by numerical simulations (Brethouwer et al. 2007; Riley and Lindborg 2008).

As already noted, this profoundly anisotropic (i.e., vertical velocities are much smaller than horizontal velocities), yet inherently three-dimensional and turbulent, flow regime is often referred to as “stratified turbulence” in the fluid dynamical literature (Lindborg 2006; Brethouwer et al. 2007). It is characterized by the development of layering whose vertical scale is set by *l*_{υ} ~ *U*/*N*, where *U* is a characteristic horizontal flow velocity. The horizontal scale *l*_{h} ≫ *l*_{υ} is set by the dissipation rate of turbulent kinetic energy *l*_{h} ~ *U*^{3}/*ε*. In this case, the horizontal Froude number *F*_{h} is given by

Scaling arguments suggest a relationship between *l*_{h} and the Ozmidov scale *l*_{O}, where

The existence of this regime, which we refer to as the LAST regime, is supported by reinterpretation of published observations by Riley and Lindborg (2008) and of idealized numerical simulations by Brethouwer et al. (2007). Since turbulent flow within the LAST regime has a horizontal power spectrum proportional to , an associated slope spectra must have positive power-law dependence on *k*_{x}, and so there exists a wavenumber *k*_{C} (i.e., length scale *l*_{C} = 1/*k*_{C}) at which there is a crossover from a slope spectra with wave-like characteristics to a slope spectra with turbulent-like characteristics.

We have considered two possible crossover models. The first is one where the observed slope spectrum is an additive combination of wave-like and turbulent-like spectra where (Klymak and Moum 2007a,b). This additive model suggests that the wavenumber *k*_{h} = 1/*l*_{h} associated with the horizontal extent of the turbulent layers *k*_{h} < *k*_{C} (i.e., the horizontal extent of layers is larger than the crossover scale) and that the existence of turbulence on scales smaller than *l*_{C} (i.e., wavenumbers *k* < *k*_{C}) does not immediately destroy wave-like behavior.

In this case, both turbulent and wave-like motions exist over a range of scales and the additive crossover will be smooth and curved. The predicted flow structure, showing both wave-like motions and turbulence patches at all horizontal scales, is illustrated in Fig. 13a. The inherently additive nature of the underlying power spectrum containing both wave-like and turbulence-like contributions is shown schematically in Fig. 13c. This additive model is based on the observation that internal wave and turbulent spectra decay as a function of wavenumber at different rates. Therefore, the crossover scale from one power-law description to another marks the scale at which one becomes more dominant. The crossover scale simply reflects a change in the balance of two physical processes acting over a range of scales, and the crossover scale itself has no particular physical significance.

Because of the central scaling assumptions of the LAST regime, the vertical scale *l*_{υ} ≪ *l*_{h} with *l*_{υ} ≫ *l*_{O}. Thus, inherently anisotropic turbulence occurs for all horizontal scales *l*_{O} ≤ *l* ≤ *l*_{h}. For vertical scales smaller than *l*_{O}, stratification is, in some sense, no longer sufficiently strong to affect turbulence. It is therefore possible for isotropic turbulence with a classical inertial range to occur for horizontal scales smaller than *l*_{O} provided that the Ozmidov scale is sufficiently large compared to the Kolmogorov microscale, *l*_{K} = (*ν*^{3}/*ε*)^{1/4}, where *ν* is the kinematic viscosity of the fluid.

This final condition for the existence of the LAST regime (i.e., *l*_{O} ≫ *l*_{K}) is equivalent to the requirement that the buoyancy Reynolds number ≫ *O*(1), where

It is debatable what constitutes an appropriately large value of for the existence of the LAST regime. Shih et al. (2005) suggest that if > *O*(100), then the system is fully energetic (i.e., its dynamics are free of viscous effects). In contrast, Bartello and Tobias (2013) showed that a −5/3 spectral dependence occurs if > *O*(10) based on very high-resolution numerical simulations.

The alternative, and our favored, onset model is illustrated in Figs. 13b and 13d. In this case, there is a pronounced change in slope at the crossover length scale *l*_{C} that separates wave-like and turbulence-like spectra. At some horizontal length scale (e.g., *l*_{h} of the layers central to the LAST regime), waves break down catastrophically and practically no wave-like dynamics survive to higher wavenumbers. Wave energy is injected into the turbulent regime at this characteristic onset scale. Conversely, little turbulence exists at scales larger than the crossover length scale. Therefore, the forward cascade of turbulence ensures that the spectrum for all wavenumbers greater than the crossover scale is completely dominated by turbulence dynamics. The slope spectrum has a +⅓ power-law dependence on horizontal wavenumber. This dependence is assumed to be associated with the LAST regime for *k*_{C} = *k*_{h} < *k*_{x} < *k*_{O} and with classical isotropic turbulence for *k*_{O} < *k*_{x} < *k*_{K} = 1/*l*_{K}. Since the power-law dependence of spectra is expected to be identical both above and below the Ozmidov scale, it is reasonable to assume that any premultiplying factors that scale spectral power will be the same on either side of the crossover. The predicted flow structure shows wave-like motions at large and intermediate scales but patches of turbulence at intermediate and smaller scales (Fig. 13b). The inherently onset nature of the underlying power spectrum, comprising a wave-like power spectrum at low wavenumbers and a turbulence-like power spectrum at high wavenumbers, is shown schematically in Fig. 13d.

For this end member, the crossover scale has a physical meaning that corresponds to the scale at which turbulence onsets and internal waves break down. The mechanism underlying such a process is probably a scale-selective physical process that leads to a catastrophic decrease of energy within the internal wave regime. Candidate processes for such a scale-selective onset include primary internal wave instabilities and nonlinear interaction within the wave field. In essence, the crossover scale represents an injection scale for the forward cascade of turbulent energy within the LAST regime, and it is reasonable to suppose that *l*_{C} corresponds to the typical horizontal scale *l*_{h} of the anisotropic and high-aspect ratio layers characteristic of this regime (Brethouwer et al. 2007). Little coherent internal wave dynamics can be expected to survive at larger wavenumbers since the wavefield breaks down because of the onset of spatially and temporally incoherent turbulent motions. Thus, a sharp crossover marks the sudden onset of stratified turbulent behavior that has limited overlap with the internal wave regime.

It is important to emphasize that the LAST regime is an idealized model for turbulence within a stratified fluid that is dynamically unaffected by rotation. The scale of the turbulent layer may be such that rotation might affect its ultimate horizontal extent. Nonetheless, the dynamics of turbulence within that layer is small enough and fast enough for rotation to be dynamically unimportant. An additional constraint is that the crossover length scale is sufficiently small and that the flow velocities are sufficiently large so the effects of rotation can be neglected. In particular, the anisotropic turbulent layers required for the LAST regime to exist are not necessarily manifestations of the low-frequency “vortical mode” with nonzero potential vorticity affected by planetary rotation (Thorpe 2005). Finally, we note that there are alternative explanations for the existence of power spectra with a power-law decay of at wavenumbers that are inconsistent with isotropic turbulence. For example, Hua et al. (2013) suggested that baroclinic instability of preexisting quasigeostrophic vortices could give rise to this spectral slope. However, our observations suggest that such instability dynamics are not necessary for the manifestation of dependence. The precise nature of the crossover between internal wave and turbulence regimes is challenging to determine by experiment or by numerical simulation because of the required range of length and time scales. Our observations provide an important constraint.

## 5. Diapycnal diffusivity

Diapycnal diffusivity *K*_{T} can be determined across the seismic image using slope spectra calculated from tracked reflections (Holbrook and Fer 2005; Sheen et al. 2009; Holbrook et al. 2013). First, slope spectra of all individual tracked reflections >640 m are calculated using the methodology described in section 3. Spectra are fitted using three power-law functions with exponents of −1, ⅓, and 2. These lines correspond to the internal wave, turbulent, and white noise regimes, respectively (Fig. 14). Second, these starting fits are only used as the basis for isolating that part of each spectrum that corresponds to turbulence. Since the power of the turbulent regime is more sensitive to energy level, we can exploit this portion of the spectra to calculate *K*_{T}. As already noted, because of the continuity that must apply between spectra associated with inertial convective isotropic turbulence below the Ozmidov scale and LAST regime turbulence above the Ozmidov scale, it is straightforward to convert into *ε* using Eq. (5). Following Osborn (1980), *ε* is converted into *K*_{T} using

where, for simplicity, Γ = 0.2. The value of *N* at any depth is given by the average profile shown in Fig. 2. In this way, we can determine the spatial distribution of *K*_{T} (Fig. 15a). Its average value is 3.1 × 10^{−5} m^{2} s^{−3}, which is broadly consistent with regional hydrographic studies (Garabato et al. 2004; Waterhouse et al. 2014). The spatial variation of *K*_{T} closely follows the geometry of the thermohaline structure. For example, reduced values of *K*_{T} occur over an eddy structure located at a range of 70 km and at a depth of 900 m. Bands of changing values of *K*_{T} crosscut the image, dipping in the opposite direction to the bathymetric slope. The apparent increase in *K*_{T} toward the sea surface is probably an artifact caused by an increase in ambient noise. We note in passing that Sheen et al. (2009) carried out similar mixing calculations based on spectral analysis of both internal wave and turbulent regimes. However, a direct data transform analysis of their processed image highlighted the drawbacks of their particular implementation of a frequency–wavenumber migration algorithm (Holbrook et al. 2013). Furthermore, Sheen et al. (2009) used a less robust form of reflection tracking that introduced spectral artifacts, especially at *k*_{x} > 10^{−2} m^{−1}. Consequently, Fig. 2b of Sheen et al. (2009) differs in several respects from Fig. 15b.

Values of *ε* and *N* can be used to calculate the variation of Ozmidov length scale *l*_{O} across the image using Eq. (1). We obtain *l*_{O} values of *O*(0.1–1) m, which agree with those previously observed (e.g., Gargett and Hendricks 1981). These values are substantially smaller than the length scales at which the spectral characteristics of turbulence (i.e., ) are observed.

Previous analyses of seismic reflection images exploited both internal wave and turbulent regimes to constrain dissipation, and hence diapycnal diffusivity, using the Osborn (1980) model. These approaches assumed a power-law exponent of −0.5 for the internal wave slope spectrum, in accordance with the GM75 model (Garrett and Munk 1975). However, competing models for the exponent of the internal wave slope spectrum exist, and values between 0 and −1 could reasonably be used. An attractive property of the onset model is that diffusivity calculations are independent of the slope chosen for the internal wave regime. To compare onset and additive values of *K*_{T}, we chose a value of −0.5 for the exponent of the internal wave regime in agreement with the GM75 model and with previous seismic oceanographic studies (Holbrook and Fer 2005; Krahmann et al. 2008; Sheen et al. 2009; Holbrook et al. 2013).

An onset model necessarily yields higher estimates for *K*_{T}. This outcome occurs for purely geometric reasons since, for a given value of *K*_{T}, an additive spectrum will always have higher amplitude than an onset spectrum. Hence, when fitting slope spectral data, a lower *K*_{T} will be required to match the amplitudes observed in the input data if the additive model is used. In areas where the signal-to-noise ratio is ≥4, the average increase in the value of *K*_{T} is by a factor of ~2 with considerable spatial variation. Note that we do not calculate *K*_{T} from the internal wave regime. Instead, we assume that this regime is well represented by a single power-law relationship. We then use either an additive or an onset model in the fitting stage. The diffusivity *K*_{T} is calculated from the turbulent component alone. This procedure sidesteps the vexed issue of equating *K*_{T} with power of the internal wave regime.

The fact that similar reflections are sometimes located above and below one another means that individual undulations are not statistically independent. This possibility could affect uncertainties at lower wavenumbers, but a detailed study is beyond the scope of this study. Calculated mixing rates, which rely on the higher wavenumber portion of spectra, are unlikely to be adversely affected by a lack of statistical independence. Other sources of uncertainty can be estimated and their effects propagated using Eqs. (5) and (10). For example, Γ and *C*_{T} have uncertainties of at least ±0.04 and ±0.05, which yield uncertainties in log_{10}*KT* of ±0.04 and ±0.08 logarithmic units (i.e., ~9% and ~20%), respectively (Moum 1996; Sreenivasan 1995). Note that the likely uncertainty in the sound speed profile used for depth conversion yields a small shift in *K*_{T} of ~±0.025 logarithmic units (i.e., ~5%).

The buoyancy frequency *N* is probably the most important source of uncertainty in this study, particularly since coeval hydrographic measurements are unavailable (Fig. 2). The mean value of *N* observed between 500 and 1500 m depth and within ±2 months of the survey month is 1.32 cph with a standard deviation of *σ* = 0.3. Over 90% of *N* measurements fall between 0.5 and 2.5 cph. If *N* ± 1*σ* is propagated through Eqs. (5) and (10), the resulting log_{10}*K*_{T} for a given *ε* changes by about −0.5 and +0.3 logarithmic units (i.e., a decrease of 70% or an increase of 100%), respectively (cf. Figs. 15b–d). This uncertainty in *K*_{T} is small compared to the observed spatial variation of *K*_{T}. Furthermore, since legacy buoyancy frequency profiles tend to have similar shapes but different magnitudes, it is likely that this uncertainty yields a static shift away from the correct values rather than variable spatial patterns. An important caveat exists for regions where thermohaline structures manifest horizontal variability. For example, the region above the eddy in Fig. 15 may have lower rates of mixing. If so, higher stratification (i.e., *N* ≈ 5 cph) caused by vertical compression of isopycnal surfaces could account for this observation. If the observed variation in is solely caused by buoyancy frequency changes and if *ε* is fixed at 10^{−10} m^{2} s^{−3}, *N* would have to vary between 0.5 and 7 cph, which is a larger range than hydrographic observations could reasonably support (Fig. 16).

One final source of uncertainty arises from fitting noisy spectra. In Figs. 14b–d, the identified turbulent subrange is fitted by systematically varying *K*_{T}. In each case, the misfit *χ*^{2} is plotted as a function of *K*_{T}. Well-defined global minima exist and, for an appropriate tolerance (e.g., twice the minimum value of *χ*^{2}), the uncertainty in *K*_{T} is no worse than one-half of an order of magnitude. The uncertainty that arises from actual identification of the turbulent subrange, which we believe to be robust, is beyond the scope of this contribution. It is important to emphasize that all of these sources of uncertainty do not affect our two principal conclusions. First, the lowest wavenumber portion of the −5/3 subrange cannot be accounted for by isotropic (i.e., Kolmogorov) turbulence but is consistent with the LAST model (Lindborg 2006). Second, a sharp onset crossover between internal wave and turbulent regimes exists.

## 6. Conclusions

We show that horizontal slope spectra obtained by tracking reflections across a two-dimensional seismic image have the expected power-law relationships. The high quality of these data, combined with autotracking methodology and spectral analysis, permit closer investigation of the crossover from internal wave to turbulent regimes for vertical displacement power spectra. This crossover occurs at horizontal length scales that are substantially larger than that those considered plausible for isotropic turbulence. Instead, it is more likely that crossover is caused by the onset of a flow regime that we have referred to as the layered anisotropic stratified turbulent (LAST) regime.

Our results suggest that crossover between regimes is rapid. In particular, we do not observe a transitional subrange that would be characteristic of an additive model in which internal waves and turbulence coexist over a range of scales. This observation suggests that there is a switch in the governing fluid dynamics from internal waves to turbulence without a significant overlap of the two regimes. A sharp transition is suggestive of an instability or nonlinear process that causes the internal wavefield to break down catastrophically so that little energy remains within the wavefield at smaller scales. This breakdown to the LAST regime occurs at a well-defined length scale that is substantially larger than the Ozmidov scale. Central to our interpretation is the existence of a scale-selective mechanism that destroys the wavefield and sets the characteristic large injection scale of the turbulent dynamics. It remains a challenge to identify this mechanism.

## Acknowledgments

M.F. is supported by the Department of Earth Sciences. Research activity of C.P.C. is supported by EPSRC Programme Grant EP/K034529/1 (Mathematical Underpinnings of Stratified Turbulence). We thank C. Bond, A. Dickinson, K. Gunn, W. S. Holbrook, J. Klymak, J. Moum, and S. Thorpe for their help. We are very grateful to J. Klymak for generously making available his MATLAB toolbox for calculating Garrett–Munk spectra.

## REFERENCES

*Geophys. Res. Lett.*,

**32**, L15604, doi:.

*Ten Chapters in Turbulence*, P. A. Davidson, Y. Kaneda, and K. R. Sreenivasan, Eds., Cambridge University Press, Cambridge, 269–317.

**22**, 192–205, doi:.

*The Turbulent Ocean.*Cambridge University Press, 439 pp.

*Seismic Data Analysis.*Society of Exploration Geophysicists, 2027 pp.

## Footnotes

Department of Earth Sciences contribution number esc.3582.