## Abstract

Motivated by observations of a strong near-inertial wave signal at the base of the semipermanent anticyclonic Cyprus Eddy during the 2010 Biogeochemistry from the Oligotrophic to the Ultraoligotrophic Mediterranean (BOUM) experiment, a numerical study is performed to investigate the role of near-inertial/eddy interactions in energy transfer out of the mixed layer. A hybrid temporal–spatial decomposition is used to split all variables into three independent components: slow (eddy) and fast (inertial oscillations + waves), which proves useful in understanding the flow dynamics. Through a detailed energy budget analysis, we find that the anticyclonic eddy acts as a catalyst in transferring wind-driven inertial energy to propagating waves. While the eddy sets the spatial scales of the waves, it does not participate in any energy exchange. Near-inertial propagation through the eddy core results in the formation of multiple critical levels with the largest accumulation of wave energy at the base of the eddy. A complementary ray-tracing analysis reveals critical-level formation when the surface-confined inertial rays originate within the negative vorticity region. In contrast, rays originating outside this region focus at the base of the eddy and can propagate at depth.

## 1. Introduction

Understanding the manner by which wind-generated inertial oscillations (IOs) leave the oceanic mixed layer continues to be a topic of interest since it represents an important pathway from a surface energy source to eventual dissipation at depth (Ferrari and Wunsch 2009). Whether most of the dissipation occurs in the mixed layer or in the stratified interior, at critical layers or through shear instability, remains an open question. IOs excited by the wind exhibit constant phases and amplitudes over wide regions, which translates to very slow vertical propagation. Yet, observations indicate that IO energy leaves the mixed layer after a wind event on relatively short time scales, suggesting that smaller horizontal scales must be produced through a dephasing mechanism. D’Asaro (1989, 1995) noted that the meridional wavenumber of inertial currents grows as *βt*, where *β* is the planetary vorticity gradient, and this mechanism can produce sufficiently small horizontal scales to trigger vigorous inertial pumping into the thermocline.

Enhanced inertial pumping also occurs in the presence of mesoscale eddies (van Meurs 1998); in this case, the role of *β* is replaced by that of vorticity gradients. When *β* and eddy anticyclonic vorticity act in concert, inertial pumping is even more pronounced (Balmforth and Young 1999). Kunze (1985) demonstrated the broadening of the inertia–gravity wave band in the presence of anticyclonic vorticity with a ray-tracing analysis. Observations have since confirmed the trapping and accumulation of near-inertial energy in anticyclonic vortices (Kunze et al. 1995; Kunze and Toole 1997; Elipot et al. 2010; Joyce et al. 2013). Kunze and Toole (1997) and Kunze and Boss (1998) found that the trapped energy is dominated by azimuthal and radial mode-1 waves. The dominance of these modes inside a barotropic anticyclonic vortex was demonstrated analytically by Kunze and Boss (1998). Llewellyn Smith (1999) used the near-inertial oscillation approximation to obtain the temporal evolution of trapped near-inertial oscillations in anticyclones. In contrast, near-inertial oscillations originating inside a cyclone are not trapped and rapidly propagate into the deeper ocean (Kunze 1985; Lee and Niiler 1998). In the absence of modulating planetary or vorticity gradients to produce small horizontal scales, very little internal wave energy radiates out of the mixed layer (Voelker et al. 2019).

Further theoretical and numerical studies have elucidated the mechanisms by which the frequency of purely horizontal IOs is modulated in the presence of background mean flow. A shift to superinertial frequencies can arise through advection, refraction by the vorticity or phase dispersion (Young and Ben Jelloul 1997; Klein and Llewellyn Smith 2001; Klein et al. 2004). Numerical simulations of wind-generated IOs in a mesoscale eddy field by Danioux et al. (2008) revealed the excitation at depth of motions with frequency 2*f*, where *f* is the Coriolis frequency. In a follow-up study, Danioux and Klein (2008) demonstrated that 2*f* waves are excited by resonant interactions between IOs and the mesoscale field. A similar explanation was given by Wagner and Young (2016) who showed, by means of an asymptotic analysis, that while near-inertial waves can extract energy from the mean flow, energy transfer to 2*f* waves occurs with the mean flow acting as a catalyst. In a study of near-inertial waves trapped in strong fronts, M. Claret et al. (2018, unpublished manuscript) noted the beginnings of an energy cascade to higher frequencies and proposed a new mechanism by which wind-forced inertial energy is transferred to superinertial frequencies at depth through resonant wave–triad interactions. Kawaguchi et al. (2020) observe the development of multiple inertial-frequency harmonics at the base of an anticyclonic eddy in the Sea of Japan following a storm. They attribute the emergence of these frequency peaks to resonant near-inertial wave interactions.

