Using a novel model and observations from the Canary Current system of a coherent eddy of 100 km diameter, we explore the interaction between a realistic internal gravity wave field and this eddy. We study wave refraction and energy transfers between the waves and the eddy induced by the eddy’s lateral and vertical shear. Waves lose energy at the eddy rim by vertical shear and gain outside of the eddy rim by horizontal shear. We find large vertical wave refraction by vertical shear at the eddy rim, where waves break and mix density, which can have wide ranging effects on the ocean’s circulation. These results are important to understand the ocean’s mixing and energy cycle, and to develop eddy and wave parameterisations.

1 Introduction↩︎

Internal gravity waves (IGW) are ubiquitous in the ocean and their spectra are observed to be quite continuous in wavenumber and frequency space, as described by the Garrett-Munk model spectrum [1], [2]. IGWs are generated by a variety of mechanisms, primarily by winds and tides [3], and also by meso-scale eddies [4]–[6] in the ocean. Although IGWs evolve over fast time scales, their propagation and associated energy transfers can occur over large spatial scales, ranging from a few meters to hundreds of kilometers. Interactions between IGWs and meso-scale eddies, which evolve at much slower time scales, are important for the energy transfers in the ocean, however, these interactions make it challenging to separate, diagnose and quantify the associated oceanic processes. This study focuses on the evolution of IGWs propagating in a coherent meso-scale eddy using a novel Internal Wave Energy Model (IWEM) based on the six-dimensional radiative transfer equation which allows to separate some of the processes as described below.

During their propagation, IGWs are refracted by their changing environment and can interact with the mean flow, topography, small-scale turbulence, other waves such as Rossby waves, or with IGWs themselves. For instance, when IGWs propagate in a vertically sheared mean flow, energy can be transferred from the waves to the mean flow and vice versa. Similarly, due to horizontal shear waves can exchange energy with the mean flow. In a critical layer, IGWs get trapped with growing vertical wavenumber (but constant horizontal wavenumber) due to the vertical shear, which eventually leads to wave breaking and thus wave dissipation, sometimes called critical layer absorption. A similar process is wave capture [7] where horizontal shear is involved and all three wavenumber components grow without bounds which may also lead to wave breaking and dissipation.

During the breaking, the IGW energy is partly converted to small-scale turbulent energy (which is the dissipated wave energy) and is transferred to smaller scales, where it is converted into both internal energy by molecular dissipation and potential energy by turbulent density mixing. This transfer to potential energy by density mixing is known to be an important driver of the large-scale ocean circulation [3], [8]. However, mixing by wave breaking is not only generated by wave dissipation by vertical and horizontal wave refraction due the mean flow. The energy transport towards waves with smaller vertical wavelengths, which are prone to breaking as a result of resonant triad wave-wave interactions within the IGW field itself, is also thought to be an important driver of wave breaking and mixing. The importance of vertical and horizontal wave refraction for wave dissipation as compared to other processes, however, remains largely unknown, and is thus the focus of the present study.

Estimates of IGW energy dissipation and turbulent mixing from observational studies [9], [10] confirm that the complexity of the vertical and horizontal structure of IGW induced mixing is not captured by a constant vertical diffusivity, as typically used in ocean models. Addressing these issues, [11] and [12] proposed a promising concept based on a prognostic spectral energy balance for IGWs, called the Internal Wave Action Model (IWAM) with six-dimensions, with three dimensions in the spatial physical space and three in the wavenumber space (in addition to the time dimension). However, no actual application was published from the IWAM model due to analytical, practical, and computational limitations owing to the need to parameterize wave generation, energy transfers, and dissipation. Consequently, [13] drastically simplified this six-dimensional problem, and by integrating the spectral energy balance over the wavenumber space an energetically consistent parameterization for diapycnal diffusivity was given [14].

For the spectral energy balance to hold, it is necessary to assume slowly varying mean flow and stratification compared to wave lengths and periods, i.e. the WKBJ approximation. Recently, [15] derived a kinetic equation similar to the spectral energy transfer equation using the Wigner transform [16] without assumptions of spatial scale separation. However, in contrast to [15], we use in the present study the spectral energy balance as in [11] (and also [17]), for reasons discussed in the concluding section. We resolve the spectral energy balance in its full form with the novel numerical model called the Internal Wave Energy Model (IWEM), to understand the interaction between the geostrophically balanced mean flow in form of an coherent eddy, with IGWs in the ocean.

In Section 2, we revisit the propagation and refraction of IGWs in a mean flow with vertical and meridional shear in the Wentzel-Kramer-Brillouin (WKB) approximation. With the novel numerical code IWEM, we describe the evolution of a continuous wave field by the spectral energy balance of IGWs which includes transport, refraction, and wave-mean flow interactions within the limits of the WKB approximation. We use the turbulent kinetic equation to infer vertical mixing related to wave dissipation by vertical and horizontal wave refraction. In Section 3, numerical experiments with IWEM based on observations of an coherent meso-scale eddy in the Canary Current System are discussed to understand the role of the horizontal and vertical gradients of the mean flow. We close with a summary and discussion of the results.

2 Energetics of IGWs↩︎

The ocean can be described as a body of water with characteristic stratification \(N\), where IGWs are superimposed on a mean current \(\boldsymbol{U}\) with (usually) large temporal and spatial scales in contrast to the shorter periods and wavelengths of the waves. Considering a horizontally and vertically sheared horizontal mean flow \(\boldsymbol{U}_h(\boldsymbol{x}_h,z)\), stratification frequency \(N(z)\), and assuming that the temporal and spatial scales of \(\boldsymbol{U}_h\) and \(N\) are large compared to the periods and wavelengths of the waves, the WKB approximation can be applied. Freely propagating waves with wave ansatz \(\exp i(\boldsymbol{k}_{h}\cdot \boldsymbol{x}_{h} + k_zz -\omega_{enc}t)\) and wave vector \(\boldsymbol{k}=(\boldsymbol{k}_h,k_z)=(k_x,k_y,k_z)\) obey the local dispersion relation for the intrinsic frequency \(\omega=\omega_{enc}-\boldsymbol{k}_h\cdot \boldsymbol{U}_h\) of IGWs \[\label{eq2467} \left(\omega_{enc} - \boldsymbol{U}_h\cdot \boldsymbol{k}_h\right)^{2}= \omega(\boldsymbol{k}_{h},k_z,\boldsymbol{x}_{h},z,t) ^{2}= \frac{N^{2}k_{h}^{2}+f^{2}k_z^{2}}{k_{h}^{2}+k_z^{2}}\tag{1}\] where \(\omega_{enc}\) is the frequency of encounter (sometimes also called absolute frequency). The intrinsic frequency \(\omega\) follows from the wave properties of IGWs propagating in a motionless background and is observed by moving along with the mean current. The frequency of encounter \(\omega_{enc}\) is noticed by a stationary observer and is shifted by \(\boldsymbol{k}_h\cdot\boldsymbol{U}_h\), the Doppler shift.

2.1 Spectral energy balance↩︎

The IGW field in the ocean is depicted here by a superposition of wave packets. The wave packets have a dominant amplitude, wave vector, and frequency, that change according to the familiar ray equations [18]. Propagating wave packets can interact with the background stratification \(N\) or mean flow \(\boldsymbol{U}_h\), they can be influenced by wave-wave interactions and other non-linear processes, and there can be forcing or dissipation mechanisms that create new wave packets or destroy the existing ones. For instance, IGWs can be excited by wind-driven vertical velocity fluctuations at the base of the surface mixed layer, or by horizontal (tidal) currents flowing over bottom topography. However, here we exclude forcing, explicit dissipation, and non-linear processes and focus on progpagation, refraction, and mean-flow interaction in spin-down simulations of a given wave field.

To describe the energetics of an ensemble of IGWs with energy density in space and time \(E(\boldsymbol{x}_h,z,t)\) and corresponding wave action density \(A=E/\omega\), we define its energy power spectral density \(\mathcal{E}\), such that \(E=\int \mathcal{E} \, d \boldsymbol{k}_h \, d k_z\), and the corresponding wave action density \(\mathcal{A}=\mathcal{E}/\omega\). From wave action conservation of IGWs follows the spectral balance \[\label{eq2463462} \partial_t \mathcal{A} + \nabla_{\boldsymbol{x}_h}\cdot (\dot{\boldsymbol{x}}_{h} \mathcal{A}) + \partial_z(\dot{z}\mathcal{A}) +\nabla_{\boldsymbol{k}_h}\cdot (\dot{\boldsymbol{k}}_{h}\mathcal{A}) + \partial_{k_z}(\dot{k_z}\mathcal{A})=\mathcal{S}\tag{2}\] The left hand side of Eq. (2 ) describes changes of \(\mathcal{A}\) due to tendency, horizontal and vertical propagation by group velocity \((\dot{\boldsymbol{x}}_h,\dot{z})=(\nabla_{\boldsymbol{k}_h}\omega_{enc},\partial_{k_z}\omega_{enc})\), and refraction in horizontal and vertical wave vector components by \((\dot{\boldsymbol{k}}_h,\dot{k}_z)=-(\nabla_{\boldsymbol{x}_h}\omega_{enc},\partial_{z}\omega_{enc})\). The term \(\mathcal{S}\) on the right hand side describes processes violating wave action conservation, such as wave-wave interaction, forcing, or dissipation. Reformulating Eq. (2 ) for energy \(\mathcal{E}= \omega\mathcal{A}\), an additional term arises on the right hand side \[\label{eq2463463} \partial_t \mathcal{E} + \nabla_{\boldsymbol{x}_h}\cdot (\dot{\boldsymbol{x}}_{h}\mathcal{E})+\partial_z (\dot{z}\mathcal{E}) +\nabla_{\boldsymbol{k}_h}\cdot (\dot{\boldsymbol{k}}_{h}\mathcal{E}) + \partial_{k_z}(\dot{k}_z\mathcal{E})=\omega\mathcal{S}+ \left( \mathcal{E}/\omega \right) \dot{\omega},\tag{3}\] where the changes in frequency are defined as: \[\label{eq4} \dot{\omega}= - \left( \dot{\boldsymbol{x}}_h- \boldsymbol{U}_h \right) \cdot (\boldsymbol{k}_h\cdot \nabla_{\boldsymbol{x}_h} \boldsymbol{U}_h) -\dot{z}(\boldsymbol{k}_h \cdot \partial_z \boldsymbol{U}_h)\tag{4}\] with \(\dot{\boldsymbol{x}}_h- \boldsymbol{U}_h\) being the instrinsic group velocity. The term \((\mathcal{E}/\omega) \dot{\omega}\) in Eq. (3 ) is a source of mechanical energy gathered along the wave path by interaction of the wave field with the background flow. This energy can be transferred in both directions, i.e. from the mean flow to the wave field or the other way around. Note that wave-wave interaction and wave forcing as part of \(\mathcal{S}\) are neglected in Eq. (3 ), while dissipation is implicitly present as discussed next.

