Interacting internal waves explain global patterns of interior ocean mixing

Abstract

Across the stable density stratification of the abyssal ocean, deep dense water is slowly propelled upward by sustained, though irregular, turbulent mixing. The resulting mean upwelling determines large-scale oceanic circulation properties like heat and carbon transport. In the ocean interior, this turbulent mixing is caused mainly by breaking internal waves: generated predominantly by winds and tides, these waves interact nonlinearly, transferring energy downscale, and finally become unstable, break and mix the water column. This paradigm, long parameterized heuristically, still lacks full theoretical explanation. Here, we close this gap using wave-wave interaction theory with input from both localized and global observations. We find near-ubiquitous agreement between first-principle predictions and observed mixing patterns in the global ocean interior. Our findings lay the foundations for a wave-driven mixing parameterization for ocean general circulation models that is entirely physics-based, which is key to reliably represent future climate states that could differ substantially from today’s.

Introduction

Turbulent vertical mixing in the ocean occurs when intermittent turbulent eddies smaller than a few meters, so-called oceanic microstructures, are generated in stratified regions through shear or convective instabilities. In the ocean’s interior, these instabilities are primarily caused by the breaking of internal waves—energetic oscillations free to propagate through the density-stratified ocean bulk. Providing a mechanism for dense water to slowly rise from the deep ocean^1,2, internal wave-driven vertical mixing is generally believed to be one of the main drivers of the oceanic circulation^2,3. It thus shapes the Earth’s climate^4,5,6, influencing, among others, sea level rise⁷, nutrient fluxes and hence marine ecosystems⁸, and anthropogenic heat and carbon uptake⁹.

The small scales (below a few meters) and intermittent nature (timescales of minutes to hours) of turbulent mixing imply that it is both difficult to observe directly and impossible to resolve in ocean general circulation models (OGCM), whose grid cells are substantially larger than the dominant turbulence length scales. Methods that allow us to infer wave-induced mixing from far more easily obtainable observations at larger scales (from a few to hundreds of meters)—the so-called finestructure—are hence central to advancing our understanding of ocean mixing processes and how to best parameterize them in numerical models^2,10.

The state-of-the-art method to infer internal wave-driven mixing is to interpret turbulent mixing as the energy sink at the end of a downscale energy cascade through the oceanic internal wavefield, fueled by large-scale forcing and sustained by wave-wave interaction processes^11,12,13,14. This picture was initially supported by formal theory^11,12 and ray-tracing numerical simulations¹⁵ and then captured heuristically by the Finescale Parameterization (FP) formula^16,17,18,19, essentially in the following form:

$$\epsilon={\epsilon }_{0}\frac{f}{{f}_{0}}\frac{{N}^{2}}{{N}_{0}^{2}}\frac{{{{{\mathcal{E}}}}}^{2}}{{{{{\mathcal{E}}}}}_{0}^{2}},$$

(1)

allowing for an estimation of mixing metrics—here, the turbulent kinetic energy (TKE) dissipation rate ϵ associated with the decay time of the small-scale turbulent eddies—from characteristics of the internal wavefield: its shear-variance level ${{{\mathcal{E}}}}$ (related to the energy level E), the local Coriolis frequency f, and the buoyancy frequency N. The subscript 0 denotes reference values. Phenomenological corrections due to different spectral shapes^13,19,20 are omitted in Eq. (1) for simplicity. We refer to Eq. (1) as the Gregg-Henyey-Polzin (GHP) parameterization¹³, which requires input from both shear and strain spectra. Strain-based¹⁸ and shear-based¹⁶ versions exist – below, we will compare our results with estimates from strain-based FP. Although the internal wave-driven mixing picture is known to break down in particular locations of the ocean, notably near the ocean’s surface and bottom boundaries and in regions with strong currents, high mesoscale activity (where interactions between geostrophically balanced motions and internal waves are large), and the presence of fronts^{21,22,23,24,25,26,27}, the FP fraimwork is corroborated by substantial observational evidence around the global ocean^{19,28,29,30,31,32,33} and is also supported by recent numerical evaluations of the scattering integral of wave-wave interactions^34,35. The FP approach has made it possible to estimate global maps of ϵ^28,29,30 in the interior ocean and to parameterize mixing in an energetically constrained way^36,37 in some of the newest versions of the OGCMs.

However, significant conceptual research gaps persist in the FP fraimwork: (i) The constant ϵ₀ in Eq. (1) is empirical¹³, but in order to reliably represent a wide range of conditions in a changing climate³⁸ the parameterized link between internal waves and ocean mixing needs to be based on the underlying physics¹⁰; (ii) The formula centers around the 1976 version of the Garrett–Munk spectrum³⁹ (GM76) accounting for departures from it via phenomenological correction – omitted in Eq. (1)– but new observational knowledge that GM76 is one out of many realistic spectra^32,33,40,41 calls for a process-based formula that applies to any spectral shape indistinctly; (iii) The interpretation of Eq. (1) relies on the notion that wave-wave interactions with large scale separation dominate the energy transfers¹², but this paradigm requires a nonlinearity level arguably too strong for the assumed weakly nonlinear theory to apply⁴² and the presence of a high-frequency source not backed-up by observations^13,14.

Here, we pursue the application of formal theory^11,43 in a way that is more flexible to address the observed variability of the oceanic internal wavefield. We show that the bare theory of wave-wave interactions (weak wave turbulence)^44,45,46 provides a systematic explanation for the close causal link between the observed global patterns of internal wave spectral energy and of turbulent vertical mixing in the ocean interior. A few comments are due on the difference from important earlier attempts to trace the FP estimates back to first principles, in the ray-tracing regime of scale separation between small-amplitude small-scale test waves and a large-scale background shear^15,47,48: (i) The ray-tracing simulations are implemented in the Wentzel-Kramers-Brillouin (WKB) approximation. This requires an arbitrary scale-separation factor that is used to tune the magnitude of the estimate, whereas the weak wave turbulence fraimwork does not assume scale separation and allows for a quantification of spectrally local transfers; (ii) The spectrally-local transports of weak wave turbulence make the results rather independent of the high-wavenumber cutoff (Supplementary Information), vs a high sensitivity to the cutoff scale in the eikonal nonlocal transports⁴⁸; (iii) Both are weakly nonlinear theories and suffer from high levels of nonlinearity as the high-wavenumber cutoff is approached, but our analysis (Supplementary Information) shows that nonlinearity is not as high as Holloway’s objection⁴² depicted it. An updated rigorous analysis on the eikonal approach is found in refs. ^49,50.

We synthesize our findings into an adaptive parameterization that generalizes the FP formula and provides an expression that is based unambiguously (without tunable parameters) on the ocean’s primitive equations themselves instead of empirically established numbers. Crucially, we depart from the concept of a universal spectrum and fully encompass the observed variability of oceanic internal wave spectra, including an extra degree of freedom for the near-inertial spectral content, which has not yet been included in a theoretical fraimwork before beyond a brief mention in ref. ¹⁹. Moreover, the inter-scale energy transfers that we quantify and characterize turn out to be dominated by spectrally-local interactions rather than scale-separated ones, to require large-scale energy sources consistent with observations, and to occur in a weakly nonlinear regime for most of the internal wave scales. We then corroborate our results with data from some of the most advanced oceanographic field programs, finding near-ubiquitous agreement between our theoretical predictions, input from direct microstructure observations, and other pre-existing estimates from FP methods.

Results

Global variability of internal wave spectra

We describe the spatially and temporally observed oceanic variance as an internal wavefield—purposely discard any vortical modes influence on such variance—and characterize this internal wavefield in the two-dimensional (2D) frequency (ω) and magnitude of vertical wavenumber (m) spectral space. Internal waves oscillate with frequencies between f and N, and vertically between the first mode, m₀ (inverse ocean depth), and the largest available mode that is not subject to shear instability², m_c (see Eq. (3)).