The theory of critical levels, defined as levels at which a wave’s phase velocity matches the velocity of the mean flow, was originally developed to explain the stalling of wave energy propagation in parallel flows (Drazin and Reid 1981; Maslowe 1986; and references therein), but it was later expanded to more general situations by Kunze (1985) who introduced the concept of a critical layer for nonparallel flows, for example, flows with curvature. In this more general situation, a critical layer is defined as the location where a wavenumber component goes to infinity and the corresponding group velocity to zero. In frontal regions where wave propagation is laterally constrained by the surrounding flow, the critical layer is closer to a critical line and in axisymmetric flows, such as the one considered in our study, the critical layer collapses to a single critical point along the axis of symmetry (see section 5). Here, we shall use the expression “critical layer” in the general sense, as defined above for nonparallel flows, in keeping with widely used terminology in the oceanographic community. In contrast, a caustic or turning point occurs where a group velocity component may approach zero, but the total wavenumber remains finite. Lamb and Shore (1992) described a horizontal critical layer where *N* → 0 laterally. In this case, a horizontal component of the wavenumber goes to infinity. Kunze (1986) reported the first observational evidence for trapping and amplification of near-inertial waves in an anticyclone, although Perkins (1976) had previously observed subinertial frequency near-inertial waves in the presence of an eddy. Kunze et al. (1995) confirmed with fine- and microstructure measurements that turbulent dissipation is a dominant sink for the trapped near-inertial energy, as hypothesized by Lueck and Osborn (1986). Whitt and Thomas (2013) used ray tracing and numerical simulations to illustrate trapping and amplification of near-inertial waves in strongly baroclinic currents with the strongest amplification occurring at sloping critical layers oriented along the tilted isopycnals of the current. Joyce et al. (2013) described observations of near-inertial waves in a Gulf Stream warm-core ring and demonstrated with a numerical simulation that the waves were likely forced by wind-driven oscillations. In contrast to Whitt and Thomas (2013), they did not find evidence of critical-layer formation; instead, the accumulation of near-inertial energy below the eddy coincided with the vertical position at which the waves attain their minimum intrinsic frequency. Danioux et al. (2015) derived a new conservation law from the Young and Ben Jelloul equation (Young and Ben Jelloul 1997) for near-inertial energy evolution and used it to demonstrate that the concentration of near-inertial energy in anticyclones arises as a direct consequence of the shrinking spatial scales of wind-driven IOs by a geostrophic eddy field. Xie and Vanneste (2015) adopted a generalized-Lagrangian-mean formulation to identify a new interaction between IOs and the geostrophic flow that they called “stimulated wave generation.” This interaction draws available potential energy from the mean flow to excite waves. Xie and Vanneste (2015) find that the kinetic energy of IOs is conserved, implying that the wave energy must necessarily come from the mean flow. However, a caveat for the existence of this mechanism is that it is limited to very small Rossby-number regimes.

The present numerical study is motivated by observations of intense near-inertial wave activity above and below the core of a semipermanent anticyclonic eddy in the eastern Mediterranean, during the Biogeochemistry from the Oligotrophic to the Ultraoligotrophic Mediterranean (BOUM) campaign (Cuypers et al. 2012). Its primary objective is to explain the dynamics behind these observations.

A numerical framework well suited for this problem is the *f*-plane Boussinesq model (Winters et al. 2004; Winters and de la Fuente 2012). Our numerical simulation is performed with fine vertical resolution throughout the domain in order to fully capture the development of the wave field inside and its fate below the eddy core.

The rest of the paper is organized as follows. Section 2 provides a brief description of field observations made during the BOUM campaign. Section 3 presents the numerical approach and analysis tools, including a new decomposition splitting the flow into eddy, inertial oscillations (horizontal wavenumber *k*_{h} = 0), and propagating (|*k*_{h}| > 0) inertia–gravity waves, which is used to analyze our simulation output. This is followed by section 4 with results of our numerical simulation, comparisons with the observations using time–depth time series, and a wavelet analysis to identify vertical and horizontal scales present in and below the eddy core. Our results are complemented with energy budget and ray-tracing analyses. We find that under weakly dissipative conditions and in the intermediate-Rossby regime, the generated wave field draws its energy entirely from the kinetic energy of the wind-generated inertial oscillations. This behavior is similar to cases of waves generated inviscidly at a front where there is no energy exchange with the mean flow (Whitt and Thomas 2013; Wagner and Young 2016; M. Claret et al. 2018, unpublished manuscript). It is also consistent with the notion that zero-frequency vortical motions can act as catalysts in facilitating energy transfer among waves of equal frequency (Lelong and Riley 1991; Bartello 1995; Wagner and Young 2016). Discussion of our results in the context of previous studies and future directions are provided in section 5.

## 2. Observations

One of the primary objectives of the BOUM field experiment was to study the possible impact of submesoscale features on the vertical transport of nutrients (Moutin and Prieur 2012). Three different Mediterranean eddies were sampled, but our present focus is on the Cyprus Eddy [labeled eddy C in Cuypers et al. (2012)] for which a strong near-inertial signal was recorded at the eddy base.

The Cyprus Eddy is a semipermanent anticyclonic eddy on the slope of the Eratosthenes Seamount in the eastern Mediterranean Sea. A 3-day station (3.75 inertial periods at 33°36′N) was performed near the eddy center whose position was determined prior to the station [see Moutin et al. (2012) for details]. The eddy was sampled down to 500 m every 3 h, with a CTD providing hydrological properties with 1-m vertical resolution and a 300-KHz lowered-broadband-acoustic-Doppler-current profiler (LADCP) providing currents with 8-m vertical resolution. A single full-depth (923 m) CTD profile was also realized.

In addition to the 3-day station, a (northwest–southeast) cross section of temperature within the eddy was sampled down to 800 m using eight XBT profiles with a vertical resolution of 2 m and horizontal resolution of ≈16 km. Multiple transects across the eddy enabled the characterization of the velocity in the upper 150 m from the ship-mounted ADCP and allowed accurate determination of the eddy vorticity (Moutin and Prieur 2012).

Observations revealed a shallow mixed layer and a largely unstratified core extending from roughly 100 to 350 m, flanked by two regions of sharp density variation and strong vertical shear in the upper layer and at the base of the eddy (Fig. 1a). A clear near-inertial wave pattern with vertical wavelength ≈125 m was observed at depths [350–600] m (Fig. 4a). A horizontal wavelength of ≈100 km was inferred from the linear inertia–gravity wave dispersion relation, using local buoyancy and Coriolis frequencies *N* and *f*, respectively.

A complete description of the BOUM experiment can be found in Moutin et al. (2012) and a detailed analysis of the microstructure measurements that motivated the present numerical study is given in Cuypers et al. (2012).

## 3. Numerical approach

The problem is posed as an initial-value problem in the numerical model “flow_solve” (Winters et al. 2004; Winters and de la Fuente 2012). Model equations are the three-dimensional nonhydrostatic Boussinesq equations on the *f* plane,

where **u**_{h} = {*u*, *υ*} is the horizontal velocity, *w* is the vertical velocity, *f* is the Coriolis frequency, *ρ*_{0} is the constant background density, $\rho \xaf\u2061(z)$ is the mean density, *ρ* is the perturbation density, *p* is the perturbation pressure, and ∇_{h} denotes the horizontal gradient operator.