2.2 Wave breaking and dissipation↩︎

As noted earlier, in this study we exclude any explicit wave dissipation. On the other hand, dissipation (and thus non-zero \(\mathcal{S}\) in Eq. (3 )) is still implicitly present in our model by a possible wave energy flux \(F_{crit}\) across the large vertical wavenumber threshold \(k_z^c>0\) \[\begin{align} \label{fcrit} F_{v}= \int \left\{ \begin{array}{ll} \left[ \dot{k_z} \mathcal{E} \right]_{k_z = k_z^c} & if ~\dot{k_z} >0 \\ 0 & else \end{array} \right\} \, d \boldsymbol{k}_h - \int \left\{ \begin{array}{ll} \left[ \dot{k_z} \mathcal{E} \right]_{k_z = -k_z^c} & if ~\dot{k_z} <0 \\ 0 & else \end{array} \right\} d \boldsymbol{k}_h \end{align}\tag{5}\] and the flux \(F_{h} = F_{h}^x +F_{h}^y\) across the large horizontal wavenumber threshold \(k_h^c>0\) \[\begin{align} \label{fcapt} F_{h}^x &=&\int \left\{ \begin{array}{ll} \left[ \dot{k_x} \mathcal{E} \right]_{k_x = k_h^c} & if ~\dot{k_x} >0 \\ 0 & else \end{array} \right\} d k_zd k_y -\int \left\{ \begin{array}{ll} \left[ \dot{k_x} \mathcal{E} \right]_{k_x = -k_h^c} & if ~\dot{k_x} < 0 \\ 0 & else \end{array} \right\} d k_zd k_y ~~~~ \end{align}\tag{6}\] and similar for \(F_{h}^y\). We set \(\pm k_z^c\) and \(\pm k_h^c\) to the maximal values of the model domain in wavenumber space. The positive definite energy fluxes \(F_{v}\) and \(F_{h}\) across the wavenumber thresholds are taken to equal the energy transfer from the waves to small-scale turbulence and thus density mixing or further dissipation to internal energy. This happens, for example, for the case of wave breaking at critical layers, where waves are refracted with growing vertical wavenumber but constant horizontal wavenumber due to the vertical shear only. Note, however, that often there is both horizontal and vertical shear, such that it is difficult to relate the flux \(F_v\) only to the special case of critical layer absorption. Further, \(N(z)\) can also drive a flux \(F_v\) even without any vertical shear. In the following, we refer to wave dissipation by vertical wave refraction for wave energy to cross the vertical wavenumber threshold \(k_z^c\) generating the flux \(F_{v}\), see Eq. (5 ), and similar for \(F_h\) from Eq. (6 ). We will test below if the vertical wavenumbers (or related Richardson number) of the wave field at the thresholds are already large (or small) enough to warrant this approach.

To treat the fate of the energy transfer \(F_{v}+F_{h}\) to small-scale turbulence, we consider the turbulent kinetic energy (TKE) equation. Under the assumption of stationarity and horizontal homogeneity and neglection of the transport terms, the TKE equation simplifies to a local balance of production terms and dissipation \(\epsilon\), i.e. the Osborn-Cox relation (see [19]) given by \[-\overline{\boldsymbol{u}_{h}^{'}w^{'}}\cdot\partial_z \overline{\boldsymbol{u}}_{h}+\overline{b^{'}w^{'}}-\epsilon = 0 \label{eq2464461}\tag{7}\] The primed terms in Eq. (7 ) stand for small-scale turbulent fluctuations, \(\overline{\boldsymbol{u}}_{h}\) denotes wave velocity, and \(\epsilon\) dissipation of TKE. The vertical turbulent buoyancy flux is given by \(\overline{b^{'}w^{'}}\), with buoyancy \(b=-g\rho/\rho_{0}\) and density \(\rho\). The shear production term \(\overline{\boldsymbol{u}_{h}^{'}w^{'}} \cdot\partial_z \overline{\boldsymbol{u}}_{h}\) represents the wave-mean flow energy exchange, i.e. the flux \(F_{v}+F_{h}\). With the diffusive parameterisation \(\overline{b^{'}w^{'}} = \kappa N^2\), it follows that \[\label{eq2464462} \kappa=\frac{\gamma}{1+\gamma}\frac{F_{v}+F_{h}}{N^2}\tag{8}\] where \(\kappa\) is the turbulent diffusivity and \(\gamma\) the relative mixing efficiency, (see [20]). The energy transfer by wave breaking \(F_{v}+F_{h}\), assumed to be equal to the shear production term \(\overline{\boldsymbol{u}_{h}^{'}w^{'}} \cdot\partial_z \overline{\boldsymbol{u}}_{h}\), can either be dissipated into heat or leads to an increase of potential energy by density mixing. The fraction of the energy dissipated into heat or converted into potential energy by density mixing is represented by the relative mixing efficiency \(\gamma\). Based on fine scale measurements by [19] and e.g. [21], the relative mixing efficiency is \(\gamma\approx 0.2\) for stratified flows such as the ocean. Typical values of local turbulent diffusivity in the interior of the ocean are \(\kappa\approx 10^{-5} \, \rm m^{2}s^{-1}\), inferred from typical observational estimates of dissipation in the range of \(\epsilon \approx 10^{-10}\) to \(10^{-9} \, \rm m^{2}s^{-3}\).

2.3 The Garrett-Munk model spectrum↩︎

The Garrett-Munk (GM) model spectrum [1], [2] is a representation of the energy density spectrum \(\mathcal{E}\) of IGWs derived from observations, and is used here as the initial condition for our simulations described below. The GM model spectrum describes a horizontally isotropic and vertically symmetric spectrum with a spectral power law of -2 in frequency and vertical wavenumber. There are several different versions of the GM model spectrum, here we use the following form \[\label{N123} \mathcal{E}(k_z,\omega,\phi,z) = \frac{E(z)}{2 \pi}A(k_z,\omega)B(\omega) ~~,~~ \mathcal{E}(\boldsymbol{k},z)= \mathcal{E}(k_z,\omega,\phi,z) \, J^{-1}\tag{9}\] with wavenumber angle \(\phi\) and \[\begin{align} A(k_z,\omega) = \frac{\tilde{A}(k_z/k_z^{*}(\omega))}{k_z^{*}(\omega)} ~~,~~ \tilde{A}(\zeta)=\frac{n_{a}}{1+\zeta^{2}} ~~,~~ B(\omega)=\frac{n_{b} f}{\omega \sqrt{\omega^{2} - f^{2}}} \\ J=\frac{\partial(\boldsymbol{k})}{\partial(k_z,\omega,\phi)} ~~, ~~ k_z^{*}(\omega)=\frac{ \sqrt{N^{2}-\omega^{2}}}{c^{*}} ~~,~~ c^{*}=\int_{-h}^{0} \frac{N(z)}{j^* \pi} \,dz \end{align}\] where \(J\) denotes the Jacobi determinant for the coordinate transformation \((k_z,\omega,\phi)\to(\boldsymbol{k})\), \(k_z ^{*}\) the bandwidth of the vertical wavenumber spectrum equivalent to the vertical mode number \(j^{*}=10\) and \(\zeta=k_z/k_z^{*}(\omega)\). The energy density spectrum \(\mathcal{E}\) is factorized with the wavenumber shape function \(A\) and the frequency shape function \(B\). Both shape functions integrate to one by choice of the parameter \(n_a\) and \(n_b\). The energy density spectrum is then normalized to the total mechanical energy of the wave field in physical space, such that \[\begin{align} E(z)=E_{0}N(z)/N(z_0) = \int_{f}^N \int_{-\infty}^\infty \int_0^{2\pi} \mathcal{E}\, d \omega d k_z d \phi \end{align}\] where \(E_0\) denotes the total energy content \(E_0\) at \(z_0\), and the factor \(N(z)/N(z_0)\), the so-called WKB-scaling. We use \(E_0 = 3 \times 10^{-3} \, \rm m^2/s^2\) and \(z_0=0 \, \rm m\) as classical parameter, but note that this is probably an underestimation of the typical wave energy [22]. While the original GM spectrum uses analytical expressions for the parameters \(n_a\) and \(n_b\), we choose them such that \(E(z) = \int_{-{\boldsymbol{ k }}_h^c}^{{\boldsymbol{ k }}_h^c} \int_{-k_z^c}^{k_z^c} \mathcal{E}\, d \boldsymbol{k}_h d k_z\) at all depths.