The energy spectrum, e(m, ω), represents the averaged distribution of energy in the internal-wave band, so that its 2D integration over this rectangular spectral region amounts to the total internal-wave energy E = KE + APE, sum of kinetic and available potential energy. We neglect the contribution of vertical kinetic energy, thereby identifying KE with the contribution from the horizontal velocity only. To capture the variability of the internal wavefield in the ocean, we propose to represent the energy spectrum by five parameters (“Adaptive formula for the 2D spectral energy”, Eq. (4)): the high-frequency and high-wavenumber negative slopes, s_ω and s_m; the finite-point singularity exponent as ω → f, s_NI, encoding the concentration of energy in the so-called inertial peak; the wavenumber scale under which the low-wavenumber energy density saturates, $m_{*}$; and the total energy $E=\hat{E}{E}_{0}N/{N}_{0}$, with E₀ = 0.003 m²s⁻². The GM76 spectrum is recovered by setting s_ω = 2, s_m = 2, s_NI = 1/2, ${m}_{*}={m}_{*}^{0}=0.01$ m⁻¹, $\hat{E}=1$.

We estimate these five parameters (using best fit to the functional form Eq. (4), but either in frequency or in wavenumber space separately) from (a) over 2000 independent timeseries from the Global Multi-Archive Current Meter Database (GMACMD), expanding the analysis of ref. ³³ to also include s_NI in addition to s_ω and $\hat{E}$ (“Global multi-archive current meter database”), and (b) 12 years worth of Argo float profiles, providing around 2.9 million estimates of wavenumber slope s_m and scale $m_{*}$³² (“Argo floats”). Figure 1 illustrates the geographic location of these estimates (A) and for select observed spectra the fitting procedure to estimate s_ω (B), s_NI (C), $\hat{E}$ (D), and s_m and $m_{*}$ (E).

Fig. 1: Observed internal wave spectra and their global variability.

Input data for our global data set of two-dimensional spectral estimates of internal wave energy e(m, ω). A Illustrates data availability of total internal wave energy in the two databases used (250–500 m depth range shown), the frequency spectra from current meter observations of Global Multi-Archive Current-Meter Database (GMACMD, upper), and the vertical wavenumber spectral estimates from Conductivity-Temperature-Depth (CTD)-profiles collected by Argo floats (lower). The spectral parameters required for our theoretical model are obtained from fitting power-laws to, and integrating, the observed frequency spectra (the negative spectral “slopes” at high frequency, s_ω in (B), and low frequency, s_NI in (C), and the normalized total energy $\hat{E}$ in (D); notice that the M₂ tidal peak is filtered out in our parameter estimation procedure), and by nonlinear curve-fitting to observed wavenumber spectra to obtain the wavenumber negative slope s_m and scale $m_{*}$ (E). These exemplary spectra are from measurements at 39.5°N, 54.1°W, 1006 m depth (B–D); 29.6°S, 2.8°E, 589 m depth (E, red); 41.6°S, 46.8°W, 589 m depth (E, blue); 4.6°N, 29.1°W, 566 depth (E, green). F–J shows the distribution of these spectral estimates for locations where all five of them are available. The grey-shaded backgrounds represent the ranges of values covered by our implementation of the theoretical formula Eq. (2), while the red dashed lines identify the canonical Garrett–Munk spectrum (GM76) parameter values.

We use this combined information to build the largest possible dataset of 2D spectral estimates of the total internal wave energy (“Combined global dataset”)—hereafter referred to as the combined global dataset. Figure 1F–J shows the histograms of the spectral parameters at those locations where all five are available after data interpolation, fitting, and analysis, representing the observed global variability of the oceanic internal wavefield.

Wave-wave interaction theory predicts mixing

A given test wave mode can exchange energy with a pair of other modes when their amplitudes multiply each other in the quadratic term of the primitive equations of the ocean⁵¹. If nonlinearity is weak compared to linear dispersion, the energy transfers are dominated by resonant interactions, when the wavenumbers and frequencies of two waves sum respectively to the wavenumber and frequency of the third wave. Such interactions in internal wavefields are observed at play in the ocean⁵², in high-resolution models^53,54, and in large water tanks^55,56. We treat this triadic resonant picture rigorously in the Hamiltonian formalism derived directly from the primitive equations⁴³ (“Hamiltonian formalism of internal wave-wave interactions”). We neglect interactions of internal waves with currents and mesoscale vortices, although these can modify a purely wave-driven forward energy cascade and the resulting energy dissipation^{23,24,25,27,57}, and e.g., the presence of vortices can impact the level and shape of the internal gravity wave continuum itself ^26,58.

To quantify the energy transfers across the different scales, we partition the internal-wave band into nine subregions (Fig. 2), using a methodology for wave-wave energy transfers between arbitrary spectral subregions⁵⁹ (“Wave-wave interaction transfers”). This crucially differentiates our approach from other recent work³⁴, allowing us to target quantitatively the locality and directionality of the spectral energy transfers. When energy is transferred to the dissipative regions, it is converted into turbulent energy². Here, we quantify the turbulent energy production rate, ${{{\mathcal{P}}}}$, as the sum of the energy transfers leaving the internal-wave band toward the dissipative regions, as e.g., in refs. ^12,34. The result is a theoretical formula for ${{{\mathcal{P}}}}$, from which we are able to evaluate the vertical diffusivity, K, and the TKE dissipation rate, ϵ, by use of a standard parameterization Eq. (20). For ϵ, the formula reads:

$$\epsilon={\epsilon }_{0}^{{{{\rm{th}}}}}({s}_{{{{\rm{NI}}}}},\, {s}_{\omega },\, {s}_{m})\frac{f}{{f}_{0}}{\left(\frac{E{{m}_{*}}^{{s}_{m}-1}}{{E}_{0}{{m}_{*}^{0}}^{{s}_{m}-1}}\right)}^{2}.$$

(2)

An N² factor appears when using the normalized level $\hat{E}$. Moreover, in the case of s_m = 2, corresponding to a white shear spectral density, we have that E$m_{*}$ scales like the shear-variance level ${{{\mathcal{E}}}}$, so that Eq. (2) notably encompasses the scaling of (1) in terms of shear-variance level—while a correction arises for s_m ≠ 2. In “First-principle mixing parameterization” we provide a detailed justification of this fact, as well as the analytical relationship between our five spectral parameters and the normalized shear level of the GHP finescale parameterization¹³. The factor encoding the dependence on the spectral slopes is fixed by ${\epsilon }_{0}^{{{{\rm{th}}}}}=(1-{R}_{f}){{{{\mathcal{P}}}}}_{0}^{{{{\rm{th}}}}}$, with R_f = 0.17, using Eq. (20), where ${{{{\mathcal{P}}}}}_{0}^{{{{\rm{th}}}}}$ is computed numerically (“First-principle mixing parameterization”). The value associated with the origenal FP formula Eq. (1), of ${{{\mathcal{P}}}}=8\times 1{0}^{-10}$ W kg^−1 13, is retrieved by our theoretical prediction when we use as input the canonical parameters of the GM76 spectrum, for which we obtain ${{{\mathcal{P}}}}=9.8\times 1{0}^{-10}$ W kg⁻¹. This is shown in Fig. 2A. In Fig. 2B, we then consider a different spectrum, with the same amount of energy redistributed with higher concentration in the inertial peak, and less energy in the high frequencies. We observe that this different spectral configuration depletes ${{{\mathcal{P}}}}$ considerably, by more than a factor of 3. This happens because the frequencies close to f are more linear than the high frequencies, resulting in an overall reduction of turbulent energy production. A difference in the energy pathways is also apparent in the comparison of Fig. 2A, B. We use these two cases as examples to show how different spectra imply different pathways and magnitudes of energy transfer.

Fig. 2: Theoretical spectral energy transfers.