Dissipation is governed by sixth-order hyperviscosity/hyperdiffusion, with coefficients *ν*_{6} = *κ*_{6}. Boundary conditions are doubly periodic in the horizontal plane and free slip, rigid lid in the vertical direction.

### a. Initial conditions

The initial condition (written in Cartesian coordinates to conform to the numerical model formulation) consists of an axisymmetric, surface-intensified anticyclone in geostrophic and hydrostatic equilibrium (Fig. 1), superimposed with a horizontally uniform zonal wind impulse of amplitude *B*_{0} that is confined to the mixed layer,

Here, the eddy center is at {*x*, *y*} = {0, 0} and the minimum vorticity is 2*A* (*A* < 0). The initial wind impulse decays exponentially at a rate −*γ* from its peak value at the surface. Away from the eddy, the wind impulse excites horizontally uniform inertial oscillations

where *B*(*t* = 0) *= B*_{0}. Creating IOs as a response to an initial condition rather than through continuous forcing avoids perturbing the eddy too strongly away from its equilibrium state and reduces the impact of the low-frequency ageostrophic secondary circulation that dominates in forced cases, obscuring the weaker wave signal of interest. Moreover, since we did not have sufficient information on the structure or intensity of high-frequency winds during the BOUM campaign, and since total energy did not decrease appreciably over the course of our simulation, external forcing was deemed unnecessary.

### b. Flow decomposition

A hybrid temporal/spatial flow decomposition for the flow variables, similar in spirit to the one used by Lelong and Kunze (2013), provides an efficient means of splitting flow variables {*u*, *υ*, *w*, *ρ*, *p*} into eddy (slow), IO, and inertia–gravity wave (IGW) components and has proven helpful in the interpretation of numerical results. The first step in the decomposition rests on the assumption of a time scale separation between the eddy turnover time and the inertial period and relies on a Reynolds decomposition to separate slow (vortex) and fast (IO + IGW) components. This separation is achieved with a running-mean window in time. The width of the averaging window is chosen such that the projection of the fast onto slow component is minimized. Here, a window spanning three inertial periods has been used to obtain a clear delineation between inertial and mean-flow time scales, ensuring that slow and fast components are orthogonal.

Applying the decomposition to *u* yields

for the slow component, where *T* denotes the inertial period. The fast velocity is then defined as the difference of total and slow components

The second step hinges on the fact that IOs lack a characteristic horizontal length scale and separates IO and IGW contributions by applying a horizontal-mean operator over the domain

The IGW velocity is defined as the residual

Remaining variables *υ*, *w*, *ρ*, and *p* are similarly decomposed.

### c. A forced wave equation

The IGW velocity accounts for all of the fast vertical velocity since the IO velocity is purely horizontal. This distinguishing feature can be used to identify wave sources at early times by formulating a wave equation for *w*, derived from the inviscid forms of (3.1),

Here, $N={\u2061(g/\rho 0)[d\rho \xaf\u2061(z)/dz]}1/2$ is the buoyancy frequency. The right-hand side (rhs) is made up of nonlinear terms *F*_{i} and includes all possible interactions but since there is no IGW signal initially, interactions involving *w* will be neglected since they arise at higher order in amplitude. The remaining *F*_{i} terms

include slow/slow and slow/IO contributions. IO/IO terms are identically zero since IO velocities are horizontally uniform. Here, *b* = −*gρ*/*ρ*_{0} denotes buoyancy. Of course, not all solutions of (3.9) exhibit wavelike behavior. Solutions also include vertical velocities associated with nonpropagating, low-frequency balanced ageostrophic circulations excited by slow/slow interactions. The only nonlinear terms on the rhs of (3.9) capable of exciting *propagating waves* are slow/IO interaction terms since they alone project onto both IGW frequency and finite horizontal wavenumber bands. As in Young and Ben Jelloul (1997) and Lelong and Kunze (2013), we anticipate that waves initially excited by nonlinear slow/IO interactions will have horizontal wavelengths matching the horizontal scale of the vortex, and near-inertial frequencies imposed by the mixed layer IOs.

Early time characteristics of the excited wave field can be obtained by computing the slow/IO rhs terms (3.12) with the initial conditions given in section 3a. The *F*_{1} does not contain any slow/IO contributions, a consequence of ∇ ⋅ **u**_{s} = 0, and *F*_{2} vanishes because of the axisymmetry of the eddy. This leaves *F*_{3} as the only nonzero term with a slow/IO projection capable of exciting inertia–gravity waves. The slow/IO contribution to *F*_{3} may be rewritten succinctly as

This expression represents the vertical gradient of the advection of eddy vertical vorticity *ζ*_{s} by inertial oscillations **u**_{IO}, which shifts the effective Coriolis frequency to lower values inside the eddy. It can be related to the refraction term in the Young and Ben Jelloul (1997) wave equation and to the pseudomomentum of Xie and Vanneste (2015). Examination of the spatial structure of (3.13) yields further insight. In particular, the vertical position of maximum forcing can be inferred. The vertical dependence of the interaction is governed by the factor

which achieves a maximum at *z* = 0 and vanishes at *z* = −*γ*/2*β*, where *β* governs the vertical extent of the eddy and *γ* is a measure of mixed layer depth.

Since **u**_{IO} decays rapidly with depth, significant wave forcing at early times can only occur near the surface. As time progresses, the downward propagation of IGWs will introduce near-inertial frequencies at depth and trigger a slow/IGW interaction. This interaction occurs at later times, when the terms involving *w* in (3.9) can no longer be neglected, and it will not be considered here. The theoretical tools developed in this section will be employed in the analysis of our numerical simulation, presented below.

## 4. Numerical results

### Simulation of the Cyprus Eddy

Observations made during the BOUM campaign are used to initialize the simulation. The background density profile (Fig. 1a) was constructed from in situ temperature and salinity measurements down to 850 m taken in the vicinity of Cyprus (Cuypers et al. 2012). Below 850 m, the density profile was extrapolated to 1670-m depth to allow downward wave propagation below the eddy without premature interference from bottom-reflected waves. The computational domain, chosen sufficiently large to allow any generated waves to propagate away from the eddy, has dimensions *L*_{x} × *L*_{y} × *L*_{z}, with *L*_{x} = *L*_{x} = 750 km and *L*_{z} = 1670 m.