2.4 Instability at wavenumber thresholds↩︎

For the assessment of the stability of the wave field, the shear generated by the total wave ensemble given by \(\mathcal{E}\) can be analysed and a critical value for the vertical wavenumber threshold \(k_z^c\) calculated. Wave instabilities can be assessed by the bulk Richardson number \(Ri_{\mathrm{Bulk}}=N^2/(\partial_z {\boldsymbol{ u }})^{2}\), where \(Ri_{\mathrm{Bulk}}<1/4\) characterises an unstable situation and \(Ri_{\mathrm{Bulk}}>1/4\) a stable one. Following [23], \(k_z^c\) can be inferred in terms of a critical Richardson number sustained by the wave field, so that waves with \(|k_z|>k_z^c\) are particularly prone to break, \[Ri_{\mathrm{Bulk}}^{-1}=\int_{-k_z^c}^{k_z^c } d k_z \int_f^N d \omega \, Ri^{-1}(k_z,\omega) =\int_{-k_z^c }^{k_z^c} dk_z \int_f^N d \omega \, \frac{\omega^{2}+f^{2}}{\omega^{2}N^{2}}k_z^{2}\mathcal{E}(k_z,\omega) \label{ri}\tag{10}\] where \(Ri(k_z, \omega)\) is the spectrum of the Richardson number, for which \(\partial_z {\boldsymbol{ u }}\) is expressed by the shear spectrum in the second step.

However, the value for \(k_z^c\) is not uniquely given in literature. Different authors use corresponding wavelengths from \(\lambda=10 \, {\rm m}\) down to \(\lambda=1 \, {\rm m}\). According to [24] and in agreement with [25], the IGW spectral shape obtained from measurements in the Northwest Atlantic Ocean comprises primarily of IGWs with vertical wavelengths \(10\, {\rm m}\leq\lambda< 100 \, {\rm m}\), finestructure and decaying internal waves between \(1 \, {\rm m}\leq\lambda< 10 \, {\rm m}\), and small scale turbulence at centimeter scales. [24] find that the observed spectrum results in Ri\(_{\mathrm{Bulk}}=1\) and suggests a cut-off vertical wavelength \(\lambda_z^c=10m\). [23] predicted Ri\(_{\mathrm{Bulk}}^{-1}(k_z^c)=0.5\) with a roll-off at \(k_z^c=0.6 \, \rm m^{-1}\) (wavelength of about 10 m), showing more agreement with observations by [26], where \(Ri^{-1}_{\mathrm{Bulk}} \approx 0.3\) [27] for spectra close to the model spectrum by [28]. Moreover, [29] found Ri\(_{\mathrm{Bulk}}\) between one and three from observations in the upper ocean, where \(Ri_{\mathrm{Bulk}}=1\) indicates high non-linearity of the IGW field and the upper bound of \(Ri_{\mathrm{Bulk}}=3\) indicated saturated situations with maximum allowable shear.

Here, \(k_z^c = 1 \, \rm m^{-1}\) (or wavelength 6.28 m) is used and a zonal and meridional wavenumber threshold of \(k_h^c = 0.5 \, \rm m^{-1}\). This choice results in a \(Ri_{\mathrm{Bulk}}^{-1}=0.9\) under consideration of a depth-invariant stratification of \(N=5 \times 10 ^{-3} s^{-1}\) and the GM model spectrum for \(\mathcal{E}\) from above. For the case of a depth-dependent exponential stratification \(N=N_1(z)\) given by Eq. (11 ), we find \(Ri_{\mathrm{Bulk}}^{-1}\approx 1.4\) at the surface, and \(Ri_{\mathrm{Bulk}}^{-1}\approx 0.2\) in the deep ocean. These predictions of IGW instability are in partial agreement with [24]. However, \(Ri_{\mathrm{Bulk}}^{-1}\) depends on the total wave energy, i.e. on the parameter \(E_0\) of the GM model spectrum. With larger \(E_0\), \(Ri_{\mathrm{Bulk}}^{-1}\) will increase and \(k_z^c\) will decrease. Note that three times larger values than used here for \(E_0\) are found by [22] from ARGO float observations, i.e. the value for \(E_0\) used here is most likely a low bias.

3 Numerical simulations↩︎

The numerical simulations in this study are performed with the Internal Wave Energy Model (IWEM), which integrates Eq. (3 ) in time. A detailed description of the model is provided in Appendix A. All numerical experiments are spin-down simulations initialized with the IGW field given by the GM model spectrum from Section 22.3, embedded in a stationary horizontally and vertically varying mean flow and a stationary, vertically varying, but horizontally constant stratification. The integration time is six days, for which we assume the mean flow and stratification to be stationary. There is no wave forcing, and dissipation only by the energy transports across the vertical and horizontal wavenumber thresholds as described in Section 22.2. The initial wave field is horizontally homogeneous, isotropic in \(\boldsymbol{ k}_h\) and symmetric in \(k_z\), but the total wave energy is vertically increasing towards the surface since the GM-spectrum is proportional to the local stratification \(N\) (see Section 22.3).

The domain for the simulation has a horizontal extent of \(240 \times 240 \, {\rm km}\) and its vertical extent is \(L_z=4\, {\rm km}\) deep. In physical space, we use \(40 \times 40 \times 60\) grid points in the zonal, meridional, and vertical directions respectively. The resolution in the horizontal direction is \(\Delta x = \Delta y= 6 \, \rm km\). The vertical resolution decreases exponentially from the surface layer with \(\Delta z =15 \, \rm m\) to the bottom layer with \(\Delta z=74 \, \rm m\). In wavenumber space we use \(60 \times 60 \times 60\) grid points, with horizontal wavenumber \(\Delta k_{x,y}=4 \times 10^{-3} \, \rm m^{-1}\) at \(k_{x,y}=0 \, \rm m^{-1}\) increasing to \(\Delta k_{x,y}=18 \times 10^{-3} \, \rm m^{-1}\) at \(\pm k_{x,y}^c\), and vertical wavenumber \(\Delta k_{z}=8 \times 10^{-3} \, \rm m^{-1}\) at \(k_{z}=0 \, \rm m^{-1}\) increasing to \(\Delta k_{z}=36 \times 10^{-3} \, \rm m^{-1}\) at \(\pm k_{z}^c\). The model time step is \(2s\). The total number of grid points is about \(20 \times 10^9\), which corresponds roughly to a global ocean model with horizontal resolution of 1 – 2 \(\rm km\). Using even finer grids, the model quickly becomes computationally very expensive such that the simulations presented here are already at the upper end of available resources to us.

The background flow in our model simulations is inspired by measurements of a coherent meso-scale eddy in the tropical northeastern Atlantic Ocean provided by the REEBUS¹ project and shown in Fig. 1 (see also Appendix B), which we think is representative for a typical coherent eddy at mid latitudes. The eddy measurements were made in the Eastern Boundary Upwelling System at the Mauritanian coast west of Africa in the Canary Current System in December 2019, between 14N to 15N and 24W to 26W and was named \(CE\_2019\_14N\_25W\) by REEBUS. The observed coherent eddy is of cyclonic type with positive zonal velocity south of the eddy center and negative north of it (Fig. 1). At the eddy rim, at a distance of around 45 km from the eddy center on either side, the velocities are the largest and reach \(0.4 \, \rm ms^{-1}\) and \(-0.6 \, \rm ms^{-1}\) at the surface. Between the eddy rim, the horizontal change of the velocity is approximately linear. Outside of the eddy, velocity decays exponentially in the horizontal. As suggested by [31] this motivates the idealized eddy structure as detailed in Appendix C, which is also shown in Fig. 1 (grey contour lines). Since the velocity observations reach only to 1000 \(\rm m\) depth, we use the first baroclinic vertical mode to extrapolate the observations down to the bottom at 4000 \(\rm m\). We also have a further simulation with an anticyclonic background flow which is simply a mirror of the cyclonic one. However, if not otherwise noted we always refer to the cyclonic eddy in our discussion of the results.

The background stratification is computed from in-situ measurements at the eddy center. A mean background stratification \(\overline{N}\) is computed by smoothing the observations with a moving mean of 150 m with a minimum of ten data points, shown in shown in Fig. 1. \(\overline{N}\) is then fitted with an exponential function \(N_{1}(z)\) given by \[\begin{align} \label{N1} N_{1}(z)&=& 2.5 \times 10^{-3} \, {\rm s^{-1} } + 12.5 \times 10^{-3} \, {\rm s^{-1} }\, \exp{\frac{z+50 \, {\rm m} }{100 \, {\rm m}} } \end{align}\tag{11}\]

