Microstructure observations and mixing parameterizations along an Atlantic transect in very weak turbulence

Microstructure measurements of shear and temperature can be used to calculate ocean turbulent dissipation rates and diffusivities. Here microstructure observations are taken along an transect in the North Atlantic, that includes observations of very weak deep ocean turbulence. In this paper we show the necessity of using the thermistor probes, instead of the more common shear probes, to calculate dissipation rates when they are smaller than $\mathrm{1}\times {\mathrm{10}}^{-\mathrm{10}}$ W kg⁻¹. Profiles of combined dissipation rates from the shear and thermistor probes are then compared to the finescale strain parameterization and Thorpe sorting method. Based on this comparison, recommendations and restrictions are suggested for applying both parameterizations in a weakly turbulent environment. The results indicate that temperature-based strain provides improved estimates of dissipation rates in the deep ocean where density gradients are small, while density-based strain provides better results otherwise. We find that Thorpe based estimates are very accurate when pre-existing knowledge of the turbulent kinetic energy dissipation rate ε is used. When this knowledge is not available, using climatological mean estimates of ε can allow for more detailed estimates of dissipation by applying the Thorpe resorting method. Finally, we employ the triple decomposition framework to get more insights in the relative roles of dianeutral and isoneutral mixing processes, and use this to calculate the dianeutral and isoneutral diffusivities. It turns out that the triple decomposition is generally not a good predictor of the isoneutral diffusivity. Overall, this paper has assessed the potential of direct observations and parameterizations of dissipation and showed that dissipation rates can be estimated quite well within a factor 5 between different methods, but it becomes difficult to achieve higher accuracy.

1 Introduction

Dianeutral mixing – mixing across layers of different water density – is caused by a wide variety of processes, such as winds, tides and geostrophic currents (MacKinnon et al., 2017; Whalen et al., 2020). It plays an important role in circulation (de Lavergne et al., 2022) – involved in the closure of the Meridional Overturning Circulation (Wunsch and Ferrari, 2004; Melet et al., 2022; Cimoli et al., 2023) –, the transport and distribution of tracers in the ocean such as heat, carbon (Tatebe et al., 2018) or contaminants and nutrients (Friedrich et al., 2011; Spingys et al., 2021), eventually impacting regional and global climate (Pradal and Gnanadesikan, 2014). Away from surface and bottom boundary layers, in the ocean interior, mixing is sustained by a field of internal waves (MacKinnon et al., 2017). The breaking of these internal waves sets off an energy cascade down to the smaller scales that eventually leads to irreversible energy dissipation and turbulent tracer mixing. The rate of energy dissipation can be recalculated into a mixing strength or diffusivity (Osborn, 1980), which can subsequently be used to study the role of mixing in the ocean and climate system through for example numerical simulations (Melet et al., 2013).

Observational studies have shown that the dianeutral mixing strength is highly variable throughout the ocean, spanning several orders of magnitude (Polzin et al., 1997; Naveira Garabato et al., 2004). Weak diffusivities of 𝒪(10⁻⁵) m² s⁻¹ have been observed in the ocean interior (Munk and Wunsch, 1998), while stronger mixing rates (𝒪(10⁻³) m² s⁻¹) are found closer to boundaries and over rough topography (Polzin et al., 1997; Waterhouse et al., 2014; Wynne-Cattanach et al., 2024). The patchy and intermittent nature of dianeutral mixing makes it a difficult quantity to measure. Meanwhile, values found in the quiescent ocean interior can be extremely low, adding difficulty to obtaining such observations because of the high instrument precision required for measuring such low values (Scheifele et al., 2018).

Obtaining observations of the turbulent dissipation rate ε, and subsequently the diffusivity, requires the resolution of the observations to be 𝒪(mm-cm) in order to resolve the velocity shear variance at the scales of the turbulent inertial subrange. This can be achieved for example with (vertical) microstructure profilers. These profilers are usually equipped with a combination of airfoil shear probes, measuring vertical velocity shear $\left(\frac{\partial u}{\partial z}\right)$, and fast-response thermistors (FP07), measuring the temperature gradient $\left(\frac{\partial T}{\partial z}\right)$. This allows for measurements of the turbulent kinetic energy (TKE) dissipation rate ε and the temperature variance dissipation rate χ, directly at the scales of dissipation. The specialized nature of this type of instruments cause these kind of observations to remain scarce, especially when covering the full ocean depth (Waterhouse et al., 2014). Consequently, alternative methods have been developed to indirectly estimate dissipation and mixing strengths from more readily available data such as ship-based Conductivity-Temperature-Depth (CTD) probes, moorings, and autonomous platforms (e.g., gliders or Argo floats). Some examples are the finescale parameterization (Gregg, 1989; Polzin et al., 1995, 2014), the Thorpe resorting method (Thorpe, 1977; Dillon, 1982), inverse modelling (Ganachaud and Wunsch, 2000; Sloyan and Rintoul, 2000; Zika et al., 2010; Groeskamp et al., 2017; Hautala, 2018; Kusters et al., 2024), and internal wave ray tracing (de Lavergne et al., 2020). However, the application of these parameterizations relies heavily on a variety of assumptions, adding uncertainty to estimates. The combination of observations and indirect estimates of dissipation and turbulent mixing strength, allows for calibration, fine-tuning and checking of the indirect estimates. The indirect estimates often provide possibilities to estimate mixing at spatial and temporal scales that are impossible for microscale observations. Hence combining both approaches greatly improves our understanding of turbulent mixing processes and our ability to include these effects in numerical simulations and other studies (Fox-Kemper et al., 2019).

The goals of this study are to (1) compare direct microstructure observations from both shear probes and thermistor probes with each other, and (2) compare these direct estimates to indirect estimates obtained using the finescale parameterization and the Thorpe resorting method, both derived from finescale temperature and salinity data collected with a profiler-mounted CTD, and finally (3) understand if the triple decomposition framework (Ferrari and Polzin, 2005; Naveira Garabato et al., 2016; Castro et al., 2024) can give insight for the interpretation of the results.

Various studies have reported on comparisons between microstructure observations and one or more parameterizations for the dissipation rate. These studies were mostly situated in either energetic locations such as in the Antarctic Circumpolar Current (ACC) (Thompson et al., 2007; Waterman et al., 2013; Frants et al., 2013; Park et al., 2014), in the proximity of rough topography (Ferron et al., 1998) or in the upper parts of the ocean (Whalen et al., 2015; Howatt et al., 2021; Roget et al., 2023), i.e. locations where sensor limitations are of lesser concern. Few studies do make a comparison between direct measurements of ε and one or more parameterizations in a low energetic environment, such as the Arctic Ocean (Fine et al., 2021; Baumann et al., 2023). However, these studies often only use shear-based observations of ε and are thus unable to consider values below the shear noise floor. Whereas for the few studies that are concerned with thermistor-based ε estimates in the deep sea (e.g., Scheifele et al., 2018; Yasuda et al., 2021), a comparison with parameterizations is not within the scope of those studies.

This study uses data from two research cruises across the subtropical North Atlantic (Fig. 1) as part of the Mixation and Nanoplastics 2 projects. For large parts of this dataset, there is weak stratification and (very) low dissipation rates, often below the noise floor of the more conventionally used shear probes, but yet within the range of the thermistor probes. Here we assess the difference between the two types of probes on the microstructure profiler, and we compare these microstructure measurements to two parameterizations – a finescale strain parameterization and a Thorpe scale-based method –, in this dataset with weak turbulence. Additionally, we explore practical improvements to such indirect methods, e.g. the use of strain based on temperature gradients for the use in the finescale parameterization and objective criteria to exclude false overturns in the Thorpe resorting method.

The paper is structured as follows. In Sect. 3, the limits of both the shear and thermistor probes will be explored, showcasing the need for using thermistor probes instead of shear probes in areas of the deep sea where turbulent dissipation rates are below $\mathrm{1}\times {\mathrm{10}}^{-\mathrm{10}}$ W kg⁻¹. In Sect. 4, the dissipation rates from the microstructure data are compared to the finescale-strain parameterization that uses strain based on temperature gradients and the buoyancy frequency from the CTD data. We showcase that in the deepsea, the temperature-based finescale parameterization is preferred over the more commonly used density based finescale parameterizations, due to higher instrumental noise in the salinity probe. In Sect. 5, the CTD-temperature data is used at smaller scales for the Thorpe resorting method. Here we find that Thorpe parameterization is able to reproduce the dissipation rates well, also in the weakly turbulent regions, but that these results are sensitive to methodological choices. We then apply the triple decomposition framework (Ferrari and Polzin, 2005; Naveira Garabato et al., 2016; Castro et al., 2024) to the collected CTD and mircostructure data (Sect. 6). This framework provides further insight into the relative roles of isoneutral and dianeutral processes in setting the observed T-S relations that impact the used parameterizations. These observations also allows for estimating the dianeutral diffusivities and mixing efficiency (Sect. 7). Finally we end with a discussion and conclusions (Sect. 8).

Figure 1 Map with the locations of the stations of the Mixation 2 project (R/V Atlantic Explorer, August 2023, red circles) southeast of Bermuda, and the stations of the Nanoplastics 2 project (R/V Pelagia, November–December 2023, orange triangles), roughly along the transect along 30° N.

2 Data

In this study we use a set of microstructure profiles collected during the Mixation and Nanoplastics 2 cruises (Fig. 1). A total of 12 full-depth microstructure profiles (VMP-6000, Rockland Scientific International Inc) with standard (pumped) CTD data (SeaBird Electronics SBE 3/4/5, 64 Hz) were collected to approximately 30 m above the seafloor. Except st. NP3, where the VMP experienced an early release approximately 700 m above the seafloor. Additional minor problems with the FP07 thermistors were encountered at some stations, such as anomalous spikes likely caused by collisions or general malfunctioning of one of the probes. These problems were taken care of in the post-processing and at any time, at least one of the probes was functioning as intended. The stations of the Mixation cruise are referred to as M1, M2, M3 and MM. The stations of the Nanoplastics cruise will be called NP1–NP8. An extensive description of all used methods is provided in the appendices.