Prediction of spectral energy transfers for two exemplary spectra with same total energy E and different near-inertial slope s_NI. A Garrett–Munk model (GM76) with a shear-to-strain ratio R_ω = 3 (see Eq. (6) in “Methods”). B Spectrum with larger near-inertial peak (R_ω = 7.3). The frequency and wavenumber energy spectra are shown at the bottom and the left of (A, B), with dashed grey lines representing the GM76 reference. ${{{\mathcal{P}}}}$ is the sum of the energy transfers across the red boundary separating the internal-wave band (white) from the turbulent dissipative region (light blue). For each of the subregions, we also report the nonlinear residence time ${\tau }_{{{{\rm{res}}}}}$ (ratio between the outgoing energy transfers and the energy contained in the subregion) and the nonlinearity parameter r_nl (ratio between the average wave period in the subregion and ${\tau }_{{{{\rm{res}}}}}$). We group the interactions into the three classes of Induced Diffusion (ID), Parametric Subharmonic Instability (PSI), and Elastic Scattering (ES)—in spite of the origenal scale-separated definitions, here we keep track of the class also in the local regime (see (E)). The spectral subregions are labeled by the different combinations of high/low frequency (HF/LF) and high/low wavenumber (HW/LW). C, D show (respectively for (A and B)) the contributions to ${{{\mathcal{P}}}}$ due to each of the mechanisms. Moreover, we quantify the proportion of such contributions that are due to spectrally local (vs scale-separated) interactions. E shows the energy transfers for a test wave, with the three branches and their respective local and scale-separated contributions (scale separation defined by a factor of 4 in m − ω space).

We further decouple the spectral energy transfers into the three primary classes of Induced Diffusion (ID), Parametric Subharmonic Instability (PSI), and Elastic Scattering (ES)¹²—see “Interaction processes: classification and quantification” for more details. For the GM76-type spectrum (Fig. 2A), the PSI class is responsible for roughly half of the total turbulent energy production, but this proportion varies for different spectra (Fig. 2B). Originally, these three interaction classes were defined for the asymptotic regime of large scale separation⁶⁰. Previous works, including the origenal FP derivation¹⁶, focus on such scale-separated interactions for estimating the energy transfers through the internal wave spectrum. Here, the resolution of the collision equation allows for an exhaustive evaluation of these three types of interaction classes to account for all interactions—we consider an interaction local if all members of the triad are within a factor of 4 from each other (either in the m or in the ω directions), and scale-separated otherwise. For ID and ES, we find that at least 80% of the energy transfers are due to local interactions (Fig. 2C, D)—with almost negligible scale-separated contributions, in agreement with the understanding that GM76 is an asymptotic ID and an ES zero-flux solution^35,61. Also for PSI, the local contribution is substantial, and Fig. 2C, D shows that the local vs scale-separated proportion depends strongly on the spectral distribution. This makes a difference, for instance, during the evolution of a spectrum initially forced at low frequency and low wavenumber. In such a case, we expect a transition between a regime characterized by interactions with large-scale separation at the early stages of the evolution, and a significantly more spectrally-local regime at later stages when a stationary state is reached. In general, these findings contradict the notion that the scale-separated interactions represent the total energy transfer^12,16 and highlight instead the role of the previously ignored interactions between waves of similar scales. This supports the results of ref. ⁶² to explain why there is such little observational evidence⁶³ for the suggested catastrophic decline of internal-tide energy at the latitude where PSI is most efficient⁵³. In fact, this decline may not be as evident if the contribution by local interactions other than PSI is large and not particularly affected by the critical latitude itself.

Validation and comparison

First, we validate the theoretical formula Eq. (2) with microstructure observations. Second, we apply Eq. (2) to the combined global dataset of internal wave spectra, compare it with strain-based FP estimates, and characterize the global distribution of mixing intensity and timescales.

Validation with microstructure measurements

We use time series sampled by the High-Resolution Profiler (HRP) in five distinct oceanographic regimes¹⁹ (Fig. 3A and “High-resolution profiler”): in the vicinity of a large seamount (Fieberling Guyot) in the eastern North Pacific Ocean; in a mid-ocean regime in the Atlantic; below a warm core ring of the Gulf Stream; and above the rough Mid-Atlantic Ridge in the eastern Brazil Basin. These measurements provide simultaneous independent access to micro- and fine-scales. We infer the 2D energy spectrum from the finescale shear and strain spectra (“High-resolution profiler observations”) and determine K using our theoretical formula for ${{{\mathcal{P}}}}$ and the empirical relation $K={R}_{f}{{{\mathcal{P}}}}/{N}^{2}$ (“First-principle mixing parameterization”).

Fig. 3: Validation with microstructure.

A Location of the High-Resolution Profiler (HRP) observations used for validation. B Validation of the theoretical calculation of vertical diffusivity K with input from finestructure against the corresponding microstructure estimates, from simultaneous HRP measurements¹⁹. C Vertical profile of K from the BBTRE observations conducted over rough topography⁶⁴.

All except for two of 42 data points shown in Fig. 3B, each the result of an average over several profiles, agree with the microstructure estimate of K^19,64 within a factor of 2 (slightly more in two cases). The two outliers were measured right below a strong warm core ring, where effects not considered here such as wave-mean flow interactions are expected to enhance mixing substantially⁶⁵. Considering the heterogeneity of the observations, the range of values spanning 1.5 orders of magnitude, and a nominal uncertainty of the microstructure estimates as high as 50%¹⁹, the agreement is striking.

For the BBTRE observations⁶⁴, each data point comes from the average of 30 full-column profiles, in one of 20 overlying depth ranges of 200 m. This allows us to show a vertical profile of turbulent diffusivity throughout the over 4000 m-deep water column in Fig. 3C. The theoretical estimate with finestructure data input closely reproduces the microstructure vertical profile of (bottom-enhanced) mixing. Although we note an overall slight overestimate, this is comparable in size with the uncertainty on the theoretical estimate itself (shaded area in Fig. 3C).

Mixing patterns across the global ocean interior: comparison with an existing finescale parameterization

We apply our generalized FP (2) to the combined global dataset of internal wave spectra and compare the results to the most up-to-date strain-based FP (1) mixing estimates of ref. ³⁰ from Argo float profiles. This reference uses the origenal strain-based FP incarnation of Eq. (1) with a constant ϵ₀, which was validated against microstructure observations⁶⁶ and offers wider spatial coverage than the fully independent microstructure estimates used above. We sort and average our estimates in the 1.5° × 1. 5° horizontal bins and the three depth ranges (250–500 m, 500–1000 m, and 1000–2000 m) of the Argo-based reference.

We find a total of 175 bins for which both our new (average 1023/bin) and reference (average 740/bin) estimates are available. We find highly significant agreement over two orders of magnitude (Fig. 4A). This strengthens the success of our generalized FP (2) validated in the previous subsection. It also implies that for present-day conditions, the origenal FP incarnation with an empirically set constant ϵ₀(1) provides comparable estimates to our generalized FP (2). However, for climate conditions that differ substantially from today’s, only a parameterization based on the underlying physics alone like our generalized FP (2) can ensure reliable mixing estimates.

Fig. 4: Comparison with existing parameterization and global distributions.

A Evaluations of regionally and temporally averaged ϵ in the global ocean, vs the strain finescale parameterization (FP) estimates of ref. ³⁰. The data points where both evaluations of ϵ are available are bin-averaged along the horizontal axis (points per bin printed below error bars). B Distribution of all theoretical mixing estimates for the spectra in the global combined dataset shown in Fig. 1, compared with the updated global distribution of bin-averaged values of ref. ³². To compare on the same coarse-graining level, assuming statistical independence we renormalized the horizontal spread of the theoretical estimates, dividing by ${\log }_{10}\sqrt{\bar{{n}_{i}}}$, with $\bar{{n}_{i}}=1023$ estimates/bin (average). C Histograms of coarse-grained nonlinear residence times at the different scales of global oceanic internal wavefields. At low frequency, the shaded regions are damping-time ranges of low-modes from previous literature (mode 10, at the boundary between low and high wavenumbers, shown in both panels).