Before we proceed to a discussion of the results of the model simulations, we check the validity of the WKB approximation, inherent to the model approach here. For the approximation to hold, it is necessary that the background flow \({\boldsymbol{ U }}({\boldsymbol{ x }},z)\) and stratification \(N(z)\) are only slowly changing compared to vertical and horizontal wave length, i.e. \(|{\partial}_z N|/ |k_z N|\ll 1\), \(|{\partial}_z {\boldsymbol{ U }}_h |/ |k_z {\boldsymbol{ U }}_h| \ll 1\), and \(|{\partial}_{{\boldsymbol{ x }}} {\boldsymbol{ U }}_h | / |{\boldsymbol{ k }}_h {\boldsymbol{ U }}_h| \ll 1\). Fig. 2 shows the maximal corresponding ratios. For both \({\partial}_z N\) and \({\partial}_z {\boldsymbol{ U }}\) the ratios are indeed very low, except for very small \(k_z\) and the zero crossing of \({\boldsymbol{ U }}\) at mid depth. The ratios for lateral shear stay well below one everywhere. We do not regard the large ratios at the smallest \(k_z\) as a problem, since here turning points are expected which do not show any singular behaviour for the wave energy, and which are also limited by the finite depth. Therefore, the WKB approximation appears to be valid for our model configuration.

3.1 Waves in a coherent eddy↩︎

We first discuss a simulation with stratification \(N_1(z)\) and the idealized cyclonic coherent eddy. The speed of the cyclone is shown in a horizontal plane close to the surface in Fig. 3 (a) and a meridional cross-section of zonal velocities of the eddy center in Fig. 3 (d). We choose to show in Fig. 3 the WKB-scaled wave energy \(\mathcal{E}\), i.e. \(N_0 \mathcal{E}/N(z)\), so that the depth-dependency of \(\mathcal{E}\) vanishes at time \(t=t_0\). We also use the so-called stretched vertical coordinates \(z^*(z)\) instead of \(z\) with \(dz^* = N(z)/N_0 \, d z\), where \(N_0=\int_{-h}^0 (N(z)/h) dz\) and \(h=4000 \, \rm m\), to better show the dominant processes close to the surface.

After six days of simulation, we see already a notable change in wave energy in Fig. 3 (b,e). Wave energy has remained relatively stable at the eddy center over the whole water column, while wave energy decreases the most at the eddy rims close to the surface (Fig. 3 b,e). Wave energy changes over time (Fig. 3 c,f) show an overall decrease in wave energy towards the surface, with maxima at the vicinity of the eddy rim and a minimum at the eddy center. The waves at the eddy center have gained a small amount of energy while IGWs at the eddy rim have lost energy throughout the water column. It will be shown below that a similar inhomogeneous distribution of wave energy within the coherent eddy can also be seen in the available observations. This energy change is the main result of this paper which we aim to understand in the following.

In a simulation with an anticyclonic eddy using the same stratification \(N_1(z)\) but simply the negative of the eddy velocities of the cyclone, the energy changes are almost identical. In both spin-down simulations, there is a total wave energy loss of \(14.7\%\) relative to the initial energy integrated over the whole domain. At the eddy rim there is a relative energy loss of \(28.8\%\) integrated over the water column, while the eddy center gains \(2.8\%\). Note, that the asymmetries in Fig. 3 are due to the finite size effect of simulating the propagation of IGWs in a radial symmetric eddy in Cartesian coordinates. The results are symmetric in the eddy core, but the discontinuity in the eddy velocities at the circular eddy rim lead to slight asymmetries from the eddy rim outwards.

3.2 Wave-mean flow interaction↩︎

We first consider the wave energy change due to wave-mean flow interaction. The wave-mean flow interaction is given by the term \(\mathcal{E}\dot{\omega}\omega^{-1}\) in Eq. (3 ) and depends thus directly on the changes in frequency, see Eq. (4 ). We split the discussion into effect of lateral and vertical shear, i.e. \[\begin{align} \mathcal{E}\dot{\omega}\omega^{-1} = \mathcal{E}\dot{\omega}_h\omega^{-1} + \mathcal{E}\dot{\omega}_z\omega^{-1} \end{align}\] with \(\dot{\omega}_h = - \left( \dot{\boldsymbol{x}}_h- \boldsymbol{U}_h \right) \cdot (\boldsymbol{k}_h\cdot \nabla_{\boldsymbol{x}_h} \boldsymbol{U})\) and \(\dot{\omega}_z = -\dot{z}(\boldsymbol{k}_h \cdot \partial_z \boldsymbol{U})\). All terms are shown in Fig. 4.

Consider the effect of lateral shear first. The horizontal intrinsic group velocity can be written as \(\dot{\boldsymbol{x}}_h- \boldsymbol{U}_h= v(k,|k_z|) (\cos \phi, \sin \phi)\), with \(\boldsymbol{k}_h = k (\cos\phi, \sin\phi)\) and appropriate function \(\nu\). Then \(\dot{\omega}_h\) becomes \[\begin{align} \dot{\omega}_h = -v k \left( \cos^2 \phi \partial_x U + \cos \phi \sin \phi ( \partial_y U + \partial_x V )+ \sin^2 \phi \partial_y V \right) \end{align}\] Since \(\tau_{s}\)=\(\partial_y U + \partial_x V =0\), and \(\partial_x U = \partial_y V =0\) within the idealized eddy (Fig. 1 c,d), it holds that \(\dot{\omega}_h=0\) there, and the wave-mean flow interaction by \(\mathcal{E}\dot{\omega}_h \omega^{-1}\) vanishes from the eddy rim towards the center, as seen in Fig. 4 a) and d). Outside the eddy rim, the waves gain energy due to \(\mathcal{E}\dot{\omega}_h \omega^{-1}\), which amounts to 6.7% of the initial wave energy when integrated over the whole model domain, with a maximum close to the eddy rim near the surface of 48.8% relative to the initial water column’s energy.

Outside the eddy rim, the waves gain energy due to \(\mathcal{E}\dot{\omega}_h \omega^{-1}\) because of a developing anisotropic spectrum, which can be understood as follows: For an initially isotropic spectrum \(\mathcal{E}/\omega = \mathcal{A}(\phi,k,|k_z|)\) with \(\partial_\phi \mathcal{A}=0\), we have \[\begin{align} \int \mathcal{A}\dot{\omega}_h d k d \phi d k_z = - \pi \left( \partial_x U + \partial_y V \right) \int \mathcal{A} v k d k d k_z = 0 \end{align}\] since \(\partial_x U + \partial_y V =0\) for the idealized eddy, \(\int_0^{2 \pi} \cos^2 \phi d \phi=\int_0^{2 \pi} \sin^2 \phi d\phi = \pi\), and \(\int_0^{2 \pi} \cos \phi \sin \phi d\phi=0\). For an isotropic spectrum without \(\phi\) dependency, as specified in the initial conditions, the wave-mean flow interaction by lateral shear is thus vanishing everywhere in the idealized eddy. Only with anisotropy in \(\mathcal{E}\), non-vanishing wave-mean flow interaction by lateral shear develops outside the eddy rim, as seen in Fig. 4 (a,d).

The horizontal anisotropy in \(\mathcal{E}\) at the eddy rim is generated in the first place by the effect of vertical refraction, energy transfer to mean flow by vertical shear, and energy dissipation by vertical shear, which is largest at the eddy rim, and discussed in detail below. Using single column simulations (excluding dimensions \(x\) and \(y\) but not \(k_x\) and \(k_y\), not shown) with \({\boldsymbol{ U }}(z)\) from the eddy rim and the same \(N(z)\) as before, i.e. excluding all horizontal processes, we see a similar horizontal anisotropy developing in \(\mathcal{E}\). The effect of vertical refraction, energy transfer to mean flow by vertical shear, and energy dissipation by vertical shear deforms the initially isotropic isolines (circles) of \(\mathcal{E}\) (or \(\mathcal{A}\)) in the \(k_x - k_y\) plane towards ellipses with a major axis roughly oriented towards the eddy center with maximal energies inside the ellipse. This is shown in Fig. 5 for the northern and eastern position at the rim.

At the eddy center, however, vertical shear and thus the effects of vertical refraction, energy transfer to mean flow, and energy dissipation are vanishing and no anisotropy is generated by them. Also the horizontal refraction \(\dot{\boldsymbol{k}}_h\) does not generate any anisotropy in \(\mathcal{E}\) at the center, since here \(\partial_x U = \partial_y V=0\) and \(\partial_y U = - \partial_x V\) (Fig. 1 c,d), and so \[\dot{\boldsymbol{k}}_h = - \nabla_{\boldsymbol{x}_h} \omega_{enc} =-\nabla_{\boldsymbol{x}_h} ( \boldsymbol{k}_h\cdot \boldsymbol{U}_h ) = ( k_y , - k_x ) \partial_y U \sim \boldsymbol{z} \times \boldsymbol{k}_h\] such that transports by \(\dot{\boldsymbol{k}}_h\) are perpendicular to \(\boldsymbol{k}_h\) with no effect on the isotropic isolines (circles) of \(\mathcal{E}\) (see also Fig. 5 a). Since the horizontal intrinsic group velocity \(\dot{{\boldsymbol{ x }}}_h - {\boldsymbol{ U }}\) does not depend on wave angle \(\phi\) and thus also does not generate any anisotropy from an initially isotropic spectrum, and since for the Doppler shift \({\boldsymbol{ U }}\) in the group velocity the same holds, \(\mathcal{E}\) stays perfectly isotropic (Fig. 5 d). Only an isotropic loss of energy due to horizontal wave propagation \(\dot{{\boldsymbol{ x }}}_h \mathcal{E}\) is visible at higher wavenumbers in Fig. 5 (d). This isotropic loss adds to the other effects also at the rim.

