Journal cover Journal topic
Ocean Science An interactive open-access journal of the European Geosciences Union
Journal topic
Ocean Sci., 14, 403–415, 2018
https://doi.org/10.5194/os-14-403-2018
Ocean Sci., 14, 403–415, 2018
https://doi.org/10.5194/os-14-403-2018

Research article 06 Jun 2018

Research article | 06 Jun 2018

# High-resolution diapycnal mixing map of the Alboran Sea thermocline from seismic reflection images

High-resolution diapycnal mixing map of the Alboran Sea thermocline from seismic reflection images
Jhon F. Mojica1,2, Valentí Sallarès2, and Berta Biescas2,3 Jhon F. Mojica et al.
• 1Center for global Sea Level Change CSLC – NYUAD, Abu Dhabi, United Arab Emirates
• 2Institute of Marine Sciences, ICM-CSIC, Barcelona, Spain
• 3Consiglio Nazionale delle Ricerche CNR-ISMAR, Bologna, Italy

Correspondence: Jhon F. Mojica (jhon.mojica@nyu.edu)

Abstract

The Alboran Sea is a dynamically active region where the salty and warm Mediterranean water first encounters the incoming milder and cooler Atlantic water. The interaction between these two water masses originates a set of sub-mesoscale structures and a complex sequence of processes that entail mixing close to the thermocline. Here we present a high-resolution map of the diapycnal diffusivity around the thermocline depth obtained using acoustic data recorded with a high-resolution multichannel seismic system. The map reveals a patchy thermocline, with spots of strong diapycnal mixing juxtaposed with areas of weaker mixing. The patch size is of a few kilometers in the horizontal scale and of 10–15 m in the vertical one. The comparison of the obtained maps with the original acoustic images shows that mixing tends to concentrate in areas where internal waves, which are ubiquitous in the surveyed area, become unstable and shear instabilities develop, enhancing energy transfer towards the turbulent regime. These results are also compared with others obtained using more conventional oceanographic probes. The values estimated based on the seismic data are within the ranges of values obtained from oceanographic data analysis, and they are also consistent with reference theoretical values. Overall, our results demonstrate that high-resolution seismic systems allow the remote quantification of mixing at the thermocline depth with unprecedented resolution.

1 Introduction