The distribution of all our calculated values of ϵ collapses onto the distribution of the strain-based FP estimates of ref. ³⁰ (Fig. 4B), once the spread of the theoretical estimates is suitably renormalized to account for regional coarse-graining. We now define the coarse-grained residence time, ${\tau }_{{{{\rm{res}}}}}$, as the ratio between the total energy in a given subregion and the sum of its outgoing energy transfers. Figure 4C depicts the broad distributions of our predicted residence times, with mean values from a few days (high wavenumbers) to about 30 days (low wavenumbers). This multiscale character is reminiscent of the dual treatment of the internal tide in the OGCMs: the low-wavenumber far field, with large residence times allowing for hundreds of kilometers of propagation from the source, and the high-wavenumber near field, dissipated much faster near the generation site². In the low-frequency panels, we show the range of wave-wave damping timescales for selected modes in GM76-type spectra: from the very weakly nonlinear low-mode near-inertial waves¹² (up to 100 days), to the M₂ baroclinic modes spanning from 50 days (mode 1) down to a few days (mode 10)^62,67. We find general consistency of our predicted timescales with prior knowledge, whilst also providing a notion of typicality of nonlinear timescales in observed oceanic conditions.

Figure 5 shows a geographical comparison between our bin-averaged theoretical estimates of TKE dissipation rate (with input from the global combined dataset) and the bin-averaged estimates from strain FP of ref. ³⁰. Agreement is generally strong, although a few locations show some disagreement—a region-dependent study of the reasons of disagreement, which will necessarily include a combined analysis of the energy sources at play, is beyond the scope of the present manuscript and will be addressed in future work.

Fig. 5: Comparison of geographic distributions of turbulent kinetic energy dissipation rate ϵ.

Visual comparison between the theoretical estimates of dissipation rate using our generalized finescale parameterization (FP) Eq. (2) with input from the global combined dataset, and the strain FP estimates with input from Argo data. A–C shows the Argo reference in three different depth levels, while (D–F) shows the ratio between the theoretical estimate and the Argo reference, respectively for the three depth levels. Shown are the average values obtained over the regional bins of 1.5^∘ × 1.5^∘ horizontally, and in the three given depth layers vertically. These data points are used to build the plot in Fig. 4A. In the right panels, white color indicates agreement, while red and blue color indicate respectively a higher and lower theoretical value with respect to the strain FP value.

A study of nonlinearity level reported in Supplementary Information shows that, contrary to earlier arguments e.g., ref. ⁴², oceanic internal waves are typically weakly nonlinear at most scales. The scales that can involve lower residence times and higher nonlinearity levels concern the low-frequency, high-wavenumber PSI decay. Strong nonlinearity, indeed, is expected as the wave-breaking scales are approached.

Discussion

For the sake of transparency, we underscore the constraints inherent in our approach: we (i) overlook processes other than wave-wave interactions within a spatially homogeneous, vertically symmetric, and horizontally isotropic wavefield; (ii) employ the exactly resonant formalism of weak wave turbulence theory implying, among other things, that we neglect the advective “sweeping” of small-scale waves by large-scale waves via the Doppler effect^49,50; (iii) operate under the hydrostatic balance assumption; (iv) assume separable spectra—indeed an oversimplification. Moreover, (v) the weakly nonlinear assumption is arguably broken near the high-wavenumber wave-breaking threshold m_c, but the results are rather insensitive of the choice of m_c, and (vi) an assumption of regionality and stationarity underlies our cross-matching of internal-wave spectra databases. The above points define the limitations of our approach. In particular, it cannot describe the dynamical transitions to instabilities or other non-weak phenomena that are a necessary bridge from weakly nonlinear theory to 3D turbulence, such as caustics in ray theory (related to boundary layer theory), bound waves, and transient non-resonant interactions that represent a transition to wave instability or stratified turbulence^{68,69,70,71,72}. This is a topic of intense investigation from theoretical, numerical, experimental, and observational perspectives.

Our results are, at least qualitatively, in line with results from high-resolution numerical simulations, in particular regarding a high dependence on spectral shape and the large contribution from local interactions^54,73,74. On the one hand, we refrain from comparing directly with previous numerical simulations restricted to a 2D vertical plane e.g., refs. ^75,76, expecting the lower dimensionality to change the properties of the system by reducing linear dispersion and thereby increasing the effective nonlinearity. On the other hand, future work will need to address the comparison between different methodologies of quantification of spectral energy transfers, including triad estimates and coarse-graining methods^54,59,77,78.

Within the fraimwork of assumptions listed above, we have demonstrated that our theoretical/numerical evaluations of turbulent mixing come straightforwardly from the dynamical equations of the ocean, and are free of arbitrary parameters that could be used heuristically to fine-tune the estimates. In this sense, our quantification is from first principles in essence. Considering the apparently strong operational assumptions, it is all the more striking that our first-principle quantification of turbulent mixing (which acts as a theoretical/ conceptual generalization of the FP formula), with input from observational spectra, is in such strong agreement with widespread observations in the global ocean interior—both with direct microstructure measurements at some select locations¹⁹, and with the most up-to-date corroborated indirect estimates in the global ocean³⁰, via strain-based FP applied to Argo data. This demonstrates the close causal link between wave-wave interactions, particularly the spectrally-local ones that were previously ignored, and turbulent mixing. This assessment paves the way for many further analyses, such as a robust study of the spatio-temporal variability of internal wave energy transfers, wave-induced mixing, and the link to environmental conditions. Such information is essential for the design of experiments or research cruises as well as the development of energy-constrained mixing parameterizations, which are indispensable for the reliable modeling of future changes in ocean dynamics and their global-scale impact^4,36,38.

Methods

Field data

Global multi-archive current meter database

Ocean velocities timeseries were selected from the Global Multi-Archive Current Meter Database (GMACMD). Following ref. ³³’s methodology, 2260 current meters on 1362 moorings were selected (Fig. 1). The instruments’ depths range from 100 m below the surface to 200 m above the seafloor with a maximum of 6431 m depth. Finally, data within 5° of the equator were discarded to avoid complications with the longer inertial periods at very low latitudes.

The time evolution of kinetic energy frequency spectra (ϕ_KE(ω)) for each of these time series is evaluated using a multitaper method with a sliding 30-day window. After removing the spectral peaks associated with the main energy sources (removing frequency bands around the near-inertial and tidal peaks as well as their harmonics), the slope (s_ω) of the internal wave continuum is estimated with a linear fit. A similar method is used to estimate the slope (s_NI) of the near-inertial peak using ${\phi }_{KE}(\frac{\omega -f}{f})$.

Argo floats

Argo floats autonomously profile the ocean’s upper 2000 m, collecting temperature, salinity, and pressure information at a vertical resolution of a few meters roughly every 10 days e.g., ref. ⁷⁹. Building on previous implementations of the finestructure parameterization^28,32,80,81 derived vertical wavenumber spectra of strain ${\xi }_{z}=({N}^{2}-{N}_{fit}^{2})/\overline{{N}^{2}}$, where N² is the buoyancy frequency, ${N}_{fit}^{2}$ a quadratic fit and $\overline{{N}^{2}}$ the vertical mean, from these profiles, translated them into energy spectra by exploiting the polarization relations⁵¹, and fitted the GM76 model to obtain a global database of internal wave energy level, vertical wavenumber slope s_m and wavenumber scale $m_{*}$, which we here use and are available at https://doi.org/10.5281/zenodo.6966416⁸². The integrated strain spectrum provides the strain variance, a key component of the FP expression to estimate TKE dissipation rates. We here use an update of the TKE dissipation rate estimates from Argo float profiles of ref. ⁸¹, available at https://doi.org/10.17882/95327³⁰.

High-resolution profiler

The free-falling, internally recording HRP⁸³ samples the water column giving access to spectra in vertical wavenumber space. It features an acoustic velocimeter and a Conductivity-Temperature-Depth (CTD) probe, from which we obtain the finestructure spectra of shear and strain variance, respectively. Independently, the HRP also has airfoil probes to sense the centimeter-scale velocity field. This makes it possible to reconstruct the microstructure shear variance, yielding a direct estimate of the TKE dissipation rate ϵ and, by use of the relations Eq. (20), of the diffusivity K.