An exactly north-south oriented ellipse at the eddy rim is still not effective in generating mean-flow interaction by \(\dot{\omega}_h\) since it generates an anomaly \(\mathcal{A}'\sim - \cos 2 \phi\) with \[\begin{align} \int \mathcal{A}' \dot{\omega}_h d \phi \sim - \int v k ( - \cos 2 \phi ) \cos \phi \sin \phi ( \partial_y U + \partial_x V ) d \phi = 0 \end{align}\] The same is true for a general orientation of the ellipse towards the eddy center. However, an anomaly with \(\sim \sin 2 \phi\) becomes effective at the northern position at the rim. In fact, if the ellipse is slightly tilted from north-south direction direction as seen in Fig. 5 (b), it has a component \(\sim \sin 2 \phi\). With \({\partial}_x V>0\) but \({\partial}_x V \ll {\partial}_y U\) the transports have indeed a tilt in positive \(\phi\) direction and an additional anomaly \(\mathcal{A} \sim -\sin 2 \phi\) develops with \[\begin{align} \int \mathcal{A}' \dot{\omega}_h d \phi \sim \int v k \sin 2 \phi \cos \phi \sin \phi ( \partial_y U + \partial_x V ) d \phi = v k ( \partial_y U + \partial_x V ) \pi /2 >0 \end{align}\] which generates gain of wave energy related to \(\dot{\omega}_h\), as seen in Fig. 5 (b). The same is true for a general orientation of the ellipse towards the eddy center, see e.g. Fig. 5 (c). The slight tilt in positive \(\phi\) direction of the orientation of the anisotropic ellipses in \(\mathcal{E}\) thus explains the energy gain outside the eddy rim by the effect of lateral shear.

The small tilt is generated by horizontal refraction \(\dot{\boldsymbol{ k}}_h\). North or south of the eddy rim is \(\partial_x U= \partial_y V=0\) and \(\partial_y U \gg \partial_x V >0\) (Fig. 1 c,d). The wave refraction is here roughly given by \[\dot{\boldsymbol{k}}_h = - \nabla_{\boldsymbol{x}_h} \omega_{enc} =-\boldsymbol{k}_h\cdot \nabla_{\boldsymbol{x}_h} \boldsymbol{U}_h = -( k_y \partial_x V, k_x \partial_y U ) \approx -(0, k_x ) \partial_y U\] which is also indicated in Fig. 5 b). This flow in the \(k_x\)-\(k_y\)-plane tilts the ellipse slightly to generate the wave energy gain.

In case of the anticyclone, (Fig. 5 e and f), lateral wave-mean flow interaction can be explained analogously, keeping in mind that the horizontal velocities rotate in opposite direction. The development of anisotropy as well the orientation of the ellipses towards the eddy center is identical to the cyclone, now however with \(\dot{\boldsymbol{k}}_h\) pointing in opposite direction, which thus generates the slight tilt opposite to the one of the cyclonic case.

Fig. 4 (b,e) shows that the effect of the vertical shear in \(\mathcal{E} \dot{\omega}/\omega\) is negative, i.e. the waves lose energy. This is because of the following reason: The sign of the vertical group velocity \(\dot{z}\) in \(\dot{\omega}_z =-\dot{z}(\boldsymbol{k}_h \cdot \partial_z \boldsymbol{U}_h)\) depends on the sign of \(k_z\). Since that sign changes when the wave is reflected at the top or bottom boundary while all other wave properties are conserved, the effect due to \(\dot{\omega}_z\) is completely reversed after reflection, without any net effect after a full cycle of reflection at surface and bottom. Only dissipation during the wave propagation can break this cycle, and in our model this dissipation is the energy flux \(F_v\) across the vertical wavenumber boundaries. When the waves propagates towards a critical layer, they always lose energy to the mean flow, thus we see only wave energy loss by the wave-mean flow interaction related to \(\dot{\omega}_z\). To understand this, consider a wave ray propagating in a vertically sheared mean flow and assume also constant \(N\) for simplicity. The wave ray has wave action \(A=E/\omega\) starting at time \(t_0\) with initially zero mean flow and approaches a critical layer at time \(t\). Because of wave action conservation \(\dot{A}=0\), and \(\omega(t) = |f| = \omega(t_0) -| U_h k_h |\) \(\to\) \(E(t)/|f|= E(t_0)/\omega(t_0)\). Since \(|f|/\omega(t_0) = |f|/(|f|+ |k_h U_h |) <1\), wave energy is lost to the mean flow during the approach of the ray into the critical layer. In the critical layer, the wave breaks and the remainder of the energy is transferred to small-scale turbulence (or generates secondary waves). In our model, the remainder of the energy propagates beyond the positive or negative vertical wavenumber boundary which we diagnose as energy flux \(F_v\), and which is available for mixing. Note that only near-inertial waves transport a significant fraction of energy into the critical layer and to turbulence.

We do not have the rather special case of a pure critical layer due to the additional horizontal shear. We, however, observe a “critical layer like” wave dissipation by vertical wave refraction, such that the above explanation transfers to our case. To justify that, Fig. 6 shows the outward flux \(\dot{k_z}\mathcal{E}\) at the \(k_z^{c}\)-threshold at the eddy rim where critical layer absorption is strongest (see below). The outward dissipative flux \(\dot{k_z}\mathcal{E}\) is largest at \(k \to 0 \, \rm m^{-1}\), i.e. most energy leaves at \(\omega \to f\), indicative of a critical layer. In Fig. 6 (b) we show the vertical group velocity \(\dot{z}\) (color), which is highest for \(k_z \to 0 \, \rm m^{-1}\) and tends to zero for large \(k_z\). The figure also shows that the intrinsic horizontal group velocity \(\partial_{k}\omega(k,k_{z})\) (or \(\dot{\boldsymbol{x}_h}-\boldsymbol{U}_h\)) also tends to zero for large \(k_z\) (contours). We can interpret all these features as "critical layer like" wave dissipation behaviour.

Towards the eddy rims, the vertical shear in the mean flow increases, thus \(\dot{\omega}_z\) strengthens, with the result that the waves lose more energy to the mean flow there. Towards the eddy center \(\partial_z \boldsymbol{U}_h \to 0\) and the vertical wave-mean flow interaction vanishes. The magnitude of the vertical wave-mean flow interaction is increasing towards the surface because of the stronger vertical shear of the mean flow. The total relative energy loss by vertical wave-mean flow interaction over the whole domain amounts to 5% of the initial wave energy, with 27.3% relative to the initial energy at the northern position at the rim.

For the case of the anticyclone, the vertical wave-mean flow interaction is again analogous to the case of the cyclone, keeping in mind, that now the vertical shear in velocities has the opposite sign to the one for the cyclone. Nevertheless, the magnitude of the vertical wave-mean flow interaction increases as well towards the surface due to the stronger vertical shear of the mean flow and leads to a loss of energy, again enhanced at the vicinity of the eddy rim.

The net effect of the wave-mean flow interaction is loss of wave energy from eddy rim to the center and gain outside of similar magnitude (i.e. the eddy would be accelerated inside and decelerated outside) such that both effects appear to cancel each other. We calculate the integrated net effect of wave-mean flow interaction as \[\begin{align} WM = \frac{ \int \int \mathcal{E}\dot{\omega}/\omega \, d ({\boldsymbol{ x }},z) d ({\boldsymbol{ k }},m)dt}{ \int \int |\mathcal{E}\dot{\omega}/\omega| \, d ({\boldsymbol{ x }},z) d ({\boldsymbol{ k }},m)dt} \end{align}\] which is \(WM=0.14\) for both cyclone and anticyclone. There is thus a small net gain of energy by the waves. Wave-mean flow interaction contributes to a total energy gain of \(1.7\%\) relative to the initial energy and integrated over the whole cyclone or anticyclone after the six days of simulation.

3.3 Wave propagation↩︎

The total energy change by wave propagation integrated over wavenumber space and the six days of simulation is shown in Fig. 7 for the horizontal plane close to the surface and the meridional cross-section through the eddy. The effect of wave propagation is as large as the mean-flow interaction and partly counterbalances it, but is more inhomogeneous in its spatial pattern and hence more difficult to understand. Horizontal and vertical propagation (not shown) contribute with similar magnitudes to the total wave propagation effect in Fig. 7. The vertically integrated relative maximal energy loss by wave propagation is \(40.6\%\) relative to the initial energy. Wave propagation contributes to a total energy loss of \(14.1\%\) relative to the initial energy and integrated over the whole cyclone or anticyclone after the six days of simulation.

The intrinsic horizontal group velocity \(\dot{ \boldsymbol{x}}_h - \boldsymbol{U}_h\) depends on the wavenumber magnitudes, but is identical for different wavenumber angles \(\phi\). For the initially isotropic wave field we thus expect no net effect of the horizontal transport \(\dot{ \boldsymbol{x}}_h \mathcal{E}\), except for a radiation out of the model domain, thus a general decrease of the wave energy. In principle, this can be seen in \(\nabla \cdot \dot{ \boldsymbol{x}} \mathcal{E}\), but at the eddy rim, however, there is predominantly wave energy gain by \(\nabla \cdot \dot{ \boldsymbol{x}} \mathcal{E}\) (not shown) which is related to the developing anisotropy at the eddy rim. Outside the eddy and towards the center there is loss of wave energy by the term, with largest magnitudes at the surface and bottom. Towards the eddy center, the effect of \(\nabla \cdot \dot{ \boldsymbol{x} } \mathcal{E}\) vanishes, since there the wave field stays isotropic.