Initial vortex strength *A* = −1.32 × 10^{−5} s^{−1}, zonal wind impulse amplitude *B* = 0.2 m s^{−1}, radial, vertical and mixed layer *e*-folding scales *α* = 8 × 10^{−10} m^{−2}, *β* = 8.2 × 10^{−6} m^{−2} and *γ* = 0.01 m^{−1} are chosen to best replicate the observations. The local Coriolis frequency *f* = 8.7 × 10^{−5} s^{−1}.

Observed and initial model vorticity and isopycnal displacements are displayed in Fig. 1. When compared with the Cyprus observations (Fig. 1a), our analytic eddy profile (Fig. 1b) is more symmetric and exhibits a shallower, more stratified core. Density gradients above and below the core are less pronounced in our initial condition, especially at the base of the eddy. As a result, simulated vertical and horizontal shears are weaker than in the observations. Both observed and simulated eddies are characterized by *O*(−*f*/3) vorticity in the core, surrounded by a ring of weaker *O*(+*f*/4) vorticity (Fig. 1).

Vertical cross sections of the fast velocity component through the eddy are shown in Figs. 2a–d at four different times. The initial impulse (not shown) is horizontally uniform zonally and is confined to the top of the domain. Within an inertial period, it gives rise to pure IO. By five inertial periods (Fig. 2a), intensification of the near-inertial signal inside the eddy is visible, and, by 20 inertial periods, downward pumping in the eddy core [the so-called chimney effect described by Lee and Niiler (1998)] is observed (Fig. 2b). Over time, the pumping becomes more vigorous, resulting in accumulation of near-inertial energy at the base of the eddy (Figs. 2c,d). Animations of the flow field clearly show upward phase propagation.

Horizontal cross sections of vertical velocity close to the surface reveal a precessing mode-1 azimuthal perturbation. The spiral waves that develop are part of a symmetry-breaking solution that arises in response to the initial axisymmetric vortex being perturbed. The precessing mode-1 perturbation acts as a nonstationary wave source that produces spiral rather than axisymmetric wave trains (Fig. 3).

A comparison of observed and simulated zonal velocities over three different time periods (19.9–23.1, 35–38.2, 69.5–72.7 days) is shown in Fig. 4. To minimize the impact of differences in observed and model stratifications, the depth is WKB-scaled (the physical depth is also included in each panel). Figures 4b–d demonstrate the progressive tilting of phase lines, indicative of shrinking vertical wavelengths. At early times (Figs. 4a,b), the simulated vertical wavelength is significantly larger than in the observations with WKB-scaled wavelength on the order of 100 m (non-WKB-scaled 400 m) as compared with the observed WKB-scaled vertical wavelength of 50 m (non-WKB-scaled 150 m). By the end of the simulation (Fig. 4d), the vertical wavelength in the region of maximum velocity becomes comparable to the observations.

A frequency–depth plot for the zonal velocity shows a clear bias toward subinertial frequencies in the eddy core (Fig. 5). Both sub and superinertial frequency bands broaden at the base, although there is relatively little energy outside the (0.9*f*, 1.1*f*) frequencies. The eddy manifests itself at low frequencies down to ≈700 m. Interestingly, 2*f* peaks, albeit very weak, are present near the surface and at the base of the eddy. In our simulation, the 2*f* peak at the surface is likely due to the interaction of outward radiating and reentrant near-inertial oscillations in the mixed layer through the periodic boundaries of our computational domain. The 2*f* peak between 500 and 600 m, on the other hand, may be a manifestation of the Danioux et al. (2008) or Wagner and Young (2016) excitation mechanisms for 2*f* waves, hinting at the development of a wave continuum to higher frequencies, as suggested by M. Claret et al. (2018, unpublished manuscript). A 2*f* peak at depth was also present in the observations (Cuypers et al. 2012).

Time–depth plots of the zonal velocity decomposed into slow, IO, and wave components through the eddy core are shown in Fig. 6 for the duration of the simulation. As expected, the slow velocity *u*_{s} is close to zero (Fig. 6a) and would vanish completely, had the eddy not been perturbed by the initial IO impulse. At long times, a faint *u*_{s} signal appears at a depth corresponding to the base of the eddy (Fig. 6a). This is an indication that the slow/fast decomposition slowly breaks down, particularly in regions where nonlinear energy exchanges occur, for example, at the base of the eddy. As anticipated from its definition, the horizontally uniform IO component remains confined to the surface layer (Fig. 6b). The IGW component, on the other hand, intensifies quickly as it travels downward (Fig. 6c). During the first seven inertial periods, upward IGW energy propagation is visible, suggesting a reflection of IO energy at the base of the mixed layer. Byun et al. (2010) observed a similar phenomenon of downward phase (upward energy) propagation within the thermostad of a semipermanent anticyclonic eddy in the southwestern East Sea/Sea of Japan and showed analytically that a reflection occurs when wave energy propagation is controlled by subinertial vertical shears rather than by buoyancy. By 15 inertial periods, a strong downward-propagating energy signal dominates. This signal intensifies over time at the eddy base. The vertical wavelength and the group velocity gradually decrease, as evidenced by the stalling of the signal between 400 and 600 m. The maximum amplitude of the IGWs occurs at that depth (Fig. 6c). Also visible in the IGW velocity at depths between 100 and 400 m are several weaker signals with decreasing vertical scales. These small-scale signals first develop above the dominant downward-energy-propagating signal. Over time, more layers form at progressively shallower depths.

A cross section of the rhs slow/IO term (3.13), computed with simulated **u**_{s} and **u**_{IO} averaged over 10–20 inertial periods through the vortex center plane (Fig. 7), confirms that early time forcing of *w* takes place near the surface, *z* = 0. The fainter signal at depth signals the presence of near-inertial waves at depth. The strong forcing signal at the surface displays short vertical scales comparable to the depth of the mixed layer and horizontal scales on the order of the eddy radius.

