Barotropic tidal oscillations over seafloor topography generate baroclinic tides that may be damped in turn via nonlinear triad interactions with internal gravity waves, fueling the ambient wave field. We derive the kinetic equations for this tidal damping and the energy transfer to the ambient wave field and compute damping times and energy transfer rates for the M2 tide and a Garrett–Munk-like ambient wave field. We show that parametric subharmonic instability (PSI) interactions are important, where the tide interacts resonantly with two background waves, each of half the tidal frequency. PSI is restricted to the latitude belt 28.8°N/S and yields under typical conditions damping times of about 20 days for tides with low vertical wavenumber. Damping times decrease with equivalent mode number j roughly as 1/j2. Outside the critical latitudes PSI is not possible, and damping times are from one to two orders of magnitude larger. The energy transfer to the ambient wave field is concentrated at half the tidal frequency ω at all latitudes within the critical latitude belt. Outside, the transfer is much smaller and peaks at ω + f and N. An estimation of the tidal spectral transfer on the global scale is hampered by insufficient knowledge of the baroclinic energy distribution over the vertical modes. Using results from a numerical circulation model with tidal forcing, we find an energy transfer from the tide to the ambient wave field of typically 0.3 TW, about half of what is currently proposed for the conversion of barotropic to baroclinic energy.
The dominant energy source of internal gravity waves in the ocean is attributed to the conversion of barotropic to baroclinic tidal energy by tidal flow over rough topography as described by Bell (1975) and Llewellyn Smith and Young (2002), among others. It is argued that this route provides about 1 TW for ocean mixing (Wunsch and Ferrari 2004), with up to 85% provided by the principal lunar semidiurnal component, the M2 tide (Egbert and Ray 2001). Recent estimates by Falahat et al. (2014) and Vic et al. (2019) point toward somewhat lower values, that is, 0.5 TW, for the first 10 modes of the M2 constituent. The tidal source has tidal frequency and mostly long wavelengths. To enable mixing by overturning of small-scale waves a transfer of baroclinic tidal energy to these small scales in the background wave continuum has to take place. Tidal energy must be supplied to the wave continuum and a downscale transfer in the wave continuum must follow, the mechanism for this route being scale-bridging nonlinear wave–wave interactions. Furthermore, bottom scattering (Müller and Xu 1992) and wave–eddy interactions (Savva and Vanneste 2018) have been addressed. While the barotropic-to-baroclinic conversion, basically a linear process, has recently attracted considerable observational and theoretical research, for example, Nycander (2005), the subsequent downscale transfer is much less investigated. Early work by Olbers (1976), McComas and Bretherton (1977), and Müller et al. (1986) and more recently by Polzin and Lvov (2011) and Eden et al. (2019) showed that wave–wave interactions indeed provide energy to the small-scale overturning waves on time scales of a few days to a week. A similar time scale has been found for wave–eddy scattering (Savva and Vanneste 2018) and scattering of low-mode tidal waves at the continental shelf, whereas bottom scattering in the open deep ocean seems to establish only an inefficient transfer (Kelly et al. 2013), which is in contrast to the role of bottom scattering in the model of Eden and Olbers (2014). On the observational side, using SSH satellite data, Zhao et al. (2016) provide us with a rather small residence time (total energy over total conversion flux) of baroclinic tidal energy of 1–1.5 days, corresponding to an e-folding scale of 400 km for propagation of mode 1 over the World Ocean. This is an extremely small estimate compared to the attenuation scales 750–3000 km found by Alford et al. (2019) from the same data. Nevertheless, a gap remains in this scenario right at the start of the transfer route from tidal energy to small-scale overturning: the damping of the baroclinic tide and the associated energy transfer to the wave continuum. This is the topic of the present study.
Several studies have investigated the interaction of a spectral peak with a background wave continuum (Olbers 1974; Pomphrey et al. 1980; Olbers and Pomphrey 1981) and more recent evaluations have been done by Hibiya et al. (1998), Hibiya et al. (2002), and Eden and Olbers (2014) (note that there was a sign error in their transfer rates) and Onuki and Hibiya (2018). The results of these studies seem quite diverse. Olbers (1974) uses a rather unrealistic box-shaped model spectrum (GM72) of the background internal wave field, while Pomphrey et al. (1980) use a more realistic spectral shape (GM76) but apply a quite severe cutoff with respect to high vertical wavenumbers that are allowed to interact with the tide. The calculation of Olbers and Pomphrey (1981) with this “amputated” background spectrum is for a latitude of 30° where the parametric subharmonic instability (PSI) interactions of the M2 tidal frequency are not possible, missing thus the most important sum interaction and thus all of the PSI process. Scaling of the results to lower latitudes, as discussed in Olbers and Pomphrey (1981) and in Olbers (1983), then only reflects the slow decay with about 100 days for mode 1 and 10 days for mode 10, arising from the difference interactions. Including PSI interactions, Onuki and Hibiya (2018) find a much faster decay of 20 days for mode 1. How this tidal energy spreads over the wave continuum has not yet been investigated but will be presented here.
Since baroclinic tides make up a major source of internal gravity wave energy in the ocean and consequently for mixing, the decay of the tide by PSI has become a paradigm of downscale transfer, in particular because the process was promoted by the numerical experiments of MacKinnon and Winters (2005). These showed an almost catastrophic breakdown of northward propagating baroclinic tides at the critical latitude 28.8° into a pair of near-inertial waves and smaller scale turbulence. A similar more detailed scenario was simulated by Gerkema et al. (2006). Increased dissipation was thus expected near this latitude but substantial attempts to find such signature in observational data have largely failed: though there were some indications of PSI in phase relations there was no catastrophic drain of energy for the tidal motion (Alford et al. 2007; MacKinnon et al. 2013). An exceptional case of increased diapycnal diffusivity south of the critical latitude is reported by Hibiya and Nagasawa (2004) and Hibiya et al. (2007) based on finestructure measurements and parameterization (Gregg 1989). Following MacKinnon and Winters (2005) and other numerical modeling studies the time scale of PSI processes is longer than 2–5 days. The shortest decay time of the lowest-mode internal tides, found by Onuki and Hibiya (2018), is 20 days, which makes the role of PSI rather limited even near the critical latitude. Such values agree with the analytical work of Young et al. (2008).
We start our investigation of baroclinic tidal decay like our predecessors with the kinetic equation for resonant weak wave–wave interactions in a random-phase approximation (see, e.g., Hasselmann 1967b; Olbers 1976; Nazarenko 2011; Olbers et al. 2012). The spectrum for the internal gravity wave field is assumed to consist of a tidal peak on top of a background continuum, assumed here to be of Garrett–Munk (GM) type (Garrett and Munk 1975; Munk 1981). The kinetic equation is separated into the energy balance of the tide and that of the background wave field (section 3). The tidal decay is treated in section 4 and the spreading of tidal energy over the continuum spectrum in section 5. Section 6 opens a global perspective, evaluating the decay rate of baroclinic tides for the World Ocean. The field of internal wave energy used in these calculations is taken from Argo float data, as described in Pollmann et al. (2017). The appendixes contain details of the GM spectral class, the coupling coefficients of the interaction triads and the method of integration of the kinetic equation. A summary and discussion follows in section 7. Because parametric subharmonic instability is a dominant triad interaction in the tidal-background interaction, a short summary of PSI physics is given in section 2.
2. PSI interactions
Triad interactions occur in systems with quadratic nonlinearities. A wave with frequency ω0 interacts resonantly with waves of frequencies ω1 and ω2 by sum interactions when the resonance condition ω0 = ω1 + ω2 is satisfied, and by difference interactions in case of ω0 = ω1 − ω2 (all frequencies are taken to be positive). The amplitudes aj, j = 0, 1, 2 of the wave components are governed by the canonical triad interaction equations (Hasselmann 1967a)
where η is a complex interaction coefficient arising from the quadratic terms in the original equations of motion of the system (e.g., the Boussinesq equations), and ± corresponds to sum/difference interactions. The triad problem in Eq. (1) has three independent integrals (among them the total energy ) and can be integrated in analytical form (see, e.g., Craik 1988), showing that energy is continuously exchanged between the wave components but not lost.
For the specific conditions of initially infinitesimal waves 1 and 2, however, the sum interaction shows an instability phenomenon. Energy is then “pumped” from component 0 (the pump wave) to components 1 and 2 that grow exponentially. As long as |a0| ≫ |a1|, |a2|, the pump wave amplitude can be regarded as constant and the now linear equations for a1 and a2 with time dependence ~e−iλt imply eigenvalues (Hasselmann 1967a)
The frequency factor ω1ω2 = ω1(ω0 − ω1) peaks at ω1 = ω2 = ω0/2. The instability of such triads is referred to as parametric subharmonic instability. While the energy transfer in the frequency domain is toward lower frequencies, the transfer in wavenumber space is toward higher wavenumbers (both horizontal and vertical). In particular, the vertical wavenumbers m1 and m2 must be opposing and much larger than m0, the wavenumber of the pump wave. The PSI transfer is thus toward the dissipation range.
Parametric subharmonic excitation occurs in many mechanical, electronic and optical oscillatory devices by means of a periodic variation of a parameter of the system (see, e.g., Adler and Breazeale 1971). The system can then oscillate with a frequency equal to one-half of the driving frequency. For water waves a prototype of PSI are internal gravity waves (IWs) because these are prone to nonlinear triad interactions (see, e.g., McEwan and Robinson 1975). IWs have the restriction that their frequency must lie between the local Coriolis f and the local Brunt–Väisälä frequency N. Consider a pump wave candidate with frequency ω0 (say an internal tidal wave) at a latitude φ where1f = 2Ω sinφ. PSI requires that ω0 > 2f such that ω1 and ω2 exist above f. This defines a critical latitude by ω0 = 4Ω sinφc. If the pump wave is at φ = φc where ω0 = 2f it must decay into two waves of frequency f, a genuine PSI process. If φ > φc, resonant sum interactions are impossible and the wave is stable. If φ < φc, a broad range of triad interactions is possible, always including a PSI process. For the M2 tidal constituent as the pump wave this is equatorward of 28.8° latitude, as shown in Fig. 1. Note that resonant difference interactions are not subject of such a restriction: for the M2 tide they may occur at all latitudes up to the natural bound of propagating tides at 74.5°.
The situation in a spectrum, the issue of the present study, is not that of an instability. In a spectrum usually all frequencies are populated with waves of finite energy and a PSI interaction is rather a specific type of triad interaction among others in the spectral domain. Any spectral component at frequency ω0 can interact with components at frequencies ω1 and ω2 and energies according to the sum resonance condition ω0 = ω1 + ω2. At ω0 = 2f the only triad partners for sum interaction, however, have frequencies ω1 = ω2 = f. At higher frequencies ω0 > 2f more triad combinations become possible, always including the PSI interaction ω1 = ω2 = ω0/2. PSI interactions, however, are expected to dominate the spectral rate at and a little above ω0 = 2f in GM-type spectra. This is not because of an instability situation but derives from the shape of the frequency distribution of such spectra: the IW energy is located overwhelmingly at near-inertial frequencies, and these waves are the triad partners for PSI in the vicinity of ω0 = 2f. At higher frequencies PSI interactions compete with difference interactions and it must be investigated whether PSI is dominant.
There is one further outstanding point to PSI interactions in a spectral framework. Important for the spectral interaction rate is the range of triad partners that can interact resonantly, loosely speaking the number of interacting triads for a given frequency ω0. For an IW of frequency ω0 and horizontal wavenumber k0 = |k0| the range of wavenumbers k1 of partner triads with frequency ω1 is generally finite (there is a low-wavenumber and a high-wavenumber limit) unless ω1 = ω0/2: then the range of k1 extends to infinity, that is, waves of given (ω0, k0) interact resonantly with waves of frequency ω0/2 in a broad range of scales, actually up to arbitrarily small scales. This property is documented in the present study. It singles out the PSI interaction process in the spectrum. However, it does not immediately imply that PSI should govern the spectral rate of change since there are always difference interactions and sum interactions of non-PSI quality possible.
3. Resonant interactions of a wave continuum with a monochromatic line
Our analysis is performed by representing the internal gravity wave field by three-dimensionally propagating waves, which requires a WKB condition on the typical wavelength and the vertical scale of the background Brunt–Väisälä frequency. This differs from Pomphrey et al. (1980) and Onuki and Hibiya (2018), who used vertically standing modes. Although a valid mathematical representation, we are skeptical of the standing mode description of baroclinic tides because the travel time (bottom to top of the ocean) is of similar magnitude as the nonlinear damping time of the tidal decay, making it hard to form physical modes. As shown in this study, WKB waves yield very similar decay times as the modal waves analyzed in Onuki and Hibiya (2018).
The energy of the line with wave vector2 and frequency , interacting with a spectral continuum, is governed by
where is the decay rate. A derivation of the decay rate starts with the triad scattering integral (or kinetic equation) for resonant interactions (see, e.g., Hasselmann 1967b; Olbers 1976; Nazarenko 2011; Olbers et al. 2012; Eden et al. 2019),
Here is the action spectrum of the wave field and ωj = Ω(Kj) is the (positive) frequency according to the dispersion relation of IWs. Furthermore, Tμ is the cross section given by
where is the triad interaction coefficient,3 given in appendix A. Wave triads interact by either sum (μ = +) or difference (μ = −) interactions, satisfying the resonance conditions K0 = K1 + μK2, ω0 = ω1 + μω2. We will use the fact that triad interactions conserve the total energy of the wave field (see, e.g., Olbers et al. 2012).
We insert for the action spectrum a superposition of a δ function at the wave vector and a smooth background spectrum,
with , such that
The derivation of the respective rates of change for and is straightforward. With Eq. (6) the kinetic equation [Eq. (4)] contains terms with the δ-function factor and terms without. Such terms must balance individually. The result is summarized by
The derivation is lengthy but straightforward. The decay rate is frequently called the Langevin rate (Pomphrey et al. 1980). The above balance of the background spectrum describes the triad transfer within it (by ) and the exchange interactions of the background with the line (the Λ term). The energy lost or gained by the line is exchanged with the background spectrum. Note that Eqs. (8) and (9) may as well be derived starting from the triad interaction equation [like Eq. (1)] by applying a perturbation analysis following the weak interaction theory (Hasselmann 1967b; Nazarenko 2011; Olbers et al. 2012). Since energy is conserved by triad interactions we have the constraint
which is valid independently on the energetic state of the line.
including in particular forcing and dissipation terms (possibly in the associated boundary conditions) and implementing the fluxes due to propagation and refraction. When evaluating the interaction terms in Eq. (11) for a certain background spectrum and set up a global view of the tidal decay we implicitly assume that this spectrum is close to observed one in the ocean and close to the solution given by the above balance.
4. Decay of the tidal line
In this section we will discuss the sum and difference interactions of the tidal decay, separating into λ+ and λ−. The sum interaction in Eq. (8) (λ+, first term) with is positive definite, it contributes to only for . For M2 this is south of 28.8° latitude, as depicted in Fig. 1. The difference interaction (λ−, second term) with contributes to for all frequencies and it can have both signs. For M2 it operates to 74.5° latitude, where the M2 tide ceases to exist as a freely propagating wave. Note that the difference interaction is the only one possible for , that is, north of 28.8° latitude.
In this study we first discuss results for M2, when f = 2Ω sinφ varies with the latitude φ and N is fixed at 5.2 × 10−3 s−1. In the global application in section 6 we use a globally varying Brunt–Väisälä frequency. In the following the line frequency is denoted by ω for brevity, with horizontal wavenumber k according to k = k(cosϕ, sinϕ). We use polar coordinates (k, ϕ) for the action and energy spectra and transform the dependence of the vertical wavenumber m to frequency ω, using the dispersion relation
Here σ is the sign of m. The background action spectrum is then written as .
For the difference interaction we obtain
and the contribution from sum interactions is given by
The first is the Jacobian of the transformation, Ji = J(ki, ωi), and the second expression arises from the elimination of the frequency δ-function by integration over ϕ1. Furthermore γ = ϕ − ϕ1. Note that (1/2)kk1|sinγ| is the area of the triangle in wavenumber space that is enclosed by the triplet k, k1, and k2 Prescribing k and ω, the rates λ± describe the decay of a line in the spectrum with vertical wavenumber m. We use |m| ≈ jπN/∫N dz to define an equivalent vertical mode number j of WKB waves.
The integration of λ− and even more so of λ+ is quite cumbersome. First, the domain of integration has a complicated shape. Figure 2 displays cosγ as function of ω1 and k1 for given tidal frequency ω and the particular latitude 18°. The boundary of the domain where resonant interactions can take place are given by the integration limits k1min and k1max, which are functions of ω1. They are defined by the bounding curves cosγ = ±1 in Fig. 2 and can be derived in analytical form. The integrand, which will be discussed in detail in section 5, is singular on the boundaries of the k1 range owing to the factor 1/|sinγ|. At the limiting wavenumbers k1min and k1max we have sinγ = 0, cosγ = ±1. The triad wave vectors are aligned in this case. The singularity goes as and , and we deal with this complication by a partly analytical integration of the square root behavior, as described in Olbers (1974).
We evaluate the decay rate and the spectral transfer for spectra of GM type (Garrett and Munk 1975; Munk 1981) that have a slope of −2 at high frequency and −s in vertical wavenumber (see appendix B). In particular examples two values of wavenumber slope will be considered, s = 2 and s = 3. The other parameters of the spectra, the total energy E and the wavenumber scale , are given in appendix B.
a. Difference interactions
For difference interactions the line interacts with two waves of higher frequency according to ω = ω1 − ω2 and K = K1 − K2. In Fig. 3 we display the cross section T− of this triad as function of ω1 and k1. It is strongly increasing with wavenumber k1. The integration domain for difference interactions has a cusp at ω1 = 2ω where the lower limit of resonant k1 extends to zero. Note that this is a PSI interaction with ω = ω2 = ω1/2, that is, the triad partner with frequency ω1 is here the pump wave. Integration problems are not expected there since the integration always has a finite range. The GM spectrum is isotropic and symmetric and λ− (as well as λ+) becomes independent of φ.
We show a calculation of λ− in Fig. 4 for three modes 1, 3, and 10. The background spectrum is GMCL (a GM-type spectrum without inertial cusp, see appendix B). Other spectra are discussed later. The resolution for the integration is 3000 points for k1/k and 2000 points for ω1. Increasing the latter to 5000 points yields almost identical results.4 The integration is thus stable and decay times between 1 (for high modes and low latitudes) and 100 days (for low modes) are obtained. Note that the slope s = 3 produces negative λ− at higher latitudes, indicating that there the line grows at the expense of the background wave continuum. It is found that the decay time 1/λ− is generally much larger than the tidal period . Only very high modes at low latitudes may decay within a tidal period, violating the assumption of weak interactions.
b. Sum interactions
For the sum interaction the line interacts with two waves with lower frequencies ω1, ω2 < ω according to ω = ω1 + ω2 and K = K1 + K2. Contrary to the difference interaction, which may yield negative λ−, the sum interaction always has a positive definite λ+. Note that the integration domain for λ+ goes to zero when ω approaches 2f and vanishes for ω < 2f (for M2 north at latitudes higher than 28.8°). Sum interactions of the tide with the internal wave continuum are only possible equatorward of this critical latitude. Continuum partners of the resonant triad are, however, not only near-inertial waves but all waves with frequencies ω1 above the local value of f and below ω − f. At a latitude close to the critical latitude this range is becoming increasingly narrow. A complication may occur with the frequency part of the GM spectrum when the critical latitude is approached: the triad partners then reach frequencies very close to f, where GM76 becomes singular (the integral is still finite). For this reason, to test for stability of the integration, we shall first take the cuspless GMCL for simplicity where the enhanced inertial peak is absent (see appendix B).
Everywhere in the integration domain, as for λ−, we have an integrable singularity on the rim where sinγ goes to zero. These singularities are treated as before for the difference interaction by analytical means. In Fig. 5 the cross section T+ for sum triads is displayed (for mode number j = 1 and both signs of m1). Unlike T− the cross section T+ generally increases with decreasing wavenumbers k1. Note that there is a PSI resonance at the frequency ω1 = ω/2. Here the test wave with frequency ω decays into two waves of half the frequency. The cross section T+ is generally small at high wavenumbers in this PSI tongue of the integration domain and the integrand does not gather much power either. Still, exactly at the PSI resonance ω1 = ω/2, the k1 integral extends to infinity and—depending on the spectral slope of the background spectrum—a singular behavior could happen. Figure 6 displays the k1 dependence of the integrand for a gridded range of ω1 values for resonant sum interactions (each curve is for a constant ω1), using GMCL as background spectrum with wavenumber slopes s = 2 and s = 3. The cross section T+ develops into a constant behavior at high k1 and the total integrand approaches a behavior . In Fig. 6 we first note the singularity from the factor 1/sinγ at both ends of the respective curves for constant ω1. Second consider the curve extending to the largest resolved k1: this corresponds to the PSI interaction ω1 = ω2 = ω/2. Apparently, the wavenumber spectrum with wavenumber slope s = 2 (GM76) results in an integral of and a spectrum with s = 3 in an integral of . The decay rate for GM76 thus does not lead to finite values unless we consider a high-wavenumber cutoff in which case it will depend on this cutoff wavenumber. On the other hand the spectrum with slope s = 3 is integrable at the PSI point. In the examples discussed below we use spectra with s = 2 and s = 3. To have comparable conditions both cases are treated with a cutoff at corresponding to a wave scale in the dissipation range.
The integral for sum interaction faces severe problems of convergence, originating from PSI interactions. We report on our heuristic method to overcome convergence problems in appendix C. We use the criteria for convergence, developed there, to compute the decay times of tidal modes 1–10 in the background of GMCL with s = 2 and s = 3, shown in Fig. 7. We find decay times which generally are well above the tidal period, ranging from a couple of 100 days for mode 1 to a few days for mode 10.
Figure 8 shows the decay time of the tide for the three modes 1, 3, and 10 where sum and difference interactions are combined appropriately. The left panel is for the cuspless spectrum GMCL, and the right panel shows the same but using the complete GM76 spectrum that includes the GM inertial peak, both with slope s = 2. Naturally, the enhanced inertial peak in GM76 leads to higher decay rates and smaller decay times. Note the sharp decrease of decay time when the critical latitude is approached from south and the sharp increase when this latitude is passed. For GM76 the decay time of M2 becomes as small as 10 days for mode 1 at the critical latitude of 28.8°. Higher modes have much smaller decay times.
The GM class of spectra (see appendix B) that we use for the background wave continuum has a number of free parameters which enter the tidal decay rate and likewise the transfer of energy to the background is parameter dependent. These quantities thus depend on the total energy E of the continuum, on the spectral wavenumber scale and the slope s in the vertical wavenumber space. Additional parameters are the Coriolis frequency f, which we have considered explicitly in the previous analysis, and the Brunt–Väisälä frequency N, which was held constant so far. The latitudinal presentation of the decay, displayed in the previous sections, must therefore be taken with caution since under realistic conditions all parameters are expected to vary with latitude, not only f.
To assess the parametric dependence of the decay rate we have chosen the latitude 20° for sum interactions and 40° for difference interactions and varied all parameters within a reasonable range. Note that the dependence on E is trivial as λl ~ E exactly. Results for the other parameters are displayed in Fig. 9. In view of the quite large range of parameter variation the decay rates change only to a minor degree, with exception of the decay for difference interactions, where the dependences on slope and wavenumber scale are significant.
5. Effect on the background spectrum
The energy lost by the tidal decay is spread across the background spectrum according to the rate of change specified in Eq. (9), written here as with
The integral over K1 has been evaluated accounting for the δ function of wavenumbers. As we did with the tidal decay rate we transform to (k, ω, ϕ, σ) coordinates, , where J0 = J(k0, ω0) is the Jacobian given in Eq. (15). We use the transformation for the spectra in Λ as well,
where and . For an isotropic background . In addition we have Ω(−K1) = Ω(K1) for IWs.
In the following analysis it is meaningful to rename the wave vector and indices: 0 in 1, 1 in 2. Using then the symmetry of the interaction coefficients we find
and we arrive at
We should bear in mind that the two contributions to Λ occur in disjunct domains in (ω1, k1) space, as demonstrated during the discussion of the integrands of the decay rates λ±. This property holds at each latitude and is further discussed in section 7. We split Λ and in their two contributions according to ± interactions.
In the following we consider the rate of change of energy . In view of the isotropy of the background spectrum it is meaningful to integrate the rate over the angle ϕ1 such that
The derivative dω2/dϕ1 is given by Eq. (15). The terms for sum and difference interactions are readily identified with the corresponding integrands of the line decay rates λ±.
Figure 10 shows the transfer rate in wavenumber–frequency space at a specific latitude for sum and difference interactions and both signs σ1 = ±. The background spectrum is GM76 in all cases. The figures reveal that sum interactions of the tide with (ω, k) clearly favor background waves at low ω1 ≃ f and low k1 ≪ k (near-inertial waves) and background waves at ω1 ≃ ω − f and k1 ≃ k, with mutually exchanged ω2 partners. Except for these extreme interaction foci the transfer is generally to small wavenumbers k1 < k. The PSI point ω1 = ω/2 does not show any increased magnitudes of the integrand in these plots, except that the wavenumber range for sum interaction extends to infinity. Only when integrated over wavenumbers k1 PSI emerges as exceptional case (see the following discussion). The difference interactions of the tide, on the other hand, favor background waves with frequency ω1 ≃ ω + f, hence ω2 being near inertial, and wavenumbers k1 > k. Here the rate spectrum is positive (transfer of tidal energy to the background). Negative spectral rates occur at low wavenumbers and frequencies around ω1 = 2ω.
Figure 11 shows the transfer rate , integrated over wavenumbers, for sum and difference interactions as function of latitude and frequency ω1. At each latitude the PSI triad interaction is outstanding in the sum transfer rate. The transfer for difference interactions is several orders of magnitude smaller and concentrated on the rims of the frequency domain at ω1 = ω + f and ω1 = N. Referring to the sketch (shown in Fig. 16 below) of the frequency–latitude domain in section 7, the background takes over energy from the tidal decay in the sum domain along the PSI ridge up to the critical latitude and in the difference domain all along the frequency boundaries at all latitudes.
6. A global perspective
We have evaluated the decay rates λ± for the World Ocean, using for the Brunt–Väisälä frequency N(z) the WOCE-Argo data (Gouretski 2018) and for the IW energy E estimates obtained from Argo data (Riser et al. 2016) by Pollmann et al. (2017). The energy estimates used here are for the range 250–500 m. The data, and thus also the calculated decay rates, are on a 1.5° × 1.5° grid. The slope s of the wavenumber spectrum and the wavenumber bandwidth , or, equivalently, the mode number bandwidth , as defined by Eq. (B4), is set to the conventional values s = 2 and [see Pollmann et al. (2017) for further discussion]. A map of the energy E used for the further calculations is shown in Fig. 12 (left).
Figure 13 displays the decay time 1/λj with for modes j = 1, 2, 3, and 10 from our calculations, shown as global maps. Remember that modes in our WKB frame refer to equivalent vertical modes where mode j has a vertical wavenumber mj = jπN/∫N dz. The most obvious feature in Fig. 13 is the strong overall decrease of the decay time with mode number, roughly by a factor of 10 from mode 1 to 3 and again from mode 3 to 10. The decay rate λj thus approximately increases with mode number as j2, as shown in Fig. 14. In this graph the decay times 1/λj are plotted against latitude and the spreading at each latitude comes about by the variation of the Brunt–Väisälä frequency and IW energy along the associated longitude circles. The rough dependence ~j2 agrees with the analysis of Olbers (1983) used in the recent studies of de Lavergne et al. (2019) and Vic et al. (2019). However, as explained on in the introduction, the early results of Olbers (1983) only apply to the slow decay by difference interactions. Second, the abrupt change from the latitude band −28.8° to 28.8°, where in addition to difference interaction also sum interaction and PSI take place, to higher latitudes, where only difference interaction takes place, is evident for all modes. In the zonal mean of the data shown in Fig. 14 the jump of the decay time at the critical latitude is from roughly 19 (6, 3, 0.7) days to 76 (30, 17, 3) days for mode 1 (2, 3, 10). Mode 10 decays very fast in the critical latitude band on time scales comparable to the tidal period. Therefore, the assumption of weak interaction breaks down for the high modes, and such interactions need to be discussed using other means than the kinetic equation [for a discussion of this problem, see Holloway (1980) and Müller et al. (1986)]. Outside the critical band, however, the theory holds for all modes considered.
The amount of energy gained by the wave continuum depends on the baroclinic tidal energy and to evaluate this globally we need to know the ocean’s tidal energy distribution. Considerable insight has been gained in recent years on the modal and regional distribution of baroclinic tidal energy. A major observational step is the mooring program conducted at the Hawaiian Ridge and reported by Zhao et al. (2010). Tidal energy was found predominantly in the first few modes, with 2 kJ m−2 in mode 1 at the Hawaiian Ridge, decreasing strongly with mode number (e.g., a factor of 10 less in mode 3). Such pattern is confirmed by Vic et al. (2018) (0.8 kJ m−2 in mode 1) at the Mid Atlantic Ridge and is also consistent with the global map in Zhao et al. (2016) based on satellite observations (note that satellite data only resolve the first mode). The Hawaiian Ridge and the Mid-Atlantic Ridge are areas known as energetically intensive with energies of the baroclinic tide that are much larger than found in other areas (see, e.g., Wunsch 1975; Noble 1975; Schott 1977; Hendry 1977) with values of about 102 J m−2 or less. A representative value can be obtained from Zhao et al. (2016) who reports 36 PJ as globally integrated mode-1 energy, equivalent to an area mean of 102 J m−2. The global ocean model STORMTIDE (Müller 2012; Müller et al. 2012) with a resolution of 1/10° and forced by the complete lunisolar tidal potential simulated 80 PJ of kinetic energy of M2 baroclinic tidal motion (Li et al. 2015). Our damping times for mode 1 are between 10 and 60 days (Fig. 8), the lower values for sum interactions, the higher for difference interactions. Such time scales imply an energy transfer of 0.1–0.02 mW m−2 for a tidal energy of 102 J m−2, as reported in Olbers and Pomphrey (1981; note that in that study only the larger decay times for difference interactions were considered). The energies from the tidal hot spots, referred to above, easily yield transfer values exceeding 1 mW m−2 for mode 1. This is significant as source for the wave continuum (for comparison: 1 TW input globally from tidal conversion amounts to 2.9 mW m−2 as a mean).
To estimate the net transfer of baroclinic tidal energy to the wave field,
we need to know the tidal energy Ej in mode j, as a function of position (x, y, z). At present such spatial information is not available on global scales. Using the STORMTIDE model data (Müller 2012) and repeating the analysis of Li et al. (2015, see, e.g., their Fig. 2; note that this procedure does not allow the separation between individual modes), we calculate the total resolved M2 horizontal kinetic energy and average onto the 1.5° × 1.5° grid of the Argo-based energy levels of the IW continuum and derived decay rates. The total mechanical energy [as function of position (x, y, z)] of the internal tide field is given by the sum of corresponding kinetic and potential energies and can be calculated from the kinetic energy by using the polarization relations of internal gravity waves (see, e.g., Olbers et al. 2012),
with ω being the M2 frequency. The global integral of in the STORMTIDE data of Müller (2012) and Li et al. (2015) on our grid yields 76.5 PJ whereas the integrated wave energy is 110.4 PJ. We integrate over the depth range 250–500 m and average. The result, denoted by in the following, is depicted in Fig. 12 (right). We use to represent the lateral distribution of baroclinic tidal energy.
The vertical distribution will be modeled by specifying weighting functions for the vertical modes j. We assume in Eq. (21) with wj representing the fraction of energy in the respective vertical mode j (with ), which will be taken independent of the position. With our limited calculations of the global pattern of the decay rates and the surprisingly little knowledge about the modal partitioning of tidal energy we restrict the sum in Eq. (21) to five modes and assume wj ~ j−n roughly following local observations (Schott 1977; Rainville and Pinkel 2006) and numerical models with tidal forcing (e.g., Cummins and Oey 1997). The exponent n is likely larger than 2. With the present data we are able to compute the transfer integrated (or averaged) over the depth range 250–500 m, using in Eq. (21). The vertical integral is approximated by multiplying by a vertical scale D which we take from the global integrals of baroclinic energy, from the STORMTIDE data of Müller (2012) and Li et al. (2015). The calculation yields D = 1590 m for the total depth.
In view of all these uncertainties we can only give a rough estimate of the net transfer. Table 1 summarizes the results for wj ~ j−n with j = 1, …, 5 and n = 2, 3, 4, showing only differences up to a factor of 2: the values for the total transfer (last column) range from 0.1 to 0.05 TW. As can be inferred from the table, the majority of the transfer comes from the sum interactions in the latitude belt from −28.8° to 28.8°, which are dominated by interacting PSI triads. Following this particular scenario for the baroclinic energy and its modal distribution the transfer is as global integral less than 0.1 TW.
The tidal energy in the numerical simulation of Müller (2012) and Li et al. (2015) may be underestimated. Considerably larger energies result in the analysis of de Lavergne et al. (2019) with 165 PJ (mode 1), 79 PJ (mode 2), 30 PJ (mode 3), 13 PJ (mode 4), and 6 PJ (mode 5), summing to 293 PJ. Upgrading the energy level of 110 PJ of our previous estimation to 293 PJ and taking the distribution of that analysis we arrive at the values in the last two rows (marked “lav”) of Table 1 with a net transfer of 0.31 TW.
Figure 15 displays the transfer computed for the modal spreading wj, following de Lavergne et al. (2019), and integrated over the water depth, as explained above. The transfer is largest in the western and central Pacific with maximum values around 10−2 W m−2. Also the western Indian Ocean and the eastern Atlantic occasionally show such a magnitude of transfer, but in large areas of the World Ocean we find transfers less than 10−6 W m−2, in particular in the Southern Ocean and the eastern Pacific. The increase of transfer at the critical PSI latitude is clearly visible.
A net transfer of 0.1–0.3 TW is generally smaller than global estimates of internal tide generation. Based on an inverse model constrained by TOPEX/Poseidon altimeter data, Egbert and Ray (2003) estimate a barotropic-to-baroclinic energy transfer at the M2 frequency of 0.78 TW in the deep ocean, which agrees well with the 0.8 TW found in the numerical studies of Niwa and Hibiya (2011), Müller (2013), and Green and Nycander (2013, depending on the tidal dissipation parameterization, values as low as 0.5 TW were found). Owing to differences in model setup, resolution, and the practical definition of energy conversion, the modeled internal tide generation can vary substantially: Buijsman et al. (2016), for example, report a global conversion of 0.53 TW for the M2, S2, N2, and K2 constituents in a 1/12.5° HYCOM simulation, which resolves the first 2–3 internal tide modes. Global estimates of baroclinic M2 tide generation obtained from linear theory (Bell 1975; Llewellyn Smith and Young 2002) amount to 1.2 TW for depths greater than 500 m (Nycander 2005, 0.6 TW for depths greater than 1000 m), and, following a vertical normal mode decomposition, 0.71 TW for the first 10 modes and depths greater than 400 m (Falahat et al. 2014; 60% of this energy was transferred into modes 1 and 2). The latter number was reduced to 0.42 TW when the correction for supercritical slopes (Melet et al. 2013) was applied, and further to 0.32 TW when only subcritical slopes were considered. This is in line with the estimate of 0.52 TW of Vic et al. (2019), who calculated the global energy conversion below 700 only to avoid the necessity of such a correction for supercriticality.
The spread of numbers illustrates the uncertainty associated with the global quantification of internal tide generation; all of them, however, are larger by a factor of 2 or more than the energy transfer from the M2 internal tide to the continuum estimated above. Most of these global conversion estimates represent more than the first two baroclinic tide modes resolved in the STORMTIDE simulation and the lack of information on, for example, modal energy partition renders our estimate of the energy transfer from the internal tide to the wave continuum somewhat rudimentary. It seems that alternative dissipation mechanisms of internal tide energy play an important role, too.
7. Summary and conclusions
The baroclinic tide in the ocean is generated by the barotropic tide oscillating over ocean topography. It has imprinted the tidal frequency and the wavenumbers of the topography undulations. Propagating upward from the bottom it interacts by nonlinear wave–wave interactions with the omnipresent internal gravity wave field in the water column. Such interactions generally lead to a damping of the baroclinic tide and, since wave–wave interactions conserve energy, a transfer of the removed tidal energy to the background wave field. The wave–wave interaction process is an important step along the route of tidal energy to small-scale internal waves which might get unstable to break and mix ocean waters, enabling henceforth the large-scale overturning circulation of the oceans.
We have analyzed the tidal decay and spreading of energy across the background wave continuum with the same mathematical tools—weak interaction theory of a random wave ensemble—as used in previous work (see references in the introduction), leading to a kinetic equation for resonant triad wave interactions with the baroclinic tidal wave as one of the three partners and two continuum waves as the others. A new aspect of our research is that the transfer of tidal energy to the wave continuum is derived as a spectral rate in wavenumber–frequency space and analyzed in detail. The interaction process is investigated for different spectral forms of the GM class, varying slope and bandwidth in wavenumber space. Using the kinetic equation, we give evidence for the first time that it is predominantly PSI which shapes the transfer patterns. The tidal decay and energy transfer to the continuum is overwhelmingly supported by triads interacting according to the sum resonance ω = ω1 + ω2, where ω is the tidal frequency and ωj, j = 1, 2 are the interacting partners from the wave continuum. At the critical latitude (28.8° for M2) these are near-inertial waves, ωj ≃ f, which have the largest energy level in the GM-type spectra. But PSI also dominates the damping process at latitudes up to the equator, always with continuum partners with ωj = ω/2, which is the classical PSI scenario. The PSI interactions in the spectral framework are, however, not promoting an instability but rather a regular spectral transfer.
We have investigated the tidal decay process and the spectral energy exchange with the wave continuum in a process-oriented analysis where the Brunt–Väisälä frequency is held fixed but the latitude is changed to spell out the effect of PSI. Figures 3–11 refer to our findings. A globally realistic analysis follows and uses the globally varying Brunt–Väisälä frequency. Global maps of the decay times and energy transfer are displayed in the figures following Fig. 13.
On the process side our findings are summarized by the following:
The time scales of decay of the baroclinic tide crucially depend on the vertical wavelength (equivalent mode number) of the tide and the spectral shape of the continuum at low frequencies (where the overwhelming amount of energy is situated), but overwhelmingly on latitude since PSI of the M2 decay is limited to 28.8°N/S and outside only the less effective difference interactions operate.
The decay time is lowest just equatorward of the critical latitude with values depending on the spectral form of the wave continuum: about 10 days for GM76 and 20 days for a spectrum with no inertial cusp. The time scale increases toward the equator, and at latitudes higher than the critical latitude time scales of 60 days are found for GM76, decreasing toward northern latitudes. These values are for mode 1 and the standard parameters of the GM spectrum. Higher modes decay much faster.
The transfer of tidal energy to the wave continuum has a very distinct pattern. Equatorward of the critical latitude the transfer is dominated by the sum interaction PSI and generally is to wavenumbers lower than the one of the baroclinic tide at all frequencies. Difference interactions play here a minor role but determine the behavior at higher latitudes. The transfer of this type of interactions is generally positive at higher wavenumbers and negative at low ones (the tide grows at the expense of the continuum).
In the frequency spectrum of the transfer for sum interactions develops a δ-function behavior concentrated at half the tidal frequency (PSI) with an amplitude that strongly increases toward the critical latitude by orders of magnitude. Difference interactions show distinct peaks at ω + f and N at all latitudes but with a level which is orders of magnitude lower than the PSI transfer. In the frequency domain the transfer to the continuum is thus overwhelmingly to ω/2 and to a lesser amount to low-frequency near-inertial and to high-frequency near-buoyancy waves. The peculiar behavior is sketched in Fig. 16: given the tidal frequency ω and resonant sum and difference interactions ω = ω1 ± ω2 with partners from the wave continuum occur in restricted and separated domains of ω1 and latitude.
An attempt of a representative global estimate of decay times and energy transfer rates is introduced, based on maps of baroclinic tidal energy from the STORMTIDE simulation (Müller 2012) and internal wave energy estimates obtained by Pollmann et al. (2017) from Argo data (Riser et al. 2016). We also use the modal energy distribution given in de Lavergne et al. (2019) for the first five modes. With energy transfers below 0.1 TW from the STORMTIDE simulation and 0.3 TW from the analysis of de Lavergne et al. (2019), these estimates are lower than estimates of barotropic-to-baroclinic energy conversion (e.g., Falahat et al. 2014; Vic et al. 2019) at least by factor of 2, which suggests that other mechanisms such as remote dissipation, scattering at topography or mesoscale eddy structures might be important contributors to the global internal M2 tide energy cycle. To support or discard this hypothesis and to clarify the role of nonlinear wave–wave interactions for baroclinic tidal energy dissipation, it is necessary to reduce the uncertainty of our estimate by quantifying how the baroclinic tidal energy spreads over the vertical modes. Such an analysis and a subsequent refinement of our estimate are not only interesting in their own right, but help to better understand global internal gravity wave energetics and to improve mixing parameterizations based thereon (e.g., Olbers and Eden 2013; Eden and Olbers 2014).
This paper is a contribution to the Collaborative Research Centre TRR 181 Energy Transfer in Atmosphere and Ocean funded by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) Projektnummer 274762653. We acknowledge the use of the STORMTIDE data to Malte Müller, Zhuhua Li, and Jin-Song von Storch.
The Triad Interaction Coefficient
is the polarization vector of the Lagrangian displacement vector.
In this study the decay rate is evaluated for GM-type energy spectra of the form
which are horizontally isotropic and vertical symmetric. Here σ = ± stands for the sign of the vertical wavenumber, σ = signm. The spectrum is normalized (sum over σ and integral over wavenumber, frequency, and horizontal wave vector angle) to the total mechanical energy E = E(z) at depth z of the wave field. The shape functions of GM are defined by
The first frequency function is the original GM form, and the second has a modified (cuspless; referred to as CL) inertial peak. For the wavenumber slope s = 2 we are facing GM76 or GM76CL, respectively. We also use this spectral form with s = 3. The wavenumber-dependent part A(k, ω) has the same shape at each frequency and
characterizes its width in horizontal wavenumber space at each frequency. A corresponding scale for the vertical wavenumber m follows from the dispersion relation of internal gravity waves,
with Brunt–Väisälä frequency N = N(z). As standard parameters of GM we assume E = 3 × 10−3 m2 s−2, N = 5.2 × 10−3 s−1 (that is 3 cph), and . That makes . The Coriolis parameter f will be considered with its latitudinal variation. The normalization constants are as follows:
Tests of Convergence
Figure C1 shows the behavior of the integrand of the decay rate after integration over k1 as function of ω1 for difference and sum interactions. It is obvious that for the difference interaction the integrand passes smoothly across the PSI point ω1 = 2ω (left panel) whereas the sum integrand is heavily dominated by the behavior at the PSI point ω1 = ω/2 (right panel). Here the integrand develops a very narrow and large peak. To evaluate λ+ we thus have tested the convergence of the integration by implementing different cutoffs at a high wavenumber, by varying the number of grid points (e.g., 1000 or 5000) for the ω1 array, and by varying the location of the last grid point before the PSI point ω1 = ω/2. The latter is done by computing Δ = ω/2 − ω1k where ω1k is the last grid point before ω/2 in the ω1 array. Then ω1k is shifted to ω1k + dΔ with 0 < d < 1, that is, closer to the PSI point. For the right-hand side of the PSI point the procedure is repeated correspondingly.
A shift of the last grid point toward the PSI point, however, does not solve the convergence problem: it just increases the trapezoidal integration element without consideration of the sharply curved behavior of the integrand between the points. Cutting off the integrand at high wavenumbers decreases the decay rate but not very drastically, even a severe cutoff at has almost no effect.
We therefore increase the resolution between ω1k and ω1k + dΔ, introducing n/2 additional equidistant points (for the left-hand side with respect to ω/2; correspondingly for the right-hand side). For a converging behavior it is sufficient to use d = 0.9999, n = 2000 with the ω1 grid of originally 1000 equidistant points to obtain a relative accuracy of less than 1%.
In our notation f is assumed positive. Translation to the Southern Hemisphere is readily done.
Wave vectors in this study are three-dimensional, K = (k, m), where k is the horizontal and m the vertical wavenumber.
Note that in a Lagrangian treatment of the triad interactions the coefficient and the cross section are invariant with respect to permutations of all three index pairs. If the kinetic equation is derived from Eulerian equations, the interaction shares this full symmetry if resonant interactions are considered (Eden et al. 2019).