The effect of the vertical propagation \(\partial_z \dot{z} \mathcal{E}\) (not shown) is close to the surface and bottom predominantly negative at the rim, and positive towards the center, while outside the eddy no coherent pattern can be seen. The net effect by the wave propagation is a wave energy loss outside the eddy rim related to wave dissipation by wave refraction (see below), and small gain or vanishing effect towards the eddy center, see Fig. 7. If it is not a loss, the effect of vertical wave propagation tends to be a redistribution of the wave energy over the eddy. In the vertical propagation effect \(\partial_z \dot{z} \mathcal{E}\), however, we see artifacts of the chosen snapshot sampling in time (half daily in our case) which decrease using higher snapshot frequencies (which do not show up in other terms), but which makes it difficult to understand the role of the term. For the case of the anticyclone, the wave propagation is again analogous to the case of the cyclone.

We use horizontal open boundaries such that waves propagating out of the domain across the horizontal boundaries lead to the largest relative energy loss. We assess the role of open boundaries by comparing to an equivalent simulation with horizontally cyclic boundaries. The relative energy loss in time and due to the effect of wave propagation relative to the initial energy integrated over the whole domain is now reduced to approximately \(1\%\). However, the effects of wave-mean flow interaction as well as wave dissipation are almost identical for both simulations. Therefore, the waves that propagate out of the domain in the simulations with horizontal open boundaries do not affect our results owith respect to the wave-mean flow interaction and dissipation.

3.4 Vertical wave refraction and wave dissipation↩︎

In this section, we discuss the role of wave dissipation by vertical refraction in more detail and its role for mixing. The vertical wave refraction is given by \[\label{kz} \dot{k_z} = - \partial_z \omega_{enc} =- \boldsymbol{k}_h\cdot \partial_z \boldsymbol{U}_h -\frac{k_{h}^{2}}{(k_{h}^{2}+k_z^{2})^{2}}\frac{N}{\omega}\partial_z N\tag{12}\] Strong vertical shear \(\partial_z \boldsymbol{U}_h\) and large \(\partial_z N\) can refract waves such that waves eventually leave the wavenumber domain across the vertical wavenumber threshold \(k_z^c\) contributing so to the energy flux \(F_v\) , which we interpret as wave dissipation by vertical refraction. \(F_v\) contributes to a total of \(2.3\%\) of energy loss relative to the initial energy and integrated over the whole cyclone or anticyclone after the six days of simulation.

The effect of the vertical shear \(\partial_z \boldsymbol{U}_h\) dominates the effect of \(\partial_z N\) in the flux \(F_v\). Using single column simulations (excluding dimensions \(x\) and \(y\) but not \(k_x\) and \(k_y\), not shown) with \({\boldsymbol{ U }}=0\) but the same \({\partial}_z N\), we see much less wave dissipation by vertical refraction (in this case critical layer absorption) and much smaller flux \(F_v\) than including \({\partial}_z {\boldsymbol{ U }}_h\). However, the effect of \(\partial_z N\) only also drives a small flux \(F_v\), reducing the total energy by 0.6 % over the six day simulation.

Since the dominant effect by \(\partial_z \boldsymbol{U}_h\) in Eq. (12 ) depends on the horizontal wavenumber vector \(\boldsymbol{k}_h\), it is useful to divide the diagnostics into different wave compartments: north/southward and west/eastward propagating waves. That means we perform a coordinate transformation \(\mathcal{E}(\boldsymbol{k}_h,k_z,z) \to \mathcal{E}(k_z,\omega,\phi,z)\) and integrate over frequency \(\omega\) and the four major cardinal direction, where north is \(45^{\circ} <\phi <135^{\circ}\), east is \(315^{\circ} <\phi <45^{\circ}\), south is \(225^{\circ} <\phi <315 ^{\circ}\) and west is \(135^{\circ} <\phi <225^{\circ}\), where \(\phi\) is the angle with respect to the zonal wavenumber (\(k_x\)) axis. We show the compartment-wise integrated energies at the northern eddy rim as function of depth and \(k_z\) in Fig. 8 (where \(k_z<0\) denotes upward and \(k_z>0\) downward propagating waves) after six days of simulation. At the northern eddy rim is \(\partial_{z} U<0\) and \(\partial_{z} V \to 0\), i.e. \(\dot{k}_z >0\) for eastward propagating waves and \(\dot{k}_z<0\) for westward propagating waves (ignoring the minor effect of \({\partial}_z N\)). Accordingly, Fig. 8 (b) shows that for the compartment of the eastward propagating waves the upward propagating branch with \(k_z<0\) is emptied, and wave energy is transported towards the downward propagating branch with \(k_z>0\), where the flux by \(\dot{k}_z\) across the wavenumber threshold \(+k_z^c\) contributes to \(F_v\). That means that east- and downward propagating waves dissipate by vertical refraction at the northern rim. The opposite is true for the westward propagating waves (Fig. 8 d), where west- and upward propagating waves dissipate by vertical refraction. The net effect is a drainage of energy in east and westward propagating waves, and no such effect in north and southward propagating waves, generating the horizontal anisotropy and the ellipse with major axis oriented toward the eddy center seen in Fig. 5. The same holds for the other positions at the rim, and for the anticylonic eddy accordingly. The effect of \({\partial}_z N\) on \(\dot{k}_z\) is minor but can also contribute to \(F_v\) which can be seen in Fig. 8 (a) and (c) where the effect of \({\partial}_z {\boldsymbol{ U }}_h\) vanishes for the north and southward compartments. Around \(1000 \, \rm m\) depth, the maximal \({\partial}_z N>\) generates the largest magnitude of \(\dot{k}_z <0\), which give rise to a small flux \(F_v\) for upward propagating waves.

Horizontal shear can also lead to wave breaking and dissipation by horizontal refraction . Due to the lateral shear, waves refract according to \(\dot{{\boldsymbol{ k }}}_h =-{\boldsymbol{ k }}_h \cdot \nabla_{{\boldsymbol{ x }}_h} {\boldsymbol{ U }}_h\) and can eventually leave the wavenumber domain across the lateral wavenumber threshold \(k_h^c\) and contribute to the energy flux \(F_h\), which we interpret as wave dissipation by horizontal refraction. This effect is, however, two orders of magnitude smaller than \(F_v\) in all our simulations, therefore we do not discuss it in more detail.

The vertical diffusivity \(\kappa\) is computed from the flux by wave dissipation \(F_v\) from Eq. (8 ), and is shown in Fig. 9 (a,b) after six days simulation. Vertical diffusivities for the case of the cyclonic and anticyclonic eddy are identical. \(\kappa\) is large near the surface with (decreasing) values of \(\kappa\approx 10^{-5}\) to \(10^{-7} \, \rm m^2s^{-1}\) during the simulation period. The diffusivity is maximal at the eddy rim close to the surface where most of the wave breaking takes place as waves get trapped by the mean shear and dissipate. \(\kappa\) is mainly related to \(F_v\) with a minor contribution of \(F_h\). \(\kappa\) is maximal at the eddy rim and decays towards the eddy center and outside of the eddy as the vertical shear of the mean flow decreases.

The simulations show much larger diffusivities of the order of \(\kappa\approx10^{-5} \, \rm m^2s^{-1}\) towards the surface at the eddy rims during the first two days of simulation (not shown). In Fig. 9 (c) we show the evolution of the fraction between the diffusivities at the northern eddy rim against the center. While at the beginning the diffusivities at the rim are overall larger than at the center, where the shallower maxima in \(\kappa\) are predominant (not shown), the lower ones strengthen with time at the rims by more than an order of magnitude. Eventually, as there is no energy fed to the system, no waves are left to break and no more energy is transported outside of the wavenumber domain, and thus diffusivities vanish.

3.5 Comparison to observations↩︎

Available observations appear to show the key feature predicted by our model. A mooring with a downward looking ADCP with \(75\, \rm kHz\) and a sampling period of \(8 \, \rm min\) was deployed for 22 days at the core of our eddy within the REEBUS project. We divide the measured horizontal velocity time series in two series of six days since then the first time series corresponds roughly to measurements at the eddy core and the second at the eddy rim. We compute the horizontal kinetic energy spectral estimate \(\mathcal{E}_h(\omega,k_z)\) for a moving window of four days for five depth-ranges of \(206 \, \rm m\), to finally integrate in vertical wavenumber.

Fig. 10 shows the observed frequency spectra of horizontal kinetic energy densities \(\mathcal{E}_h\) at the eddy center and the rim. They contain both elevated energy levels at tidal frequencies (\(O_1\), \(M_2\) and harmonics), but also all other (wave) frequencies between the local \(f\) and \(N\) are populated with energy, akin to the GM model spectrum. At all such frequencies except for the tidal frequencies, the energy seems indeed to be higher for the eddy center than at the rim. Without more observations from other eddies, it is difficult to say if the difference of the energy level in the observations is statistically significant or just by chance. We made an attempt to proof that the spectral estimate at the eddy center is indeed statistically significantly larger compared to the eddy rim by performing a \(\chi^2\)-goodness of fit² of the fraction of the kinetic energy at the center and rim. With a simple minimum-finding algorithm we compute an expected fraction of \(1.2\) that significantly fits (\(p>0.05\)) our hypothesis.