#### 1) Spatial scales

Given the localized nature of the wave signal, we perform a wavelet analysis with the Generalized Morse wavelets (Lilly and Olhede 2012; Lilly 2017) to identify IGW wavelengths. Horizontal wavelengths, averaged over days 40–50, are examined at two different depths: 300 m corresponding to a position within the eddy core, and 500-m depth located below the eddy core (Fig. 8). At 300-m depth, dominant horizontal wavelengths are in the range 40–70 km and are radially confined in the eddy core. Smaller wavelengths of *O*(10) km occupy a narrower band about *x* = 0, the eddy center. The waves are relatively weak compared to those below the core. At 500-m depth, waves are much stronger and occupy a broader area about *x* = 0. The dominant wavelengths are 50–100 km. Weaker waves with smaller wavelengths are also visible at the edges of the vortex (Fig. 8b). Overall, the range of horizontal scales of the IGW field is constrained by the vortex dimensions.

Wavelets are also used to find vertical wavelengths in the center of the eddy (Fig. 9). Our analysis reveals the presence of waves with large vertical wavelengths (*λ*_{z} > 250 m) from the surface down to 800-m depth. A weak tongue, centered around 300-m depth, contains smaller vertical wavelengths (*λ*_{z} ≤ 100 m) corresponding to the short horizontal wavelengths seen in Fig. 8a. The presence of this diffuse peak in the core suggests shrinking vertical wavelengths, reminiscent of critical-level behavior. Peaks centered at the surface and around 500 m exhibit wavelengths between 100 and 200 m. The dominant peak at the surface corresponds to the position of maximum forcing predicted from analysis of (3.9).

#### 2) Energetics

The mean fast (IO + IGW) kinetic energy density ⟨*E*_{K}⟩ = ⟨*u*^{2} + *υ*^{2} + *w*^{2}⟩/2 evolves in time as

Here, {*u*, *υ*, *w*, *p*} denote fast variables (omitting subscripts for simplicity); *U*_{x,y,z} and *V*_{x,y,z} denote the (slow) vortex spatial velocity derivatives; h.o.t denotes the higher-order triple-correlation IGW/IGW/IGW terms, which are neglected; and *D*_{ν} denotes the kinetic energy dissipation rate. The ⟨*E*_{K}⟩_{t} is shown in Fig. 10a. Figures 10b–f represent rhs terms of (4.1), with positive or negative regions respectively indicating where (IO + IGW) kinetic energy is increasing or decreasing. The dominant rhs source/sink term is the energy-flux divergence −∇ ⋅ ⟨*pu*⟩, which depletes kinetic energy in the mixed layer and acts as a source beneath the eddy core (Fig. 10d). Buoyancy flux ⟨*bw*⟩ is nearly canceled by the fast/slow vertical shear correlations −(⟨*uw*⟩*U*_{z} + ⟨*υw*⟩*V*_{z}) (Figs. 10b,e). The buoyancy flux is also responsible for conversion of kinetic energy to potential energy. Normal Reynolds-stress deformation terms −(⟨*uu*⟩*U*_{x} + ⟨*υυ*⟩*V*_{y}) (Fig. 10f) act as a weak source and sink and dominate over off-diagonal Reynolds shear-stress terms −⟨*uυ*⟩(*U*_{y} + *V*_{x}) (Fig. 10c), which drain kinetic energy out of the eddy region. The dissipation term is very weak and is not shown.

Mean fast available potential energy ⟨*E*_{P}⟩ = ⟨*b*^{2}⟩/(2*N*^{2}) satisfies