All (in-situ) temperature and conductivity data from the CTD systems were converted to Conservative Temperature (Θ) and Absolute Salinity (S_A) using the GSW software toolbox (IOC et al., 2010; McDougall and Barker, 2011) and Neutral Density γⁿ (Jackett and McDougall, 1997). Conservative Temperature Θ (°C) is proportional to potential enthalpy (by the constant heat capacity factor $c_p^{\mathrm{0}}$ (J kg⁻¹ K⁻¹), representing the heat content per unit mass of seawater (McDougall, 2003; Graham and McDougall, 2013). Absolute Salinity S_A is designed to approximate the ratio between the mass of dissolved material and the mass of seawater (g kg⁻¹) (Wright et al., 2011; McDougall et al., 2012). It is measured on the reference-composition salinity scale (Millero et al., 2008).

In this manuscript various methods to estimate the dissipation rate are being used. To avoid confusion in the comparison of the estimates for ε we will use the following subscripts: _µU will be used to indicate variables linked to the microstructure shear data, similarly, _µT will be used for the data belonging to the microstructure thermistor data. The subscript _FS will be used for the finescale strain method and _TP will be used for the Thorpe resorting.

3 Energy dissipation rate from microstructure measurements

The microstructure measurements made at each station (Fig. 2) consist of observations made by two types of sensors, the airfoil shear-probes and the fast-response thermistors. Both type of sensors measured at a rate of 512 Hz and provide a way to obtain the turbulent dissipation rate of kinetic energy ε and thermal variance χ_Θ. ε values were obtained from shear as well as temperature measurements. ε_µU and χ_Θ are obtained from spectral integration of the microstructure velocity and temperature data (Kraichnan, 1968; Nasmyth, 1970; Oakey, 1982). ε_µT follows from fitting the temperature gradient spectrum to a theoretical Kraichnan spectrum and the fitted Batchelor wavenumber k_B (Kraichnan, 1968; Ruddick et al., 2000). A full description of the processing of the microstructure data and the applied quality control checks is provided in Appendix A. We present results from both methods and compare them, showing that for a large part of the deep ocean, shear-based estimates may overestimate dissipation rates because the actual dissipation is lower than their noise floor. Such observation motivates us to investigate the range of turbulence levels that can be resolved by each type of sensors, and discuss how shear and temperature-based ε estimates may be combined to obtain a single profile of energy dissipation rates.

Figure 2 Dissipation rates obtained from the microstructure shear probes (in blue) and the microstructure thermistor probes (in red) for each station. The background shading marks the 95 % confidence interval. A 5-point running mean has been applied to increase readability.

3.1 Determining the noise floor for shear-probes

For most of the ocean, e.g. more energetic environments and the upper ocean, there is in general a good agreement between energy dissipation rates measured by the shear probes (ε_µU) and thermistors (ε_µT), usually within a factor of 2 (Oakey, 1982; Scheifele et al., 2018; Piccolroaz et al., 2021; Yasuda et al., 2021). In less energetic environments, e.g. the abyssal ocean over smooth topography, there is often a larger discrepancy between the two types of probes (Fig. 2). These discrepancies are particularly noticeable in the deep parts of station NP3, below 2000 m, where ε_µU saturates at $\sim {\mathrm{10}}^{-\mathrm{10}}$ W kg⁻¹ whilst ε_T reaches values as low as $\sim {\mathrm{10}}^{-\mathrm{12}}$ W kg⁻¹.

Such differences are clearly reflected in the ε probability distributions. While the dissipation rates from the thermistors have the expected nearly-symmetric and approximately log-normal distribution found for ocean turbulence (Cael and Mashayek, 2021; Gregg et al., 1993), the distribution of the shear dissipation rates is more heavily skewed (Fig. 3a, c). The skewness of the distribution of ε_µU is here attributed to the noise floor of the shear probes. In a low-energetic environment as this, the presence of the noise floor can bias the estimated dissipation towards larger values, by removing the left tail of the distribution. Often, the noise floor for shear measurements is reported to be 𝒪(10⁻¹⁰) W kg⁻¹ (e.g., Fer et al., 2014) and is attributed to contamination of the shear spectrum by electronic noise of the instrument. The exact magnitude of the noise floor varies not only between instruments, platforms and manufacturers, but also between cruises due to for example small differences in instrument set-up (which is why for Fig. 3 the distributions are plotted separately for each cruise).

Figure 3 (a) Distributions of the ε_µU and ε_µT data from the Nanoplastics II stations. (b) Cross-validation between ε_µU and ε_µT for the Nanoplastics II data. Panels (c), (d) are the same as (a), (b) but for the Mixation II data. Panels (b), (d) the solid black line notes the 1:1 ratio. Other black lines mark a factor 2 and 3 difference between the sensors. Shear probe noise floors according to Table 1 are marked in red. A clear break in the agreement can be seen around the level of the noise floor of ε_µU.

We now explore various methods for defining the noise floor of shear-derived ε_µU for each dataset. Scheifele et al. (2018) suggests using the mode of ε_µU as an indicator for the noise floor, whereas Piccolroaz et al. (2021) suggests taking the lower 5th percentile. The noise floors derived from these methods for both cruise datasets are summarized in Table 1 and marked in Fig. 3b, d. However, because of the large difference between the estimates of dissipation rates from the fast-response thermistor and the shear probes for station NP3 (below 2000 m, Fig. 2), we have a unique opportunity to provide another estimate of the noise floor for shear probes. The open-ocean turbulent dissipation rates are so low in this region, that the estimates of dissipation by the shear probe are about an order of magnitude larger than the dissipation rates measured by the thermistors, and are therefore reflective of the noise in the instrument. We can thus estimate the noise floor empirically by taking the average and standard deviation over part of such noise-dominated region, where the disagreement between the ε_µU and ε_µT is largest and consistent (i.e. st. NP3 between 2000–4200 dbar). The (arithmetic) mean over this section is $\mathrm{4.5}\times {\mathrm{10}}^{-\mathrm{11}}\pm \mathrm{1.9}\times {\mathrm{10}}^{-\mathrm{11}}$ W kg⁻¹. Table 1 shows that the empirical estimate is close to the methods used by Scheifele et al. (2018) and Piccolroaz et al. (2021), ranging 1.4–$\mathrm{6.3}\times {\mathrm{10}}^{-\mathrm{11}}$ W kg⁻¹. It is noted that the Nanoplastics stations have higher noise levels than the Mixation dataset, even though the same instrument was used and that this empirical noise estimate is rather conservative for the Mixation cruise, given the solid red line in Fig. 3d. However, the noise floor of a sensor is not a hard threshold, it represents the range where the signal is dominated by an increasing amount of noise. Instead, to be sure the signal contains very little noise, we suggest to set the threshold for discarding shear-based dissipation rates at the mean plus two standard deviations, which equals $\mathrm{8.3}\times {\mathrm{10}}^{-\mathrm{11}}$ W kg⁻¹. When also accounting for other environmental and instrument-specific factors influencing uncertainties mentioned above, a value of $\mathrm{1}\times {\mathrm{10}}^{-\mathrm{10}}$ W kg⁻¹, as often found in literature, would equal the mean plus 3 standard deviations and would put a user absolutely at the save side of knowing the results are not dominated by noise.

Scheifele et al. (2018)Piccolroaz et al. (2021)

Table 1 Noise floor of ε_µU by method and split between datasets. Values are expressed in 10⁻¹¹ W kg⁻¹.

3.2 Determining limits for the fast-response thermistors

The thermal variance dissipation rate and the Batchelor wave number were obtained by fitting the observed microscale temperature-gradient spectrum to a theoretical model spectrum (Kraichnan, 1968; Ruddick et al., 2000) (see Appendix A for a more extensive description).

Then, ε_µT is derived from the fitted Batchelor wave number (${\mathit{\epsilon }}_{\mathrm{\mu }\mathrm{T}}=\mathit{\nu }{\mathit{\kappa }}_{\mathrm{T}}^{\mathrm{2}}{k}_{\mathrm{B}}^{\mathrm{4}}$), which signals the position of the spectral roll-off at high wave-numbers (Ruddick et al., 2000). Examination of the ε profiles and probability density function suggest that the noise level of the thermistor, which is known to be much lower than that of the shear probe (Kocsis et al., 1999; Scheifele et al., 2018), has not significantly impacted our observations. In our processing routines we discard any data segment where the average spectral levels fall below a factor of 1.3 (Piccolroaz et al., 2021) of the well-characterised sensor noise curve (following RSI, see Technical Note 40 available at http://www.rocklandscientific.com, last access: 23 October 2025). Only 0.49 % of data-segments were discarded following that criterion. For the retained data, the distributions of ε_µT do not show a strong deviation from the expected log-normal distribution (Fig. 3a, c), as found for the shear sensor, suggesting that these estimates are still reliable and above the noise floor.

Of bigger concern for the thermistor-based estimates are an upper limit due to the limited time-response of the thermistor (5–10 ms), which smooths out high-frequency temperature fluctuations. Whilst the limited, speed-dependent time-response is corrected for using a double-pole transfer function (Vachon and Lueck, 1984), a harder high-wavenumber cut-off is imposed by the fact that an anti-aliasing filter is applied to the microstructure temperature data. The filter is applied at $\mathrm{1}/\mathrm{5}$th of the sampling frequency (512 Hz), at ∼100 Hz The filter causes a rapid roll-off of the temperature gradient spectrum, similar to the physical roll-off occurring around the Batchelor wave number k_B. The anti-aliasing will cause the observed spectra to display an unphysical roll-off at a wavenumber of $k_{\mathrm{A}}=f_{\mathrm{A}}/W=(\mathrm{0.2}\times \mathrm{512})/\mathrm{0.8}=\mathrm{128}$ cpm, where W=0.8 m s⁻¹ is the profiling speed. When k_B is smaller than k_A, the spectrum will have two roll-offs, of which one is physical at the Batchelor wave number, and the second is instrumental at k_A≈128 cpm. In this case, the fitting algorithm will capture the physical roll-off as it should. For the case when k_B>k_A, the fitting algorithm will fit to the instrumental roll-off, rather than the physical roll-off and the resulting k_B estimate will underestimate the actual k_B. Assuming that the highest k_B that can be correctly used is k_B≈k_A, this would result in an upper limit for the thermistor estimates at ${\mathit{\epsilon }}_{\mathrm{\mu }\mathrm{T},max}=\mathit{\nu }\times {\mathit{\kappa }}_{\mathrm{T}}^{\mathrm{2}}\times (\mathrm{2}\mathit{\pi }{k}_{\mathrm{B}}{)}^{\mathrm{4}}=\mathcal{O}(\mathrm{1}\times {\mathrm{10}}^{-\mathrm{8}})$ W kg⁻¹. With κ_T being the thermal molecular diffusivity.

3.3 Combining shear-based and thermistor-based dissipation rates

The discussion above shows that the shear probes perform well above the noise floor of $\mathrm{8.3}\times {\mathrm{10}}^{-\mathrm{11}}$ W kg⁻¹, whereas the thermistor probes provide good estimates of the dissipation rate below $\mathrm{1}\times {\mathrm{10}}^{-\mathrm{8}}$ W kg⁻¹ (Fig. 3). The temperature-based dissipation rates ε_µT show that for large parts of the open ocean, the actual dissipation rate can be up to two order magnitude smaller than the noise floor of shear probes. This highlights that awareness of the sensor limitations is important for using the right sensor to get reliable results. Something that is often overlooked and suggests that shear probes and thermistor probes should be combined to get the optimal result in regions of weak turbulence.

For combining both instruments, we note that the shear probes provide the most direct method for obtaining ε, while estimating the dissipation rate from the thermistor probes requires more assumptions and is more indirect. Hence, when possible, the shear-based estimates are preferred over thermistor-based estimates. We will therefore use shear-probe estimates above the empirically determined threshold of $\mathrm{8.3}\times {\mathrm{10}}^{-\mathrm{11}}$ W kg⁻¹. If estimates are below this, thermistor-based dissipation rates are used instead. This results in a combined profile of dissipation rates, referred to as ε_vmp, for each station. We will use ε_vmp when comparing the finescale and Thorpe parameterizations with microstructure data in the remainder of the manuscript.

4 Finescale parameterization

We here apply the finescale parameterization to estimate ε. The finescale parameterization assumes that the observed variance in the strain or shear spectra at the fine scale 𝒪(10–100) m are predominantly caused by the breaking of internal gravity waves, and that the dissipation rate of turbulent kinetic energy follows from nonlinear wave-wave interactions that transfer energy to smaller scales, that ultimately break into turbulence (Thompson et al., 2007; Whalen et al., 2015; Dematteis et al., 2024). Such estimates are obtained using vertical gradients of properties such as Θ, that can be obtained from standard CTD data. This method can be either applied to shear, strain or shear and strain data (Gregg, 1989; Polzin et al., 1995, 2014). Since large-scale shear data is absent in this dataset, we consider two calculations of the finescale strain-only method, with the strain rate based on the buoyancy frequency N² and vertical temperature gradients Θ_z. The dissipation rate ε_FS is obtained via (Whalen et al., 2015)

$$\begin{array}{}\text{(1)}& {\mathit{\epsilon }}_{\mathrm{FS}}={\mathit{\epsilon }}_{\mathrm{0}}{\displaystyle \frac{\stackrel{\mathrm{‾}}{N^{\mathrm{2}}}}{N_{\mathrm{0}}^{\mathrm{2}}}}{\displaystyle \frac{\langle {\mathit{\xi }}_z^{\mathrm{2}}{\rangle }^{\mathrm{2}}}{\langle {\mathit{\xi }}_{z\mathrm{GM}}^{\mathrm{2}}{\rangle }^{\mathrm{2}}}}h\left(R_{\mathit{\omega }}\right)L\left(f,N\right).\end{array}$$

Here $\langle {\mathit{\xi }}_z^{\mathrm{2}}{\rangle }^{\mathrm{2}}$ is the strain variance and $\stackrel{\mathrm{‾}}{N^{\mathrm{2}}}$ the segment-avaraged buoyancy frequency, obtained via the method of adiabatic levelling (Bray and Fofonoff, 1981). The strain variance is compared to a theoretical Garrett-Munk (GM) reference spectrum $\langle {\mathit{\xi }}_{\mathrm{GM}}^{\mathrm{2}}{\rangle }^{\mathrm{2}}$ (Garrett and Munk, 1975; Cairns and Williams, 1976). Furthermore, ${\mathit{\epsilon }}_{\mathrm{0}}=\mathrm{6.73}\times {\mathrm{10}}^{-\mathrm{10}}$ m² s⁻² and $N_{\mathrm{0}}=\mathrm{5.24}\times {\mathrm{10}}^{-\mathrm{3}}$ rad s⁻¹ are constants (Whalen et al., 2015), and L(f,N) and h(R_ω) are correction functions for latitudinal dependences and the shear-to-strain ratio (Polzin et al., 1995; Gregg et al., 2003). The reader is referred to Appendix B for a full description of these functions.

4.1 The shear-to-strain ratio R_ω

In order to get estimates using the finescale strain parameterization, in the absence of shear data, a choice needs to be made for the shear-to-strain ratio R_ω. Different studies have used available shear data to calculate this value and reported a large range of values for R_ω, both in space and time, ranging from 1–50 (Waterman et al., 2013; Chinn et al., 2016; Fine et al., 2021), while Kunze et al. (2006) reports an average oceanic value of R_ω=7. R_ω for the GM model spectrum is R_ω=3, which corresponds to h(R_ω)=1 (Eq. B4). It can be argued that the best choice to gain optimal results would be to use a variable R_ω (e.g Sun et al., 2024), as the dissipation rates (and thus the error) are directly proportional to the choice of R_ω. Though in practice, often is chosen for either R_ω=3 or R_ω=7. With our dataset a direct comparison with microstructure data can be made and so the better choice for R_ω can be assessed. To this end, the results for three different choices of R_ω are compared to the ratio ${\mathit{\epsilon }}_{\mathrm{FS}}/{\mathit{\epsilon }}_{\mathrm{vmp}}$, for both calculations of the strain rate (Fig. 4). Estimates of ε_vmp have been averaged over the same bins as the finescale method. Arguments can be made for either choice for R_ω (3 or 7). Though for both methods, R_ω=3 gives the best fit overall, as indicated by the median of the estimates (middle orange line in the figure). Which for both methods falls within a factor of 2. Lastly, it can be seen that the estimates ε_FS,Θ tend to overestimate the estimates ε_vmp, whereas the estimates of ${\mathit{\epsilon }}_{\mathrm{FS},N^{\mathrm{2}}}$ contain both over- and underestimates.

Figure 4 The ratios between ε_FS and ε_vmp, for different values of R_ω. With (a) strain based on Θ_z and (b) strain based on N². ε_vmp has been averaged over the same segments as the segments of the finescale method. Color coding of the dots indicates average pressure of finescale estimate. The black lines indicate a factor 2 (dashed) and 3 (dotted) respectively, and the orange lines mark, from bottom to top, the 5, 25, 50, 75 and 95th percentiles. Jitter of the dots is to increase readability.

4.2 Dissipation rates from the finescale parameterization

In Fig. 5, the ratio $\frac{{\mathit{\epsilon }}_{\mathrm{FS}}}{{\mathit{\epsilon }}_{\mathrm{vmp}}}$ is plotted per station. Here the upper 100 dbar is omitted, as well as the two bins closest to the bottom, due to boundary effects (note, the surface mixed layer was shallower than 100 dbar for all profiles). Note that station NP5 is located within a fracture zone canyon with an observed bottom boundary layer from 3800 dbar and deeper, such that processes other than internal wave breaking may dominate the dissipation (Kunze et al., 2002; Polzin et al., 2014), and results should be treated with caution.

In general the two different methods (temperature-strain and density-strain) tend to reproduce similar patterns of overestimating or underestimating ε_vmp, though there are some areas with distinctly different patterns. The first area is at the depth of the Mediterranean Outflow Waters (st. NP1–NP3, 800–1500 dbar), where both methods overestimate the dissipation rate, and in particular the temperature-based estimates. Another region is at station M2 between 1000–2000 dbar, where large overestimates are observed for the temperature-based method. The finescale method with strain calculations based on vertical temperature gradients is known to be sensitive for finestructure of water-mass variability (Kunze et al., 2006; Thompson et al., 2007). Finestructure watermass variability in the form of interleaving patterns are observed in the $T/S$ diagram (Fig. B1). The depths where these patterns are observed correspond to the locations where temperature-based finescale estimates (severely) overestimate the measured dissipation rates ε_vmp. This finestructure can be caused by double-diffusive processes (e.g. salt-fingering, Laurent and Schmitt, 1999) or isopycnal stirring of tracers (e.g. Ferrari and Polzin, 2005), instead of internal waves as assumed in the method. This decouples the temperature strain from the density strain (N²) making it unsuitable for estimating dissipation (Ferris et al., 2022).

Another important region of discrepancy are the deep abyssal regions, where especially the finescale method based on the buoyancy frequency underestimates the dissipation rate. For small Θ and S_A gradients, the strain rates based on N² become sensitive to salinity spiking (Kunze et al., 2006; Thompson et al., 2007; Whalen et al., 2015). Changes in salinity as measured by the conductivity sensor are then dominated by sensor noise. Temperature sensors have higher accuracy than conductivity sensors and thus do not experience these problems in the regions as seen in (Fig. 5a, see also Sect. 5). Based on these results, we suggest to discard segments for which the maximum salinity difference (over the segment) is smaller than 0.008 g kg⁻¹ and maximum temperature difference is smaller than 0.08 °C (marked by yellow dots in Fig. 5b). Note this is empirically estimated for this dataset, so slightly different values might be needed for other areas.

An alternative explanation for the underestimates can be found in the choice for R_ω. If the actual R_ω is higher than the chosen value of R_ω=3, it would cause an underestimate in the dissipation rate. One of the causes for the higher R_ω can be the interaction of flow with topography and lee wave generation (Ferris et al., 2022). Though most of these profiles are taken over abyssal plains with relatively smooth topography, we consider this factor of limited influence to the underestimates.

Overall, the finescale estimates tend to approximate the dissipation rate reasonably well when compared to ε_vmp, also in this weak turbulent regime. By taking R_ω=3, we find that 75 % of ${\mathit{\epsilon }}_{\mathrm{FS},N^{\mathrm{2}}}$ is within a factor 5 of ε_vmp. Similarly, 79 % of ε_FS,Θ is within a factor 5 of ε_vmp. These results illustrated the caveats of the two different methods of calculating the strain rate. Though both methods perform comparably well in terms of accurately estimating the dissipation rate, the finescale method with strain based on the buoyancy frequency is preferred for the upper parts of the watercolumn, but for abyssal segments and areas with weak (salinity) gradients, the finescale method with strain based on temperature gradients appears to perform better.

Figure 5 Ratios of the finescale method and the microstructure data (averaged over the same bins). The surface mixed layer and the bottom two estimates of the finescale method are omitted due to boundary effects. With for (a) finescale strain based on Θ_z and (b) finescale strain based on N². Yellow dots mark the estimates flagged by the threshold stated in the text.

5 Thorpe Resorting

The size of vertical overturns of density, where higher density is found above lower density, can be associated with the dissipation of turbulent kinetic energy ε (Thorpe, 1977; Dillon, 1982). The size of the overturns is characterised by the Thorpe length scale L_T, usually defined as the root-mean-square (rms) of the Thorpe displacements: the shortest distance a fluid parcel needs to move in the vertical direction, in order to be at the location where it is stably stratified according to the stable background conditions. The Thorpe scale is assumed to be proportional to the Ozmidov scale $L_{\mathrm{O}}=\sqrt{\mathit{\epsilon }{N}^{-\mathrm{3}}}$, which represents largest scales of isotropic motions that are possible before stratification becomes a limiting factor (Dougherty, 1961; Ozmidov, 1965), leaving:

$$\begin{array}{}\text{(2)}& {\mathit{\epsilon }}_{\mathrm{TP}}=a^{\mathrm{2}}\stackrel{\mathrm{‾}}{L_{\mathrm{T}}^{\mathrm{2}}{N}^{\mathrm{3}}},\end{array}$$

where $a=\frac{L_{\mathrm{O}}}{L_{\mathrm{T}}}$. The method is most commonly applied to temperature overturns, because of the higher accuracy of temperature measurements compared to salinity and thus density measurements. This however, requires additional criteria to establish whether the temperature overturns actually represent density overturns. Regardless, the Thorpe resorting method has been used in a range of studies, from moored temperature sensors (van Haren and Gostiaux, 2012; van Haren, 2020; van Haren et al., 2020), (ship-based) CTD profiles (Ferron et al., 1998; Gargett and Garner, 2008; van Haren and Gostiaux, 2014) and numerical modelling (Smyth et al., 2001; Taylor et al., 2019), amongst others.

Here we calculate ε_TP from the high resolution CTD (64 Hz) mounted on the microstructure profiler, and compare this to ε_vmp (see Sect. 3). First we will aim to provide a “best” estimate of ε_TP, using all available data, including the often unavailable ε_vmp, which is used as benchmark. In doing so we also investigate what extra criteria are needed in order to allow for the use of temperature overturns as a proxy for density overturns. Then we provide estimates of ε_TP that exclude the use of measured ε_vmp, as this information is often not available. Instead we here examine the use of recently available climatological estimates of ε (de Lavergne et al., 2019, 2020) to find that this may improve estimates of ε_TP. In addition we explore possibilities for different values of the proportionality constant “a”.

5.1 Defining overturns and the Thorpe length scale

The raw temperature data from the CTD is first bin-averaged to the highest resolution (smallest vertical bin-width) for which individual measurements are distinguishable from another. This depends on the sampling frequency, the profiler fall speed and the sensor response times and is 0.08 m for this paper (Appendix C1). For each sample (measurement) in this profile, its vertical displacement is calculated as the difference between its current location and its location in the stably sorted profile. Overturns follow from the vertical integral of these displacements. An overturn is defined as the region between two zero crossings of the integral (Mater et al., 2015). Subsequently, the Thorpe length scale for each overturn is taken as the rms of all vertical displacements within each overturn (Appendix C2).

5.2 Choosing a value for “a”

ε_TP is proportional to $a=\frac{L_{\mathrm{O}}}{L_{\mathrm{T}}}$, which has previously been determined to be a=0.8 based on arithmetic averages over many overturns and a wide distribution using temperature sensors (Dillon, 1982). Other studies, using geometric means, mostly found values in the range between 0.8 and unity, but values between 0.4 and 1.3 have also been reported (Itsweire, 1984; Crawford, 1986; Ferron et al., 1998; Smyth et al., 2001; Mater et al., 2015). Here we will use two estimates of a for comparison that are, (1) the a=0.8 found by Dillon (1982), and (2) the geometric mean a over all overturns, calculated for each profile separately. The latter estimate is only possible because of the available estimates of ε_vmp that allows us to calculate a for each overturn. Over the 12 profiles, a varies between 0.80 and 1.33 with a mean of $a=\mathrm{1.0}\pm \mathrm{0.2}$. Note that using an arithmetic mean of a instead, is sensitive to high outliers and consequently overestimates ε_TP with respect to ε_vmp. Using these two choices for a, combined with different sets of criteria that determine if individual overturns are eligible to be used to calculate ε_TP, we will define 4 separate estimates of ε_TP per profile, that are mutually compared.

5.3 Criteria for selecting overturns

Not all detected overturns are suitable for estimating ε_TP for reasons discussed below. These reasons will be used to define four criteria for calculating ε_TP. The effect of these criteria on the detected overturns of an example station are shown in Fig. C1. Roughly 30 % of the detected overturns is flagged by one or more of the criteria.

5.3.1 Criterion 1: Thorpe-Scales Vs Overturn Size

When the size of the Thorpe scale L_T is larger than the vertical size of the overturn Δz_overturn,l, it implies that the fluid parcel needs to travel further than the size of the overturn itself in order to reach stably stratified waters. We consider this impossible and indicative of noise, spikes or other undefined problems. Therefore we only use overturns for which;

$$\begin{array}{}\text{(3)}& L_{\mathrm{T}}<\mathrm{\Delta }{z}_{\mathrm{overturn}}.\end{array}$$

5.3.2 Criterion 2: Measurement Accuracy vs Overturn Length

By dividing the precision of the temperature sensor by the vertical gradient of temperature, a minimal vertical length scale $L_{\mathrm{\Theta }}^{\mathrm{error}}$ is determined for which two observations are significantly different from one another. Overturns need to be larger than this minimum length scale. Details regarding the calculation of this length scale can be found in Appendix C3. We will only consider overturns for which,

$$\begin{array}{}\text{(4)}& L_{\mathrm{\Theta }}^{\mathrm{error}}<\mathrm{\Delta }{z}_{\mathrm{overturn}}.\end{array}$$

5.3.3 Criterion 3: A watermass criterion using the stability ratio

Observed temperature overturns have to be representable for density overturns. Situations can occur where a registered temperature inversion is exactly compensated by a salinity inversion, such that there is no density overturn, in which case the temperature inversion does not represent turbulence. In previous studies methods were used in which they would examine temperature-salinity relationships (Galbraith and Kelley, 1996; Mater et al., 2015) or the goodness of fit of a linear relationship between temperature and density (Gargett and Garner, 2008; van Haren and Gostiaux, 2014), among other methods. Results sometimes required tedious manual examination of the results or hard-to define thresholds for excluding data. Here we use two related, but different criteria to address this problem with the intention to simplify watermass based selection criteria. First, we take the Pearson correlation coefficient ${\mathrm{CC}}_{\mathrm{\Theta },{\mathit{\rho }}_{\mathrm{loc}}}\left(l\right)$ between Θ and locally referenced potential density ρ_loc. Assuming that turbulent overturns work against their local background stratification, the correlation coefficient is calculated over a vertical range of twice the overturn size, with a minimum of 10 m. A coefficient larger than 0.8 or smaller than −0.8, is considered a very good correlation (Asuero et al., 2006). A caveat to this criteria is that salinity might be more important in changing density than temperature, while they both work in the same direction, such that a false correlation between temperature and density could be obtained. Hence, in addition we will use that the stability ratio,

$$\begin{array}{}\text{(5)}& R_{\mathit{\rho }}={\displaystyle \frac{\mathit{\alpha }\frac{\partial \mathrm{\Theta }}{\partial z}}{\mathit{\beta }\frac{\partial S_{\mathrm{A}}}{\partial z}}},\end{array}$$

indicates the influence of vertical temperature and salinity gradients on the vertical density gradient. When |R_ρ|>1, temperature dominates the changes in density. So for the temperature overturns to be a proxy for density overturns, temperature changes should at the very least be dominant over salinity changes. We combine both criteria to only include overturns for which

$$\begin{array}{}\text{(6)}& {\mathrm{CC}}_{\mathrm{\Theta },{\mathit{\rho }}_{\mathrm{loc}}}\left(l\right)<-\mathrm{0.8}\mathrm{and}\left|R_{\mathit{\rho }}\right|>\mathrm{1}\end{array}$$

Further details for calculating and applying this criteria per overturn are provided in Appendix C. Note that a caveat of using temperature as a proxy for density overturns is that, when temperature is not dominant in setting the density changes, some overturns might be unaccounted for.

5.3.4 Criterion 4: Thorpe scale Vs Ozmidov Scale

The fundamental assumption for using the Thorpe scale, is that there is a statistical agreement (over a large number of observations) that the Thorpe scale L_T is approximately equal to the Ozmidov scale L_O (Dougherty, 1961; Ozmidov, 1965). When both length scales are not comparable in size, we argue that Eq. (2) does not hold and therefore we only use overturns for which:

$$\begin{array}{}\text{(7)}& {\displaystyle \frac{\mathrm{1}}{\mathrm{5}}}{L}_{\mathrm{O}}\left(l\right)<L_{\mathrm{T}}\left(l\right)<\mathrm{5}\times L_{\mathrm{O}}\left(l\right)\end{array}$$

Although the factor 5 is somewhat arbitrary, it turns out that this criterion is very important for improving dissipation estimates using Thorpe scales. The stricter the criterion, the fewer dissipation estimates there are, but the looser the criterion, tends to overestimate the dissipation rates (see results below). A factor 5 appears to be a reasonable balance between these two effects.

5.4 Different ways for estimating dissipation using overturns

In this study, we compare the following four methods to estimate the dissipation rate ε_TP from overturns:

• Method 1 – ${\mathit{\epsilon }}_{\mathrm{TP}}^{\mathrm{geo}}$: Here a is the profile-specific geometrical mean, and all criteria are used. Hence, we use all available information.

• Method 2 – ${\mathit{\epsilon }}_{\mathrm{TP}}^{\mathrm{Dil}}$: Here a=0.8, and all criteria are used.

• Method 3 – ${\mathit{\epsilon }}_{\mathrm{TP}}^{\mathrm{Dilx}\mathrm{4}}$: Here a=0.8, and criterion 4 is not used. This is because ε_vmp is often not available, when CTD measurements are available.

• Method 4 – ${\mathit{\epsilon }}_{\mathrm{TP}}^{\mathrm{Dil},\mathrm{clim}}$: Here a=0.8, and all criteria are used. However, criterion 4 is based on climatological parameterized estimates of ε.

Comparing these four estimates allows us to evaluate how well the Thorpe resorting method reproduces microstructure-based dissipation rates. This comparison also helps assess whether the Thorpe method can estimate low dissipation rates in the deep sea (i.e., below the shear probe noise floor), and how the results are affected when criterion 4 is omitted or replaced with climatological values. The latter is of interest as these are globally available, whereas ε_vmp is usually not. Note that we do not compare a version of method 4 where the geometric mean of a is calculated using the climatological value of ε. This is because the estimates for ε_clim tend to be larger (Fig. 6) than the instantaneous estimate from the VMP, meaning that L_O and thus a and ε_TP will be overestimated and thus has no additional value. This also has implications for the application of criterion 4 that we did not further refine.

5.5 Results

To compare values obtained in different ways over different vertical depth ranges, we will average ε_vmp and all ε_TP(l) estimates (per overturn), into vertical bins of 25 m. This implicitly assumes that the mean ε_TP(l) over the available overturns within this depth range, is representative of the total dissipation in that depth range. We then assess the Thorpe method by taking the geometric mean and standard deviation of the ratio ${\mathit{\epsilon }}_{\mathrm{TP}}/{\mathit{\epsilon }}_{\mathrm{vmp}}$, to indicate how well the estimates match. However, due to the nature of the geometric mean, it cannot directly be compared with the factor 2 uncertainty that is known for microstructure measurements. Also, outliers on either end can partially cancel out and are downweighted, so the results of the geometric mean may underestimate the severity of outliers. In addition to the standard deviation, we therefore add another metric that assesses the percentage of ε_TP that is within a factor 2, 5 or 10 from ε_vmp. This is more indicative of outliers. The results of these metrics are listed in Table 2, for all four methods. An additional cross-validation of the estimates against ε_vmp is shown in Fig. C2.

Table 2 Overview of the performance of methods 1–4 for all data. The “average ratio” is based on the geometric mean of the ratio ${\mathit{\epsilon }}_{\mathrm{TP}}/{\mathit{\epsilon }}_{\mathrm{vmp}}$, ± the geometric standard deviation. The last three columns indicate the percentage of estimates that falls within the respective factor compared to ε_vmp.

5.5.1 Comparing Method 1 (Geometric mean a) with Method 2 (a=0.8)

Method 1 is closer to ε_vmp (factor 1.1) than Method 2 (factor 1.6). However, Method 1 has less estimates within a factor 2 (48 % vs 54 %, respectively) or a factor 5 (87 % vs 92 %, respectively) than Method 2. This indicates that the Method 2 estimates have a narrower spread, but based on the ratio, there is a systematic overestimation as well. Overestimation of the microstructre estimates can be caused by a too coarse sampling resolution of the CTD instrumentation (Sheehan et al., 2023).

The a factor is a quantity that varies for every overturn. Taking the geometric average over all available data or taking a=0.8 results in dissipation estimates that perform comparable well, as these two a factors are not that different. We have not been able to further identify factors that have a strong influence in determining the magnitude of a, as to make a good estimate what a to use.

5.5.2 Comparing results when omitting Criterion 4

Both methods 1 & 2 use a direct calculation of L_O in criterion 4 which requires a-priori knowledge of ε. Fully omitting criterion 4, and taking a=0.8, is done by method 3. The impact of these choices is clearly visible in Table 2. Where the average ratio and standard deviation is much larger for method 3 compared to the other methods. Method 3 also has 6 times more outliers beyond a factor 10, than Method 1 or 2, and we assessed that these can sometimes be hundreds of times larger, i.e. large spikes. In short, within the setup of this study, we conclude that criterion 4 that assures that the Ozmidov scale and Thorpe scale are of similar order of magnitude, is of particular importance for obtaining Thorpe-based dissipation estimates that are comparable to direct observations.

Method 4 uses climatological dissipation rates ε_clim (de Lavergne et al., 2019, 2020) for calculating L_O and Criterion 4. Generally ε_clim slightly overestimates the measured dissipation rates ε_vmp (Figs. 6 and C2), but overall follows the broad structure of ε_vmp over multiple orders of magnitude. The climatological estimates do in some occasions show larger differences, e.g. differences in depth (e.g. st. NP5 and MM) or have a larger overestimate (e.g. st. NP2 and NP3).

Method 4 significantly improves the dissipation rates compared to Method 3, where criterion 4 is omitted. The estimates of Method 4 are less accurate and have a larger spread than the estimates of Methods 1 and 2 (Table 2), however, the given results seem acceptable. Especially considering that within these errorbars, the results tend to follow the large scale patterns of ε_vmp over multiple orders of magnitude (Figs. 6 and C2). Using a climatological estimate of the dissipation rate thus provides a way to improve dissipation estimates using the Thorpe resorting method when other dissipation rates are not available.

5.5.3 Summary from Thorpe estimates

In this study, we compared dissipation rates estimated from temperature overturns ε_TP to ε_vmp. With 8 cm vertical resolution of temperature we found 𝒪(10³) overturns for each profile of which about 30 %–40 % is flagged and removed per profile (mostly by criterion 3 and 4 if applied), based on criteria determining the eligibility of using the overturn for calculating the dissipation rate. For the remaining overturns, dissipation rates ε_TP are calculated, averaged in 25 m bins and compared against ε_vmp (Figs. 6; Note that while the comparison between methods is done on the 25 m-binned data, the figure shows the estimates in half overlapping 200 m bins for increased readability. The 25 m binned estimates for three example stations are shown in Fig. C2). Using all available information, roughly 50 % of the estimates ε_TP fall within a factor 2 of ε_vmp, and follow the structure of ε_vmp quite accurately over multiple orders of magnitude. This indicates the potential for using Thorpe overturns for estimating dissipation rates. However, these results relies heavily on the use of criterion 4, which omits overturns for which the Ozmidov scale is more than a factor 5 different from the Thorpe scale. Without applying this criterion, the results are prone to large outliers and it is difficult, to distinguish between outliers and accurate estimates. However, criterion 4 requires a-priori knowledge of ε (here obtained using ε_vmp) that usually is not available. When instead using globally available climatological estimates of dissipation, reasonable results are still obtained 81 % of the estimates being within a factor 5 of ε_vmp. We also found that for this particular setup, the performance of both methods 1 & 2 was not much different. However, for a different vertical bin-size (now 8 cm) and with other criteria, this may be different.

Figure 6 Dissipation estimates ε_TP for all different methods (ε_geo: Method 1, ε_Dil: Method 2, ε_Dilx4: Method 3 and ε_Dil,clim: Method 4, as defined in Sect. 5.4) are compared to ε_vmp and the used climatological estimates ε_clim. For plotting purposes all data is averaged over the same half overlapping 200 m bins.

6 The Triple Decomposition Framework

The data from the microstructure profiler provides direct observations of dissipation rate (ε) and temperature variance (χ_Θ). Using the analytical triple decomposition framework (Stern, 1967; Joyce, 1977; Davis, 1994; Garrett, 2001; Ferrari and Polzin, 2005), the microstructure data can be used to qualitatively assess the relative contributions of isoneutral and dianeutral mixing to tracer variance dissipation (Ferrari and Polzin, 2005; Naveira Garabato et al., 2016; Merrifield et al., 2016; Orúe-Echevarría et al., 2023; Castro et al., 2024). This can shed a light on the dominant processes that set the observed Θ and S_A variability. A complete description of the framework can be found in Appendix D.

Where a standard Reynolds decomposition separates variables into a mean and a fluctuating component (for any tracer C, $C=\stackrel{\mathrm{‾}}{C}+C^{\prime }$), the triple decomposition framework variables are separated into a mean or background component (C^m), an mesoscale eddy component (C^e) and a microscale turbulent component (C^t). One can then write $C=C^{\mathrm{m}}+C^{\mathrm{e}}+C^{\mathrm{t}}$ (Joyce, 1977; Davis, 1994; Garrett, 2001). This separation of scales assumes a spectral gap between the different scales, though it is not clear whether such gap really exists in the ocean (Davis, 1994; Van Haren et al., 1994). A spectral gap is required for the spatially divergent components of the Reynolds stresses. Here we used a simplified approach, ignoring these terms under the assumption of non-divergence amongst others (Appendix D). This assumption is a source of uncertainty for the method. Nevertheless, the triple decomposition is way to qualitively assess the scale transformations of T-S variance between scales (Ferrari and Polzin, 2005).

The triple decomposition framework implies a balance between the production of tracer variance at the different scales and the dissipation of tracer variance at the microscale. This creates two pathways how tracer variances cascades from the scales where it is produced down to the scales where it is dissipated (Garrett, 2001):

$$\begin{array}{}\text{(8)}& \textcolor{red}{}\end{array}$$

In (1), tracer variance is produced by mesoscale stirring. This creates finescale structures that are transferred to the microscale by isotropic turbulence (2) after which it is dissipated (4). In (3), microscale turbulence acts directly on background tracer gradients and variance gets directly transferred from the large scales to the dissipation scales (4). The tracer dissipation itself (end of path) is directly observed by the VMP (χ_Θ). The total (observed) dissipation results from the combination of downscale variance transfer by small-scale and mesoscale turbulence. At the turbulent scales Θ^t, the dianeutral production term $P_{{\mathrm{\Theta }}^{\mathrm{2}}}^{⟂}$ can be approximated by using a gradient-flux approximation (Osborn and Cox, 1972) and be written in terms of a dianeutral diffusivity D (see also Appendix D),

$$\begin{array}{}\text{(9)}& P_{{\mathrm{\Theta }}^{\mathrm{2}}}^{⟂}=\mathrm{2}〈{\mathit{u}}^{\prime }{\mathrm{\Theta }}^{\mathrm{t}}〉\cdot {\mathrm{\nabla }}_{⟂}{\mathrm{\Theta }}^m=-\mathrm{2}D|{\mathrm{\nabla }}_{⟂}{\mathrm{\Theta }}^m{|}^{\mathrm{2}}\end{array}$$

The angled brackets indicate an averaging scale that is large compared to the scales within, but smaller than one scale larger in the framework, here ∇_⟂ is the dianeutral gradient operator. If an assessment can be made of the dianeutral production at the microscale, the contribution by mesoscale stirring Θ^e follows from,

$$\begin{array}{}\text{(10)}& P_{{\mathrm{\Theta }}^{\mathrm{2}}}^{\parallel }={\mathit{\chi }}_{\mathrm{\Theta }}-P_{{\mathrm{\Theta }}^{\mathrm{2}}}^{⟂}={\mathit{\chi }}_{\mathrm{\Theta }}-\mathrm{2}D{\left({\mathrm{\nabla }}_{⟂}{\mathrm{\Theta }}^m\right)}^{\mathrm{2}}.\end{array}$$

To obtain $P_{{\mathrm{\Theta }}^{\mathrm{2}}}^{⟂}$ using Eq. (9), it is needed to obtain Θ^m from measured CTD data ($\stackrel{\mathrm{‾}}{\mathrm{\Theta }}$). Here an overline refers to data as measured by the CTD instrumentation averaged in 𝒪(1m) vertical resolution. Formally, to get Θ^m one should be averaging over appropriate long term time and length scales (Ferrari and Polzin, 2005). Only few datasets are large enough and with appropriate spacings that such averaging is possible, e.g. those of NATRE (Ferrari and Polzin, 2005) and DIMES (Naveira Garabato et al., 2016). For our single profiles we instead use the approach of Castro et al. (2024) that approximates background gradients $\left(\frac{\partial {\mathrm{\Theta }}^m}{\partial z}\right)$ using polynomal fits against density, to remove density-compensated temperature intrusions produced from eddy stirring (see Appendix D). Dianeutral diffusivities have been calculated following Eq. (12).

Similar as for Θ, the triple decomposition framework can also be applied with S_A as the tracer. Calculating the dianeutral production term for salinity $P_{S^{\mathrm{2}}}^{⟂}$ from the CTD data is done the same way as was done for temperature. However, the used profiler was not equipped with a microstructure salinity sensor, so haline variance χ_S has not been measured. To have an estimate of χ_S, we follow the approach of Castro et al. (2024) and assume,

$$\begin{array}{}\text{(11)}& {\mathit{\chi }}_S\approx \mathrm{2}D{\left({\displaystyle \frac{\partial {\stackrel{\mathrm{‾}}{S}}_{\mathrm{A}}}{\partial z}}\right)}^{\mathrm{2}}.\end{array}$$

The relative roles between isoneutral and dianeutral processes are considered based on the ratios between the dianeutral production and variance dissipation ($P_C^{⟂}/{\mathit{\chi }}_C$, Fig. 7).

The ratios show a high degree of variability for both temperature and salinity. For the majority of the data, diapycnal production dominates over eddy stirring, especially for temperature. Though there is a signal of more prominent eddy stirring for stations NP1–NP3 between 1000–1500 dbar, which coincides with the influence of Mediterranean Outflow Water. Also a stronger stirring signal is noted around 2000 dbar at station M2. In terms of salinity, a signal of isopycnal stirring dominating over diapycnal production around 1000 dbar extends over the entirety of the dataset. For the stations in the Mediterranean Outflow and from the Mixation dataset (st. M1–M3) the dominance of isopycnal stirring extends to larger depths.

The ratios between the diapycnal production and χ (Fig. 7) show a similar pattern as to where the finescale estimates overestimate the microstructure derived dissipation rates ε_vmp (Fig. 5), and also where the Thorpe estimates fail to produce estimates (Fig. 6). Those issues arise where isoneutral stirring processes are of similar or greater importance than dianeutral mixing processes (Fig. 7), and explain increased density-compensated T-S variability at the finescale. Alternatively, double diffusive processes can cause a similar effect as isoneutral stirring. Though almost the entire dataset is susceptible to double diffusion (R_ρ>1), no developed thermohaline staircases are observed in the individual Θ and S_A profiles, but the $T/S$ diagram (Fig. B1) shows increased interleaving patterns for the areas where a smaller ratios $P_C^{⟂}/{\mathit{\chi }}_C$ are found, thus supporting the idea that isoneutral stirring processes at least dominate over double diffusive processes. Note that the finescale parameterization assumes that the observed strain is predominantly caused by the non-linear wave-wave interactions (Dematteis et al., 2024). Increased temperature variability will thus cause an overestimate of the strain rate and consequently of the dissipation rate. This effect will be smaller with strain rates based on density variability (N²), as density compensation between T-S variability will reduce the effects of increased T-S variability due to isopycnal stirring. As such, these locations are also indicative of regions where temperature is not suitable as a proxy for density, as used in the Thorpe method, with perhaps most notably, the Med Outflow around stations NP1–NP3. Although these general patterns of the triple decomposition thus possibly indicate where the Thorpe or finescale method can be applied, a clear metric for rejections of those estimates is yet inconclusive and requires further research.

Figure 7 The ratios between $P_C^{⟂}/{\mathit{\chi }}_C$ expressed in  %. With (a) the tracer being temperature and (b) salinity. Ratios have been bin-averaged into 50 m bins for readability.

7 Diffusivities

Microstructure data is most commonly used to obtain estimates of the turbulent dissipation rate ε and of the dianeutral diffusivity D. For this we use the Osborn relation (Osborn, 1980),

$$\begin{array}{}\text{(12)}& D_{\mathit{\rho }}=\mathrm{\Gamma }{\displaystyle \frac{{\mathit{\epsilon }}_{\mathrm{vmp}}}{{\stackrel{\mathrm{‾}}{N}}^{\mathrm{2}}}}.\end{array}$$

Here, D is the dianeutral diffusivity (m² s⁻¹), ${\stackrel{\mathrm{‾}}{N}}^{\mathrm{2}}$ is the segment averaged buoyancy frequency obtained by the method of adiabatic levelling (Bray and Fofonoff, 1981). Γ is the mixing coefficient and often assumed to be 0.2 (Osborn, 1980; Gregg et al., 2018). Here we calculate Γ based on the flux Richardson number R_f (Fernández Castro et al., 2022), using the microstructure data as,

$$\begin{array}{}\text{(13)}& R_{\mathrm{f}}={\displaystyle \frac{D_{\mathrm{\Theta }}{N}^{\mathrm{2}}}{{\mathit{\epsilon }}_{\mathrm{vmp}}+D_{\mathrm{\Theta }}{N}^{\mathrm{2}}}},\text{where}{D}_{\mathrm{\Theta }}={\displaystyle \frac{{\mathit{\chi }}_{\mathrm{\Theta }}}{\mathrm{2}{\left(\frac{\partial \stackrel{\mathrm{‾}}{\mathrm{\Theta }}}{\partial z}\right)}^{\mathrm{2}}}}.\end{array}$$

Here D_Θ is the temperature diffusivity based on the Osborn-Cox model (Osborn and Cox, 1972), and the effective mixing coefficient is related to the flux Richardson number as (Gregg et al., 2018),

$$\begin{array}{}\text{(14)}& \mathrm{\Gamma }={\displaystyle \frac{R_{\mathrm{f}}}{\mathrm{1}-R_{\mathrm{f}}}}.\end{array}$$

For the calculation of Γ, R_f has been smoothed with a 30-point running mean. Most profiles show an increase in R_f between 1000–2000 dbar, leading to a corresponding rise in Γ, and resulting in differences in diffusivities of up to an order of magnitude (Fig. 8). Areas with increased R_f can be indicative of caveats in the underlying models and assumptions, and can partially be caused by non-turbulence effects such as isopycnal stirring or double diffusion. The resulting diffusivities range from $\mathcal{O}({\mathrm{10}}^{-\mathrm{6}}-{\mathrm{10}}^{-\mathrm{2}})$ m² s⁻¹, with increasing diffusivities closer to the bottom (Fig. 8). Enhanced diffusivities extending to mid-depths are observed for the stations at or close to the Mid-Atlantic Ridge (st. NP4–NP6).

Figure 8 The calculated flux Richardson number R_f (in red) and mixing coefficient Γ (in grey) for each station. Also shown are the dianeutral diffusivities calculated with an assumed Γ=0.2, shown as D_ρ. And the diffusivity calculated using Eq. (12), with the variable Γ, shown as D_Γ. The canonical Γ=0.2 is marked by the dashed black line. A 50-point running average has been applied on all diffusivities for readability.

Additionally, the triple decomposition framework can also be used to get estimates of the isoneutral diffusivity based on the microstructure data. However, due to the small isoneutral gradients of temperature and salinity in the abyss, the final estimates appear to be noisy and the results inconclusive. The estimates of the isoneutral diffusivity can be found in Appendix D2 for the interested reader.

8 Discussion and conclusions

In this study we analysed a set of full-depth microstructure profiles from the North Atlantic (Fig. 1). First, dissipation estimates from velocity shear probes and the fast-response thermistors were compared (Fig. 2). The thermistor-based estimates generally showed good agreement with shear-based estimates for moderate levels of turbulence, but they were sometimes up to 2 orders magnitude lower than the shear-based dissipation rates for low turbulence levels ($\mathit{\epsilon }\lesssim {\mathrm{10}}^{-\mathrm{10}}$ W kg⁻¹). The data suggested that it was unlikely that the noise floor of the thermistor-based estimates was reached, while observations in the weakly turbulent deep ocean were frequently below the noise floor of the shear-based estimates. This could lead to substantial differences when shear-based dissipation rates are used for e.g. flux calculations, as is common practice. Therefore we combine both products to obtain ε profiles that adequately capture the range of turbulence levels in our dataset (Fig. 9). Using a coincidental unique profile where dissipation rates from both shear and thermistor based where entirely separated over a large part of the water column (Fig. 2, station NP3), we could empirically determine the noise floor for the shear-based method to be at $\mathrm{8.3}\times {\mathrm{10}}^{-\mathrm{11}}$ W kg⁻¹. Hence, considering that the upper limit of the thermistor based method is $\sim \mathrm{1}\times {\mathrm{10}}^{-\mathrm{8}}$ W kg⁻¹, we combined both estimates by using thermistor-based dissipation estimates below $\mathrm{8.3}\times {\mathrm{10}}^{-\mathrm{11}}$ W kg⁻¹, and shear-based estimates above that. Note that, shear-based estimates are preferred over thermistor-based estimates, when far from the noise floor, as these are arguably more direct and require fewer assumptions. To conclude, when one expects to measure weak turbulence, we recommend using both shear and thermistor probes in measurements to be able to cover the entire turbulent range.

Figure 9 The dissipation rates ε_vmp, ε_TP and ε_FS plotted for each station. The finescale estimates are combined based on the threshold recommended in Sect. 4. The Thorpe estimates plotted correspond to ε_Dil,clim of Fig. 6 and method 4 of Sect. 5. All methods are averaged over the same segments as the finescale method.

Second, the microstructure ε estimates were compared to both the finescale parameterization and the Thorpe resorting method, obtained by using the CTD data. The finescale parameterization has been applied in its strain-only form based on strain rates calculated from both the buoyancy frequency and vertical temperature gradients. In general, both methods tend to follow the microstructure estimates reasonably well, over multiple orders of magnitude. However, the temperature-based estimates severely overestimated the measured dissipation rates at the depth of the Mediterranean Outflow Water and stations M2, MM. The triple decomposition framework supported the hypothesis that an increased isopycnal stirring signal was cause for finestructure watermass variability not associated with internal waves, leading to the mismatch. In these locations the method based on the buoyancy frequency performed better. However, the density-based method was affected by salinity spiking for abyssal parts of the profiles due to the small salinity gradients in the segments. Overall we therefore recommend using the temperature-based estimates over the density-based estimates for segments where the maximum salinity difference over a segment is smaller than 0.008 g kg⁻¹ and the maximum temperature difference is smaller than 0.08 °C. In general, the agreement in our estimates for ${\mathit{\epsilon }}_{\mathrm{FS},N^{\mathrm{2}}}$ compares well with other comparisons between the finescale method and microstructure data (e.g. Whalen et al., 2015; Baumann et al., 2023). Only few studies have reported on the use of strain based on temperature data (Thompson et al., 2007). In this study, the limitations of this method in upper ocean applications, due to the sensitivity are clearly seen. However, it is also shown that the strength of this method lies in the deeper parts of the watercolumn, where this temperature-based method outperforms the estimates based on buoyancy frequency strain.

For the Thorpe resorting method, the analysis focussed on the selection criteria to distinguish between false and actual overturns. One criterion aimed to provide an objective way for determining if a detected temperature overturn is a good proxy for a density overturns. This criterion is based on the stability ratio R_ρ and a correlation between density and temperature variations (Sect. 5). As the Thorpe method is based on the assumption that the Thorpe lengthscale is proportional to the Ozmidov lengthscale, we tested a second criterion that requires these length scales to be within a factor 5 of each other. The criterion turned out to be critical for removing outliers and providing accurate estimates of the dissipation rate. The caveat here is that this criterion needs a-priori estimates of the dissipation rates for the calculation of the Ozmidov scale. To make this criterion more generally applicable we tested the use of stationary climatologically averaged estimated of dissipation (de Lavergne et al., 2020). It turned out that applying Thorpe based dissipation estimates to CTD data under these criteria, provides very acceptable results (Fig. 9). Various other studies have previously reported on the comparison between the Thorpe method and microstructure data. Overestimating the dissipation rates by the Thorpe method is often seen and is attributed to either the sampling resolution (e.g. Howatt et al., 2021; Sheehan et al., 2023) and the relationship between L_T and L_O (e.g. Ferron et al., 1998; Mater et al., 2015). While sampling resolution depends on the instrumentation used, and can be difficult to improve. With microstructure data, an optimal factor a can be found for the relationship between L_T and L_O, but when only coarse data (e.g., CTD data) is available, the use of our criterion 4 and the climatological dissipation estimates is a useful way to improve the Thorpe based estimates, also for areas with weak dissipation rates.

It is noted that the data in this study is relatively sparse. Stations are far apart and only snapshots in time. It is known that for a good comparison between direct microstructure observations and indirect parameterizations (both finescale and Thorpe parameterizations) over larger spatiotemporal scales, the microstructure data should be averaged over equal spatiotemporal scales (Whalen, 2021). The presented results here have first been (vertically) averaged to equal scales, but optimally also horizontally or temporal averaging would be applied if a closer spacing of the stations was available. Adding more profiles will undoubtedly lead to better comparisons between the direct microstructure observations and the indirect parameterizations, but it seems unlikely that it will change the key conclusions, as the focus of this paper is on individual differences between methods rather than the assimilation of a large quantity of microstructure profiles (see also Sect. 1).

In conclusion, both parameterizations approximate the dissipation rates reasonably well, with 79 % (for ε_FS,Θ), 75 % (for ${\mathit{\epsilon }}_{\mathrm{FS},N^{\mathrm{2}}}$) and 81 % (for ${\mathit{\epsilon }}_{\mathrm{TP}}^{\mathrm{clim}}$) being within a factor 5 of the measured ε_vmp. Further analyses using the triple decomposition framework (Fig. 7) shows the influence of mesoscale stirring, producing finescale density-compensated TS variance, where also the Thorpe and finescale parameterizations either failed or showed strong overestimates and help understand why the underlying assumptions broke down.

Finally, by using the Osborn relation the dianeutral diffusivities were calculated using ε_vmp. An additional attempt to obtain estimates of the isoneutral diffusivity from combining the dianeutral diffusivities with the triple decomposition framework remained inconclusive.

This work has compared and contrasted several direct observations of dissipation rates with parameterizations (Fig. 9). This helps understand where the use of different parameterizations may or may not be appropriate, and shows the degree to which turbulence can be estimated. Obtaining estimates within a factor 2, even between shear probes, between thermistor probes and between both shear and thermistor probes, is not trivial. In that light, parameterizations that can estimate dissipation rates within a factor 2–5 from measured dissipation rates, could arguably be considered successful. Further improving both observation and parameterizations will allow us to improve global ocean mixing estimates in space and time.

Appendix A: Processing of microstructure data

In this appendix, the procedures taken for processing both the microstructure shear and microstructure thermistor data are described.

A1 Microstructure shear data

The microstructure profiler VMP-6000 of Rockland Scientific International (RSI) was equipped with two orthogonally mounted airfoil shear probes, sampling at a rate of 512 Hz. This provided independent estimates of both $\frac{\partial u}{\partial z}$ and $\frac{\partial v}{\partial z}$. Under the assumption of isotropic turbulence, each component can be independently used for the estimation of ε_µU according to (Oakey, 1982; Piccolroaz et al., 2021; Lueck et al., 2024),

$$\begin{array}{}\text{(A1)}& {\mathit{\epsilon }}_{\mathrm{\mu }\mathrm{U}}={\displaystyle \frac{\mathrm{15}}{\mathrm{2}}}\mathit{\nu }〈{\left({\displaystyle \frac{\partial u^{\prime }}{\partial z}}\right)}^{\mathrm{2}}〉.\end{array}$$

With ν being the kinematic viscosity and $〈{\left(\frac{\partial u^{\prime }}{\partial z}\right)}^{\mathrm{2}}〉$ the shear variance, $\frac{\partial u^{\prime }}{\partial z}$ the horizontal shear fluctuations. (^′) denotes a turbulent quantity. The processing of the raw microstructure data to ε_µU was done using software packages provided by RSI (Zissou Premium, v1.0). This processing involved despiking of the data, highpass filtering and the removal of coherent noise using a Goodman filter (Goodman et al., 2006). Estimates of ε_µU were obtained by fitting to a theoretical Nasmyth spectrum (Nasmyth, 1970; Oakey, 1982). Example spectra are shown in Fig. A1a.

Figure A1 (a) Two examples of the observed shear spectra. With the points used in the integration marked by the circles. Background shows the theoretical Nasmyth spectrum for various levels of ε. Vertical dotted line indicates the upper integration limit. (b) Corresponding example microstructure temperature spectra. Points used for the spectral integration are marked by the circles. Background shows the theoretical Kraichnan spectrum for various levels of χ with the Batchelor wavenumber k_B fixed at 33 cpm. Note, the Kraichnan spectra for k_B fixed at 113 cpm are omitted for clarity. Vertical dotted line shows the Batchelor wavenumber, used for the calculation of ε_µT.

Furthermore, the estimates were quality controlled according to the ATOMIX best practices (Lueck et al., 2024). These quality controls check for (i) a poor Figure of Merit, (ii) a large fraction spikes (iii) the ratio between estimates between the two probes, (iv) the number of iterations of the de-spiking routine, and (v) the amount of variance resolved. Most data was removed by the Figure of Merit criterion, which is a measure of how close the measured spectrum resembles a model spectrum (Lueck et al., 2024). In the limit of very weak turbulence, one or more assumptions underlying the spectra might not be valid, resulting in a poor resemblance and thus rejection of the estimate based on the Figure of Merit. Most of the rejected estimates originate in regions where the estimates are close to the instrument noise floor. Estimates that did not pass the quality control criteria were removed from further calculations. The estimates from the two different probes were averaged in case both passed the quality control, otherwise the single estimate was taken.

The 95 % confidence interval for dissipation estimates obtained by the spectral integration method is given by (Lueck, 2022; Lueck et al., 2024),

$$\begin{array}{}\text{(A2)}& {\mathrm{CI}}_{\mathrm{95}\mathrm{\%},\mathrm{U}}={\mathit{\epsilon }}_{\mathrm{\mu }\mathrm{U}}\exp (\pm \mathrm{1.96}{\mathit{\sigma }}_{\ln \mathit{\epsilon }}).\end{array}$$

σ_ln ε is obtained from (Lueck et al., 2024),

$$\begin{array}{}\text{(A3)}& {\mathit{\sigma }}_{\ln \mathit{\epsilon }}^{\mathrm{2}}={\displaystyle \frac{\mathrm{5.5}}{\mathrm{1}+{\left({\widehat{L}}_f/\mathrm{4}\right)}^{\mathrm{7}/\mathrm{9}}}},{\widehat{L}}_f={\displaystyle \frac{l_{\mathit{\epsilon }}}{L_{\mathrm{K}}}}{V}_{\mathrm{f}}^{\mathrm{3}/\mathrm{4}}\end{array}$$

where ${\mathit{\sigma }}_{\ln \mathit{\epsilon }}^{\mathrm{2}}$ is the variance of the dissipation estimate, L_f a non-dimensional datalength, l_ε the length of the data segment, L_K the Kolmogorov length and V_f the fraction of shear variance that is resolved by the spectral integration. A confidence interval for dissipation estimates obtained from fitting the inertial subrange follows a slightly different procedure (Lueck et al., 2024), though that does not apply to this dataset, as this generally applies to dissipation rates larger than ${\mathit{\epsilon }}_{\mathrm{\mu }\mathrm{U}}\ge \mathcal{O}({\mathrm{10}}^{-\mathrm{5}}$ W kg⁻¹) (Lueck, 2022).

A2 Microstructure thermistor data

Simultaneously to the shear probes, the two FP07 fast-response thermistors, sampling at 512 Hz, provide a way to calculate the dissipation of temperature variance χ_Θ. Here we will mostly follow the procedure as described by Piccolroaz et al. (2021). Assuming isotropy, χ_Θ can be determined from,

$$\begin{array}{}\text{(A4)}& {\mathit{\chi }}_{\mathrm{\Theta }}=\mathrm{6}{\mathit{\kappa }}_{\mathrm{T}}〈{\left({\displaystyle \frac{\partial T^{\prime }}{\partial z}}\right)}^{\mathrm{2}}〉=\mathrm{6}{\mathit{\kappa }}_{\mathrm{T}}\underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}{\mathrm{\Psi }}_T\left(k\right)\mathrm{d}k\end{array}$$

Here κ_T is the molecular thermal diffusivity and $〈{\left(\frac{\partial T^{\prime }}{\partial z}\right)}^{\mathrm{2}}〉$ the temperature gradient variance, with $\frac{\partial T^{\prime }}{\partial z}$ the turbulent temperature gradients. The temperature gradient variance equals the area under the one-dimensional temperature gradient wavenumber spectrum Ψ_T(k) (Piccolroaz et al., 2021). The measured spectrum is fitted to a theoretical Kraichnan spectrum (Kraichnan, 1968) using the Maximum Likelihood Estimation method of Ruddick et al. (2000). An example spectrum is shown in Fig. A1b. This method provides an estimate of χ_Θ, and a way to get ε_µT with an uncertainty estimate on the fitted Batchelor wave number k_B. Here ε_µT depends on the roll-off of the temperature gradient spectrum towards the Batchelor wave number (Ruddick et al., 2000; Luketina and Imberger, 2001). Quality control criteria were applied according to the methods of Piccolroaz et al. (2021); Ruddick et al. (2000) and were based on:

• the signal-to-noise ratio of the integrated spectrum and the instrument noise spectrum should be larger than 1.3.

• the mean absolute deviation (MAD) over the integrated wavenumber range should be $\mathrm{MAD}\le \mathrm{2}(\mathrm{2}/d{)}^{\mathrm{1}/\mathrm{2}}$, with d the number of degrees of freedom of the measured spectrum. The MAD criterion is comparable to the Figure of Merit criterion of the shear data.

• a log likelyhood ratio between a power-law fit and the theoretical spectrum should be LR >2.

Additionally, if estimates from two probes are within a factor of 2.8, they are averaged, otherwise the smallest of the two estimates was taken.

The MLE fitting procedure provides a standard deviation for the fitted Batchelor wave number k_B. The dissipation rate ε_µT is proportional to the Batchelor lengthscale, which is the inverse of the Batchelor wave number:

$$\begin{array}{}\text{(A5)}& {\mathit{\epsilon }}_{\mathrm{\mu }\mathrm{T}}=\mathit{\nu }{\mathit{\kappa }}_{\mathrm{T}}^{\mathrm{2}}{k}_{\mathrm{B}}^{\mathrm{4}}.\end{array}$$

With ν the kinematic viscosity and κ_T being the thermal molecular diffusivity. Through error propagation, the standard deviation of ε_µT follows from,

$$\begin{array}{}\text{(A6)}& {\mathit{\sigma }}_{\mathit{\epsilon },\mathrm{T}}=\left|{\displaystyle \frac{\partial {\mathit{\epsilon }}_{\mathrm{\mu }\mathrm{T}}}{\partial k_{\mathrm{B}}}}\right|\cdot {\mathit{\sigma }}_{k_{\mathrm{B}}}=\mathrm{4}\mathit{\nu }{\mathit{\kappa }}_{\mathrm{T}}^{\mathrm{2}}{k}_{\mathrm{B}}^{\mathrm{3}}\cdot {\mathit{\sigma }}_{k_{\mathrm{B}}}.\end{array}$$

To account for the approximate log-normal distribution of ε_µT, the variance of those estimates is given by,

$$\begin{array}{}\text{(A7)}& {\mathit{\sigma }}_{\ln \mathit{\epsilon },\mathrm{T}}^{\mathrm{2}}=\ln \left(\mathrm{1}+{\displaystyle \frac{{\mathit{\sigma }}_{\mathit{\epsilon },\mathrm{T}}^{\mathrm{2}}}{{\mathit{\epsilon }}_{\mathrm{\mu }\mathrm{T}}^{\mathrm{2}}}}\right).\end{array}$$

Finally, the 95 % confidence interval is given by,

$$\begin{array}{}\text{(A8)}& {\mathrm{CI}}_{\mathrm{95}\mathrm{\%},\mathrm{T}}={\mathit{\epsilon }}_{\mathrm{\mu }\mathrm{T}}\exp \left(\pm \mathrm{1.96}{\mathit{\sigma }}_{\ln \mathit{\epsilon },\mathrm{T}}\right).\end{array}$$

Appendix B: The finescale parameterization

The finescale strain parameterization is applied using Eq. (1). Here a more extensive description of the formula involved is given.

Figure B1 A $T/S$ diagram of the different stations. Background contours indicating σ₂ potential density are provided for reference.

B1 Obtaining the observed and model strain variance

As discussed in Sect. 4, the observed strain rate can be calculated by using either the buoyancy frequency or vertical temperature gradients (Kunze et al., 2006). Here we used both:

$$\begin{array}{}\text{(B1)}& {\mathit{\xi }}_{z,N^{\mathrm{2}}}={\displaystyle \frac{N^{\mathrm{2}}-N_{\mathrm{fit}}^{\mathrm{2}}}{\stackrel{\mathrm{‾}}{N_{\mathrm{fit}}^{\mathrm{2}}}}},\mathrm{and}{\mathit{\xi }}_{z,\mathrm{\Theta }}={\displaystyle \frac{{\mathrm{\Theta }}_z-{\mathrm{\Theta }}_{z,\mathrm{fit}}}{\stackrel{\mathrm{‾}}{{\mathrm{\Theta }}_{z,\mathrm{fit}}}}}.\end{array}$$

where Θ_z is the vertical derivative of Θ over the segment of interest. $N_{\mathrm{fit}}^{\mathrm{2}}$ and Θ_z,fit are a quadratic fit to the corresponding variable over the segment. $\stackrel{\mathrm{‾}}{N_{\mathrm{fit}}^{\mathrm{2}}}$ and $\stackrel{\mathrm{‾}}{{\mathrm{\Theta }}_{z,\mathrm{fit}}}$ are the mean of the quadratic fits. In this application, we have used vertical 200 m half overlapping segments, starting at a depth of 100 m, to avoid the surface mixed layer (the mixed layer is within the 50–80 m range). Also the two bottom most bins are omitted for boundary effects caused by the bottom boundary layer, where most internal wave breaking occurs.

The strain segments are first detrended and windowed using a Hanning window before being Fourier transformed to obtain the strain spectrum S_str of the segment. Each segment is integrated between wavelengths of 10–100 m to get the strain variance $\langle {\mathit{\xi }}_z^{\mathrm{2}}\rangle$. Each segment was integrated between 100 m wavelength and an upper limit being the highest possible wavelength between 10 and 40 m, while keeping $\langle {\mathit{\xi }}_z^{\mathrm{2}}\rangle \le \mathrm{0.2}$ as to avoid oversaturation of the spectrum.

The GM model spectrum (Garrett and Munk, 1975), following the GM76 version (Cairns and Williams, 1976), is calculated as (Gregg and Kunze, 1991),

$$\begin{array}{}\text{(B2)}& {\mathit{\xi }}_{z\mathrm{GM}}={\displaystyle \frac{E{b}^{\mathrm{3}}}{\mathrm{2}{j}_{*}\mathit{\pi }}}{\left({\displaystyle \frac{N_{\mathrm{0}}}{\stackrel{\mathrm{‾}}{N}}}\right)}^{\mathrm{2}}{\displaystyle \frac{\mathrm{1}}{{\left(\mathrm{1}+k_z/k_{z*}\right)}^{\mathrm{2}}}}.\end{array}$$

with,

$$\begin{array}{}\text{(B3)}& k_{z*}={\displaystyle \frac{\mathit{\pi }{j}_{*}}{b}}{\displaystyle \frac{N}{N_{\mathrm{0}}}}.\end{array}$$

The used constants are $j_{*}=\mathrm{3}$, b=1300 and $E=\mathrm{6.3}\times {\mathrm{10}}^{-\mathrm{5}}$. $\langle {\mathit{\xi }}_{z\mathrm{GM}}^{\mathrm{2}}\rangle$ is obtained by integrating the GM model spectrum between the same wavelengths as $\langle {\mathit{\xi }}_z^{\mathrm{2}}\rangle$.

B2 Functions for latitudinal dependence and R_ω

The function h(R_ω) describes the dependency of the ratio between shear and strain and is given by,

$$\begin{array}{}\text{(B4)}& h\left(R_{\mathit{\omega }}\right)={\displaystyle \frac{\mathrm{1}}{\mathrm{6}\sqrt{\mathrm{2}}}}{\displaystyle \frac{R_{\mathit{\omega }}\left(R_{\mathit{\omega }}+\mathrm{1}\right)}{\sqrt{R_{\mathit{\omega }}-\mathrm{1}}}}.\end{array}$$

The choice for a value of the shear-to-strain ratio R_ω in absence of shear data is discussed in Sect. 4. By choosing R_ω=3, Eq. (B4) reduces to 1.

The last term of Eq. (1) is a function to account for the latitudinal dependence of internal waves (Gregg et al., 2003). This latitudinal correction is,

$$\begin{array}{}\text{(B5)}& L\left(f,N\right)={\displaystyle \frac{f\mathrm{arccosh}\left(\frac{\stackrel{\mathrm{‾}}{N}}{f}\right)}{f_{\mathrm{30}}\mathrm{arccosh}\left(\frac{N_{\mathrm{0}}}{f_{\mathrm{30}}}\right)}}.\end{array}$$

Here f is the Coriolis parameter and f₃₀ the Coriolis parameter at a latitude of 30^∘.

B3 Additional figures

A $T/S$ diagram is provided to support the analysis in the main text.

Appendix C: The Thorpe resorting method

This appendix is concerned with the calculation of temperature overturns and the associated Thorpe lengthscale L_T (Thorpe, 1977) from the CTD sensors mounted on the VMP.

C1 Defining vertical averaging scale (δz) based on sensor characteristics

Below we determine a minimum vertical averaging length scale for averaging the measured data from the CTD, which depends on the noise that is dictated by the specifics of the used instrumentation and it interaction with the environment.

Figure C1 The overturn criteria applied on all detected overturns of a station. (a) Criterion 1: Thorpe length vs overturn size, (b) Criterion 2: measurement accuracy vs overturn size, (c) Criterion 3: stability ratio and correlation coefficient and (d) Criterion 4: relationship between Thorpe and Ozmidov scales. Red dots are flagged and will be removed from further calculations.

C1.1 Thinnest Overturns

The thinnest overturns that can be detected depend on the sampling frequency of the sensors and the fallspeed of the profiler. The used SBE3 temperature sensor and pressure sensor operated at a frequency of f_SBE=64 Hz. Then, the highest resolution sampled, L_z (m), is the fallspeed of the profiler v_vmp divided by the sampling frequency f_SBE. The maximum of L_z over the watercolumn is taken and multiplied by a factor 5 as a sort of “Nyquist frequency” (see Galbraith and Kelley, 1996), such that the thinnest detectable overturn size L_thinnest (m) is:

$$\begin{array}{}\text{(C1)}& L_z=\max \left[{\displaystyle \frac{\left(v_{\mathrm{vmp}}\right(z)}{f_{\mathrm{SBE}}}}\right]⟶L_{\mathrm{thinnest}}=\mathrm{5}{L}_z.\end{array}$$

We find L_z≈0.015 m and thus L_thinnest≈0.075 m.

Figure C2 Cross-validation of the different methods, as defined in Sect. 5.4, against the corresponding values of ε_vmp.

Figure C3 Dissipation estimates ε_TP for all different methods (ε_geo: Method 1, ε_Dil: Method 2, ε_Dilx4: Method 3 and ε_Dil,clim: Method 4, as defined in Sect. 5.4) are compared to ε_vmp and the used climatological estimates ε_clim, shown for three example stations.

C1.2 Response time of sensors

The overturns that can be detected also depend on the response time of the sensors. The response time of the temperature sensor is $t_{\mathit{\theta }}^{\mathrm{response}}=\mathrm{0.065}$ s (Seabird Scientific, 2025) and can be used to indicate over what vertical length scale the measurements are obtained (Galbraith and Kelley, 1996). The response time is recalculated into a vertical length scale L^response using:

$$\begin{array}{}\text{(C2)}& L_{\mathrm{\Theta }}^{\mathrm{response}}=\max \left[t_{\mathrm{\Theta }}^{\mathrm{response}}{v}_{\mathrm{vmp}}\left(z\right)\right]\approx \mathrm{0.06}\mathrm{m}.\end{array}$$

Note that $L_{\mathrm{\Theta }}^{\mathrm{response}}$ changes depending on the vertical speed of the VMP, which may differ somewhat per profile depending on stratification and the exact mass of the weights added on. However, based on L_thinnes, $L_{\mathrm{\Theta }}^{\mathrm{response}}$ we choose a vertical averaging scale of δ_z=0.08 m.

C2 Defining overturns and overturn-averaged values

C2.1 Defining overturns

The raw CTD data is bin-averaged to the averaging scale as defined above. Vertical overturns are defined based on this measured profile Θ(z) and on Θ_sorted. Here Θ_sorted is the sorted Θ(z) in descending order, representing a stably stratified background condition. For each Θ(z) one can find its exact equivalent in the profile Θ_sorted(z_sorted). The vertical displacement of each sample is given by the vertical difference between a value of Θ in both the measured and sorted profiles:

$$\begin{array}{}\text{(C3)}& \mathit{\delta }z=z_{\mathrm{sorted}}-z.\end{array}$$

This gives a profile of the vertical displacement of each sample. Overturns are defined by the zero crossings of the vertical integral of this displacement profile (Thorpe, 1977; Mater et al., 2015). For each overturn, the Thorpe length scale is given as the root-mean-square of the vertical displacements Δz of all samples within that overturn:

$$\begin{array}{}\text{(C4)}& L_{\mathrm{T}}\left(l\right)=\sqrt{{\displaystyle \frac{\mathrm{1}}{M_l}}\sum _{m=\mathrm{1}}^{M_l}(\mathit{\delta }{z}_m{)}^{\mathrm{2}}}.\end{array}$$

where l indicated an overturn and $m=\mathrm{1}:M_l$ the index of the samples within the overturn. M_l thus depends on the size of the overturn and can vary between overturns.

The Ozmidov scale L_Oz(l) and dissipation rate ε_TP(l) per overturn are calculated using the overturn-averaged values of N² and ε_vmp. Both N² and ε_vmp are first interpolated to the averaging scale of the CTD data, before being averaged for each overturn. The same procedure is followed for calculating the stability ratio per overturn.

C2.2 Calculating overturn-averaged variables

To calculate N² we first average the CTD data of Θ and S_A and p into 1 m vertical bins. This generally already averages over about 80 % of the overturns (not shown). Subsequently, we stabilize the profile over 50 m segments using the method from Barker and McDougall (2017) embedded in the TEOS-10 GSW-software (McDougall and Barker, 2011; IOC et al., 2010). We then interpolate the stabilized N² to the original vertical grid (of 8 cm) and only use the middle 30 m of the 50 m section. This is repeated to cover the whole vertical profile. Using the 1 m-binned data, the vertical gradients of Θ and S_A are calculated and subsequently interpolated onto the 8 cm grid. These are then used to calculate the stabillity ratio and averaged over an overturn.

C3 The measurement accuracy used for criterion 2

The accuracy or precision of the measurements have an influence on the vertical scale over which measurements are significantly distinguishable, and thus a minimal vertical length overturns need to have. Using the precision of $t_{\mathit{\delta }}=\mathrm{5}\times {\mathrm{10}}^{-\mathrm{5}}$ K for SBE3 temperature sensor (Seabird Scientific, 2025), we find

$$\begin{array}{}\text{(C5)}& {\mathit{\delta }}_{\mathrm{\Theta }}={\displaystyle \frac{\max \left[\mathrm{\Theta }\left(S_{\mathrm{A}},t+t_{\mathit{\delta }},p\right)-\mathrm{\Theta }\left(S_{\mathrm{A}},t-t_{\mathit{\delta }},p\right)\right]}{\sqrt{N_{\mathrm{O},l}}}}.\end{array}$$

Here δ_Θ is the error estimate on each Θ measurement over an overturn. We divide this error by the background gradients to obtain a vertical length scale which one must exceed to be able to distinguish one measurement from the other:

$$\begin{array}{}\text{(C6)}& L_{\mathrm{\Theta }}^{\mathrm{error}}={\mathit{\delta }}_{\mathrm{\Theta }}\times {\left({\stackrel{\mathrm{‾}}{{\displaystyle \frac{\partial \mathrm{\Theta }}{\partial z}}}}^l\right)}^{-\mathrm{1}}.\end{array}$$

The resulting criteria mean that the size of the overturns needs to be larger than the length scale associated with this error, which is 𝒪(10⁻⁴–10⁻²) m (Fig. C1).

C4 Additional figures

Based on the criteria defined in Sect. 5.2, part of the detected overturns are deemed unsuitable for estimating ε_TP. The effect of the criteria is shown in Fig. C1. Here station NP1 is taken as an example, but other stations show similar patterns.

Table 2 discussed the accuracy of the methods 1–4. In Fig. C2 a cross-validation of the estimates against ε_vmp are shown in support of Table 2.

Figure 6 showed the estimates of ε_TP for all four of the discussed methods, but with the estimates binned in halfoverlapping 200 m bins for readability. The estimates at the 25 m resolution as was used in the analysis, is shown for three example stations in Fig. C3.

Appendix D: The Triple Decomposition framework

Here we will discuss the tracer variance budget and how background profiles are obtained from an individual profile.

Figure D1 Estimated isopycnal diffusivities using the triple decomposition framework, based on both temperature and salinity data. Estimates are compared to climatological estimates of K_SG of Groeskamp et al. (2020). Depths for which $\left|{\mathrm{\nabla }}_{\parallel }{\mathrm{\Theta }}^m\right|\le {\mathrm{10}}^{-\mathrm{7}}$ are shaded in gray.

Assuming stationarity, homogeneity and non-divergence, the simplified tracer variance budgets, connecting the large-scale to the mesoscale and the mesoscale to the microscale, become (Garrett, 2001; Ferrari and Polzin, 2005)

$$\begin{array}{}\text{(D1a)}& {\displaystyle }\langle {\mathit{u}}^e{\mathrm{\Theta }}^e\rangle \cdot \mathrm{\nabla }{\mathrm{\Theta }}^m-\langle \stackrel{\mathrm{̃}}{{\mathit{u}}^t{\mathrm{\Theta }}^{\mathrm{t}}}\cdot \mathrm{\nabla }{\mathrm{\Theta }}^e\rangle & {\displaystyle }=\mathrm{0}\text{(D1b)}& {\displaystyle }\langle {\mathit{u}}^t{\mathrm{\Theta }}^{\mathrm{t}}\rangle \cdot \mathrm{\nabla }{\mathrm{\Theta }}^m+\langle \stackrel{\mathrm{̃}}{{\mathit{u}}^t{\mathrm{\Theta }}^{\mathrm{t}}}\cdot \mathrm{\nabla }{\mathrm{\Theta }}^e\rangle & {\displaystyle }=-{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\langle \mathit{\chi }\rangle .\end{array}$$

In the context of the triple decomposition framework, we will use angled brackets 〈⋅〉 to denote an average over a scale large compared to the scales within the brackets, but small compared to one scale larger in the framework. Summing these two equations results in,

$$\begin{array}{}\text{(D2)}& \underset{P_{{\mathrm{\Theta }}^{\mathrm{2}}}^{⟂}}{\underbrace{\mathrm{2}\langle {\mathit{u}}^t{\mathrm{\Theta }}^{\mathrm{t}}\rangle \cdot {\mathrm{\nabla }}_{⟂}{\mathrm{\Theta }}^m}}+\underset{P_{{\mathrm{\Theta }}^{\mathrm{2}}}^{\parallel }}{\underbrace{\mathrm{2}\langle {\mathit{u}}^e{\mathrm{\Theta }}^e\rangle \cdot {\mathrm{\nabla }}_{\parallel }{\mathrm{\Theta }}^m}}=-{\mathit{\chi }}_{\mathrm{\Theta }}.\end{array}$$

Here ∇_⟂ and ∇_∥ are the across and along isopycnal gradient operators. Equation (D2) states that the sum of the variance production on the microscale ($P_{{\mathrm{\Theta }}^{\mathrm{2}}}^{⟂}$) and the variance production on the mesoscale ($P_{{\mathrm{\Theta }}^{\mathrm{2}}}^{\parallel }$) are balanced by the variance dissipation (χ_Θ) as measured by the microstructure profiler. Using a gradient-flux approximation (Osborn and Cox, 1972), the production terms of Eq. (D2) can be written in terms of the isoneutral and dianeutral diffusivities;

$$\begin{array}{}\text{(D3)}& \mathrm{2}K|{\mathrm{\nabla }}_{\parallel }{\mathrm{\Theta }}^m{|}^{\mathrm{2}}+\mathrm{2}D|{\mathrm{\nabla }}_{⟂}{\mathrm{\Theta }}^m{|}^{\mathrm{2}}={\mathit{\chi }}_{\mathrm{\Theta }},\end{array}$$

where K is the isoneutral diffusivity and D the dianeutral diffusivity. This approach assumed that dianeutral tracer gradients are much larger than isoneutral gradients and that the isoneutral slopes are approximately flat.

D1 Assessing background gradients from CTD data

In order to calculate the dianeutral production term, it is needed to obtain dianeutral gradients of the background scale temperature distribution from single CTD profiles. For this, a rationale similar to Castro et al. (2024) is followed. As the profiles are too long to have an appropriate fit over the entire profile, a fourth order polynomal is fitted to $\stackrel{\mathrm{‾}}{\mathrm{\Theta }}\left(z\right)$ over a 200 dbar window. At the window midpoints, the vertical derivative of the of fitted polynomal is taken and stored for the background profile. For the upper and lower ends of the profile, within half a window size of boundary, the derivative of the first (cq. last) fitted polynomal is taken.

D2 Isoneutral diffusivities from the triple decomposition framework

The triple decomposition framework can be leveraged to obtain a first order estimate of the isoneutral diffusivity K (m² s⁻¹) from the microstructure data. By assuming a flux-gradient approximation, the isopycnal production term $P_{{\mathrm{\Theta }}^{\mathrm{2}}}^{\parallel }$ can be written as,

$$\begin{array}{}\text{(D4)}& P_{{\mathrm{\Theta }}^{\mathrm{2}}}^{\parallel }=\mathrm{2}\langle {\mathit{u}}^e{\mathrm{\Theta }}^e\rangle \cdot {\mathrm{\nabla }}_{\parallel }{\mathrm{\Theta }}^m=-\mathrm{2}K|{\mathrm{\nabla }}_{\parallel }{\mathrm{\Theta }}^m{|}^{\mathrm{2}},\end{array}$$

where ∇_∥ indicate the along isopycnal tracer gradients. Combining Eq. (D2) with the Eq. (12), for the dianeutral diffusivity and the expression for the isoneutral production above, results in an expression for the isoneutral diffusivity K (m² s⁻¹);

$$\begin{array}{}\text{(D5)}& K={\displaystyle \frac{{\mathit{\chi }}_{\mathrm{\Theta }}-\mathrm{2}D|{\mathrm{\nabla }}_{⟂}{\mathrm{\Theta }}^m{|}^{\mathrm{2}}}{\mathrm{2}|{\mathrm{\nabla }}_{\parallel }{\mathrm{\Theta }}^m{|}^{\mathrm{2}}}}.\end{array}$$

For the isoneutral tracer gradients (∇_∥Θ^m) we apply the VENM method of Groeskamp et al. (2019) to World Ocean Atlas 2018 1^∘ resolution hydrological dataset (WOA18, Locarnini et al., 2019; Zweng et al., 2019), providing the appropriate resolution for the isoneutral background gradients. The difference term, ${\mathit{\chi }}_{\mathrm{\Theta }}-\mathrm{2}D|{\mathrm{\nabla }}_{⟂}{\mathrm{\Theta }}^m{|}^{\mathrm{2}}$, is calculated at CTD-resolution, before being interpolated to WOA18 depths and divided by the isoneutral gradients over the range available (which is not always full depth, e.g. st. M2 in Fig. D1). When ${\mathit{\chi }}_{\mathrm{\Theta }}/\mathrm{2}$ and $D|{\mathrm{\nabla }}_{⟂}{\mathrm{\Theta }}^m{|}^{\mathrm{2}}$ are of similar magnitude (at CTD resolution), the latter can occasionally be larger than $\mathit{\chi }/\mathrm{2}$. This would lead to a negative K, such that we set the minimum to 0. After repeating this exercise for salinity data, this results in two sets of mesoscale diffusivity estimates K (Fig. D1).

Most estimates of K are centered between 𝒪(10²−10⁴), but there is a large degree of variation over 5 orders magnitude. Hence, the resulting K can at best be considered a first order estimate of K, and are of similar magnitude and variabillity as estimated by Spingys et al. (2021) and Naveira Garabato et al. (2016) using a similar method. We find that K is surface-intensified and roughly follows the same pattern as the climatological estimates of Groeskamp et al. (2020) in the upper 2000 dbar. Below that, the two estimates deviate substantially. Although the estimates of Groeskamp et al. (2020) are based on the first surface modes that are forced to decay to zero at the bottom by construction (LaCasce and Groeskamp, 2020), the large estimates of K observed below 2500 dbar (especially St. NP4–NP8) are, when compared to previous diffusivity estimates (see e.g. Abernathey et al., 2021, Table 9.1, and references therein), very high. In fact, they may result from (very) small isoneutral background temperature gradients. This is emphasized by the empirical threshold of $\left|{\mathrm{\nabla }}_{\parallel }{\mathrm{\Theta }}^m\right|\le {\mathrm{10}}^{-\mathrm{7}}$ K m⁻¹, marked by the gray shading in Fig. D1, providing a rough estimate (based on Groeskamp et al., 2019) where the estimated gradients are very small and subject to noise such that estimates of K are easily overestimated. Although some previous studies suggest a reasonable comparison between the isoneutral production and the sum of the dianeutral production and dissipation terms (Cherian et al., 2024), we here conclude that these estimates of the isoneutral diffusivity K are not very accurate. They show variations over multiple orders of magnitude, with extremely enhanced values in the deep ocean due to small isoneutral gradients. Even in the upper part of the water column, where these estimates seem to be of some value, they vary over many orders of magnitude of small ranges of depth, possibly indicating some limitation of the triple decomposition method to estimate mesoscale diffusivites. The noise in the estimates can in part be caused by the fact that these are estimates based on single profiles. It would help to average the estimates over appropriate time and length scales that fit the mixing processes associated with the isoneutral diffusivity. We conclude that the reliability of isoneutral diffusion estimates may be limited in regions of limited stirring signal, such as the ones sampled in this study. A more thorough quantification of signal-to-noise levels in the decomposition of variance dissipation rates, and in the estimation of isoneutral gradients, would be required to fully assess the reliability of the derived isoneutral diffusivities, which is beyond the scope of this study. However, the method is likely to produce better results in areas with stronger tracer gradients, where the finescale stirring signal is more prominent, such as the Antarctic Circumpolar Current (Naveira Garabato et al., 2016; Merrifield et al., 2016; Orúe-Echevarría et al., 2023) or subpolar regions (Castro et al., 2024).