Diapycnal diffusivity (kρ) around the thermocline plays a major role in controlling the strength and pattern of the ocean circulation because it determines heat and salt heterogeneity at different spatial scales. This process usually occurs in a vertically stratified regime, affecting adjacent layers with the same density but different temperature and salinity (Stewart, 2008). In terms of processes, mixing in the ocean can be separated into two categories. One corresponds to internal wave (IW) breaking, which produces turbulent motion and changes the density stratification, while the second concerns the development of high-frequency dynamic instabilities that are formed due to shear (Gregg, 1987; D'Asaro and Lien, 2000). As the spatial scale decreases, mixing leads to an unbalanced pressure field that eventually results in a collapse and dispersion of mixing waters through isopycnals (Thorpe, 2005). The value of kρ depends on the buoyancy frequency (N) and the dissipation rate (ε) as indicated by the so-called Osborn (1980) relationship:

$\begin{array}{}\text{(1)}& {k}_{\mathit{\rho }}=\mathrm{\Gamma }\mathit{\epsilon }/{N}^{\mathrm{2}}.\end{array}$

This value, where Γ= 0.2 is the empirically defined mixing efficiency (Osborn and Cox, 1972), corresponds to the mixing among isopycnal layers in the thermocline. The global mean kρ value is of the order of 10−4 m2 s−1 (Munk and Wunsch, 1998), which corresponds to the value required to keep overturning in the thermocline. It has been shown that if kρ < 10−5 m2 s−1, the energy is not enough to generate mixing (Gregg, 1989).

In a conservative flow, ε might present small variation due to dissipated heat through turbulent motions, but in the presence of strong shear ε tends to increase (Thorpe, 2005), reaching a maximum value close to the Kolmogorov scale (Gargett and Holloway, 1984). Good knowledge of its behavior provides important clues on available energy and its transfer among spatial scales.

The loss rate of kinetic energy in the turbulent motion is commonly expressed as

$\begin{array}{}\text{(2)}& & \mathit{\epsilon }=\left(\frac{\mathit{\nu }}{\mathrm{2}}\right)〈{S}_{ij}{S}_{ij}〉,\text{(3)}& & {S}_{ij}=\left(\frac{\partial {u}_{i}}{\partial {x}_{j}}+\frac{\partial {u}_{j}}{\partial {x}_{i}}\right),\end{array}$

where ν= 1.064 × 10−6 m2 s−1 is the kinematic viscosity and the tensor Sij is a function of the velocity components in the three orthogonal directions (Thorpe, 2005). Conventional in situ techniques such as vertical microstructure turbulence profilers (VMPs) or MicroRiders provide the most accurate measures of kρ, but in just one dimension. In general, although measures are accurate in the vertical dimension, sampling in the horizontal direction is much poorer, particularly in the  103–101 m range (Klymak and Moum, 2007a, b). Since this is the range of scales at which the transition between isotropic internal wave and anisotropic turbulent motions (i.e., mixing) occurs, the observational evidence of mixing patterns and the understanding of the underlying physical mechanisms are rather limited so far. Overall, direct measures and observations are too few to create a global mixing map with the required resolution to feed the models with proper values of dissipation rates (Smyth et al., 2011). In turn this makes it difficult to integrate mixing into large-scale models of ocean dynamics. Mixing effects are simulated instead through the incorporation of eddy diffusivity coefficients, which are tuned ad hoc to match the large-scale distribution of ocean observables. While this approach allows the proper reproduction of regional spatial–temporal patterns, it severely hampers the long-term predictive capability of ocean dynamics and, in turn, that of climate models. Thus improving our knowledge of the short-term and small-scale mixing mechanisms and integrating them into large-scale models remain outstanding challenges.

To overcome this issue, remote-sensing techniques have recently started to be used (e.g., Gibson et al., 2007). One of these alternative techniques is multichannel seismics (MCS), an acoustic method providing quasi-synoptic images of the thermohaline boundaries in the ocean interior to full ocean depth, with a lateral resolution of O(101 m) (Holbrook et al., 2003). Several recent works have demonstrated that it is actually possible to map kρ using measures of the horizontal wave number (kx) spectra of the vertical displacements of thermohaline boundaries imaged with MCS acquisition systems (Sheen et al., 2009; Holbrook et al., 2013; Fortin et al., 2016). However, these studies use conventional, relatively low-resolution systems with source energy concentrating below  50 Hz. In addition, due to the long wavelength source wavelet, conventional MCS systems are not well suited to image the shallowest ocean layers (i.e., < 200 m), but rather deeper water levels (> 400 m depth). At these depth levels, changes in the internal structure are usually less marked than at shallower levels, especially around the thermocline. In a recent work, it has been shown that portable, high-resolution MCS (HR-MCS) systems, which use a smaller energy but higher-frequency source (> 150 Hz), allow imaging of the thermohaline structure as shallow as  30 m with a lateral resolution of 12–15 m and a vertical resolution of 1–2 m (Sallarès et al., 2016). This resolution is three- to fourfold better than that of conventional MCS systems that have been used to image deeper ocean levels. Therefore, it has the potential to image sub-mesoscale structures and processes that affect the thermocline at scales of kilometers to tens of meters, allowing coverage of the existing observational gap. Despite their potential, HR-MCS systems have never been used to date to quantify diapycnal mixing at the thermocline depth.

Figure 1(a) Bathymetric map of the Alboran Sea and location of the data used in the study. HR-MCS profile acquired during the IMPULS 2006 experiment (black line labeled 3) from expendable bathythermograph (XBT) profilers (red circles) and an expendable conductivity–temperature–depth (XCTD) probe (green circle). Acoustic Doppler current profiler (ADCP) lines (black dotted line). Geostrophic velocity for 17 May 2006 (gray arrows). (b) Temperature–salinity diagram from the XCTD probe. σ0 is the potential density in kilograms per cubic meter. Color scale indicates depth.

Here we use the abovementioned method of extracting kρ(x,z) maps from MCS images, but applied for the first time to HR-MCS data acquired in the Alboran Sea (westernmost Mediterranean). The method to calculate diapycnal mixing maps from the horizontal wave number spectra of vertical reflector displacements is based on that proposed by Sheen et al. (2009) and Holbrook et al. (2013). The result is a high-resolution mixing map of the ocean at the thermocline depth (30–110 m) along a 35 km long transect (Fig. 1a). This method could be applied in other regions where the shallow water levels are sufficiently stratified to allow the recording of the energy reflected at variations of acoustic impedance (density × sound speed contrasts among neighboring water layers).

The rest of the paper is structured as follows: in Sect. 2 we present the hydrographic context and the observations; then we describe the acquisition system and the method applied to estimate kρ from the seismic data. The results are described in Sect. 3, whereas the discussion about the imaged structures and their likely causes is presented in Sect.  4. Finally, Sect. 5 summarizes the main conclusions.

2 Data and methodology

The Alboran Sea is characterized by the continuous exchange between Mediterranean Water (MW) and Atlantic Water (AW) through the Strait of Gibraltar. This exchange concentrates near the surface (between  30 m and  200 m), where the shallow, moderately salty, and cold incoming AW (< 50 m) interacts with the deeper, warmer, saltier, and more stable outgoing MW, producing another water mass known as Modified Atlantic Water (MAW) (Chioua et al., 2013). In this framework, internal waves, strong horizontal shear instability, and prominent thermohaline stratification are generated. These particular features reflect the complex dynamic setting of the area, with kinetic energy being transferred among isopycnals from large to small scales, leading eventually to overturning, isotropic turbulence, and irreversible mixing.

The data set used in this work, which includes collocated seismic and oceanographic measurements, was collected onboard the Spanish R/V Hesperides in the framework of the IMPULS 2006 experiment. Here, we concentrate our analysis on one of the seismic profiles, acquired using a portable HR-MCS system (IMPULS-3 profile). The acquisition started on 16 May at 23:43 UTC and finished on 17 May at 04:00 UTC. In total, some 4 h were required to record a 38 km long profile. The acquisition system consisted of a 4.75 L source with a peak frequency at 150–190 Hz. As mentioned above, the corresponding size of the Fresnel zone, a proxy of the horizontal resolution, is 12–15 m depending on the target depth. The streamer was 300 m long and had 48 channels, with a group spacing of 6.25 m. The shot interval was 15 m, giving a common midpoint (CMP) gather fold of 6. The locations of the different data are displayed in Fig. 1a.

Figure 2Depth-converted high-resolution multichannel seismic profile, with the tracked reflectors superimposed (blue lines). (See Fig. 1a for location).

This profile was first processed and then used to estimate the average kx energy spectra of vertical displacements of the imaged reflectors (i.e., the acoustic images of thermohaline boundaries). A total of 68 reflectors rather homogeneously distributed throughout the surveyed area, with lengths of 1.5–21 km and a signal-to-noise ratio higher than 8 within the frequency range of 40–240 Hz (Sallares et al., 2016), were tracked and used for the analysis (Fig. 2). Vertical profiles of temperature and pressure were recorded simultaneously with the seismic acquisition using four expendable bathythermographs (XBTs), whereas the salinity and buoyancy profiles were derived from an expendable conductivity–temperature–depth (XCTD) probe dropped 3 days after the seismic acquisition. Water current profiles recorded with an Acoustic Doppler current profiler (ADCP) in the same season as the seismic experiment but 4 years later (see location in Fig. 1a) have also been used.

The IMPULS-3 profile shown in Fig. 2 reveals a number of laterally coherent seismic reflectors that are assumed to follow isopycnals. Biescas et al. (2014) showed that this assumption is valid in regions free of salinity–temperature-compensating intrusions, which is a reasonable approximation for the Alboran Sea. The analysis of the obtained kx spectra has allowed the identification of three subranges that control dynamics around the thermocline depth (Sallares et al., 2016). At scales larger than the horizontal buoyancy wavelength (lN 90 m), motions are dominated by the IW field (IW field subrange). Then the spectra rolls off, reflecting the presence of shear instabilities of probably the Kelvin–Helmholtz (KH) type, which appear to collapse at a scale of  30 m (transitional, or instability-dominated subrange), giving rise to turbulence at even smaller scales (turbulent subrange) (Fig. S3 in the Supplement). A more detailed description of these ranges and their scales of influence is presented in Sallares et al. (2016). In the present work, we use energy levels at the turbulent subrange, obtained from the kx spectral analysis of the tracked reflectors within small analyzing windows, to estimate the lateral and vertical variations in ε and kρ along the whole profile.

Since our data set does not include direct measurements of turbulence, we use the XCTD and ADCP data to estimate a vertical profile of kρ based on Gregg's (1989) model, hereafter referred to as Gregg89. This model assumes that energy dissipation in the thermocline takes place through IW energy transfer by wave–wave interaction. As the relative local change in buoyancy frequency, which relates the level of stretching and squeezing of isopycnals, is small in the surveyed area ( 1), we assume that the assumptions and approximations of the model are valid; thus it is not necessary to consider alternative ones such as the one presented in Waterman et al. (2013). The Gregg89 model is commonly applied in the midlatitude thermocline, linking shear current at different depths. The simplest way to obtain average dissipation rates over large space and timescales is through

$\begin{array}{}\text{(4)}& & \mathit{\epsilon }=\mathrm{7}×{\mathrm{10}}^{-\mathrm{10}}{N}^{\mathrm{2}}/{N}_{\mathrm{0}}^{\mathrm{2}}<{S}_{\mathrm{10}}^{\mathrm{4}}/{S}_{\mathrm{GM}}^{\mathrm{4}}>,\text{(5)}& & {S}_{\mathrm{10}}^{\mathrm{4}}=\mathrm{4.22}{\left[{\left(\mathrm{\Delta }U/\mathrm{\Delta }z\right)}^{\mathrm{2}}+{\left(\mathrm{\Delta }V/\mathrm{\Delta }z\right)}^{\mathrm{2}}\right]}^{\mathrm{2}},\text{(6)}& & {S}_{\mathrm{GM}}^{\mathrm{4}}=\mathrm{2}{\left[\left(\mathrm{3}\mathit{\pi }/\mathrm{2}\right){j}_{x}{E}_{\mathrm{GM}}b{N}_{\mathrm{0}}^{\mathrm{2}}{k}_{x}^{c}{\left(N/{N}_{\mathrm{0}}\right)}^{\mathrm{2}}\right]}^{\mathrm{2}},\end{array}$

where N0= 5.2 × 10−3 s−1 is the reference buoyancy frequency, S10 is the shear variance calculated from the meridional (V) and zonal (U) velocity variations as a function of depth (z), SGM is the variance for the Garrett–Munk model (Gregg, 1989), jx is a mode number, EGM is the Garrett–Munk energy density, b is the scale depth of the thermocline, c is the spectrum slope, and kxis the horizontal wave number.

Alternatively, the model proposed by Batchelor (1959), hereafter referred to as Batchelor59, estimates kρ as a function of the energy transfer from large to small scales in the turbulent regime. This model assumes that the energy exchange from mechanical to caloric due to N and ε can be approximated as

$\begin{array}{}\text{(7)}& {\mathit{\phi }}_{\mathit{\varsigma }}^{T}=\left(\frac{\mathrm{4}\mathit{\pi }\mathrm{\Gamma }}{{N}^{\mathrm{2}}}\right){C}_{\mathrm{T}}{\mathit{\epsilon }}_{\mathrm{T}}^{\mathrm{2}/\mathrm{3}}{\left(\mathrm{2}\mathit{\pi }{k}_{x}\right)}^{-\mathrm{5}/\mathrm{3}},\end{array}$

where φς is the energy spectrum of the isopycnals' vertical displacement measured in the turbulent subrange, and CT is a proportionality constant (Sreenivasan, 1996). We apply this model to estimate the mixing rates over the seismic profiles, applying a method proposed and described in previous works (e.g., Sheen et al., 2009; Holbrook et al., 2013). The main steps of the approach and the specifics of our work are described below.

Figure 3kρ(x,z) map obtained along the seismic profile indicated in Fig. 1, following the procedure explained in the text. White colored areas correspond to poorly sampled areas, with too few data to properly calculate kρ.

1. Select a local window larger than the resolution of the data but smaller that the entire seismic transect. The point is to select the smallest possible window that allows proper calculation of the reflector displacement spectra. We tested different window sizes and we found that the smallest ones that allow the production of robust results are 1200 m wide × 15 m high. Results with larger windows are comparable in terms of amplitude and shape of the imaged features, but structures and boundaries are better defined with this window size (Fig. S1). Smaller windows contain too few reflectors and produce abundant artefacts. Longer tracks are cut into shorter segments to fit inside the window. As explained in Sallares et al. (2016) and shown in Fig. S2, this does not affect the spectral values at the spatial scale range analyzed.

2. Compute the energy level of the displacement spectra in the turbulent subrange for all the reflectors inside the window, and calculate the average spectrum. The spectral subranges observed in the combined spectrum of the 68 reflectors (Fig. S3), which are also observed in most individual windows and reflectors, are used as a reference to select the scale range to compute the spectral amplitudes. In the case of the turbulent subrange, it is 13–30 m. The mean number of reflectors that fit inside the 1200 m wide × 15 m high windows is three, ranging from two to four depending on the imaged area.

3. Apply the Batchelor59 model (Eq. 7) to estimate ε, using the turbulent energy level (i.e., the average φς between 13 and 30 m) computed within the window, with Γ=0.2, CT= 0.3, and N calculated according to depth. Finally, we apply the Osborn80 relationship (Eq. 1) using the ε(x,z) values obtained above to compute ${k}_{\mathit{\rho }}\left(x,z\right).$

4. As only a few tracked reflectors are included inside each window, variances can affect the calculated diffusivities. To mitigate this effect, we slide the window in small steps (30 m in the horizontal and 3 m in the vertical direction), assigning the average value of the spectral amplitude to each local window. The fact that we incorporate a few new data at each step produces a smoothly varying map with a resolution that is similar to the window size ( 1000 × 10 m), instead of the one with sharp boundaries that is obtained without using overlapping windows (Fig. S1).

In summary, on the one hand we apply the Gregg89 model to obtain a vertical kρ(z) profile using the XCTD and ADCP data, and on the other hand we apply Batchelor59 to obtain a kρ(x,z) map using the vertical displacement spectra of the tracked reflectors. The results obtained using both methods and models are then compared to check if they are consistent, and to gain confidence in the HR-MCS methodology. We then analyze and discuss the high-resolution 2-D maps in terms of mixing.

Figure 4(a) Depth profile of ε(z) and kρ(z) obtained from XCTD and ADCP data and applying the Gregg89 model. The blue dotted line is the pelagic diffusivity in the ocean (kρ 10−5 m2 s−1), the green dotted line is the global average for overturning (kρ 10−4 m2 s−1), the red dotted line is the average vertical profile from XCTD and ADCP data (kρ 10−3.0 m2 s−1), and the gray area is the incidence range from MCS data (kρ 10−2.7 m2 s−1)(b) Turner angle showing ranges; the blue dotted lines show where the water column is unstable to diffusivity (Tu < 45) and stability (45 < Tu < 45) and prone to salt fingering (Tu > 45). (c) Buoyancy profile calculated with the XCTD data.

3 Results

The procedure described above allowed the production of a smoothly varying kρ(x,z) map that covers the whole profile (Fig. 3). The goal is being able to identify features and processes occurring in the transition between the internal wave and the turbulence subranges, such as the intensity and scales of variability in the mixing patches, the location and size of the mixing patches, and their potential relationship with oceanographic features such as IWs or wave instabilities. For this, we also use the vertical kρ(z) profile obtained from the XCTD and ADCP (Fig. 4a).

## 3.1 Probe-based kρ(z) profile

As we mentioned above, to have a reference value to compare with the MCS-based kρ(x,z) maps, we have first calculated a kρ(z) profile for shallow waters (< 200 m) using the XCTD and ADCP data and applying the Gregg89 model (Eqs. 4–6). To do this we have used ADCP measures averaged within 10 m depth bins. By doing this, we obtain an average value for the shear variance of ${S}_{\mathrm{10}}^{\mathrm{4}}=$ 0.28 s−4, whereas the reference value of the shear variance obtained from the Garrett–Munk model (Gregg, 1989) is ${S}_{\mathrm{GM}}^{\mathrm{4}}=$ 0.013 s−4. This gives an average dissipation rate < ε >   1.3 × 10−8 Wkg−1, and an average diapycnal diffusivity < kρ >  10−3.0 m2 s−1 for the targeted depth range (Fig. 4a). The kρ(z) profile obtained from the XCTD and ADCP is also shown in Fig. 4a, together with the global averages for overturning (< kρ >  10−4 m2 s−1) as well as the average pelagic diffusivity in the ocean (< kρ >  10−5 m2 s−1).

We obtain minimum values of the mixing rate at 50–55, 68–73, and 100–125 m. The absolute minimum of kρ= 10−5.2 m2 s−1 is obtained at  115 m, whereas the maximum is 10−2.1 m2 s−1 at  15 m. This gives a range of variation of 10−3.1 m2 s−1. Deeper than this, mixing variability is smaller. The Turner angle and buoyancy frequency (Fig. 4b) indicate that the region is mostly stable with a slight tendency to double diffusion (Tu  45).

It is worth noting that, at this specific location, the average vertical ε(z) and kρ(z) values are 1 order of magnitude higher that the global average ones. The higher values probably reflect the effect of overturning in the thermocline. While probe-based measurements are well suited to investigate mixing variability in the vertical dimension, they do not provide information on the variability in the horizontal dimension with a comparable level of detail. As explained above, to do this we have used estimations of ε and kρ based on the HR-MCS data, but applying the Batchelor59 model (Eq. 7) instead.

Figure 5High-resolution kρ(x,z) map overlapped with the HR-MCS image. Squares indicate location of some of the 1200 m × 15 m windows analyzed. They have been selected as examples of high-dissipation (windows W1–W3) and low-dissipation (windows W4–W6) areas. The color code of the squares is the same as for reflector spectra in Fig. 6, so that colors coincide with those of displacement spectra within the corresponding window.

Figure 6Average horizontal spectrum of the vertical displacements of reflectors inside windows W1–W6 (see location and color code in Fig. 5). (a) Spectra of individual reflectors in “high-diffusivity” areas (thin dotted lines), average within windows W1 (red solid line), W2 (yellow solid line), and W3 (orange solid line), and average of the three “high-diffusivity” windows (thick solid black line). (b) Spectrum of individual reflectors in “low-diffusivity” areas (thin dotted lines), average within windows W4 (magenta solid line), W5 (green solid line), and W6 (blue solid line), and average of the three “low-diffusivity” windows (thick solid black line). The reference lines are the theoretical slopes corresponding to the GM79 model for the internal wave subrange (brown dotted line), Kelvin–Helmholtz instabilities for the transitional subrange (dark blue dotted line), and Batchelor59 model for turbulence (dark green dotted line). Legend: values of diapycnal diffusivity using spectral values at the turbulent subrange within each of the analyzed windows (same color code as for windows W1–W6).

## 3.2 High-resolution multichannel seismic-based kρ(x,z) map

The kρ(x,z) map displayed in Fig. 3 has average values of < ε >  6.5 × 10−9 Wkg−1 and < kρ >  10−2.7 m2 s−1. These values are within the range of values obtained from the XCTD and ADCP data but, at the same time, they are over an order of magnitude higher than the global ocean reference value of kρ 10−4.0 m2 s−1 (Fig. 4a). Figure 5 displays the kρ(x,z) map superimposed with the HR-MCS data. It is worth noting that the range of horizontal variability is similar to that observed in the vertical dimension, although there is no direct visual correspondence between the kρ anomalies and IWs. The range of variability is over 3 orders of magnitude, locally reaching an extreme value of kρ 10−1.5 m2 s−1 at a depth of  55 m and at 16 km along the line, and a minimum value of ${k}_{\mathit{\rho }}\approx {\mathrm{10}}^{-\mathrm{4.5}}$ m2 s−1 at  95 m depth and 20 km along the line, which is close to the global oceanic average.

To better illustrate the procedure followed to generate the maps, several examples of kρ values obtained in “high” and “low” mixing areas, and the corresponding window average of the computed displacement spectra, are shown in Fig. 6. Numerous patches with kρ values exceeding 10−2 m2 s−1, with a characteristic size of 1–2 km in the horizontal dimension and  10 m in the vertical, are found throughout the whole section (i.e., the yellowish patches in Figs. 3 and 5). Not only the average depth value but also the vertical size of the anomalies, as well as the range of kρvariation, are in agreement with the probe-based values (Fig. 4). The contribution of the high kρ patches to the local average value is therefore outstanding, raising it from a background average value of  10−4 m2 s−1 to  10−2.7 m2 s−1.

Figure 7(a) High-resolution kρ(x,z) map overlapped with the HR-MCS image. Red lines indicate the locations of the two horizons analyzed (H1 and H2). They have been selected as examples of reflectors crossing higher mixing areas. (b) Diapycnal mixing obtained along H1 (see details of calculation in the text). (c) Signal filtered at wavelength ranges of the IW subrange (3000–100 m), (d) the transitional subrange (100–30 m), and (e) the turbulent subrange (30–13 m). The dashed red lines identify the “breaking point” referred to in the text. (f, g, h, i) Same as (b, c, d, e) but for H2.

## 3.3 Correspondence between mixing hotspots and imaged oceanographic features

To discuss the possible origin, or nature, of the mixing hotspots identified in the kρ(x,z) map (Fig. 3), we have visually compared the lateral variation in diapycnal diffusivity, with the structures imaged at the different subranges, along several individual reflectors. The analyzed reflectors have been selected as examples of high-diffusivity (H1 and H2 in Fig. 7) and low-diffusivity (H3 and H4 in Fig. 8) areas. To calculate kρ(x) along each horizon we have computed the spectral energy within a 1.2 km wide window moving laterally 30 m at each step along the whole reflector. To analyze the features that contribute to the energy spectrum in the different scales, and to compare them in turn with kρ(x), the horizons have been filtered at the scale ranges attributed to the IW (3000–100 m), transitional (100–30 m), and turbulent (30–13 m) subranges, respectively. As a reference, the local horizontal buoyancy wavelength estimated from the XCTD data is lN 90 m (Sallares et al., 2016). Although no general conclusions should be extracted from the analysis of a few individual reflectors, they show some relevant trends and correspondences to be taken into account when interpreting the results. In this sense, a clear trend that is observed in the displacement spectra is the systematic steep slope obtained at the transitional subrange between IWs and turbulence (Fig. S3). As explained in Sallares et al. (2016), this slope is consistent with numerical estimates for the evolutionary stage of the vortex sheet linked to shear instabilities (Waite, 2011), and it likely reflects the loss of energy in the wave field due to dissipation (e.g., Samodurov et al., 1995).

Figure 8(a) High-resolution kρ(x,z) map overlapped with the HR-MCS image. Blue lines indicate the locations of the two horizons analyzed (H3 and H4). They have been selected as examples of reflectors crossing lower mixing areas. (b) Diapycnal mixing obtained along H3 (see details of calculation in the text). (c) Signal filtered at wavelength ranges of the IW subrange (3000–100 m), (d) the transitional subrange (100–30 m), and (e) the turbulent subrange (30–13 m). The green circles identify the “diffusion peaks” referred to in the text. (f, g, h, i) Same as (b, c, d, e) but for H4.

Regarding horizons crossing high-dissipation areas (H1 and H2 in Fig. 7), a striking feature is the correspondence between the amplitude of the vertical displacements imaged in the transitional subrange and the variation in kρ. Hence, a variation in the amplitude of the features observed in the transitional subrange, at  34.7 km along horizon H1, and at  12.4 km along H2 (red lines in Fig. 7a), coincides with a decrease in kρ. In the case of H1, the average kρ value to the left of this point is 10−2.5 m2 s−1, while right of this point, it is 10−3.0 m2 s−1, whereas in the case of H2, the average kρ value to the left of this point is 10−4.1 m2 s−1, while right of this point, it is 10−2.9 m2 s−1. Although most of these values are higher than the average global value for meridional overturning circulation, the highest local average values are obtained in the region where the clearest, largest amplitude features, possibly representing KH billows (Sallares et al., 2016), are imaged. While we can identify a visual correspondence of high mixing values and the largest-amplitude features imaged in the transitional subrange, no direct correspondence is found with specific IWs, which are ubiquitous all along the profile.

For the low-dissipation areas (Fig. 8), we have selected H3, which is located at  35 m depth and has a length of  3.5 km (17.5–21 km along profile), and H4, located at  95 m depth and with a length of  4.0 km (32.5–36.5 km along profile). They were selected because their location coincides with a relatively weak mixing area, according to the kρ(x,z) map (Fig. 8a). As in the previous case, we have first calculated kρ using the spectral energy values within a 1.2 km wide window, moving laterally 30 m at each step along the whole reflector length. In the case of H3, the average value for the whole horizon is kρ 10−4.2 m2 s−1, considerably lower than in H1 but close to the global average value, whereas for H4, it is kρ 10−4.1 m2 s−1. In this case, we have identified some peaks at the transitional subrange that coincide with local highs in kρ(x). The peaks in H3 (at 18.4, 19.3, 19.9, and 20.4 km) and in H4 (at 33.5, 34.9, 35.3, and 35.8 km) display higher kρ values than the global ocean average. In particular, we obtain kρ 10−3.1 m2 s−1, kρ10−3.3 m2 s−1, kρ 10−3.4 m2 s−1, and kρ 10−3.5 m2 s−1 for each peak in H3, and kρ 10−3.0 m2 s−1, kρ 10−3.2 m2 s−1, kρ 10−3.1 m2 s−1, and kρ 10−2.9 m2 s−1 for each peak in H4. There are four more peaks at 19.5 and 20.1 km in H3 and 34 and 36.3  km in H4 that show no visual correspondence with structures in the transitional subrange but with larger amplitude features in the turbulent subrange; thus we hypothesize that they could be related to smaller-scale turbulent processes. The segments with no direct visual correspondence with kρ peaks, instead, display average kρ 10−4.3 m2 s−1, which is close to the global ocean average (Fig. 4a).

4 Discussion

The spatial variability observed along isopycnals based on the spectral analysis of the seismic data allows the identification of a number of local features at different evolutionary stages. These features are the manifestation of relevant oceanographic processes and structures, such as internal waves in the IW subrange, hydrodynamic instabilities in the transitional subrange, and turbulence at smaller scales. These processes are likely responsible for disruption of fine structure in the seismic image and the high-amplitude variability or the abrupt fading of some reflectors.

The large variations observed in the kρ(z) profile (Fig. 4), together with the slight tendency to double diffusion identified in the Turner angle, suggest that the system is prone to being affected by advection processes (e.g., Kunze and Sanford, 1996). Mixing appears to concentrate within the MAW, where the shear is strongest in the study area, and not deeper than > 110 m, where there is no significant shear and the system is weakly stratified. The shear-to-strain ratio calculated applying the Gregg89 model $\left({S}_{\mathrm{10}}^{\mathrm{4}}/{S}_{\mathrm{GM}}^{\mathrm{4}}=\mathrm{21}\right)$ indicates that the energy in the IW field is higher than that of the GM model, for which ${S}_{\mathrm{10}}^{\mathrm{4}}/{S}_{\mathrm{GM}}^{\mathrm{4}}$ is usually  3. We can therefore make the assumption that the energy is distributed in the whole inertial range where the water structures are stable (e.g., Munk, 1981). Similar results were obtained by Holbrook et al. (2013), who registered a shear-to-strain ratio of 17. The IWs can therefore be considered as an energy distributor from anisotropic to isotropic motions. The kρ value obtained from XCTD and ADCP using the Gregg89 model is kρ 10−3.0 m2 s−1, whereas we obtain kρ 10−2.7 m2 s−1 using MCS data and the Batchelor59 model. The ranges of variation in the two cases are also comparable, from maximum values of 10−2.2 and 10−1.5 m2 s−1 to minimum values of 10−5.4 and 10−5.7 m2 s−1 for the two methods. These similar values obtained based on different models and using independent techniques are well above the global average, suggesting that the downward energy cascade to small scales is highly efficient in the surveyed area.

We find no clear correlation between local kρ variations and the presence of individual IWs, which are clearly imaged and display a rather homogeneous distribution all along the line. Conversely, we find some hints of a direct relationship between changes in the amplitude of vertical isopycnal displacements and variations in kρ (Figs. 7 and 8), or between local peaks in the amplitude of vertical displacement and high kρ, in the transitional domain. Our interpretation is that IW-induced mixing is probably not efficient enough to sustain the overturning in the target area (Figs. 7 and 8). Instead, the correspondence between high-amplitude features in the shear-instability-dominated transitional domain with kρ suggests that the energy transfer between IWs and turbulence is probably enhanced by shear instabilities. The lack of clear correlation between IWs and turbulence agrees with the assumption of Klymak and Moum (2007a), suggesting a weak dependence of mixing rates on IW energy. Our observations indicate that the processes of IW destabilization and breaking appear to be important to efficiently allow the transfer of energy towards smaller scales and enhance mixing. In the case of the Alboran thermocline, the main mechanism could well be the development of shear instabilities, but other processes such as the interaction of IWs with rough bathymetry could also play an equivalent role in other regions and settings (Dickinson et al., 2017).

The mixing hotspots identified in the kρ(x,z) map (Fig. 3) likely represent a significant source of regional diapycnal mixing at the boundary layer between the MAW and the MW (30–200 m), which is subject to vertical stratification and shear values of 3.2 × 10−3 s−1. The mixing and energy transfer between these two water masses constitutes the main energy source of the region. In contrast to other regions (e.g., the Gulf of Mexico as in Dickinson et al., 2017), the smooth and relatively deep seafloor along the profile (> 800 m on average; Fig. S4c) suggests a small contribution of interaction with bathymetry in the generation of the mixing hotspots. Given that the MAW–MW boundary layer is subject to shear (Fig. S4a, b), and taking into account the visual correspondence between the location of the largest amplitude features in the transitional domain and high mixing values along individual reflectors, we hypothesize the existence of a direct link between mixing hotspots and IW shear instabilities. This could explain the high mixing variability throughout the surveyed area despite the ubiquitous presence of IWs.

The kρ values along H1 exceed the global average for overturning along most of the reflector (< kρ >  10−2.5 m2 s−1), with lower values only at some specific locations (Fig. 7b). The spatial correspondence between high diffusivity values and the presence of large-amplitude features at the transitional subrange, interpreted to correspond to KH-like shear instabilities (Sallares et al., 2016), is conceptually equivalent to the mechanism proposed by Gregg (1987), where mixing at the transitional subrange occurs principally at vortex sheets through wave instability. As an example, we hypothesize that in horizon H1, the presence of a vortex sheet left of  35 km produces the high mixing values, whereas to the right there is no vortex sheet and the ocean is more stable. The correspondence between diapycnal peaks and wave instabilities in horizon H3 suggests a similar situation (Fig. 8). Our results indicating a patchy ocean interior coincide with those previously presented by Sheen at al. (2009), Fortin et al. (2016), and Dickinson et al. (2017), but our work allows the extension of the conclusions to smaller-scale processes and shallower ocean levels. Additionally, we identify the development of shear instabilities as a likely relevant mechanism driving the downward energy cascade between IWs and turbulence at the thermocline depth that should be taken into consideration in ocean dynamic models.

5 Conclusions

We have used acoustic images obtained with a high-resolution MCS system to produce a 2-D diapycnal mixing map around the thermocline of the Alboran Sea. Our results confirm a high level of diapycnal variability and the presence of marked mixing hotspots in the water column. The kρ(x,z) map obtained by applying the Batchelor59 model to the seismic data has a strong variability with values ranging between < kρ >  10−1.5 m2 s−1, in the brightest hotspots, and < kρ >  10−3.3 m2 s−1, in the background. The obtained values are high enough to account for overturning at thermocline depths. The mixing hotspots have a characteristic size of 10–15 m in the vertical dimension and 1–2 km in the horizontal one, although there are also some smaller-sized ones. They are located at different depths within the thermohaline layer, although they appear to concentrate in highly sheared regions. The comparable values obtained with independent XCTD- and ADCP-based measures confirm that HR-MCS is a useful technique to study processes and structures occurring at the sub-mesoscale, which are difficult to capture and characterize by other means.

We investigate the relationship between mixing variability and ocean dynamics at different spatial scales by analyzing the spectral amplitudes along four seismic horizons in the internal waves and transitional, or instability-dominated, subranges. On the one hand, we found no clear correspondence between the location of the mixing patches and the location and amplitude of individual IWs, which are imaged all along the surveyed area. Conversely, a visual correspondence exists between the location of high-amplitude isopycnal vertical displacements at the instability-dominated transitional subrange and high mixing values in different reflectors, suggesting a causal relationship between both features. We interpret the development of shear instabilities as a mechanism that locally enhances downscale energy transfer between IWs and turbulence. Areas displaying the most vigorous instabilities coincide with the highest estimated diapycnal mixing values, which are well above the global average value for meridional overturning. This observation suggests that the energy transfer from anisotropic to isotropic scales is highly efficient at thermocline depths within the studied area.

Overall, our study shows that the HR-MCS technique can be used to study sub-mesoscale structures and processes at the thermocline level, provided that the stratification is strong enough to produce acoustic reflectivity that can be recorded by the system. The high-resolution 2-D maps produced from the seismic reflectivity could help improving estimates of the parameters to be incorporated in numerical models of ocean dynamics.

Data availability
Data availability.

Our research data can be accessed by request at the email vsallares@icm.csic.es.

Appendix A: Parameters used in text
Appendix B: Buoyancy Reynolds number

Gargett (1988) use an index to know if the system is isotropic or not, allowing us to determine if the buoyancy flux is vigorous enough to generate turbulence and therefore a high mixing level (Thorpe, 2005). The index depends on kinematic viscosity and is called the buoyancy Reynolds number:

$\begin{array}{}\text{(B1)}& {R}_{B}=\mathit{\epsilon }/\mathit{\nu }{N}^{\mathrm{2}}.\end{array}$

The mean kinematic viscosity in the ocean is ν= 1 × 10−6 m2 s−1. Some properties of the inertial subrange are consistent with isotropy for values of RB < O(102). To consider anisotropy and avoid serious underestimates of mixing, Smyth and Moum (2000) propose RB > 200 as safe for high mixing levels due to free viscous effects. For our sub-mesoscale regime, RB= 3200, a value that is compatible with the calculated mixing levels.

Supplement
Supplement.

Competing interests
Competing interests.

The authors declare that they have no conflict of interest.

Acknowledgements
Acknowledgements.

This work has been fulfilled in the framework of projects POSEIDON (ref: CTM2010-25169) and FRAME (ref: CTM2015-71766-R), both funded by the Spanish Ministry of Economy and Competitiveness (MINECO). The seismic and oceanographic data were acquired in the framework of the IMPULS survey (ref: 2003-05996-MAR), also from MINECO, and the SAGAS survey (ref: CTM2005-08071-C03-02/MAR-SAGAS). This work was also supported by the Grup de Recerca de la Generalitat de Catalunya: Barcelona Center for Subsurface Imaging (2014SGR940). Helpful comments were provided by Josep Lluis Pelegrí, Miguel Bruno, numerous colleagues at the Barcelona-CSI, and Diana Francis and David M. Holland from the Center for Global Sea Level Change (CSLC) – NYUAD, Abu Dhabi, UAE.

Edited by: John M. Huthnance
Reviewed by: two anonymous referees

References

Batchelor, G. K.: Small-scale variation of convected quantities like temperature in turbulent fluid, Fluid Mech., 5, 113–133, 1959.

Biescas B., Ruddick, B. R., Nedimovic, M. R., Sallarès, V., Bornstein, G ., and Mojica, J. F.: Recovery of temperature, salinity, and potential density from ocean reflectivity, J. Geophys. Res. Oceans, 119, 3171–3184, https://doi.org/10.1002/2013JC009662, 2014.

Chioua, J., Bruno, M., Vazquez, A., Reyes, M., Gomiz, J., Mañanes, R., Alvarez, O., Gonzalez, C., Lopez, L., and Gomez-Enri, J.: Internal waves in the strait of Gibraltar and their role in the vertical mixing processes within the bay of Algeciras, Estuarine, Costal and Shelf Science, Elsevier, 126, 70–86, 2013.

D'Asaro, E. A. and Lien, R.-C.: The Wave-turbulence transition for stratified flows, J. Phys. Oceanogr., 30, 1669–1678, https://doi.org/10.1175/1520-0485(2000)030<1669:TWTTFS>2.0.CO;2, 2000.

Dickinson, A., White, N. J., and Caulfield, C. P.: Spatial variation of diapycnal diffusivity estimated from seismic imaging of internal wave field, Gulf of Mexico, J. Geophys. Res.-Oceans, 122, 9827–9854. https://doi.org/10.1002/2017JC013352, 2017.

Fortin, W. F. J., Holbrook, W. S., and Schmitt, R. W.: Mapping turbulent diffusivity associated with oceanic internal lee waves offshore Costa Rica, Ocean Sci., 12, 601–612, https://doi.org/10.5194/os-12-601-2016, 2016.

Gargett, A.: The scaling of turbulence in the presence of stable stratification, J. Geophys. Res., 93, 5021–5036, 1988.

Gargett, A. and Holloway, G.: Dissipation and diffusion by internal wave breaking, J. Mar. Res., 42, 15–27, 1984.

Garrett, C. and Munk, W.: Internal waves in the ocean, Ann. Rev. Fluid Mech., 11, 339–369, 1979.

Gibson, C., Keeler, R., Bondur, V., Leung, P., Prandke, H., and Vithanage, D.: Submerged turbulence detection with optical satellites, Coastal Ocean Remote Sensing Conf., Paper 6680-33, San Diego, CA, 2007.

Gregg, M. C.: Diapycnal mixing in the thermocline: A review, J. Geophys. Res., 92, 5249–5286, 1987.

Gregg, M. C.: Scaling turbulent dissipation in the thermocline, J. Geophys. Res., 94, 9686–9698, 1989.

Holbrook, W. and Fer, I.: Ocean internal wave spectra inferred from seismic reflection transects, Geophys. Res. Lett., 32, L15604, https://doi.org/10.1029/2005GL023733, 2005.

Holbrook, W. S., Paramo, P., Pearse, S., and Schmitt R.: Thermohaline fine structure in an oceanographic front from seismic reflection profiling, Science, 301, 821–824, 2003.

Holbrook, S., Fer, I., Schmitt, R., Lizarralde, D., Klymak, J., Helfrich, C., and Kubichek, R.: Estimating oceanic turbulence dissipation from seismic images, J. Atmos. Ocean. Tech., 30, 1767–1788, https://doi.org/10.1175/JTECH-D-12-00140.1, 2013.

Klymak, J. M. and Moum, J. N.: Oceanic isopycnal slope spectra. Part a: Internal waves, J. Phys. Oceanogr., 37, 1215–1231, 2007a.

Klymak, J. M. and Moum, J. N.: Oceanic isopycnal slope spectra. Part b: Turbulence, J. Phys. Oceanogr., 37, 1232–1244, 2007b.

Kunze, E. and Sandford, T. B.: Abyssal Mixing: Where it is not, J. Phys. Oceanogr., 26, 2286–2296, 1996.

Munk, W. H.: A survey of internal waves and small-scale processes, Evolution of Physical Oceanography, edited by: Warren, B. A. and Wunsch, C., 264–291, MIT Press, Cambridge, Mass, 1981.

Munk, W. and Wunsch, C.: Abyssal recipes II: energetics of tidal and wind mixing, Deep-Sea Res. Pt. I, 45, 1977–2010, 1998.

Müller, P. and Pujalet, R.: Internal gravity waves and small scale turbulence, proceeding, “Aha Huliko” a Hawaiian winter workshop, Hawaiian institute of geophysics, Special Publications, 299 pp., 1984.

Osborn, T. R.: Estimates of the local rate of vertical diffusion from dissipation measurements, J. Phys. Oceanogr., 10, 83–89, 1980.

Osborn, T. and Cox, C. S.: Oceanic fine structure, Geophys. Fluid Dyn., 3, 321–345, 1972.

Sallarès, V., Mojica, J. F., Biescas, B., Klaeschen, D., and Gracia, E.: Characterization of the submesoscale energy cascade in the Alboran Sea thermocline from spectral analysis of high-resolution MCS data, Geophys. Res. Lett., 43, 6461–6468, https://doi.org/10.1002/2016GL069782, 2016.

Samodurov, A. S., Lubitsky, A. A., and Panteleev, N. A.: Contribution of breaking internal waves to structure formation, energy dissipation, and vertical diffusion in the ocean, Phys. Oceanogr., 6, 177–190, 1995.

Sheen, L. K., White, J. N., and Hobbs, R. W.: Estimating mixing rates seismic images of oceanic structure, Geophys. Res. Lett., 36, L00D04, https://doi.org/10.1029/2009GL040106, 2009.

Smyth, W. D. and Moum, J. N.: Anisotropy of turbulence in stably stratified mixing layers, Phys. Fluids, 12, 1343–1362, 2000.

Smyth, W. D., Moum, J. N., and Nash J. D.: Narrowband oscillations in the upper equatorial ocean. Part II:Properties of shear instabilities, J. Phys. Oceanogr., 41, 412–428, 2011.

Sreenivasan, K.: The passive scalar spectrum and the Obukhov-Corrsin constant, Phys. Fluids, 8, 189–196, 1996.

Stewart, R.: Introduction to Physical Oceanography, Department of Oceanography, Texas A&M University, 2008.

Thorpe, S. A.: Experiments on the instability of stratified shear flows: miscible fluids, J. Fluid Mech., 46, 299–319, 1971.

Thorpe, S. A.: The turbulent ocean, Cambridge University Press, Cambridge, 2005.

Waite, M. L.: Stratified turbulence at the buoyancy scale, Phys. Fluids, 23, 066602, https://doi.org/10.1063/1.3599699, 2011.

Waterman, S., Naveira Garabato, A. C., and Polzin, K. L.: Onternal waves and turbulence in the Antarctic Circumpolar Current, J. Phys. Oceanogr., 43, 259–282, https://doi.org/10.1175/JPO-D-11-0194.1, 2013.