where *B* = −*gρ*_{s}/*ρ*_{0} and *N*^{2} = *B*_{z}. Again, h.o.t denotes neglected higher-order correlations and *D*_{κ} is the potential energy dissipation. The sum of the middle two rhs terms in (4.2) nearly balances the vertical buoyancy flux except near the surface, as shown in Fig. 11. This implies that *u* ⋅ ∇*B* = 0 and that, away from the surface, fluid parcel motion lies primarily along isopycnals. This behavior is consistent with the statement that waves of minimum frequency *σ*_{min} have wavelength aspect ratio −*k*_{x}/*k*_{z} = *s*, where *s* is the slope of the isopycnals (Whitt and Thomas 2013; Joyce et al. 2013). Near the surface, there is a net transfer from kinetic energy to potential energy (Figs. 10a and 11a). Initially, fast energy is entirely kinetic since IOs do not have potential energy. The conversion of fast IO kinetic energy to potential energy is consistent with the generation of IGWs.

We now turn to examination of slow, IO, and IGW energies. While expressing any flow variable *ϕ* = {*u*, *υ*, *w*, *b*, *p*} as the sum of slow, IO, and IGW components is exact in the sense that

the decomposition is not orthogonal unless the cross-term contributions [underbraced terms in (4.3)] to the time-averaged variance ⟨*ϕ*^{2}⟩

are negligible. In practice, this requires extracting the slow component with a running mean defined over a temporal window larger than an inertial period and shorter than the simulation duration, to ensure a time scale separation between eddy and fast IO/wave components. Minimization of the three cross terms in (4.3), taking into consideration the constraints outlined above, was used to determine the optimal window of three inertial periods for slow/fast variable separation in our analysis.

The change in energy for the slow, IO, and IGW components as a function of depth is shown in Fig. 12a. The eddy energy loss is negligible throughout the water column. IO energy at the surface provides energy for the waves. Significant wave energy remains in the surface layer but some makes its way through the eddy core to form a maximum at the base of the eddy centered at about 500 m. Depth-integrated normalized energies over the course of 72 inertial periods are shown in Fig. 12b. The sum of the three normalized component energies is equal to the sum of the nondecomposed kinetic and available potential energies, indicating that the decomposition is orthogonal. The simulation is nearly inviscid, with total energy decreasing by less than 3%. Wave energy grows entirely at the expense of IO energy. The good agreement between the sum of slow, IO, wave energies, and the total energy computed with nondecomposed variables confirms the utility of the decomposition. This analysis is complemented in the next section with a ray-tracing approach.

## 5. Ray-tracing analysis

To better interpret the simulation results and the observations, we use ray tracing to compute the characteristics of near-inertial oscillations (NIOs) interacting with the eddy. Since pure inertial oscillations do not propagate, the ray tracing is initialized with a slightly *off-inertial* signal, hence the need to introduce new terminology.

Ray tracing makes use of several important assumptions. First, we use the WKB approximation, which assumes that the NIOs evolve in a slowly varying medium. This assumption is not fulfilled in our case since the eddy and the NIOs have similar horizontal *O*(10 km) and vertical *O*(100 m) scales. Yet, as shown by several studies, the WKB approximation remains useful outside its strict regime of validity and can shed some insight into the physics of our problem (Kunze 1985; Sartelet 2003; Chavanne et al. 2010; Sheen et al. 2015; Cuypers et al. 2017). We use the eikonal equations derived by Kunze (1985) that take into account the interaction of NIOs with a mean flow in the limit of weak Ro and high geostrophic Richardson number Ri_{g} = *N*^{2}/[(∂*u*_{s}/∂*z*)^{2} + (∂*u*_{s}/∂*z*)^{2}], conditions satisfied by the Cyprus Eddy (Ro ≃ 0.2; Ri_{g} ≥ 10). Since we consider the propagation of a single NIO at a time, wave–wave interactions are not considered. Moreover, we assume that the eddy field is unaffected by the NIO propagation. The energy budget shown in the previous section demonstrates that there is indeed very little energy variation of the slow eddy component. In the following discussion, we assume that the eddy is represented by (3.2) (minus the IO contribution).

The eikonal equations that determine the position **x**(*x*, *y*, *z*) and wavenumbers **k**(*k*, *l*, *m*) of the NIO energy rays (Lighthill 1965) are

with *ω*_{o} being the observed (Eulerian) frequency that is conserved along the rays in a stationary background; *ω*_{o} is related to *ω*_{i}, the NIO intrinsic frequency, through the Doppler shift relationship *ω*_{i} = *ω*_{o} − **k** ⋅ **U**, where **U** = (*u*_{s}, *υ*_{s}) is the velocity of the background (eddy) field; *ω*_{i} satisfies the extended dispersion relationship (Kunze 1985) for NIOs in the presence of a nonhomogeneous background

Here, *f*_{e} = *f + ζ*/2 is the effective inertial frequency taking into account the additional background rotation induced by the eddy, and

is the effective buoyancy frequency felt by the NIOs, taking into account the (weak) background baroclinicity. The minimum frequency above which propagation is possible is *ω*_{i} = *f*_{e}. This frequency is reached when $Ne2=0$, that is when the ray trajectory (and particle orbit) parallels isopycnals as noted by Kunze (1985) and Whitt and Thomas (2013). Therefore, NIO rays are expected to tilt following the isopycnal slope, as observed in our simulations and confirmed using energetics arguments in the previous section. This behavior was also noted by Joyce et al. (2013). Yet, it is worth noting as well that $Ne2$ can sometimes become negative in regions where isopycnals are steep, allowing propagation for *ω*_{i} < *f*_{e}.

To illustrate the ray propagation in space, we first consider two typical cases. In the first case, NIOs are initialized outside the negative vorticity region of the eddy with a slightly superinertial observed frequency *f*_{e} ≤ *f* < *ω*_{o} in order to ensure propagation. We choose *ω*_{o} = 1.04*f*, which is in the range of extension of the inertial peak near the surface (Fig. 5). In the second case, we investigate the potential trapping of NIOs in the negative vorticity core of the eddy with NIOs initialized within the core with a slightly subinertial frequency *f*_{e} < *ω*_{o} < *f*. We fix *ω*_{o} = 0.95*f*, again falling within the range of extension of the inertial spectral peak. Note that in each case, we choose to fix the observed rather than the intrinsic frequency because it is the quantity characterized from the usual Eulerian spectrum and it is conserved along the ray path. We initialize each ray with a horizontal wavelength of 100 km based on the dominant wavelength found in the simulations near the surface (Fig. 3) and we fix an initial propagation direction along *x* (*l* = 0). The superinertial rays initialized outside of the eddy have an initial position (*x*_{0}, *y*_{0}) along the *x* = −60 km line with *y*_{0} varying between −70 and 70 km, relative to the vortex center. As shown in Fig. 13a, when *y*_{0} is negative the intrinsic frequency experiences a red shift −**k** ⋅ **U**, *ω*_{i} < *ω*_{o}. The NIO then behaves more inertially, as evidenced by initial trajectories closer to the horizontal. In contrast, rays with initially positive *y*_{0} experience a red shift of *ω*_{i} and behave more superinertially. As a consequence, their propagation toward the bottom is more rapid. NIO ray trajectories are also clearly influenced by the horizontal advection of the eddy (represented somewhat arbitrarily in Fig. 13 by the isosurface *ζ* = −0.04*f*) and a part of the ray initialized close to *y*_{0} = 0 experiences a spiraling trajectory reminiscent of the spiral pattern observed in the simulations (Fig. 2). Those rays also focus horizontally near the base of the eddy around *z* = 400 m, suggesting an increase of NIO energy density at that depth. Near the base of the eddy, a clear decrease of the vertical wavelength down to *λ*_{z} = 400 m is also observed. Yet, all the rays eventually escape the eddy and see a rapid increase of their wavelength as they propagate outside the eddy.