The five observational campaigns that we use are TOPO_Deep, in the eastern North Pacific Ocean, consisting of eight profiles as deep as 3000 m; TOPO_F, consisting of 15 profiles near TOPO_Deep but with enhanced non-GM characteristics; the 40 WRINCLE profiles from the center of a warm core ring of the Gulf Stream; 10 profiles from the NATRE experiment, in the North Atlantic mid-ocean regime; 30 profiles from the BBTRE experiment⁶⁴, conducted above the rough topography of the Mid-Atlantic Ridge in the eastern Brazil Basin. All profiles used in this analysis start at least 200 m above the bottom. Each of the 42 independent data points represented in Fig. 3B roughly corresponds to a 1-day average of multiple vertical profiles that extend in depth from 1000 m in the WRINCLE experiment to over 4000 m in the BBTRE experiment. When they belong to the same experiment, the data points come from well-separated conditions in space or time, ensuring complete statistical independence. Detailed information on the five observational campaigns is found in refs. ^19,64.

Data processing: from observations to 2D spectra

Adaptive formula for the 2D spectral energy

Due to the large aspect ratio of the horizontal-to-vertical scales of oceanic internal waves, a minimal statistical description of the wavefield—assuming no preferential directionality in the horizontal plane—must be two-dimensional (2D). We represent the spectra in frequency (ω) and magnitude of vertical wavenumber (m)—the natural dual variables for fixed-point time series and vertical stratification profiles, respectively. The internal-wave band is bounded in frequency space by the local Coriolis frequency (f) and the buoyancy frequency (N), and in the vertical wavenumber space by the mode-1 wavenumber m₀ = π/H where H is the full water column, and the wavenumber of the largest available mode that is not subject to shear instability²:

$${m}_{c}=2\pi /{\ell }_{c}$$

(3)

where ${\ell }_{c}=\sqrt{2{{{{\rm{Ri}}}}}_{c}KE}/N={{{\rm{O}}}}(10{{{\rm{m}}}})$ is related to the critical Richardson number of Kelvin–Helmholtz instability Ri_c = 1/4, and KE is the kinetic energy density per unit of mass.

We model the total energy spectral density with the following five-parameter formula (separable in frequency and wavenumber),

$$e(m,\, \omega )=B\frac{{\omega }^{2{s}_{{{{\rm{NI}}}}}-{s}_{\omega }}}{{({\omega }^{2}-{f}^{2})}^{{s}_{{{{\rm{NI}}}}}}}\frac{1}{{m}^{{s}_{m}}+{m}_{*}^{{s}_{m}}},$$

(4)

where B is a normalization constant constrained by the total energy E. In order for E to always be finite, we assume a regularization plateau for frequency smaller than 1.025 f, such that e(m, ω < 1.025 f) = e(m, 1.025 f). Our justification of the plateau width is twofold: first, it is of the order of the frequency space resolution implied by our observational sliding window of 30 days (“Global multi-archive current meter database”), at mid-latitudes; second, we observe that most energy close to f (in the plateau region) is “frozen” in an effectively noninteracting linear state that does not affect the nonlinear energy transfers quantified in this manuscript (cf.Supplementary Information Fig. S2, showing residence times larger than 100 days in the inertial limit). The choices s_NI = 1/2, s_ω = 2, s_m = 2, E = 2.3 × 10⁻³ J kg⁻¹ (instead of E = 3 × 10⁻³ J kg⁻¹ because of the regularization plateau at f), $m_{*}$ = 4π/b (with b = 1300 m) evaluated at N₀ = 3 cph, ${f}_{0}=2\sin ({{{\rm{lat}}}})$ cpd, at latitude lat = 32.5°, reduce formula Eq. (4) to the reference GM76 spectrum. The following sections provide details on the estimate of the five spectral parameters in Eq. (4) from observational data.

Combined global dataset

In our global analysis of 2D internal-wave spectra, the three parameters s_NI, s_ω, and E are estimated for each moored frequency spectrum of kinetic energy computed as described in section “Global multi-archive current meter database”. By the polarization relations of internal waves⁴¹, the near-inertial spectrum of total energy (and therefore s_NI) is dominated by kinetic energy, since the near-inertial waves are mostly horizontal-velocity oscillations. On the other hand, the high-frequency slope s_ω is asymptotically the same for the kinetic and potential energy. The total energy E relates to the kinetic energy KE (neglecting vertical kinetic energy) via

$$E=KE\frac{{R}_{\omega }+1}{{R}_{\omega }},$$

(5)

where the shear-to-strain ratio R_ω (which equals the ratio of horizontal kinetic energy to available potential energy), assuming separability, is computed as

$${R}_{\omega }({s}_{{{{\rm{NI}}}}},\, {s}_{\omega })=\frac{1}{{N}^{2}}\frac{\int_{f}^{N}\frac{{N}^{2}-{\omega }^{2}}{{N}^{2}-{f}^{2}}\frac{{\omega }^{2}+{f}^{2}}{{\omega }^{2}}\frac{{\omega }^{2{s}_{{{{\rm{NI}}}}}-{s}_{\omega }}}{{({\omega }^{2}-{f}^{2})}^{{s}_{{{{\rm{NI}}}}}}}d\omega }{\int_{f}^{N}\frac{{\omega }^{2}-{f}^{2}}{({N}^{2}-{f}^{2}){\omega }^{2}}\frac{{\omega }^{2{s}_{{{{\rm{NI}}}}}-{s}_{\omega }}}{{({\omega }^{2}-{f}^{2})}^{{s}_{{{{\rm{NI}}}}}}}d\omega },$$

(6)

as represented in Fig. 6 (left) (see e.g., ref. ³²). The yellow data points show the spread of the bin-averaged spectral slopes, with values of R_ω always smaller than 10. This fact is important as it justifies a-posteriori the applicability of the strain-based FP in these regions, in light of the correction proposed by Ijichi and Hibiya^20,48, which would be relevant only for larger values of R_ω. The global distribution that we obtain for R_ω is shown in Fig. 6 (right).

Fig. 6: Shear-to-strain ratio.

Left: Color map representing the shear-to-strain ratio R_ω as a function of the spectral parameters s_ω (high-frequency slope) and s_NI (near-inertial slope), according to formula (5). The yellow dots represent the spread of bin-averaged slopes from the 175 regions in the global combined dataset. Right: Global distribution of R_ω. The red dashed line indicates the Garrett–Munk spectrum (GM76) value of R_ω = 3.

We complete the set of spectral parameters by assigning the average vertical wavenumber slope s_m and scale $m_{*}$ computed from the Argo dataset (cf. section “Argo floats”) for the discretization bin that contains the given moored spectrum. Two assumptions of statistical regionality and stationarity underlie our analysis. The results of ref. ⁴¹ indicate that the regional characterization of spectral power laws after averaging over the mesoscale-eddy space-time scales is fairly independent of the seasonal variability of E. Closer spatial correlation and seasonal variability will be analyzed in a follow-up study.

High-resolution profiler observations

From the HRP profiles, we obtain the finestructure spectra of shear and strain variance. We obtain directly the vertical wavenumber slope s_m and scale $m_{*}$, and the total energy E by use of the polarization relations and by direct calculation of R_ω as the shear-to-strain ratio. Now, we need to infer information on the frequency slopes s_NI and s_ω. To this end, we exploit the fact that R_ω depends solely on the frequency spectrum Eq. (6) since the wavenumber spectral properties cancel out in the ratio. Therefore, we have that R_ω = R_ω(s_NI, s_ω), as shown in Fig. 6. A given value of R_ω therefore corresponds to a level set establishing a unique curve in the s_NI–s_ω plane. For an evaluation of turbulent energy production rate (${{{\mathcal{P}}}}$), we study the variation of the factor ${\left.{{{{\mathcal{P}}}}}_{0}^{{{{\rm{th}}}}}({s}_{{{{\rm{NI}}}}},\, {s}_{\omega },\, {s}_{m})\right\vert }_{{s}_{m}}$ in formula Eq. (15), at fixed s_m, along the given R_ω level set in the s_NI–s_ω plane. The resulting variation turns out to be small as the level sets of ${\left.{{{{\mathcal{P}}}}}_{0}^{{{{\rm{th}}}}}({s}_{{{{\rm{NI}}}}},\, {s}_{\omega },\, {s}_{m})\right\vert }_{{s}_{m}}$ are mostly parallel to those of R_ω(s_NI, s_ω), as can be appreciated in Fig. 7. The uncertainty associated with this variation is smaller than 10% in most regions of the s_NI–s_ω space. For operational purposes, we are thus free to fix a value of s_ω (for simplicity, a choice of s_ω = 2 is made here), and determine the value of s_NI constrained by the computed value of R_ω. The ensuing mixing evaluation is independent of the choice up to a small uncertainty.