However, the difference in the observations clearly has a smaller magnitude than in the model. The corresponding spectra from the model – after a coordinate transformation of \(\mathcal{E}_h(\boldsymbol{k})\) to frequency space \(\mathcal{E}_h(\omega,k_z) =\frac{1}{2} \frac{N^2-\omega ^2}{N^2 -f^2} \frac{\omega^2+f^2}{\omega^2}\mathcal{E}(\omega,k_z)\) – at \(t=0\) and \(t=5 \, \rm d\) for the eddy center and rim are also shown in Fig. 10. Since there is no forcing in the model, and since we initialized with the generic GM model spectrum, it is possible that the larger difference in energy levels at eddy rim and center in the model would decrease over time similar to the observations, including appropriate forcing and other components neglected so far in Eq. (3 ). We also compare the spectra for single column and eddy simulations, where the additional horizontal shear in the eddy simulations leads to elevated energies for the center and rim.

We also observe prominent elevated IGW energies at the frequency \(\omega=2.5 \times 10^{-3} \, \rm s^{-1}\) (or \(log_{10}\omega f^{-1}\approx 1.9\)) in the observations both at rim and center. This frequency is in fact the smallest value to which the Brunt-Väisälä frequency converges and therefore suggest that it is a distinct prevalent turning point. Although observations show a peak in energy at this frequency, our simulations do not model this peak, and so we conclude, that it may be caused by other wave processes.

Lastly, we compare the slopes of the observations (where the tidal signal has been filtered out) and model results fitted with a linear regression. The slopes fitted up to the frequency \(\omega=2.5 \times 10^{-3} \, \rm s^{-1}\) show a very good agreement with a mean slope of \(s=-1.2\) (\(R^{2}=0.9\)), where the single column simulations show the biggest difference due to the missing horizontal wave processes. Recently, [32] also approximate a stationary energy spectrum that satisfies a power law with slope \(s=-1\) from simulations based on the ray-tracing equations. The fits show steeper slopes for higher frequencies that range from \(s=-2.1\) to \(-2.4\).

4 Summary and discussion↩︎

A novel numerical approach based on the radiative transfer equation is presented to study the propagation and refraction of IGWs in a meso-scale eddy in the ocean. In its full complexity the model includes six dimensions, three in physical and three in wavenumber space, in addition to the time dimension. The simulations presented here are the first ones, to our knowledge, to comprehensively describe the evolution of the full IGW spectrum inside a meso-scale coherent eddy, within the WKBJ-approximation, and detail the effects of wave propagation, reflection, refraction, interaction with the mean flow, and breaking related to wave dissipation by vertical and horizontal refraction.

Stratification and mean flow for the simulation are motivated from the observations of the coherent cyclonic eddy structure \(CE\_2019\_14N\_25W\) in the Canary Current System, which we think is representative for a typical coherent eddy at mid-latitudes. Mean flow and stratification are assumed to be stationary during a six day simulation in which the wave field evolves and interacts with the background field, but no additional wave forcing is applied. We use the GM model spectrum as the initial condition. A reflection boundary condition at the bottom and surface is used and open boundaries are implemented in the lateral direction, where energy freely flows out of the model domain, but no energy flows in.

In wavenumber space, we use also open boundaries at the wavenumber thresholds \(k_z^c=\pm 1 \, \rm m^{-1}\) vertically and \(k_h^c=\pm 0.5 \, \rm m^{-1}\) horizontally, and interpret the fluxes across \(k_z^c\) and \(k_h^c\) towards smaller wavelength as wave breaking and dissipation by vertical and horizontal refraction, respectively. For given \(\mathcal{E}\) from the GM model spectrum, \(k_z^c\) is large enough for waves to be considered as unstable with Ri\(_{\mathrm{Bulk}}^{-1}\geq 1\) and thus are prone to breaking at \(k_z^c\). The vertical wavenumber threshold \(k_z^c\) is not uniquely given in literature, although a slightly larger vertical wavelength of \(\lambda_c=10 \, \rm m\) (i.e. \(k_z^c=0.63 \, \rm m^{-1}\)) has been reported previously (cf. [24], [25]). We test the impact of the choice \(k_z^c=0.63 \, \rm m^{-1}\) on our results using single column simulations initialised with identical horizontal wavenumber space, wave energy (the GM model), mean flow corresponding to the northern eddy rim and stratification, but different \(k_z^c\). Reducing the vertical wavenumber cutoff from \(k_z^{c}=1 \, \rm m^{-1}\) to \(k_z^{c}=0.63\, \rm m^{-1}\) leads to a \(1.4\) times larger wave energy loss after 6 days, by wave-mean flow interaction and dissipation. Furthermore, the resulting diffusivities are enhanced by an averaged value of \(1.5\) times. Note that only constructively superimposing waves are considered here, thus a condition Ri\(_{\mathrm{Bulk}}^{-1}>4\) may overestimate the breaking. However, Ri\(_{\mathrm{Bulk}}^{-1}\) depends on the total wave energy that we use as the initial condition, i.e. on the parameter \(E_0\) of the GM model spectrum. With larger \(E_0\), Ri\(_{\mathrm{Bulk}}^{-1}\) will increase and \(k_z^c\) will decrease. Note that three times larger values than used here for \(E_0\) are found by [22] from ARGO float observations, i.e. the value for \(E_0\) used here is most likely a low bias.

There is no known analogous treatise for \(k_h^c\) as far as we know. Although we find that fluxes across \(k_h^c\) are of two orders of magnitude smaller than fluxes across \(k_z^c\), this may differ for cases with strong horizontally sheared mean flows. Thus, a proper assessment of horizontal wavenumber thresholds may thus be required in such cases, but here we found the effects to be minor.

The key results from this study are as follows:

• Wave energy significantly decreases at the eddy rim, while the eddy center shows an increase in the wave energy. We made an attempt to find such a difference in wave energy in the available observations, which tend to show also larger wave energy in the eddy center, although with a much weaker difference in the amplitudes.

• Waves gain energy from the mean flow outside of the eddy rim by lateral shear, while from eddy rim towards the center the effect vanishes. By vertical shear, the waves always lose energy to the mean flow, and this effect is largest at the eddy rim and towards the surface. The net effect is to accelerate the eddy inside and to decelerate it outside its rim, with net deceleration overall.

• Lateral shear generates wave energy gain only due to a developing horizontal anisotropy outside the eddy rim. This anisotropy is given by an ellipse with maximal energy towards its center and its major axis oriented towards the eddy center, but slightly tilted clockwise. The ellipse is generated in the first place by the effect of vertical refraction by vertical shear. The tilt is important for the effective wave energy gain and generated by horizontal refraction.

• Vertical shear extracts wave energy by propagation into critical layers, without such dissipation wave-mean flow interaction by vertical shear has no net effect because of vertical reflection.

• The effect of lateral wave propagation is as large as the mean-flow interaction and partly counterbalancing it, but is more inhomogeneous in its pattern and more difficult to understand.

• The effect of vertical propagation tends to be a redistribution of energy, but at the eddy rim a loss of wave energy because of transports towards critical layers.

• Wave breaking related to dissipation by vertical refraction due to vertical shear is two orders of magnitude larger than the effect of horizontal shear (dissipation by horizontal refraction).

• Wave breaking and dissipation by vertical refraction occurs predominantly at the eddy rims close to the surface where vertical shear and changes in \(N\) are largest, and yields overall 2.3 % energy loss over the integration period. The effect of \({\partial}_z {\boldsymbol{ U }}_h\) dominates but \({\partial}_z N\) can also drive a minor energy loss by dissipation by vertical refraction.

• Related vertical diffusivities are also largest at the eddy rims close to the surface and smallest at the eddy center and decay with depth, with maximal values ranging from \(\kappa = 10^{-5} \, \rm m^2s^{-1}\) at the start to \(\kappa = 10^{-7} \, \rm m^2s^{-1}\) at the end of the simulation. These mixing rates appear low but note that they depend on our choice of \(E_0\) of the GM model spectrum we use as initial condition, which is likely to have a low bias.

• The anticyclonic eddy behaves similar to the cyclonic one with respect to the results mentioned above, except for the slight tilt of the horizontal anisotropy which is anticlockwise, opposite to the cyclonic eddy.

Dissipation measurements within the coherent eddy in the Canary Current system taken during the observational campaign will be used in the future to compare with the predicted enhanced mixing at the eddy rim, but are not yet available to us. At the rim of a coherent meso-scale eddy in the Lofoten Basin, [33] indeed find enhanced dissipation rates by turbulent mixing and relate that to internal wave breaking. However, they also find enhanced dissipation at the eddy center, in contradiction to our model simulation. This difference might be related to different energy sources such as the instability of the background shear which [33] also find, or to the many limitations of our model approach.

The following limitations are inherent to our approach. Wave propagation and refraction lead to numerical grid dispersion errors which are getting smaller using finer grid resolution, but fine grids quickly become computationallvery expensive for our model. The simulations presented here are already at the upper end of available resources to us such that using finer grids are not possible at the moment. However, since the spectra and spatial distributions we simulate here remain continuous without large gradients on the small scales where the numerical errors are largest, we regard the effect of numerical grid dispersion errors as minor. The issue could be tested using some kind of self-refining or unstructured grids (triangular domain decomposition may have better properties than the rectangular grid used here).