For the second case, subinertial rays initialized within the eddy all remain trapped within the negative vorticity core. Indeed, as they propagate outward and toward the region of higher *f*_{e} they reach a turning point where *ω*_{i} = *f*_{e} and then reflect horizontally. As they propagate downward, the region in which propagation is possible vanishes and the NIOs spiral downward toward a critical point where their vertical wavelength collapses. Note that this point is never reached as it requires an infinite integration time. In real observations, wave breaking will also occur when some critical vertical wavelength is reached, and this process may rather manifest itself in observations as “layers” (Kunze et al. 1995). This process has been previously described by Kunze (1985) and Whitt and Thomas (2013) in the case of a baroclinic front, studied in numerical simulations by Lee and Niiler (1998) and observed by Kunze et al. (1995) in an anticyclonic eddy. It is worth noting that here, in fact, we observe some limited regions where $Ne2$ becomes slightly negative, allowing limited propagation in the region where *ω*_{i} < *f*_{e}, the so-called anomalously low frequency range. For this regime, we observe that ray characteristics have the same slope sign upon horizontal reflection in contrast to what is observed in the usual *f*_{e} < *ω*_{i} range. In the latter range, an incident ray propagating downward results in an upward-propagating reflected ray. This behavior is similar to what was described by Whitt and Thomas (2013) for this frequency range although their simulations differ in that the set of equations they used is not restricted to weak baroclinicity.

The initial position and resulting Doppler shift also have a strong impact on the trapped ray trajectory. Indeed, NIOs with an initial position *y*_{0} < 0 have a blue-shifted intrinsic frequency, propagate faster at depth, and establish a “deep” critical layer around 450 m (Fig. 13b) whereas the rays with *y*_{0} > 0 have red-shifted intrinsic frequency, have slower propagation at depth and establish a critical layer at shallower 200–300-m depth. The formation of multiple critical layers is present in the simulations where several layers of short vertical NIOs form between 250- and 500-m depth (Fig. 6c).

To investigate more systematically the behavior of the subinertial rays trapped within the eddy, we computed the maximum depth reached by NIO rays in the space parameter of the initial conditions *ω*_{0} ∈ [0.8*f* 1.0*f*] and *y*_{0} ∈ [−5050] km. The initial horizontal wavelength and initial *x* position are kept constant at *λ*_{h} = 100 km and *x*_{0} = 0 km. For the sake of simplicity, we kept the initial NIO propagation direction along *x* > 0 but the results can of course be generalized to any radial direction owing to the azimuthal symmetry of the eddy field. The trajectories were integrated over a large number (300) of inertial periods. This bound on time integration is somewhat arbitrary but actually, the rays did not show any significant propagation after approximately 100 inertial periods. We observe first that propagation is inhibited in a significant fraction of parameter space, that is the waves are evanescent because their initial intrinsic frequency is below the effective inertial frequency (*ω*_{i} = *ω*_{o} − **k** ⋅ **U** < *f*_{e}) (Fig. 14). Naturally, the possible range of initial positions for which propagation is possible decreases continuously with the initial observed frequency *ω*_{o}. The effective inertial frequency is most subinertial at the center of the eddy (*x* = 0, *y* = 0), favoring propagation. In the meantime, the Doppler shift effect is positive and favors propagation for negative initial positions along *y*_{0} and tends to inhibit it for positive initial positions along *y*_{0}. The initial position *y*_{0} allowing the widest range of initial observed frequency (for a NIO propagating along the *x* axis) is therefore shifted from the center of the eddy toward *y*_{0} = −15 km where maximum NIO flux may be expected.

To illustrate the propagation of subinertial NIO energy with depth in the numerical simulation, we show in Fig. 15 the time–depth plot of the modulus of wave velocity |*u*_{IGW}| at the center of the eddy, bandpass filtered in two frequency ranges [0.9*f* 1.05*f*] and [0.8*f* 0.9*f*] (limiting frequencies indicated by vertical lines on Fig. 5). We compare this propagation with the ray-tracing trajectories in (*z*, *t*) space, originating from *y*_{0} = −15 km. A large fraction of the trajectories originating from an initial position *y*_{0} = −15 km stall between 400- and 600-m depth, which corresponds to the region of maximum velocity amplitude in the simulation. We also observe a decreasing trend of the rays’ vertical wavelength with time (indicated with color) as well as an increasing wavelength with stalling depth showing a qualitative agreement with the simulations runs. For instance, the rays stalling around 550-m depth (maximum of wave velocity amplitude in the simulation) show a decrease of their wavelength from 300 to 100 m between days 20 and 70, in good agreement with the simulations where a decrease in the same range is present (Fig. 4 and section 4a). These results strongly suggest that a critical layer formation at the base of the eddy is taking place with a progressive decrease of the vertical wavelength and group velocity over time. It can also be noted that the rays with the lower initial *ω*_{0} tend to be trapped at shallower depth and generally show smaller wavelengths of a few tens of meters. There is much weaker energy in the lower frequency range [0.8*f* 0.9*f*] shown in Fig. 15a, (note the different scaling of the gray color bar). In that range, the maximum in the velocity amplitude appears at a constant depth between 200- and 300-m depth. The ray tracing is consistent with a rapid stalling of downward propagation around 200-m depth and continuous decrease of the wavelength with time.

Ray tracing illustrates the fact that the three-dimensional character of the flow is of fundamental importance because it strongly impacts the Doppler shift and, eventually, the ray trajectories. This was previously pointed out by Kunze (1985) who noted the importance of the angle of attack for rays approaching a two-dimensional front. The a priori determination of the surfaces on which critical layers will form, given only the initial intrinsic frequency and the stationary background, is not possible except for the very specific case where **k** ⋅ **U** = 0 such that *ω*_{i} = *ω*_{0}, which imposes a fixed initial direction of propagation normal to the flow. Using the reasonable assumption of a stationary background flow, ray tracing has reproduced many of the features of the numerical simulations, including identification of the depths at which NIO energy increases and the associated NIO vertical wavelengths time evolution.

## 6. Discussion and conclusions

Our numerical simulation has reproduced the salient features of observations made in the vicinity of the Cyprus Eddy during the BOUM campaign. Observed and simulated vertical wavelengths become comparable at large times. We have used an analytic axisymmetric Gaussian velocity profile to model the Cyprus Eddy, and while observed vorticity and mean stratification were used to constrain the analytical model, the vertical structure of the Cyprus Eddy could not be exactly replicated with our analytic profile. Nonetheless, our simulation has demonstrated the development of a rich spectrum of vertical scales inside and below the eddy.