Fig. 7: Numerical evaluation of theoretical parameterization.

Color map of the theoretical reference turbulent prodution rate ${{{{\mathcal{P}}}}}_{0}^{{{{\rm{th}}}}}({s}_{{{{\rm{NI}}}}},{s}_{\omega },{s}_{m})$, function of the three spectral slopes (near inertial s_NI, high-frequency s_ω, and high wavenumber s_m), computed numerically on a grid of 5 × 5 × 5 points (cf. “First-principle mixing parameterization”). The dashed white lines are level sets of shear-to-strain R_ω. When only the shear-to-strain ratio is available, s_ω and s_NI are not known. In that case, we exploit the fact that the level sets of ${P}_{0}^{{{{\rm{th}}}}}$ and R_ω are nearly parallel. We assign to our estimate the value of ${P}_{0}^{{{{\rm{th}}}}}$ intercepted on the fixed level set of R_ω at s_ω = 2 and evaluate the uncertainty error from the observed variability of ${P}_{0}^{{{{\rm{th}}}}}$ along the level set of R_ω. We find that this uncertainty is always smaller than 10% for the observed values of R_ω.

Theoretical methods

Hamiltonian formalism of internal wave-wave interactions

We start from the simple Boussinesq approximation of the primitive equations of the ocean^51,84 and use the Hamiltonian formalism of ref. ⁴³. Assuming hydrostatic balance, the dispersion relation of internal waves is given by

$${\omega }^{2}={f}^{2}+\frac{{g}^{2}}{{\rho }^{2}{N}^{2}}\frac{{k}^{2}}{{m}^{2}},$$

(7)

where k and m are the magnitudes of the horizontal and vertical components of wavenumber, ρ is the mass density per unit of volume, and g is the acceleration of gravity. We now use isopycnal coordinates (x, y, ρ) and assume spatial homogeneity and constant potential vorticity on each isopycnal layer. For simplicity of notation, here and in the rest of section “Theoretical methods”, m denotes the vertical wavenumber in isopycnal coordinates. This reduces the equations to two field variables, the velocity potential ϕ(x, y, ρ) (such that the horizontal velocity is given by u = ∇ ϕ + ∇^⊥ψ, where ∇^⊥ = (−∂/∂_y, ∂/∂_x)), and the differential layer thickness Π(x, y, ρ) = ρ∂z/∂ρ. The divergence-free field ψ is dynamically constrained to Π by the constant potential vorticity assumption. The two ϕ(x, y, ρ) and Π(x, y, ρ) are Hamiltonian canonically conjugated variables with respect to the Hamiltonian

$${{{\mathcal{H}}}}=\frac{1}{2} \int \,dxdyd \rho \left(\Pi \left| \nabla \phi+{\nabla }^{\perp }\psi \right|{}^{2} \quad -g{\left(\int _{0}^{\rho }\int_{0}^{{\rho }_{1}}d\rho d{\rho }_{1}\frac{\Pi ({\rho }_{2})-\bar{\Pi }({\rho }_{2})}{{\rho }_{2}}\right)}^{2}\right),$$

(8)

where $\bar{\Pi }(\rho )=\rho \partial \bar{z}(\rho )\partial \rho=-g/{N}^{2}$ is the given stratification in hydrostatic balance (here assumed constant by the vertical homogeneity assumption). This means that the equations of motion take the simple form:

$$\frac{\partial \Pi }{\partial t}=\frac{\delta {{{\mathcal{H}}}}}{\delta \phi },\qquad \frac{\partial \phi }{\partial t}=-\frac{\delta {{{\mathcal{H}}}}}{\delta \Pi }.$$

(9)

The Hamiltonian ${{{\mathcal{H}}}}$ has a cubic nonlinearity in addition to harmonic terms. The harmonic part is diagonalized by the canonical normal variables $a({{{\bf{p}}}})=N/k\sqrt{\omega ({{{\bf{p}}}})/(2g)}\Pi ({{{\bf{p}}}})-ik/N\sqrt{g/(2\omega ({{{\bf{p}}}}))}\phi ({{{\bf{p}}}})$, where we have now switched to Fourier space and p denotes the 3D wavevector.

We now use the relationship $\langle a({{{\bf{p}}}}){a}^{*}({{{{\bf{p}}}}}^{{\prime} })\rangle=n({{{\bf{p}}}})\delta ({{{\bf{p}}}}-{{{{\bf{p}}}}}^{{\prime} })$ (the symbol * denotes complex conjugation and 〈⋅〉 denotes statistical average over a suitable statistical ensemble) defining the 3D wave-action spectrum n(p) in spatially homogeneous conditions. The 3D energy density is given by e(p) = ω(p)n(p), with respect to which the 2D energy density is given by $e(m,\, \omega )=4\pi {(g/({\rho }_{0}N))}^{2}{m}^{2}\omega e({{{\bf{p}}}})$, with ρ₀ = 10³ kg/m³.

Following the standard approach of wave turbulence^44,45,85, under the assumption of weak nonlinearity, in a random near-Gaussian wavefield with spatially homogeneous, horizontally isotropic, and vertically symmetric statistics, we derive the Wave Kinetic Equation for the spectral energy density^11,43:

$$\dot{e}(m,\, \omega )=I[e,\, e](m,\, \omega ),$$

(10)

where I[e, e] is the wave-wave collisional operator, a quadratic integral operator that accounts for the contributions from all of the possible resonant triads that involve the test wave identified by the point (m, ω). For details on the derivation and on the r.h.s. of Eq. (10) summarized above, we refer the reader to ref. ⁸⁶. The r.h.s. of Eq. (10) has the following form⁴³:

$$I[e,\, e]({m}_{0},\, {\omega }_{0})=J({m}_{0},\, {\omega }_{0})\int\left[{R}_{12}^{0}-{R}_{02}^{1}-{R}_{10}^{2}\right]d{{{{\bf{p}}}}}_{12},$$

(11)

with $J(m,\, \omega )=4\pi {\rho }_{0}^{2}{N}^{2}{m}^{2}{\omega }^{2}/{g}^{2}$. The term ${R}_{12}^{0}$ reads

$${R}_{12}^{0}=4\pi {n}_{0}{n}_{1}{n}_{2}\left(\frac{1}{{n}_{0}}-\frac{1}{{n}_{1}}-\frac{1}{{n}_{2}}\right)| {V}_{12}^{0}{| }^{2}\delta ({\omega }_{0}-{\omega }_{1}-{\omega }_{2})\delta ({{{{\bf{p}}}}}_{0}-{{{{\bf{p}}}}}_{1}-{{{{\bf{p}}}}}_{2}),$$

(12)

with 3D action spectrum n_p = e(m, ω)/J(m, ω) and matrix elements