Because of the assumed stationarity of the mean flow and stratification, the wave field might continuously build up energy at locations where the wave-mean flow interaction is strong despite the large energy sinks by wave-mean flow interaction or wave breaking. The missing implementation of the mean flow’s reaction related to the wave energy changes by wave-mean flow interaction as well as changes in stratification caused by induced mixing due to wave breaking, can lead to an unrealistic representation of the evolution of the wave field. However, without implementation of an interactive mean flow and stratification these effects are difficult to assess and left for future studies. Further, the mean flow is not in balance with the background density field in our simulations, but we regard this effect as minor since the refraction by changing density should be dominated by the background vertical gradient \(N^2(z)\). Both issues could be assessed by comparing with a numerical simulation in physical space of waves interacting with a balanced eddy flow. There, however, other processes as wave-wave interaction and wave scattering at the mean flow and stratification are also present, but which have been neglected in the present model simulations.

Recently, [15] derived (for \(N=const\)) a kinetic equation using the Wigner transform [16] similar to Eq. (3 ), which avoids the WKBJ assumption of slowly varying background conditions. In contrast, this assumption is necessary in our approach using the spectral energy balance Eq. (3 ), but note that we tested the WKBJ assumption in Fig. 2 without finding significant violations. The kinetic equation by [15] could in principle also be used instead of the spectral energy balance Eq. (3 ) (although a realistic \(N(z)\) from observations as used here could not be used). In the kinetic equation by [15] the transport terms in wavenumber space and the energy transfer to or from the mean flow in Eq. (3 ) are replaced by a scattering term involving all other wavenumbers. Practically, this scattering integral over all wavenumbers would drastically slow down a parallel integration since for each time step a sum over the (parallelized) wavenumber space needs to be performed. Conceptually, it is not obvious how to implement dissipation in the kinetic equation by [15] since the fluxes across wavenumber boundaries are missing there and are replaced by the scattering integral. In contrast, Eq. (3 ) is easily parallelized and dissipation is naturally implemented as described above. A possible modification of IWEM to the kinetic equation by [15] is therefore postponed for now to a later study.

Wave forcing by tidal flow at the bottom or oscillatory Ekman pumping at the surface appears simple to implement to our model. Here, we ignore these source terms and consider idealized spin-down simulations only, because of the other limitations which we also currently not resolve. Another possible forcing which is more difficult to include is the secondary wave generation during the wave breaking process. Currently, we assume that all the wave energy propagating to smaller wavelength across the thresholds \(k_z^c\) and \(k_h^c\) is converted entirely to small-scale turbulence. However, a certain amount of the energy will certainly generate other waves at other wavelengths during the breaking process, but this process is largely unknown and needs further attention.

Another limitation of our study is that as initial condition we use the generic GM model spectrum, but inside the eddy the spectrum is changing such that it is likely that in a stationary state the wave field will differ from GM inside the eddy. We found that even in single column simulations (not shown) with vanishing mean flow but depth-dependent stratification, the initial GM model spectrum is not stationary but evolves in time due to vertical wave propagation and refraction. This effect is related to the scaling of the GM model spectrum by \(N\), i.e. by the term \(E(z) = E_0 N(z)/N(0)\) in Eq. (9 ). The GM model spectrum only remains stationary for \(N=const\) in our model. Similar to [34], different power laws of \(N\) for the formulation (larger and smaller then one) have been tested for stationarity of the GM model spectrum, and it was found that other descriptions even increase the non-stationarity. Thus, the description of the GM model spectrum by the scaling \(\mathcal{E}\sim N(z)/N_0\) appears to be a reasonable zero order approximation, but may not be the best representation [1]. In fact, the observed GM model spectrum might be far from stationary such that relatively large wave sources, transports and sinks are needed to explain the observations. It is likely that for the GM model spectrum with its power law of \(-2\) the downscale energy transfers by non-linear wave-wave interactions also play an important role. Such transports are also not considered here and including a full evaluation of the scattering integral as in [35] appears in principle possible but computationally expensive. Some kind of approximation as in [36] might help in this respect and is a subject of future work.

We thank Oliver Bühler and another anonymous reviewer for their comments and suggestions on this manuscript. This paper is a contribution to subprojects W6 and W2 of the Collaborative Research Centre TRR 181 “Energy Transfers in Atmosphere and Ocean” funded by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under project number 274762653. This paper is also associated with the research project REEBUS (Role of Eddies in the Carbon Pump of Eastern Boundary Upwelling Systems) funded by the Bundesministerium für Bildung und Forschung (BMBF, German Federal Ministry of Education and Research). We thank M. Dengler and J. Karstensen from GEOMAR, Kiel for providing the observational data and collaboration in REEBUS.

The novel numerical spectral Internal Wave Energy Model (IWEM) can be accessed via https://github.com/ceden/IWEM.git. The observational data used is taken from the data portal of the Helmholtz Centre for Ocean Research Kiel - GEOMAR, with restricted access.

The Internal Wave Energy Model (IWEM) integrates the six-dimensional Radiative Transfer Equation Eq. (3 ) in time, and simulates the interaction of the IGW field with the background medium. The model describes the propagation and refraction of IGWs in a variable background stratification and mean flow and the energy exchange of the wave field with the mean flow.

Using a staggered C-grid discretisation, wave energy \(E\), stratification \(N\) and mean flow velocities \(U,V,W\) take values at the center of the grid boxes, while zonal, meridional and vertical group velocities are located at the eastern, northern and upper side of the grid boxes. For the wave refraction in \(k_x,k_y,k_z\)-direction the discretisation is analogous. We use an equidistant grid in the meridional direction, a surface refined non-equidistant grid in vertical direction, and also a non-equidistant grid in wavenumber space, with finer resolution for small wavenumber components. For parallelisation, each processor unit solves Eq. (3 ) on a subdomain of the model.

At the surface and bottom, waves are reflected back, so that upward propagating waves with \(k_z<0\) propagate downward with \(k_z>0\) after reflection and vice versa. All other boundaries are open, so that energy can freely flow out of the domain, but no energy can propagate into the domain. For the computation of fluxes a necessary amount of numerical diffusion is introduced, so that energies remain positive definite. All fluxes are therefore computed with a second order advection scheme with superbee flux limiter and forward time-stepping. The numerical model code can be accessed at https://github.com/ceden/IWEM.git.

The observed zonal velocities in the meridional eddy cross section reach down to 1000 m depth. This is insufficient for a complete representation of the evolution of the IGW field, for which we need the mean flow over the whole depth. Moreover, the observed velocities may include IGW signals which we need to remove. To approximate the mean flow to the whole depth and to smooth the observations, a projection on the vertical modes is used. The orthogonal vertical eigenfunctions \(\Phi_n\) result from the vertical eigenvalue problem \[\label{phi} \frac{\mathrm{d} }{\mathrm{d} z}f^{2}N^{-2}\frac{\mathrm{d} \Phi _{n}}{\mathrm{d} z}+f^{2}c_{n}^{-2}\Phi_{n}=0 \;\;with \;\;\frac{\mathrm{d} \Phi_n}{\mathrm{d} z} =0 \;\; at \;\;z=0, -h\tag{13}\] with approximate solution \[\Phi_{n}\approx \sqrt{\frac{2N}{h\widetilde{N}}}cos \left ( \int_{-h}^{z} N/c_{n} \mathrm{d} {z}'\right ), \;\;c_{n}\approx h\widetilde{N}/(n\pi), \;\;\widetilde{N}=\frac{1}{h}\int_{-h}^{0}N\mathrm{d} z\] Velocity \(U\) is then given by \(U=\sum_{n=0}^{\infty }u_{n}\Phi_{n}(z)\) with \(u_{n}=\int_{-h}^{0}U \Phi_{n}\mathrm{d} z\). Since \(U\) is not available for the whole depth down to 4000 m, the amplitudes \(u_n\) for each mode are fitted for the available observations in the upper 1000 m of the ocean. The eigenfunctions are calculated for \(\overline{N}\) and the fitted profiles \(N_1(z)\) and \(N_2(z)\). We carry out the fit with the Python function \(scipy.optimize.curve\_fit\) first for the barotropic mode, subtract then the resulting velocities from the observed ones and follow this procedure iteratively up to the 10\(^{th}\) baroclinic mode, to obtain the mean flow. This is done for each meridional position along the eddy cross section. The barotropic and the first three baroclinic modes are the most powerful ones. The fit of these four modes results in a Pearson coefficient of \(\rho=0.8\), hence the fit and the observed velocities correlate very well. The velocities resulting from the sum of the first 4 modes represent hence most of the mean flow and are used for the eddy simulations.

Following, for example [37] and [31], eddies captured in the tropical Atlantic can be well described by the following idealized structure: exponential decay outside of the eddy and a linear transition in-between eddy rims. The eddy rim is located at a radius \(r=\sqrt{x_{rim}^2+y_{rim}^2}=45 \, {\rm km}\) from the eddy center at \((x,y)_{center}=(120,120)\, {\rm km}\). For the vertical direction, the first baroclinic mode \(\Phi_1(z)\) from Appendix B is used. \[\begin{align} \label{velos} U(x,y,z)&=-0.2\sqrt{L_x^{2}+L_y^{2}} \Phi_1(z) R(x,y) \sin(\tan^{-1} y/x) \\ V(x,y,z)&=0.2\sqrt{L_x^{2}+L_y^{2}} \Phi_1(z) R(x,y) \cos(\tan^{-1} y/x) \\ R(x,y)&=\left\{ \begin{array}{ll} \sqrt{x^2+y^2} & for~ \sqrt{x^2+y^2} < r \\ r \exp\left( -\left( \sqrt{x^2+y^2} -r\right)/2r \right) & for ~ \sqrt{x^2+y^2} > r \end{array} \right. \end{align}\tag{14}\] where \(L_x\) and \(L_y\) denote the horizontal extent of the domain.

References↩︎