We extended the standard Reynolds-averaging decomposition into slow/fast components by splitting the fast flow into nonpropagating inertial oscillations (IOs) and propagating inertia–gravity waves (IGWs) and this proved useful in interpreting simulation results. In particular, we have demonstrated that in this weakly dissipative regime, the wind-driven inertial oscillations provide the energy source for the IGWs while the eddy acts as a catalyzer: It does not exchange energy with IO and IGW components but nonetheless plays an essential role by providing inertial oscillations with finite horizontal scales needed to trigger vertical propagation. The lack of energy exchange between the fast (IO + IGW) component and the eddy is consistent with Asselin and Young (2020) who find that fast (IO + IGW) feedback onto the eddy field depends on the wind strength but remains weak even in relatively strong wind. In our regime, where wind forcing produces IO amplitudes of *O*(10) cm s^{−1}, these authors report no discernible energy exchange between the two components. The eddy may also be instrumental in the presence of a 2*f* peak at depth. A ray-tracing analysis has identified vertical wavelengths of the excited internal waves and demonstrated the occurrence of two distinct dynamical mechanisms depending on whether the rays originate inside or outside the eddy. The formation of critical layers at different depths inside the eddy core is due to rays originating at the surface and inside the eddy. In contrast, rays that are initially outside the eddy core are not associated with critical layers, focus at the base of the eddy, and eventually propagate out. The IOs in our simulation are horizontally uniform and both the presence of critical levels inside and at the base of the eddy core as well as focusing at the base are expected. The accumulation of near-inertial energy at the bottom of the eddy is indicative of a critical level. It is unclear whether this behavior was present in the observations, but it is not possible to draw any conclusions with data confined to a three-day window and without any knowledge of the wind event that produced the IOs. As we have seen, critical levels develop over very long times. Monthlong observations in an anticyclonic eddy in the Sea of Japan show near-inertial energy accumulation at several critical depths, in agreement with ray-tracing analysis (Kawaguchi et al. 2020). In more realistic conditions, it is highly likely that shear instability would occur and lead to turbulent dissipation as vertical wavelengths shrink. Near-inertial accumulation can be explained by the existence of a vertical position below which propagating waves become evanescent. Mathematically, this level delimits hyperbolic (characteristic of propagating waves) from elliptic (evanescent, nonpropagating) solutions and corresponds to the position of the minimum allowed wave frequency, as derived in appendix B of Joyce et al. (2013). This position is independent of the strength of the inertial oscillations and is entirely determined by geometry of the vortex and the corresponding geostrophic Richardson number, as predicted by our ray-tracing analysis. Our simulations show that critical levels develop inside the eddy core and, over time, appear at progressively shallower depths. This behavior is also exhibited in our ray-tracing analysis, which demonstrates that shallow critical levels develop over longer time scales. The ray tracing also illustrates the fact that the three-dimensional character of the flow is of fundamental importance because it strongly impacts the Doppler shift and eventually the ray trajectories, as previously pointed out by Kunze (1985).

Future directions will address the eventual dissipation of IGW energy, which was absent in the present study. Indeed, critical levels within the core may provide preferential sites for energy dissipation, but an accurate assessment of dissipative impacts will require implementing forcing and sponge layers, longer simulations and even finer vertical resolution, and seeding the flow with a broadband spectrum of internal waves. Other potential areas of research will explore the impact of different wind configurations (e.g., variability in space and time), the possible existence of a wave energy cascade below the eddy, as first reported by M. Claret et al. (2018, unpublished manuscript) and possibly seen in our simulation, and cases where the eddy is embedded in a time-evolving background flow.

## Acknowledgments

This study was initiated during author M.-P. Lelong’s visit to LOCEAN and was funded by the LABEX-IPSL Invited Scientist Program. Discussions with Eric Kunze, Mariona Claret, Jeffrey Early, and Kraig Winters, as well as pertinent comments by an anonymous reviewer, greatly improved the paper. We also thank Kraig Winters for usage and generous support of his Boussinesq code “flow_solve” (Winters and de la Fuente 2012). Additional funding for Lelong was provided by NSF-OCE1031002 and NSF-OCE1851572.

## REFERENCES

*J. Phys. Oceanogr.*

*J. Mar. Res.*

*J. Atmos. Sci.*

*Geophys. Res. Lett.*

*J. Phys. Oceanogr.*

*Biogeosciences*

*Geophys. Res. Lett.*

*f*frequency

*J. Phys. Oceanogr.*

*J. Phys. Oceanogr.*

*J. Fluid Mech.*

*β*effect

*J. Geophys. Res.*

*J. Phys. Oceanogr.*

*Hydrodynamic Stability*

*J. Geophys. Res.*

*Annu. Rev. Fluid Mech.*

*J. Geophys. Res. Oceans*

*Prog. Oceanogr.*

*J. Mar. Res.*

*Quart. J. Roy. Meteor. Soc.*

*J. Phys. Oceanogr.*

*J. Phys. Oceanogr.*

*J. Phys. Oceanogr.*

*J. Phys. Oceanogr.*

*J. Phys. Oceanogr.*

*J. Phys. Oceanogr.*

*J. Geophys. Res.*

*J. Fluid Mech.*

*J. Fluid Mech.*

*IMA J. Appl. Math.*

*IEEE Trans. Signal Process.*

*J. Phys. Oceanogr.*

*J. Geophys. Res.*

*Annu. Rev. Fluid Mech.*

*Biogeosciences*

*Biogeosciences*

*Deep-Sea Res. Oceanogr. Abstr.*

*J. Atmos. Sci.*

*Geophys. Res. Lett.*

*J. Phys. Oceanogr.*

*Dyn. Atmos. Oceans*

*J. Fluid Mech.*

*J. Phys. Oceanogr.*

*Ocean Modell.*

*J. Atmos. Oceanic Technol.*

*J. Fluid Mech.*

*J. Mar. Res.*