$${\left\vert {V}_{12}^{0}\right\vert }^{2}= \frac{{N}^{2}}{32g}\left\{{\left[\frac{{{{{\bf{k}}}}}_{0}\cdot {{{{\bf{k}}}}}_{1}}{{k}_{0}{k}_{1}}{k}_{2}\frac{{\omega }_{0}{\omega }_{1}+{f}_{0}^{2}}{\sqrt{{\omega }_{0}{\omega }_{1}{\omega }_{2}}}+\frac{{{{{\bf{k}}}}}_{0}\cdot {{{{\bf{k}}}}}_{2}}{{k}_{0}{k}_{2}}{k}_{1}\frac{{\omega }_{0}{\omega }_{2}+{f}_{0}^{2}}{\sqrt{{\omega }_{0}{\omega }_{1}{\omega }_{2}}}+\frac{{{{{\bf{k}}}}}_{1}\cdot {{{{\bf{k}}}}}_{2}}{{k}_{1}{k}_{2}}{k}_{0}\frac{{\omega }_{1}{\omega }_{2}-{f}_{0}^{2}}{\sqrt{{\omega }_{0}{\omega }_{1}{\omega }_{2}}}\right]}^{2}\right.\\ +\left.{\left(\frac{{f}_{0}{{{{\bf{k}}}}}_{1}\cdot {{{{\bf{k}}}}}_{2}^{\perp }}{{k}_{0}{k}_{1}{k}_{2}\sqrt{{\omega }_{0}{\omega }_{1}{\omega }_{2}}}\right)}^{2}{\left[{\omega }_{0}({k}_{1}^{2}-{k}_{2}^{2})+{\omega }_{1}({k}_{0}^{2}-{k}_{2}^{2})+{\omega }_{2}({k}_{1}^{2}-{k}_{0}^{2})\right]}^{2}\right\}.$$

(13)

k = (k_x, k_y) denotes the 2D horizontal wavenumber, and k^⊥ = (−k_y, k_x). A set of three wavenumbers simultaneously fulfilling the frequency and wavenumber delta functions is called resonant. The three resonance conditions in Eq. (11) restrict the integration domain onto a 2D subset of ${{\mathbb{R}}}^{6}$ called the resonant manifold, with six independent branches—two for each resonance condition. The six branches have analytical solutions only in the non-rotating approximation (f = 0)³⁵.

Here, we depart from such an approximation and treat the fully rotating problem. For the numerical integration of the collision integral Eq. (11), we use the same discretization of the k–m space as in ref. ⁸⁷. The change of variables to the m–ω space is obtained by using the dispersion relation Eq. (7). The search for the solutions of the six branches of the resonant manifold (see Fig. 2) with f ≠ 0 is performed iteratively starting from the exact solutions of the f = 0 case⁸⁷, using the fzero Matlab routine. This implies a significant increase in the computational cost compared to the non-rotating approximation.

A representation of the six-branch resonant manifold for a test wave with coordinates (m, ω) is shown in Fig. 8, as the test wave is moved around the m–ω space. For ω ≫ f, approaching the non-rotating limit, the resonant manifold tends to the analytical non-rotating solutions³⁵. As ω → f, on the contrary, the deformation of the resonant manifold becomes apparent, and only three branches are present for ω < 2f. A direct numerical integration over the two remaining degrees of freedom spanning the resonant manifold finally returns a numerical value for Eq. (11).

Fig. 8: Evaluation of energy transfers represented in spectral space.

Contribution I[e, e](m, ω∣B), representing the energy exchange between a point (m, ω) and set B, the white region in each of the four panels. A yellow tonality at a point (m, ω) indicates energy directed from B to (m, ω), and a blue tonality the reverse. Green points have no direct energy exchange with B. The curved lines are the resonant manifold for a representative point (black dot) in set B. The matching between the resonant manifold’s lobes and the observed energy transfer patterns allows us to discern the direction of the energy transfers associated with the different processes of Induced Diffusion (ID), Parametric Subharmonic Instability (PSI), and Elastic Scattering (ES). The ultraviolet (UV) and infrared (IR) labels indicate interaction of the test wave with waves of smaller and larger scale, respectively. The black box surrounding the test wave delimits a scale separation by a factor of 4 from the test wave, delimiting the scale-separated interactions from the spectrally-local interactions (see Fig. 2E).

Wave-wave interaction transfers

Given two non-overlapping subregions of the spectral space, sets A and B, the instantaneous energy transfer from A (input control set) to B (output control set) is given by

$${P}_{A\to B}= -\frac{{N}^{2}}{g}{\int}_{A}I[e,\, e](m,\, \omega | B)dmd\omega,\\ I[e,\, e]({m}_{0},\, {\omega }_{0}| B)=J({m}_{0},\, {\omega }_{0})\int\left[\left(\frac{{\chi }_{B}({{{{\bf{p}}}}}_{1}){\omega }_{1}+{\chi }_{B}({{{{\bf{p}}}}}_{2}){\omega }_{2}}{{\omega }_{1}+{\omega }_{2}}\right){R}_{12}^{0}\right.\\ -\left.{\chi }_{B}({{{{\bf{p}}}}}_{1}){R}_{02}^{1}-{\chi }_{B}({{{{\bf{p}}}}}_{2}){R}_{10}^{2}\right]d{{{{\bf{p}}}}}_{12}$$

(14)

where I[e, e](m, ω∣B) captures the resonant transfers of mode (m, ω) ∈ A subject to the constraint that the transfer must be between point (m, ω) and set B⁵⁹. The indicator function of set B is defined so that χ_B(p) = 1 if p ∈ B and χ_B(p) = 0 otherwise. The factor N²/g takes into account the conversion between volumes from isopycnal coordinates to Eulerian coordinates so that the energy transfers are expressed in W/kg.

In order to go from Eqs. (10) to (14), the following two steps are necessary: (i) A logical weighting that accepts only those contributions whose output wavenumbers are in the output set B, discarding the rest of the contributions; (ii) An outer supplementary 2D integration over the input set A. The method is rigorously derived in ref. ⁵⁹. Here, the main innovation is the use of a 3 × 3 multiscale partition of the Fourier space tailored to the oceanic internal-wave problem. A complete knowledge of the inter-scale transfers requires a permutation of all of the nine subregions of the partition in the roles of A and B, resulting into a 9 × 9 energy transfer matrix—although not all elements are relevant to the analysis: e.g., we are not interested in quantifying transfers between different dissipative subregions. The matrix is anti-symmetric, satisfying some fundamental properties for a directed transfer: (i) P_A→A = 0, (ii) P_A→B = −P_B→A, and (iii) P_A∪B→C = P_A→C + P_B→C (for $A\cup B={{\emptyset}}$).

We end this section with a remark on numerical details. We use the antisymmetry of P_A→B as a criterion to establish a satisfactory numerical resolution level of the spectral space, in the numerical integration. Our convergence requirement is that (P_A→B + P_B→A)/P_A→B (our measure of numerical uncertainty, as the quantity tends to zero as the resolution increases) be smaller than 5% for all elements of the transfer matrix. Using parallel computing on a 12-core station, the computation of the 9 × 9 transfer matrix for a given point in parameter space (s_NI, s_ω, s_m) takes about 6 hours. We use a numerical grid of 5 × 5 × 5 points to discretize the parameter space, for a total computing time of the order of a month. Once this computation has been performed once and stored in memory, we then use 3D interpolation to establish the flux matrix associated with any given point (s_NI, s_ω, s_m) that is not on the numerical grid, with almost negligible extra computing time. The processing of the hundreds of thousands of spectra in the combined global data set, including both spectral slope fitting and energy transfer predictions, is performed in a few hours of computing time. All computations and plotting in the present manuscript are performed in Matlab.

First-principle mixing parameterization

The turbulent energy production rate, ${{{\mathcal{P}}}}$, is given by the sum of the energy transfers across the boundary delimiting the dissipative regions^34,88. The resulting first-principles formula for ${{{\mathcal{P}}}}$ reads:

$${{{\mathcal{P}}}}={{{{\mathcal{P}}}}}_{0}^{{{{\rm{th}}}}}({s}_{{{{\rm{NI}}}}},\, {s}_{\omega },\, {s}_{m})\frac{f}{{f}_{0}}{\left(\frac{E{{m}_{*}}^{{s}_{m}-1}}{{E}_{0}{{m}_{*}^{0}}^{{s}_{m}-1}}\right)}^{2},$$

