Spectral energy transfers by internal gravity wave–wave interactions for given empirical energy spectra are evaluated numerically from the kinetic equation that is derived from the assumption of weak interactions. Wave spectrum parameters, such as bandwidth, spectral slope, and Coriolis frequency f, are varied, as is the spectral resolution. In agreement with previous studies, we find in all cases a forward energy cascade toward smaller vertical and horizontal wavelengths. Energy sinks due to the transfers are predominantly at frequencies between 2f and 3f. While the mechanism of the energy transfer differs partly from findings of previous studies, a parameterization for internal wave dissipation—which is used in the fine structure parameterization to estimate dissipation and mixing rates from observations—agrees well with the numerical evaluation of the energy transfers. We also find a dependency of the energy transfers on the spectral slope, offering the possibility to decrease the bias of the fine structure parameterization by improving the knowledge about the spatial variations of this (and other) spectral parameter.
The energy spectra of superinertial fluctuations in the ocean show over wide regions a similar spectral shape that is referred to as the Garrett–Munk (GM) spectrum (e.g., Garrett and Munk 1975; Cairns and Williams 1976; Munk 1981). The spectral slopes in vertical and horizontal wavenumber and frequency tend to be close to 2 for large wavenumbers and frequencies, but can also vary regionally (Polzin and Lvov 2011). It can also be shown that those superinertial fluctuations seen in the ocean are consistent with the polarization relation of internal gravity waves (e.g., Müller et al. 1978). The universal GM spectrum is thought to be generated by the energy transfer in wavenumber space because of nonlinear interactions of gravity waves: If some forcing process such as tidal flow over topography generates waves at a certain wavenumber and frequency, nonlinearities can cause an energy transfer to waves of different wavenumbers and frequencies, such that in steady state—when dissipation of the waves balances the forcing—all possible wavenumbers and frequencies of internal gravity waves are populated with energy at a different level and a certain spectral shape results, as seen in the observations.
Energy transfers by nonlinear wave–wave interactions can be very complicated, and in general it appears to be impossible to predict them. Under the assumption of slowly changing wave amplitudes, however, it is possible to predict the rate of change of wave energy by nonlinear wave–wave interactions for a given spectrum using the so-called scattering integral or kinetic equation (Hasselmann 1966). The kinetic equation for internal gravity waves was evaluated by, for example, Olbers (1976), McComas and Bretherton (1977), and Pomphrey et al. (1980) for different versions of the GM spectrum, under different approximations, and using different numerical methods. Assuming a steady balance between forcing and dissipation of the wave field, the energy transfers calculated by those studies would allow us to localize the regions of forcing and dissipation in wavenumber space. The results and interpretations of the different studies share some common features, but unfortunately they are also in parts contradictory, which might be due to the different GM-like spectra, the different assumptions, or the different numerical methods used in those studies. Furthermore, Lvov et al. (2012) recently put doubts on the validity of the weak-interaction assumption and thus on the studies we referred to above.
The dissipation of internal gravity waves is thought to be related to breaking of the waves by gravitational or shear instability, that is, to an energy transfer from the wave field to small-scale turbulence, which ultimately leads to density mixing and an increase in potential energy in the stratified interior of the ocean. Such a potential energy increase can then drive the large-scale circulation of the ocean with potentially important effects on Earth’s climate. In particular, the global meridional overturning of the ocean is driven—in addition to the winds over the Southern Ocean—to a substantial part by this interior mixing (e.g., Talley et al. 2003). Knowing the forcing and dissipation of the internal wave field in the ocean, therefore, would help to understand and to predict the large-scale ocean circulation and climate. A parameterization for the dissipation of internal waves can be inferred from the earlier studies on the kinetic equation (Olbers 1976; McComas and Müller 1981), which forms the basis of the so-called fine structure parameterization (e.g., Gregg 1989; Polzin et al. 2014). This parameterization is often used to estimate dissipation and mixing rates from observed density and velocity profiles (e.g., Kunze et al. 2006; Whalen et al. 2012), and much of our knowledge about the global distribution of turbulent mixing in the ocean comes from it.
Given the discrepancies between the previous studies on the kinetic equation for internal gravity waves and the relevance for our view of global mixing, we revisit the issue in this study by numerically recalculating the kinetic equation for GM-like spectra based on a new analytical derivation. While most of the previous studies focus on the resonant, that is, most important for infinitely long time scales, interactions only, Lvov et al. (2012) discuss the effect of using also the nonresonant interactions in an approximate form. Here, we include all resonant and nonresonant interactions with greater accuracy. Further, we use a numerical method with an error in energy conservation as low as the numerical precision allows. The error level is reported to be on the order of 1% in Olbers (1976) using a rather coarse numerical method, while the other studies do not give information on that. We test the validity of the weak-interaction assumption, use different forms of GM-like spectra, and investigate the dependency of the results on spectral parameters, resolution, and domain size. We also compare the results of the evaluation of the kinetic equation with the fine structure parameterization (e.g., Gregg 1989; Polzin et al. 2014), which is used to estimate dissipation rates and mixing coefficients from observations (e.g., Kunze et al. 2006; Whalen et al. 2012). Models of the internal wave energy propagation (Olbers and Eden 2013) are also using the same parameterization for the dissipation of internal wave energy.
In section 2, we specify the (GM-like) energy spectrum and the kinetic equation we use in this study, while the numerical methods used to evaluate the kinetic equation are presented in section 3. In section 4, we present the results of the evaluation using different numerical methods, test the validity of the weak-interaction assumption, discuss the effect of different spectral resolution and the dependency on spectral parameters and Coriolis frequency on the evaluation in comparison to the fine structure parameterization of wave dissipation, and investigate the mechanism of the energy transfers. In the last section, we summarize and discuss the results of the evaluation of the kinetic equation. In the appendix, we give details on the interaction coefficients of the kinetic equation.
2. GM-like spectra and kinetic equation
There are various versions of GM-like spectra, differing in details of their analytical form; see, for example, Polzin and Lvov (2011) for a discussion. Here, we use the general wave spectrum considered also by Müller et al. (1978) and vary parameters in the expression. The energy density of the internal gravity wave field in the ocean is expressed as,
are dimensionless functions with the normalizations and (the latter for ) given by and . The parameter is given by with , and with the modal bandwidth in the range 3–15. While the total energy level may vary in physical space, the spectral shape is fixed by the functions A and B and the parameters (or ), s, and r. The spectral shape depends on the horizontal wavenumber magnitude k and the frequency ω but not on the horizontal wavenumber vector angle ϕ and follows for large k and ω a and power law, respectively. Note that GM76 (Cairns and Williams 1976; Munk 1981) uses with , while GM75 (Garrett and Munk 1975) uses and with .
Since it is more convenient for the numerical evaluation of the nonlinear wave–wave interaction to work in (k, m) space, where k denotes the horizontal wavenumber vector and m the vertical wavenumber, we transform the energy density from to space, which yields,
with the squared internal wave frequency . If not otherwise noted, we set in this study the parameters , , , and (corresponding to a vertical wavelength for small ω and a bandwidth of for a depth of 3000 m) and , that is, GM76.
The spectral energy density given by Eq. (1) or (3) can be used to evaluate the kinetic equation, which describes the energy transfers due to weak nonlinear wave–wave interactions in the internal gravity wave field. A kinetic equation can be derived for many dynamical systems. Often, a rather complicated Lagrangian or Hamiltonian description of the equations of motions forms the starting point as in, for example, Hasselmann (1966), but the normal Eulerian form of the equations can also be used. A Fourier ansatz allows us to identify the linear-wave solutions of the system under consideration—in our case, the Navier–Stokes equations for an incompressible fluid under rotation featuring internal gravity waves—but the ansatz also generates from the nonlinear terms of the equations of motion an integral over wavenumbers, coupling the rate of change of a particular wave amplitude to the amplitudes of other waves. To simplify and to better understand these complicated wave–wave interactions, the weak-interaction assumption (Hasselmann 1966; Nazarenko 2011) can be invoked. It implies that the amplitudes of the waves are only weakly changing by the nonlinearities, a state that is also sometimes called wave turbulence. It is then possible to transform the integral originating from the nonlinear terms into the so-called scattering integral or kinetic equation for the (statistical) second moment of the amplitudes, which is proportional to the wave energy. As seen below, the kinetic equation still involves integrals in wavenumber, but allows us to predict the energy transfers by the nonlinear terms in the equations of motions for a given energy spectrum—in our case, Eq. (1) or (3).
The kinetic equation considered here is given by
with and the notation , for the energy spectrum and (positive) internal gravity wave frequencies. The interaction coefficients C we use here and the time response function can be found in the appendix. The function depends on the argument ω and time (and is complex as the coefficients C) but for time it follows that . The interactions for which the argument ω of the time response function Δ vanishes become, therefore, resonant and will dominate over all other nonresonant interactions with nonvanishing arguments on long time scales .
The derivation of the particular version of the kinetic equation [Eq. (4)] is discussed in detail in Eden et al. (2019). We use only the first two terms of their Eq. (25), which describe sum and difference interactions, and exclude the last term of their Eq. (25), which is never resonant and in general small in all examples we tested. There are several other versions of the kinetic equation for internal gravity waves in the literature as discussed by Lvov et al. (2012), with different interaction coefficients and using different approximations. The coefficients C given in the appendix can be related to the ones given in Olbers (1974, 1976) and Müller and Olbers (1975a) for the resonant interactions but differ in general for the nonresonant interactions.
3. Methods to evaluate the kinetic equation
We apply here the method by Eden et al. (2019) to calculate all resonant and nonresonant interactions from Eq. (4) (method 1), a simplified version of that method (method 2), and a method to calculate only the resonant interactions (method 3). For the first two methods, it is necessary to specify a time scale Δt for which the time response function Δ in the kinetic equation is evaluated. Here, we use Δt = 10 days if not noted otherwise. The effect of this choice is discussed below.
a. Method 1
For method 1, the kinetic equation is evaluated numerically as in Eden et al. (2019) on an equidistant grid in K space with an uneven number of grid points in all (positive and negative) directions including the point . It is found that the wavenumber grid needs to be of this shape to allow for exactly matching triads, for which all three wavenumbers can be exchanged. When that exchange was not possible, we were not able to exactly conserve energy in the calculation of the kinetic equation. Energy conservation is checked by evaluating , and it is given by zero for the first method within the numerical precision error.
b. Method 2
Since the GM-like spectrum [Eq. (1)] is isotropic in horizontal wavenumber k, the resulting energy transfer from Eq. (4) shares the same property, that is, and where ϕ denotes the horizontal wavenumber vector angle. It is in principle therefore sufficient to evaluate Eq. (4) for a certain direction in only, say for . The energy transfer should then conserve energy by itself, that is, it holds that,
On the discrete grid in space, this property is not exactly satisfied anymore. Evaluating Eq. (4) for a certain direction in k only generates a small relative error of a few per mill in the total energy conservation, but the numerical calculation is of course much less cost intensive. We will use both approaches here, the exact three-dimensional calculation (method 1) and the approximate one for a certain direction in k only (method 2).
c. Method 3
The delta function in wavenumber in Eq. (4) can be used to resolve the integral over , but still, a three-dimensional integral needs to be solved numerically for each point in K space [or for the points of a certain slice of it, for example, ], which is numerically cost intensive. Concentrating on the resonant interactions only, that is, in Eq. (4), a further simplification of the calculation of the kinetic equation becomes possible. The kinetic equation for time simplifies to
with wave action spectrum and the scattering cross section,
For given , and , the condition yields an implicit condition on , that is, for the sum interaction [first term in Eq. (6)], and similar for the difference interaction [second term in Eq. (6)]. All integrals in showing up in Eq. (6) can then be resolved according to
with , which we find numerically from the condition with a standard root-finding algorithm, which converges after a few iterations; g′ is identical for sum and difference interactions. We find that energy conservation using this method is also only approximate, and we use it (for a slice in K space) to reduce the computational costs further as a third method (method 3). Note that method 3 is similar to the evaluation method of Eq. (6) described by Olbers (1974, 1976).
4. Numerical evaluation of the energy transfers
Figure 1 shows for the spectrum given by Eq. (1) and using method 1. The domain for the calculation has an extent of 25 km × 25 km in the horizontal and 2 km in the vertical. The grid resolution is about 100 m in the horizontal and about 8 m in the vertical, with the corresponding constant grid spacing in wavenumber space. For method 1, all interactions in all positive and negative directions of the wavenumber vector are calculated. This requires us to compute about 1014 triad interactions for this grid, which takes about 2 days of computing (real) time using 100 nodes with 36 CPUs each of the supercomputing facilities of the DKRZ in Hamburg, Germany.
In Fig. 1, the energy transfer is shown integrated over wavenumber angle ϕ as the so-called content spectrum, that is, , which we choose to better visualize the content of energy transfer in logarithmic wavenumber space: The volume under the surface shown in the figures is directly proportional to the spectral content of the corresponding function. The figure shows only for positive m since it is vertically symmetric. There is energy loss predominantly for waves with ω between exactly 2f and about 3f over a wide range of k and m, while there is energy gain predominantly at at larger k and m, but also for .
Figure 2 shows the energy transfer on a larger grid but using method 2. The domain for the calculation also has a larger lateral extent; the grid resolution is about 90 m in the horizontal and about 2 m in the vertical. For waves with , the energy loss extends now over the whole wavenumber space. The energy gain is predominantly at large k and m, but it is difficult to see if the energy gain takes place predominantly for or , since is in general even noisier than for method 1.
Figure 3 shows the effect of using different times to evaluate the function in the kinetic equation. For the function converges to , and only resonant interactions are possible, whereas for finite times also, nonresonant interactions contribute to . Truly resonant interactions, for which holds, are hardly possible on the discrete three-dimensional wavenumber grid, which we use to calculate the kinetic equation for methods 1 and 2. Correspondingly, the effect of using larger is that results of the calculation become noisier, while for smaller the results become smoother. We expect that the noise in the calculation for large would vanish if we used a better resolution in wavenumber space, since then more triad interactions on the discrete grid would approach the resonant condition , but we have not tested this assumption. However, the principal character of the evaluation does not change for large even on the grid shown in Fig. 1, and it is also similar to the evaluation for , as shown next. Therefore, we can conclude that the choice of finite or infinite plays no large role for the evaluation.
Figure 4 shows the energy transfer using method 3, where only the resonant interactions for infinite times are calculated, for a grid with about 40-m horizontal resolution and 2-m vertical resolution. The relative error in energy conservation is 0.6%. The energy transfer is now much smoother than when including also nonresonant interactions for finite times. For , there is another region with , and for even larger ω, large positive can be observed, as already indicated in Figs. 1 and 2. However, in this case, turnover time scales of the interaction are becoming shorter than the wave periods—as demonstrated in the next section—such that the weak-interaction assumption breaks down, and the results should not be interpreted literally.
We tested the impact of different grid sizes. The large energy transfers for tend to show up with high horizontal and vertical resolution, while the transfers for smaller ω remain insensitive to the grid resolution and domain size. We believe that the large transfers for are caused by a finite box effect: The spectrum and thus the energy transfers are horizontally isotropic, but the rectangular grid that we use is anisotropic. Therefore, certain triads with members in the corners of the rectangular grid in phase space cannot be realized in all lateral directions. The effect is then what can be seen in Fig. 4, which is reminiscent of wiping out a ramp feature in energy density for large wavenumbers, as seen in, for example, McComas and Bretherton (1977). A solution to the problem would be to use an infinite grid in the horizontal, which seems impossible, or a finite but horizontally isotropic grid, that is, a cylinder in physical space. But to treat this arrangement, one would need to use Hankel transforms in the horizontal instead of simple Fourier transforms, and this is beyond our scope. We therefore believe that the specific pattern for is artificial, but we cannot answer the question of what would happen with a better-suited grid.
a. The normalized Boltzmann rate
The underlying assumption of the kinetic equation is that wave amplitudes are only slowly changing, and it is necessary to test if this is really the case for the wave spectra of Eq. (1). This is in particular relevant if we consider only the resonant interactions that dominate after long times. The Boltzmann rate is given by and represents an inverse turnover time scale of the energy transfer due to wave–wave interactions. Following Nazarenko (2011) and Lvov et al. (2012), λ is normalized by the wave period, that is, is the normalized Boltzmann rate. If , the turnover time scale of the wave–wave interaction is large compared to the wave period and the weak-interaction assumption is valid, but if , the weak-interaction assumption is violated and the validity of the kinetic equation is in question.
Figure 5 shows for the energy transfers in Figs. 2 and 4. For both cases, the normalized Boltzmann rate is indeed very small for most regions in wavenumber space. Only for large m and k and for do the rates become large. This implies that the weak-interaction assumption is valid over most parts of the wavenumber domain shown here, except for large m, k, and large ω. The Boltzmann rate is also similar to the one shown in Lvov et al. (2012, their Fig. 5). It depends on the total wave energy level as . We have used here the canonical value of suggested by the earlier studies on the internal gravity wave spectrum in the ocean cited above. However, comparing this value with the global estimates by Pollmann et al. (2017) of wave energy from Argo floats shows that this value is rather high, which would indicate lower values of in most regions of the ocean. On the other hand, even if and the formal assumption of Hasselmann’s weak-interaction theory breaks down, that does not mean that there are no energy transfers or that the results of the evaluation are completely worthless, as they still may indicate qualitatively the nature of the wave–wave interactions. Direct numerical simulations (with their own problems) might prove useful to study the nature of the interactions for with strongly nonlinear behavior.
b. Energy transfer in vertical wavenumber
Figure 6a shows integrated over all horizontal wavenumbers as a function of vertical wavenumber m using different grids and method 3. Regions for which the normalized Boltzmann rate is becoming too high are excluded, that is, is set to zero for the integral over k, with the rationale explained above. Including leads to an increase of the magnitude of the minimum in all curves (not shown) but has no effect for larger m.
For all estimates of energy transfer in this study, it holds that there is energy loss for vertical wavelengths larger than about 50–70 m and energy gain for smaller wavelengths, indicative of a downward energy cascade. The same is true for the energy loss and gain as functions of horizontal wavenumber (not shown) but with a separation wavelength of about 600 m. However, the separation wavelength and the region of maximum energy gain depends on the grid resolution. The higher the resolution, the smaller the separation wavelength. This is an indication that the integral in Eq. (6) may not converge to a finite value for , since the maximum energy gain also increases with increasing grid resolution. This effect seems to be less prominent in the GM75 spectrum (Fig. 6b), for which the spectral slopes of in k and m are slightly steeper (−2.5 instead of −2 for GM76). Here, the separation wavelength also decreases with increasing resolution, but the maximum energy gain decreases with increasing grid resolution.
c. Parameter dependency of the energy transfers
We can compare the negative energy transfer and its dependency on parameters with the fine structure parameterization (e.g., Gregg 1989; Polzin et al. 2014). Olbers (1976) and McComas and Müller (1981) suggest the form for the dissipation of internal wave energy, which is also found by Henyey et al. (1986; but with slightly modified f dependency relevant only for small ). This expression for the dissipation of internal wave energy is used in the fine structure parameterization to obtain an observational estimate of the dissipation of small-scale turbulent kinetic energy (TKE) and turbulent mixing rates, assuming a constant ratio between production of TKE by internal wave dissipation and decrease of TKE because of the transfer to potential energy (Osborn 1980). Figures 7a and 7b compare D (using ) with integrated over the frequency range for different f and , showing indeed a surprisingly good agreement. Similar results are obtained using other regions of the (k, m) plane or for , but in these cases, different values for are needed. Olbers (1976) has shown that the above scaling is retrieved by analytical means.
Henyey et al. (1986) and McComas and Müller (1981) use the GM76 spectrum with spectral slopes, while Olbers (1976) uses the GM75 spectrum with a different slope. Here, we see that varying the spectral slope also generates different dissipation rates as shown in Fig. 7c. The extended parameterized form,
fits well using . The range of slopes shown in Fig. 7c roughly corresponds to what is reported by Polzin and Lvov (2011) for different regions of the ocean. Using a fixed value for the slope r thus yields a corresponding error (factor of 4–5) in existing fine structure estimates of the dissipation rate. The error could in principle be reduced by applying Eq. (9) instead.
Besides the amplitude, the different spectral slopes do not change the character of the energy transfer drastically, since there is still negative predominantly for and energy gain at large k and m. The main effect is that the region of energy loss in shifts to smaller m for larger r. The same effect as varying r can be seen by varying , that is, the region of energy loss in shifts to smaller m for smaller . For different f, the vertical separation wavelength does not change.
d. Mechanism of the energy transfers
The first term in Eq. (4) in the integral is called the sum interaction, and the second term is called the difference interaction. Both are shown in Fig. 8, demonstrating that the sum interaction is the dominant one for (as long as ω0 does not become too large), while the difference interaction is the only active one for . Similar results are seen varying the parameters as in the previous section. A resonant sum interaction is only possible for , which explains the sharp change in at 2f. Apparently, the sum interaction starts here to dominate the difference interaction, a feature usually attributed to parametric subharmonic instability (PSI; McComas and Bretherton 1977). PSI is an instability of a wave of frequency ω0, which decays in a weakly nonlinear system at the expense of waves with infinitesimally small amplitudes and half the frequency, that is, . The situation for a spectrum where all frequencies are populated with finite energy is not an instability but rather a specific type of triad interaction among other resonant triads in the spectral domain. At , however, the only triad partners have frequencies of . At higher frequencies , more triad combinations become possible, always including the PSI interaction . The fact that PSI interactions dominate at and a little above in GM-type spectra is caused by the shape of the frequency distribution of such spectra: most of the energy is located at near-inertial frequencies, and these waves are the triad partners for PSI at . At higher frequencies, the PSI interaction competes with difference interaction and is therefore no longer dominant.
PSI is identified by McComas and Müller (1981) as the dominant energy transfer mechanism toward smaller m (besides induced diffusion at higher frequencies). We approximately isolate the interactions by PSI by integrating in Eq. (6) only the negative sum interactions for smaller than 10% of and only the positive difference interactions for smaller than 10% of . Figure 8 shows this estimate and demonstrates that PSI indeed dominates for ω0 close to 2f and leads to an energy transfer toward small m. Using less (more) than 10% as a threshold for the estimate sharpens (smoothens) the maximum of at . However, PSI accounts only for a small fraction of the total transfer, and for , other, more general triad interactions seem to dominate the energy loss in the frequency band . On the other hand, we expect that a more energetic near-inertial peak in the spectrum would increase PSI and thus the increase the transfers out of 2f–3f to frequencies below 2f. However, we do not expect that the principal pattern of the transfer (transfer out of 2f–3f to lower and higher frequencies) would change in a significant way.
5. Summary and discussion
Observed spectra of internal gravity waves in the ocean show over wide regions a similar spectral shape that is called the Garrett–Munk (GM) spectrum. We have evaluated numerically the kinetic equation for the (weak) interaction of internal gravity waves in GM-like spectra using three different methods. The first method considers all resonant and nonresonant interactions in three-dimensional wavenumber space with energy conservation at numerical precision, while the second (numerically cheaper) method exploits the horizontal isotropy of the observed wave spectrum but has a larger error in energy conservation. The third (even cheaper) method considers only resonant interactions among the waves that dominate over the nonresonant ones on long time scales. All methods show an energy transfer from regions in wavenumber space with to lower and higher frequencies and to higher vertical wavenumbers, for a wide range of domain sizes and parameters of the GM-like spectra. Since the principal character of the evaluations does not depend on the choice of time , we can conclude that nonresonant interactions that are active only for finite are of less importance.
In general, the time scales for the wave–wave interactions stay large compared to the wave periods such that the underlying weak-interaction assumption used to derive the kinetic equation appears valid for most interactions. This holds both for resonant and nonresonant interactions. However, for and for very large m, the interactions become too fast and the assumption is violated. A similar behavior is reported by Lvov et al. (2012) using approximate nonresonant interactions, but here we argue that the interactions in the wave field leading to the most important energy transfers can be well described by the weak-interaction theory. To our present knowledge, the effect of too-fast interaction time scales at large ω and m can only be assessed by direct numerical simulations. However, such simulations and comparisons to the kinetic equation have their own limitations; examples are numerical grid dispersion errors and the unavoidable numerical damping. We aim to discuss this in a follow-up study, and first results suggest that there is indeed no significant effect of the too-fast interaction time scales; that is, the energy transfers in the direct numerical simulations agree well with the predictions from the evaluation of the kinetic equation and remain small for large ω and m.
One advantage of our method compared to the previous ones by Olbers (1976), McComas (1977), and Pomphrey et al. (1980) is that we do not focus only on the resonant interactions for , but also assess the effects of nonresonant interaction at finite time scales . In Lvov et al. (2012), the effects have been only approximately incorporated. We find for smooth, GM-like spectra given by Eq. (1) that a finite only results in a noisier transfer spectrum that provides no additional insight relative to the pure resonant limit. For a spectrum with a more complex structure, that is, with asymmetry in m, anisotropy in k, and narrowband peaks at tidal or other frequencies, this might be different and demands further work.
In all cases we considered, we find a forward energy cascade toward large m and k within the wave field with energy gain for vertical wavelengths smaller than 50–100 m and loss for larger wavelengths. For the spectral slope r = 2 (as in the GM76 spectrum) the maximum energy gain increases with resolution pointing toward a possible divergence of the integrals in the kinetic equation for infinitively large m, which is not seen for slope r = 2.5 (as in the GM75 spectrum). However, this issue might be more of academic nature, since variability on wavelengths smaller than a meter would be rather interpreted as turbulence and not as waves in the interior ocean. On the other hand, there might well be wave–wave interactions where one or two partners of the interacting triad are given by small-scale turbulence, a process similar to scattering at topography (e.g., Müller and Xu 1992) but beyond the scope of this study. In any case, a stationary spectrum as largely observed points toward dissipation of wave energy in regions of energy gain and forcing in regions of energy loss.
Previous studies (Olbers 1976; McComas and Müller 1981) developed a parameterization for the energy transfer in the internal wave field, which is used today to estimate dissipation rates of gravity waves, small-scale turbulence, and mixing coefficients (e.g., Gregg et al. 2003)—the so-called fine structure parameterization. We found good agreement of our calculations with the parameterization by the previous studies, which is surprising since the energy transfers also disagree to some extent with the previous studies: Parametric subharmonic instability (PSI) is an important process but not the only active one and is only responsible for a small fraction of the total energy transfer. Furthermore, there appears to be an equal share of net energy transfer from to larger and smaller ω0 and not only directed to smaller ω0 as for PSI, as argued by McComas and Müller (1981). The other process that was discussed by McComas and Müller (1981)—induced diffusion—involves interaction at large ω0 where the weak-interaction assumption is most likely violated. Despite these disagreements, we find good agreement of the calculated net energy transfer with the parameterization, which is derived for PSI and induced diffusion only.
We propose here to extend the fine structure parameterization with a functional dependency on the spectral slope parameter given by Eq. (9), which we find from an empirical fit to the results of the evaluation of the kinetic equation varying the slope parameter r. Figure 9 shows an observational estimate of the distribution of spectral slopes in the global ocean. Density profiles from the Argo dataset are used to calculate vertical wavenumber spectra from 200-m segments of strain , where Nfit is a smooth quadratic fit to the buoyancy frequency profile calculated from the Argo data and its vertical mean. Wave energy spectra are obtained from using the polarization relation of internal gravity waves. For , the dependency on m in Eq. (1) can be approximated as , and the slopes r shown in Figs. 9a and 9b are estimated from linear fits to from each spectrum.
More sophisticated spectral fits, as proposed in, for example, Müller et al. (1978) using all available data, are clearly in demand, but our first rough estimate shows already a large range of the slope parameter in the internal wave spectra. On average, the slope r decreases with depth and with latitude: in the lower depth range, that is, between 1000 and 2000 m, the variability is weakest with a mean slope of 1.60 ± 2.46, while the highest values of r as well as the strongest variability are observed in the 200–500-m-depth range with a mean slope of 2.31 ± 3.25. Also the decrease of r with latitude is especially pronounced between 200 and 500 m, where 50% (75%) of the slope estimates exceeding twice the global average for this depth range occur equatorward of 23° (37°). Among the three ocean basins, the Atlantic features the smallest slopes and weakest variabilities in any depth range, as well as the strongest concentration of high slope estimates in the (sub)tropics above 500-m depth. From these variations and the dependency in Eq. (9), changes of the dissipation rates and thus the mixing coefficients of a factor of 4–5 can be expected compared to estimates using a fixed slope parameter.
Figure 9b shows a rather complex spatial structure in the spectral slope similar to the total energy content shown in Fig. 2 of Pollmann et al. (2017). Such spatial variations in slope, total energy content, and in the bandwidth parameter m* are also discussed in Polzin and Lvov (2011), where a similar variation in the spectral law for the frequency dependency is found, which we have not varied here. We believe that more research should be devoted to observe and to understand the spatial variations of the internal wave spectrum—including the role of spectral peaks in the continuous spectrum—and its consequences for the nonlinear energy transfers implied by the kinetic equation.
The effects of forcing and dissipation and the effect of energy transports in physical space and phase space (refraction) can be added to Eq. (4) to form a complete energy balance of internal waves in the ocean. Other effects like scattering at topography, at the balanced flow, at small-scale turbulence, and at surface waves also need consideration, but methods to do so are, at hand, using weak-interaction theory; see, for example, Müller and Olbers (1975b), Müller (1977), Olbers and Herterich (1979), Müller and Xu (1992), and more recently Savva and Vanneste (2018). We propose to quantify the role of these effects to balance the nonlinear energy transfer and its spatial variations given by Eq. (4) in order to understand the internal gravity wave field in the ocean and its role for mixing.
This paper is a contribution to the Collaborative Research Centre TRR 181 “Energy Transfer in Atmosphere and Ocean” funded by the Deutsche Forschungsgemeinschaft [German Research Foundation (DFG)]—funding number 274762653. The authors are grateful to two anonymous reviewers for helpful comments. The Argo data were collected and made freely available by the International Argo Program and the national programs that contribute to it (http://www.argo.ucsd.edu; http://argo.jcommops.org). The Argo Program is part of the Global Ocean Observing System.
Interaction Coefficients of the Kinetic Equation
The interaction coefficients C are derived in the way detailed in Eden et al. (2019) from the Navier–Stokes equations for an incompressible fluid under constant rotation and are given by
with , with the left and right eigenvectors and of the linear system matrix,
with the notation , similar for P and , and with the normalization,
The eigenvectors are orthogonal to each other and normalized to one, and their components of and correspond to horizontal and vertical velocity and buoyancy. The interaction coefficients are symmetric in the last two arguments and .
The vector denotes counterclockwise rotation of the vector k by ; that is, for . The time response function is given by
The interactions become resonant for .