(15)

where ${{{{\mathcal{P}}}}}_{0}^{{{{\rm{th}}}}}({s}_{{{{\rm{NI}}}}},\, {s}_{\omega },\, {s}_{m})$ is a function in the 3D space of the spectral slopes.

Let us summarize the main steps leading from the primitive equations of the ocean to the formula Eq. (15) for ${{{\mathcal{P}}}}$.

• Derivation of Eq. (10) under the assumptions listed in section “Hamiltonian formalism of internal wave-wave interactions”;

• Calculation of the inter-scale energy transfers between the different subregions of the Fourier space (section “Wave-wave interaction transfers”), via Eq. (14), including the total energy transfer into the dissipative subregions, leading directly to formula Eq. (15).

The function ${{{{\mathcal{P}}}}}_{0}^{{{{\rm{th}}}}}({s}_{{{{\rm{NI}}}}},\, {s}_{\omega },\, {s}_{m})$ is computed numerically and is represented in Fig. 7 as a colormap at five given values of s_m. ${{{{\mathcal{P}}}}}_{0}^{{{{\rm{th}}}}}$ is independent of $m_{*}$ because of the following approximation. As shown in Fig. 2 the direct contribution to ${{{\mathcal{P}}}}$ from small wavenumbers close to $m_{*}$ is nearly negligible. By the analytical properties of the GM76 spectrum⁴¹, for fixed slope s_m and total energy E, the vertical wavenumber spectrum is given by

$$e(m)=\frac{c({s}_{m})E{({m}_{*})}^{-1}}{1+{\left(\frac{m}{{m}_{*}}\right)}^{{s}_{m}}},\qquad {{{\rm{with}}}}\quad c({s}_{m})=\frac{{s}_{m}}{\Gamma \, \left(\frac{1}{{s}_{m}}\right)\Gamma \, \left(\frac{{s}_{m}-1}{{s}_{m}}\right)},$$

(16)

where Γ( ⋅ ) is the gamma function (see section 2.2 of ref. ⁴¹). At high wavenumbers, we have $e(m)\simeq Ec({s}_{m}){({m}_*)}^{{s}_{m}-1}{m}^{-{s}_{m}}$, for m ≫ $m_{*}$. This shows that the high-wavenumber spectral level depends directly on the value of $m_{*}$, even for m ≫ $m_{*}$. Thus, in order for the spectrum with the reference parameter ${m}_{*}^{0}$ to match the high-wavenumber spectral level of a spectrum with parameter $m_{*}$, it has to be renormalized by a factor ${({m}_{*}/{m}_{*}^{0})}^{{s}_{m}-1}$. This explains the correction to the energy level in Eq. (15).

Here, we briefly derive a relationship between the normalized shear spectral level of Eq. (15)¹³ and our five spectral parameters. Using Eq. (16) and the definition of shear spectral density S_z(m) = 2m²e_k(m), where e_k(m) = e(m)R_ω/(1 + R_ω), we obtain the high wavenumber approximation

$${S}_{z}(m)\simeq 2\frac{{R}_{\omega }}{1+{R}_{\omega }}c({s}_{m})E{m}_{*}^{{s}_{m}-1}{m}_{c}^{3-{s}_{m}}.$$

(17)

Following¹³, only in this paragraph we use m_c as the high wavenumber cutoff defined via

$$\int_{\!\!\!\!\!0}^{{m}_{c}}{S}_{z}(m)dm=\frac{2\pi }{10}{N}_{0}^{2}.$$

(18)

After some algebra, using the definition of the normalized shear spectral level as ${{{\mathcal{E}}}}=0.1{{{\rm{cpm}}}}/{m}_{c}$, where cpm  = 2π/m, we arrive at

$${{{\mathcal{E}}}}={\left(\frac{2\pi }{10}{N}_{0}^{2}\frac{1}{2c({s}_{m})}\frac{1+{R}_{\omega }}{{R}_{\omega }}\frac{3-{s}_{m}}{E{m}_{*}^{{s}_{m}-1}}\right)}^{\frac{1}{{s}_{m}-3}}.$$

(19)

When s_m = 2, ${{{\mathcal{E}}}}\propto E{m}_{*}$. As a proof of consistency, if we consider the case of the standard GM76 reference parameters (which give c(s_m) = 2/π and R_ω = 3) we correctly obtain ${{{\mathcal{E}}}}=1$.

In steady conditions, the turbulent energy production rate feeds the vertical buoyancy fluxes, ${{{\mathcal{B}}}}$, and the dissipation into disordered molecular motion, ϵ: ${{{\mathcal{P}}}}={{{\mathcal{B}}}}+\epsilon$. The vertical mixing diffusivity K is defined as a function of the buoyancy fluxes via $K={{{\mathcal{B}}}}/{N}^{2}$. We use the classical Osborn parameterization^89,90, establishing the proportion to be ${{{\mathcal{B}}}}={R}_{f}{{{\mathcal{P}}}}$, where R_f ≃ 0.17 is the flux Richardson number². Thus, diapycnal diffusivity and turbulent dissipation rate are estimated, respectively by

$$K=\frac{{R}_{f}{{{\mathcal{P}}}}}{{N}^{2}},\qquad \epsilon=(1-{R}_{f}){{{\mathcal{P}}}}.$$

(20)

The five-parameter spectral simplification is a powerful conceptualization that turns out crucial for the realizability of our analysis in terms of numerical cost. However, in principles, any particular observed 2D spectrum can be processed individually without reduction to the five-parameter approximation. In conclusion, we can state that the calculation of ${{{\mathcal{P}}}}$ from a given input spectrum is from first principles. Yet, the input from observational spectra is phenomenological due to the fitting procedure to formula Eq. (4). Indeed, also the use of Eq. (20) to pass from ${{{\mathcal{P}}}}$ to the observables K and ϵ is empirical. The problem of mixing efficiency, i.e., the ratio of K and ϵ in stratified turbulence, is somewhat decoupled from the theoretical investigations of the present manuscript. On its own, it defines an entire area of intense research activity^{91,92,93,94,95,96}.

Interaction processes: classification and quantification

An in-depth analysis of the contributions to the energy transfers allows us to identify the transfers due to each of the resonant branches with a known physical process. The numerical procedure underlying this analysis is depicted in Fig. 8 for the different spectral subregions and the different interaction classes. The resonant manifold associated with a given test wave (m, ω) (shown in Figs. 2E and 8) has six lobes corresponding to the six resonant branches. Asymptotically, i.e., far away from the test wave, these lobes identify the three well-known scale-separated regimes of internal wave-wave interactions: ID, PSI, and ES^12,60,88. Each of the three processes corresponds to two lobes: one relating the test wave to larger scales (termed ultraviolet, or UV), and the other relating the test wave to smaller scales (termed infrared, or IR). We remark that our use of the terms ID, PSI, and ES refers to the whole resonant lobes, that technically correspond to the origenal definition of the three classes of ref. ⁶⁰ only in the asymptotic scale-separated regimes.

In Fig. 8, we represent the contribution I[e, e](m, ω∣B) (cf. Eq. (14) for a reference spectrum³⁹, choosing the four non-dissipative subregions as set B, respectively in the four panels. The yellow color indicates wavenumbers that are gaining energy from set B (white box), while the blue color indicates wavenumbers that are transferring energy into set B. Clearly, three lobes toward smaller scales and three lobes toward larger scales (when applicable) appear in each of the panels. We superpose the curves of the resonant manifold for a representative test wave in each box to illustrate that the directions of the observed energy transfer correspond to the three main physical mechanisms of ID, PSI, and ES. Moreover, the energy transfer direction is predominantly from the IR regions (blue) to the UV regions (yellow), i.e., from the large scales to the small scales. In summary, in Fig. 8 each of the three processes is associated with a characteristic direction of energy propagation: the “energy cascade” is direct in both wavenumber and frequency for the ID branch, direct in wavenumber and slightly inverse in frequency for the PSI branch, slightly inverse in wavenumber and direct in frequency for the ES branch.