# Interpreting Observed Interactions between
Near-Inertial Waves and Mesoscale Eddies

###### Abstract

The evolution of wind-generated near-inertial waves (NIWs) is known to be
influenced by the mesoscale eddy field, yet it remains a challenge to
disentangle the effects of this interaction in observations. NIWs are often
modeled using a slab mixed-layer model with no horizontal structure. Here, the
theoretical model of Young and Ben Jelloul, which describes the evolution of
NIWs in the presence of a slowly-evolving mesoscale eddy field, is compared to
observations from a mooring array in the Northeast Atlantic Ocean. The model
captures the evolution of both the observed NIW amplitude and phase much more
accurately than the slab mixed-layer model, and it allows attributing the
evolution to specific physical processes. The model reveals that differences
in NIW amplitude across the mooring array are caused by refraction of NIWs
into anticyclones. Advection and wave dispersion also make non-negligible
contributions to the observed wave evolution. Stimulated generation, a process
by which mesoscale kinetic energy acts as a source of NIW potential energy, is
estimated to be 20 in the region of the mooring array. This is two orders of
magnitude smaller than the global average input to mesoscale kinetic energy
and likely not an important contribution to the mesoscale kinetic energy
budget in this region.

##  1 Introduction

Near-inertial waves (NIWs), internal waves with a frequency close to the
inertial frequency $f$, are resonantly excited by atmospheric winds exerting a
stress on the ocean’s surface. It has long been recognized that these waves
can interact with mesoscale eddies and that this interaction may be important
in the life cycle of wind-generated NIWs. Observational evidence of
NIW–mesoscale interactions is accumulating, but interpreting the observed NIW
evolution in the presence of mesoscale eddies remains challenging. Here, we
employ the theoretical framework of Young and Ben Jelloul (1997, from hereon
YBJ) to identify and interpret such interactions in mooring observations in
the Northeast Atlantic Ocean.

NIWs are associated with a large vertical shear, which can result in shear
instabilities and vertical mixing (for a review, see Alford et al., 2016).
These shear instabilities are a key mechanism by which atmospheric storms can
cause the surface mixed layer to deepen. Jochum et al. (2013) showed that an
improved representation of NIWs in a climate model led to a deepening of the
mixed layer on average, which in turn resulted in significant changes in sea
surface temperatures, winds, and precipitation.

The amount of NIW shear in the upper ocean, and as a consequence the amount of
mixed layer deepening, depends on both the energy input into the NIW band and
how rapidly NIWs propagate to depth. The vertical propagation originally
presented a paradox. The vertical group velocity of NIWs varies as
$\kappa^{2}$, where $\kappa$ is the horizontal wavenumber. The atmospheric
storms that generate NIWs are typically $O$(1000) in size and generate NIWs
with a similar scale. Estimates of the vertical group velocity based on this
horizontal scale are much slower than observed (D’Asaro et al., 1995). A
reduction in the horizontal scale of the NIWs is required to obtain a group
velocity that matches observations. There are two possible mechanisms by which
this is thought to occur: (i) meridional variations in $f$ cause an increase
in the meridional wavenumber ($\beta$ refraction; Gill, 1984), or (ii)
interactions with the mesoscale circulation can imprint structure onto the
wave field at the scale of mesoscale eddies ($\zeta$ refraction; YBJ; Kunze,
1985). The latter process causes concentration of NIW energy into anticyclones
and subsequent downward propagation in “drainpipes” (e.g., Asselin and Young,
2020).

There is a long history of NIW observations from in situ measurements. The
Ocean Storms Experiment (D’Asaro et al., 1995) was a groundbreaking
observational campaign to study the life cycle of wind-generated NIWs. The
experiment tracked the response of the surface ocean to strong wind forcing in
a region with weak mesoscale eddies. It was found that $\beta$ refraction was
the dominant process driving the observed evolution of the meridional
wavenumber of the waves. The process of $\zeta$ refraction was found to have
much less of an effect than expected (D’Asaro, 1995), which YBJ argued was due
to strong dispersion in a weak eddy field. More recently, Thomas et al. (2020)
used ship-based observations of NIWs in a strong dipole vortex to estimate the
rate of change in time of the NIW horizontal wavenumbers. Here, the expected
change in horizontal wavenumbers from $\zeta$-refraction was consistent with
the data. The observations further showed a NIW beam at depth, indicating that
the scale reduction had allowed the NIWs to propagate vertically out of the
mixed layer. Essink et al. (2022) studied typhoon-forced NIWs in the Kuroshio
Current. They observed the trapping of NIWs in a region of anticyclonic
vorticity, followed by downward propagation. They also measured the vertical
structure of turbulence and showed how this was influenced by NIW dynamics. Yu
et al. (2022) investigated the interaction of NIWs and mesoscale eddies in
observational data of NIWs from the Northeast Atlantic (the same data as we
will use below). They showed that regions of elevated NIW kinetic energy are
statistically associated with mesoscale anticyclones and that the submesoscale
vorticity exerted little control on the horizontal concentration of NIWs and
the subsequent propagation to depth.

It has also been proposed that NIWs not only react to the presence of
mesoscale eddies but feed back on the eddies and affect their evolution.
Approximately 80$\%$ of the ocean’s kinetic energy exists as mesoscale motions
(Ferrari and Wunsch, 2009). The geostrophic constraint on mesoscale eddies
traps energy at large scales, and it is not entirely clear how the energy
input into mesoscale motion is balanced by dissipation (Müller et al., 2011).
A number of mechanisms by which mesoscale eddies lose energy are known,
including dissipation in bottom boundary layers (e.g., Arbic and Flierl,
2004), the generation of dissipative lee waves (e.g., Nikurashin et al.,
2013), energy loss near western boundaries (e.g., Zhai et al., 2010), and the
top drag arising from the current dependence of the wind stress (e.g., Dewar
and Flierl, 1987; Renault et al., 2016). The extraction of energy from
mesoscale eddies by NIWs presents another possibility (Xie and Vanneste, 2015;
Rocha et al., 2018). Given the great importance of the mesoscale to the
transport of heat and carbon (e.g., Jayne and Marotzke, 2002; Gnanadesikan et
al., 2015), even small changes (see discussion in Asselin and Young, 2020) to
the mesoscale eddy field caused by NIWs may be significant to the impact of
the ocean on climate.

To understand the propagation of NIWs and by extension the role they play in
upper-ocean mixing, we need to understand the dynamics governing their
evolution. Slab mixed layer models (Pollard and Millard, 1970, from hereon PM)
are a commonly used tool to model NIW evolution (e.g. D’Asaro, 1985; Alford,
2001, 2020; Guan et al., 2014). PM showed that a slab model can reproduce key
features of the NIW evolution observed by moorings. With the horizontal NIW
velocity denoted by $(u,v)^{\mathrm{T}}$, the PM model can be written as:

$\displaystyle\partial_{t}u-fv=-ru+\frac{\tau^{x}}{\rho_{w}H_{m}},$ (1a)
$\displaystyle\partial_{t}v+fu=-rv+\frac{\tau^{y}}{\rho_{w}H_{m}},$ (1b)

where $r^{-1}$ is a decay timescale, $\rho_{w}$ is the density of water,
$H_{m}$ is the mixed layer depth, and $\bm{\tau}$ is the wind stress. This
model does not explicitly represent any physical process affecting the
evolution of NIWs other than the wind forcing and Coriolis effect. All other
processes are subsumed in the linear drag, which can be thought of as a
parameterization primarily for the propagation of NIWs out of the mixed layer.
Observations indicate that this process is slow relative to the inertial
frequency (e.g., D’Asaro et al., 1995), and so we require $r\ll f$. Other than
this restriction, $r$ is a tunable parameter.

Despite its successes in capturing some aspects of the observed NIW signal,
the PM model cannot help us to understand propagation of NIWs out of the mixed
layer because by assumption the NIW field remains horizontally uniform. The
model captures neither $\beta$ refraction nor the interaction with mesoscale
eddies.

A more general framework to understand the evolution of NIWs was devised by
YBJ by assuming a time scale separation between the fast waves and the slowly
evolving mesoscale flow. In the YBJ framework, the horizontal NIW velocity is
first complexified (written as $u+iv$). Since NIWs have a frequency close to
$f$, it is convenient to write this complexified velocity as
$u+iv=e^{-ift}\partial_{z}M(x,y,z,t)$. The function $\partial_{z}M$ describes
the slow evolution of the envelope that modulates the NIW phase and amplitude.
The YBJ equation describes how $\partial_{z}M$ evolves in the presence of
geostrophic mesoscale eddies. On an $f$-plane, the equation reads

$\partial_{zzt}M+J(\psi,\partial_{zz}M)+\frac{iN^{2}}{2f}\nabla^{2}M+\frac{i\zeta}{2}\partial_{zz}M=\partial_{zz}\mathcal{F}-\nu\nabla^{4}\partial_{zz}M,$
(2)

where $\psi$ is the geostrophic streamfunction, $\zeta=\nabla^{2}\psi$ is the
geostrophic vorticity, $N^{2}$ is the stratification, $\nabla^{2}$ is the
horizontal Laplacian operator, $\mathcal{F}$ is a forcing term that represents
the momentum flux due to the surface wind stress, and $\nu$ is a hyper-
diffusivity included for numerical stability (see Section 3). In (2) and
throughout this paper, we assume the streamfunction to be barotropic, although
baroclinicity in the geostrophic eddy field can also be taken into account.
Unlike the PM model, the YBJ equation does not have a tuneable parameter.

The second term in (2) represents advection of NIWs by the mesoscale flow. The
third term represents changes in the NIW field due to dispersion. The presence
of the dispersion term means that an initially localized wave packet will
spread out as time progresses. The fourth term is responsible for the process
of $\zeta$ refraction. This term sets into motion the imprinting of mesoscale
structure onto an initially horizontally uniform wave field. This
$\zeta$-refraction shifts the phase of the NIWs, which we can see by
neglecting all other terms in the YBJ equation:

$\partial_{zzt}M=-\frac{i\zeta}{2}\partial_{zz}M.$ (3)

Assuming a steady vorticity field yields solutions of the form
$\partial_{zz}M=C(x,y,z)e^{-i\zeta t/2}$. Spatial heterogeneities in the
mesoscale vorticity will result in spatial heterogeneities in the NIW phase.
The dispersion term in the full YBJ equation (2) then acts on these phase
gradients and fluxes energy into anticyclones and out of cyclones. This can be
seen from the YBJ kinetic energy budget (Appendix A):

$\partial_{t}\mathcal{K}+J(\psi,\mathcal{K})+\nabla\cdot\mathbf{F}_{\mathcal{K}}+\partial_{z}G_{\mathcal{K}}=\gamma_{\mathcal{K}}^{F}+d_{\mathcal{K}},$
(4)

where $\mathcal{K}$ is the NIW kinetic energy density,
$\mathbf{F}_{\mathcal{K}}$ and $G_{\mathcal{K}}$ are the horizontal and
vertical energy fluxes due to dispersion, $\gamma_{\mathcal{K}}^{F}$ is the
energy input by the forcing, and $d_{\mathcal{K}}$ is the dissipation due to
the hyperviscosity term. The horizontal energy flux $\mathbf{F}_{\mathcal{K}}$
can be expressed in terms of gradients of the phase $\Theta$ of $M$ (Rocha et
al., 2018):

$\mathbf{F}_{\mathcal{K}}=\frac{N^{2}|M|^{2}}{2f}\nabla\Theta.$ (5)

The spatial heterogeneities that the mesoscale vorticity $\zeta$ induces in
the NIW phase thus cause a transfer of NIW energy in the horizontal. Once this
$\zeta$ refraction has imprinted the horizontal structure of eddies onto an
initially uniform NIW field, advection will also act on the resulting
gradients and stir the NIW field. Dispersion remains important and helps waves
escape straining regions (Rocha et al., 2018). The scale reduction also
accelerates the propagation of NIWs out of the mixed layer. This is encoded
into the YBJ dynamics through the vertical flux $G_{\mathcal{K}}$, which
requires horizontal gradients in the NIW field to produce energy fluxes to
depth (Appendix A).

The energy input into the NIW band by the winds $\gamma_{\mathcal{K}}^{F}$ is
also known as the NIW wind work. The difference between the NIW wind work and
the NIW kinetic energy that propagates out of the surface layer is the energy
that is available for NIW mixing in the surface ocean. Both processes
therefore influence the mixed layer depth. There has been extensive effort to
estimate the NIW wind work. It can be calculated directly from concurrent
observations of winds and NIW surface currents, but this is possible in a few
locations only. Slab models have been used to obtain global estimates of the
NIW wind work (Alford, 2001) but these are suspected to be overestimates. This
is primarily attributed to the fact that these models poorly represent the
various processes that cause NIWs to leave the mixed layer. This can include
the vertical propagation discussed above but the employed models often also do
not represent dissipation of NIWs by conversion to turbulent kinetic energy,
which would reduce the projection of the wind stress onto the waves
(Plueddemann and Farrar, 2006; Alford, 2020). High-resolution ocean models can
also be used to estimate the wind work, but here problems can arise from
limitations in the reanalysis data used to force the models. For example,
Flexas et al. (2019) showed that the NIW wind work was poorly represented in a
high-resolution ocean model due to the wind forcing missing variability on
scales less than 6 hours and 15 km. The power input is larger in high-
resolution coupled simulations but still substantially lower than
observational estimates from slab models (von Storch and Lüschow, 2023).

The mesoscale can act as a source of wave potential energy in a process known
as stimulated generation. In unforced and inviscid YBJ dynamics, the wave
kinetic energy is conserved but NIWs can gain or lose potential energy through
interactions with the mesoscale. In the YBJ framework, the domain-integrated
potential energy is not conserved (Appendix B):

$\partial_{t}\langle\mathcal{P}\rangle=\Gamma_{\mathcal{P}}^{R}+\Gamma_{\mathcal{P}}^{A}+\mathcal{D}_{\mathcal{P}},$
(6)

where $\mathcal{P}$ is the potential energy density,
$\Gamma_{\mathcal{P}}^{R}$ is the production of potential energy by
refraction, $\Gamma_{\mathcal{P}}^{A}$ is the production of potential energy
by advection, and $\mathcal{D}_{\mathcal{P}}$ is the potential energy
dissipation by hyperviscosity. The formation of horizontal structure in an
initially uniform NIW causes an increase in NIW potential energy via
$\Gamma_{\mathcal{P}}^{R}$ and $\Gamma_{\mathcal{P}}^{A}$. In the first phase
of interaction, $\Gamma_{\mathcal{P}}^{R}$ dominates as the waves have the
mesoscale structure imprinted on them. At later times,
$\Gamma_{\mathcal{P}}^{A}$ becomes the most important term. An extension of
YBJ describing the coupled evolution of NIWs and a quasi-geostrophic mesoscale
field was derived by Xie and Vanneste (2015) (see also Wagner and Young,
2016). This extension shows that the refractive and advective sources of
potential energy to the NIWs appear as sinks in the mesoscale energy budget.
This allows us to interpret $\Gamma_{\mathcal{P}}^{R}$ and
$\Gamma_{\mathcal{P}}^{A}$ as energy transfers from the mesoscale, despite the
fact that we impose the mesoscale field in the simulations described below.
Theoretical and numerical studies (Xie and Vanneste, 2015; Rocha et al., 2018;
Asselin and Young, 2020) have investigated the process of stimulated
generation, yet its importance in the real ocean remains poorly constrained.

Several attempts have been made to reconcile available observations with our
theoretical understanding of NIW–mesoscale interactions. Work on the
NIW–mesoscale interaction prior to YBJ had been based on ray tracing theory,
which additionally assumes that the waves have much shorter spatial scales
than the background mesoscale flow (Kunze, 1985). Ray tracing predicts the NIW
frequency to be shifted by $\frac{1}{2}\zeta$. This prediction, however,
applies only in regions of the ocean where the WKB limit of ray tracing is
appropriate. YBJ argued that this was not the case for the Ocean Storms
Experiment, suggesting that the observations were taken in a region where the
waves are instead in the so-called strong-dispersion limit. In this limit,
refraction of the large-scale wave field is strongly opposed by dispersion,
and the frequency shift is much smaller than predicted by WKB theory. This
provides a compelling potential explanation of D’Asaro’s (1995) observation
that the NIW frequency shift was at least five times smaller than
$\frac{1}{2}\zeta$ during Ocean Storms.

This interpretation of the Ocean Storms Experiment was pursued further by
Balmforth et al. (1998), who ran spin-down simulations of the YBJ equation.
NIWs were initialized in the mixed layer and evolved in the presence of an
idealized, barotropic mesoscale eddy field. Qualitative comparisons between
the simulations and observations showed that YBJ dynamics were not
inconsistent with the observed time for NIWs to escape the mixed layer.
Balmforth and Young (1999) showed that including the $\beta$-effect improved
the agreement with observations. Because an idealized eddy field was used,
however, no quantitative conclusions could be drawn about the ability of YBJ
to capture the observed evolution.

More recently, Asselin and Young (2020) investigated the fate of NIWs as they
propagate into a baroclinic mesoscale eddy field using numerical simulations
of the extended YBJ system that also accounts for the effect of the waves on
the mesoscale. They observed the initially horizontally uniform NIWs
undergoing scale reduction by $\zeta$-refraction and then propagate downwards
in anticyclones. The vertical wave propagation terminated at depth due to the
decay of the baroclinic vorticity away from the surface. For strong NIWs, they
also found that the mesoscale eddy field was weakened due to stimulated
generation.

In this study, we use the YBJ framework to interpret the observed evolution of
NIWs in the Northeast Atlantic Ocean. We use an array of nine moorings to
capture some of the mesoscale variations in the NIW field. The YBJ framework
allows us to attribute the observed NIW evolution to a set of well-defined
physical processes. We integrate the three-dimensional YBJ equation using
observational inputs for the wind forcing, mesoscale streamfunction, and
stratification, and we compare these simulations to simpler slab models. We
show that the YBJ model offers significant improvements in modeling NIW
evolution, without the need for any tuning. We use the YBJ energy budgets to
provide a dynamical interpretation of spatial and temporal variations in the
NIW field and quantify the relative importance of the various physical
processes involved.

##  2 Observations

###  2.1 NIW Data

Figure 1: (a) Location of the OSMOSIS mooring array (yellow square) overlaid
on a map of the bathymetry. (b) Geometry of the mooring array that consists of
a central mooring (C) surrounded by an inner and an outer group of moorings.
The inner moorings are labeled with an “I” followed by the cardinal direction
of the mooring in relation to the central mooring. The outer moorings are
labeled in a similar way but with an “O” identifier. (c) Time-series of NIW
amplitude extracted from the OSMOSIS mooring observations from September 2012
to September 2013. The solid line indicates the observations at the central
mooring while the orange ribbon represents the range across the mooring array.
The gray shading indicates the five NIW events discussed in text. To
illustrate the process of extracting the NIW signal, the inset shows the raw
velocity (top) and NIW velocity (bottom) during event 5. The amplitude of the
envelope modulating the NIW velocity is $|\partial_{z}M|$.

We extract observations of NIWs from data collected for the Ocean Surface
Mixing, Ocean Submesoscale Interaction Study (OSMOSIS; Buckingham et al.,
2016). As part of the study, measurements were taken from nine moorings
anchored over the Porcupine Abyssal Plain in the northeastern part of the
Atlantic Ocean from September 2012 to September 2013. The moorings were all
distributed less than 10 km from a central mooring located at
($48.6875^{\circ}$N, $-16.1875^{\circ}$E) (Fig. 1a,b). In this study, we focus
on observations from the central mooring and the four outermost moorings. Each
mooring took measurements of the horizontal velocity once every ten minutes
using an acoustic current meter (ACM). On the central mooring, there were 13
ACMs spaced nominally between 50 m and 500 m depth. On the outer moorings
there were five ACMs nominally spaced over the same depth range. The deepest
ACM failed on the outer northeast (ONE) mooring, and so measurements are only
available to a depth of 361 m.

The NIW signal is extracted using a Gaussian filter in frequency space. At
this latitude, the near-inertial frequency (1.1e-4) is close to the frequency
of the M2 tidal constituent (1.4e-4), and so we choose the filter width of
1.2e-5 to exclude the sizable tidal motion from the filtered time series. The
conclusions we draw below are not sensitive to the filter width so long as the
filter does not include the M2 peak. We identify five events in the year-long
time series where strong NIWs were forced relatively coherently across the
mooring array (Fig. 1c).

The actual depth of the instruments varied in time due to drag exerted by the
flow on the moorings. Over the course of the year, there were several times
where the moorings experienced knockdown by up to 200 (Callies et al., 2020).
We flag these events if the knockdown on the central mooring is more than 10
when averaged with a 1-day running mean. Outside of these knockdown events the
variations in mooring depth about the nominal depth is $\sim$2\. The
interpretation of the observed NIW signal during these knockdown events is
complicated because the filter we use to extract the NIW signal does not
commute with evaluating the signal along the trajectory of the ACMs. We do not
attempt to explain the observations during the knockdown periods.

###  2.2 Stratification and Mixed Layer Depth

We need the stratification and mixed layer depth as observational inputs to
the numerical YBJ simulations. The stratification influences the dispersion
term in the YBJ equation and the mixed layer depth informs us about what depth
range to force the waves over. Two ocean gliders sampled across the mooring
array during the OSMOSIS study period (Damerell et al., 2016; Thompson et al.,
2016). The gliders measured temperature and salinity in the upper 200 m of the
water column approximately once every two hours. Following Damerell et al.
(2020), we calculate the mixed layer depth $H_{m}$ from the glider
observations as the depth at which either the potential density $\rho$ or the
temperature $T$ has changed more than a given threshold from the value at 5
depth. These thresholds are $\Delta\rho=\qty{.03}{\kilo\per\cubed}$ for
density and $\Delta T=\qty{0.2}{}$ for temperature. The final mixed layer
depth is taken as the shallowest of the two estimates.

The stratification $N^{2}$ is calculated using:

$N^{2}=-\frac{g}{\rho}\frac{\partial\rho}{\partial z},$ (7)

where $g$ is the acceleration due to gravity. The stratification is then
averaged over a given event and the glider trajectories to result in a single
spatially and temporally averaged profile used for simulations.

###  2.3 Wind Data

For the wind forcing, we use the European Centre for Medium-Range Weather
Forecasting ERA-5 reanalysis (Hersbach et al., 2018). We extract time series
of the 10 m zonal ($u_{w}$) and meridional ($v_{w}$) winds with hourly
resolution at the grid point that contained the mooring array. Following
Pollard and Millard (1970) we convert this to a wind stress using a bulk
aerodynamic drag formulation. In terms of the complexified wind velocity
($\mathcal{U}_{w}=u_{w}+iv_{w}$), the complexified wind stress $\tau_{w}$ is
given by:

$\tau_{w}=\rho_{a}C_{D}|\mathcal{U}_{w}|\mathcal{U}_{w},$ (8)

where $\rho_{a}$ is the density of air and $C_{D}$ is the drag coefficient,
which we calculate using the speed-dependent formulation of Large and Pond
(1981).

###  2.4 Altimetry

To characterize the mesoscale eddy field, we use observations of the sea
surface height (SSH) from the Data Unification and Altimeter Combination
System’s (DUACS) delayed-time (DT) 2018 release Taburet et al. (2019). The
DUACS DT2018 SSH maps are provided at a nominal $(1/4)^{\circ}$ and daily
resolution. We convert these measurements to a geostrophic streamfunction
using $\psi=gh/f$, where $h$ is the SSH and $f$ is the latitude-dependent
Coriolis parameter.

##  3 Models

###  3.1 The PM Model

We begin by writing the PM model in the language of YBJ. Adding (1a) to
$i\times\eqref{eq:PM2}$ and multiplying by $e^{-ift}$ yields

$\partial_{tz}M=-r\partial_{z}M+\frac{\tilde{\tau}}{\rho_{w}H_{m}},$ (9)

where $\tilde{\tau}=e^{-ift}\tau_{w}$ is the back-rotated, complexified wind-
stress. For each NIW event, we solve the PM equation using the reanalysis wind
stress and the mixed layer depth from the gliders. In order to better
understand the role of $\zeta$ refraction in the life cycle of wind-forced
NIWs, we perform a second run of the PM model (denoted by PM+$\zeta$)
augmented by the refraction term from the YBJ equation:

$\partial_{tz}M=-r\partial_{z}M-\frac{i\zeta}{2}\partial_{z}M+\frac{\tilde{\tau}}{\rho_{w}H_{m}}.$
(10)

This shifts the NIW phase and hence changes the phasing of the NIWs relative
to the wind. It does not, however, capture the horizontal energy transfers
induced by $\zeta$ refraction because these transfers require horizontal
structure in the wave field, which is absent from the horizontally uniform PM
model.

For each event, we use a constant mixed layer depth that is an average of the
time-varying mixed layer depth over duration of the event. Since the events
are relatively short, the error in this approximation is minimal for most
events. We make this choice of a constant mixed layer depth to avoid
complications that arise otherwise, especially in the YBJ model discussed
below.

The parameter $r$ is intended to account for all of the processes that
decrease the wave amplitude in the region of interest. This is primarily
thought of as vertical propagation of NIWs out of the mixed layer, but other
processes such as advection and dissipation may also cause the NIW amplitude
to decrease. Treating all of these processes as a Rayleigh drag term with a
single decay parameter represents a drastic simplification in the PM model. It
also introduces a free parameter. Previous studies have chosen $r$ such that
the resulting solutions fit observations as closely as possible. The original
PM paper used both $r^{-1}=\qty{4}{}$ and $r^{-1}=\qty{8}{}$. Alford (2001)
used a damping that would correspond to about $r^{-1}=\qty{4}{}$ at the
latitude of the OSMOSIS mooring array. Yu et al. (2022) used
$r^{-1}=\qty{16.7}{}$ to estimate the NIW wind work during the OSMOSIS
experiment using the PM model. We speculate that they had to use very weak
damping because their wind data was taken from ERA-interim reanalysis, which
has 6-hourly analysis steps with forecasts used to increase the time
resolution to 3 hours, which may suppress the wind power at frequencies
important for NIW generation. A similarly weak damping is likely unsuitable
here, given that ERA-5 reanalysis winds have more power in the near-inertial
band (see discussion in Flexas et al., 2019). Nevertheless, we vary $r^{-1}$
between 4;16.

For the PM+$\zeta$ model, we calculate the vorticity using the streamfunction
from altimetry. The data processing is as described below for the YBJ model,
and we select the vorticity value at the grid point nearest to the center of
the mooring array.

We initialize the simulations with no waves and then allow the model to spin
up before the main forcing period for each event. We choose the initial time
to be when the observations show relatively little waves. This is done by eye.
These initial times are followed by a strong forcing event and so the NIW
signal is dominated by the newly generated waves, implying that the error from
using an initial condition with no waves is relatively small. This discussion
also applies to the YBJ simulations we run. We integrate both of the models
(PM and PM+$\zeta$) using a Crank–Nicolson scheme.

###  3.2 The YBJ Model

The YBJ equation is a three-dimensional partial differential equation, making
it substantially more computationally expensive to integrate than the PM
model. We solve the YBJ equation using the pseudospectral solver Dedalus
(Burns et al., 2020) with a mixed explicit and implicit diagonal RK2 scheme.
As discussed above, we start the simulations with no waves. We use a domain
that is 400 x 400 in the horizontal (centered on the mooring array) and 4 km
deep. Each dimension is discretized with 128 modes. The vertical dimension is
finite and represented using Chebyshev polynomials. The horizontal dimensions
are made periodic and represented using Fourier modes. The stratification and
the wind forcing are taken to be horizontally uniform. To construct the
mesoscale streamfunction used in the simulations, we take $\psi$ from
observations and calculate the vorticity by taking the finite-difference
Laplacian of $\psi$ on the sphere. We then interpolate this vorticity onto the
Cartesian simulation grid centered on the location of the central mooring. We
apply a taper to the vorticity, such that it goes to zero on the domain
boundaries. We then invert this tapered vorticity field for the streamfunction
on the periodic domain, making the resulting streamfunction periodic as well.
Our conclusions are neither sensitive to the form nor width of the taper. For
the results shown below, we use the following form of the taper:

$\mathcal{T}(x,y)=\cos^{2}{\frac{\pi x}{L_{x}}}\cos^{2}{\frac{\pi y}{L_{y}}},$
(11)

where $L_{x}$ and $L_{y}$ are the $x$ and $y$ lengths of the domain,
respectively.

The YBJ equation can be formulated on a $\beta$-plane by making the
substitution $\zeta/2\rightarrow\zeta/2+\beta y$. We perform all of our
analysis on an $f$-plane (i.e., $\beta=0$) for two reasons: (i) We focus our
analysis on a region in the Northeast Atlantic Ocean where $\beta L\ll\zeta/2$
with $L$ a typical meridional scale of the waves. This relative scaling varies
regionally in the ocean. For example, Thomas et al. (2020) studied a region
that was similarly dominated by the vorticity and found that the $\zeta$
refraction process was more important to the NIW evolution there than $\beta$
refraction, whereas $\beta$ refraction appeared to be important for Ocean
Storms (D’Asaro et al., 1995). (ii) The $\beta$ term adds difficulty in
simulating the YBJ equation numerically (see Balmforth and Young, 1999).

For the forcing, we need to prescribe a vertical profile that determines at
what depths the momentum is deposited. We specify the profile for the body
force $\partial_{z}\mathcal{F}$ such that it is constant in the forcing layer
and then decays rapidly to zero below:

$\partial_{z}\mathcal{F}=\frac{\alpha\tilde{\tau}}{H_{m}\left[\alpha+\ln(2\cosh\alpha)\right]}\left[1+\tanh{\left(\alpha\frac{z+H_{m}}{H_{m}}\right)}\right],$
(12)

where $\alpha$ is a parameter determining the steepness with which the body
force falls off below the forcing layer. We use $\alpha=2$ in all the
simulations presented below. The prefactor ensures that
$\mathcal{F}=\tilde{\tau}$ at $z=0$.

When the winds blow on the ocean surface, they generate turbulence that mixes
momentum downwards. If a mixed layer already exists and is not too deep, the
momentum input from the wind will be rapidly homogenized within the mixed
layer (Pollard and Millard, 1970; Kato and Phillips, 1969). For some events at
the OSMOSIS site, the mixed layer was up to a few hundred meters deep. These
deep mixed layers are likely the result of convection driven by buoyancy
forcing rather than the mechanical wind forcing (Thompson et al., 2016). In
these cases, it is unlikely that the momentum is uniform across the mixed
layer, especially if the buoyancy forcing has ceased. The depth structure of
NIWs obtained from the OSMOSIS mooring confirms this picture. The waves are
initially forced over a layer that is thinner than the mixed layer before they
propagate to depth. To avoid forcing over an unrealistically large depth, we
cap the forcing layer at 80m. This value is guided by the observations and
represents an average depth over which waves are forced when the mixed layer
is deep. We discuss possible ways to improve this representation below.

In order to solve the YBJ equation, we must further specify vertical boundary
conditions. The requirement that the vertical velocity is zero on the top and
bottom of the domain translates to the requirement that $M$ be horizontally
uniform at the top and bottom boundaries (Young and Ben Jelloul, 1997). Since
$M$ is determined by $\partial_{z}M$ up to an arbitrary horizontal function,
without loss of generality we specify $M(x,y,-H,t)=0$. The top boundary
condition can be found by vertically integrating the YBJ equation twice and
requiring that $M$ be horizontally uniform at $z=0$ (Appendix B). The result
is

$\partial_{t}M(x,y,0,t)=\mathcal{F},$ (13)

which we integrate a priori and then use as a Dirichlet boundary condition on
$M$.

##  4 Results

###  4.1 Case Study: Fall Event

Figure 2: (a) Snapshot of mesoscale vorticity field in the inner quarter of
the simulation domain on 2012-10-18 which was the peak of event 1. The stars
show the locations of three specific moorings. The red star denotes the ONE
mooring which is in a region of cyclonic vorticity, the black star denotes the
central mooring and the blue star denotes the OSW mooring which is in an anti-
cyclonic region. (b) Time series of the vorticity during event 1 at each of
the three moorings denoted by the stars above. During the main forcing period
the vorticity changes sign across the mooring array while at later times the
entire mooring array transitions to being in a region of cyclonic vorticity.

We begin with a detailed analysis of the simulation results for event 1, which
occurred in the fall (Fig. 1c). This event is by far the most energetic NIW
event observed throughout the year. The main forcing for event 1 occurred when
the mooring array straddled a dipole in the mesoscale vorticity (Fig. 2),
making it a good candidate to see the effect of $\zeta$ refraction.

####  4.1.1 NIW Amplitude and Phase

Figure 3: (a) NIW amplitude in the PM model compared to observations at the
central mooring (dashed black line). The PM model was run with a range of
difference values for the damping parameter $r$. Solid lines indicate the
results of the PM model without refraction and dotted lines indicate the
results with the refractive term added. The gray shaded region is a period of
mooring knockdown. (b) As in (a) but for the NIW phase. (c) NIW amplitude in
the YBJ simulation (solid lines) compared to observations (dashed lines) at
the central mooring (black) as well as the ONE (red) and OSW (blue) moorings.
(d) As in (c) but for the NIW phase.

The observed NIW amplitude begins increasing around Oct. 14 and reaches a
maximum value on Oct. 20 before returning to background levels by the end of
our simulation period on Nov. 10 (Fig. 3a). There are large variations in the
peak amplitude over the mooring array, although these differences disappear by
the end of the event (Fig. 3c). The phase increases sharply near the beginning
of the forcing period as it aligns with the wind (Fig. 3b). Following this, it
slowly increases for most of the event and is relatively homogeneous across
the mooring array except for a period following the initial forcing (Fig. 3d).

The PM model has trouble capturing the NIW evolution during this event (Fig.
3a,b). Using a value of the damping parameter $r^{-1}=\qty{4}{}$ results in a
peak amplitude close to that of the observed NIW peak amplitude in the central
mooring, but such a relatively strong damping causes the NIW amplitude to drop
off much too quickly compared to the observations. A more realistic amplitude
decay is achieved when weaker damping is applied (especially with
$r^{-1}=\qty{16}{}$), but then the peak NIW amplitude is overestimated
substantially. The simulated phase bears little resemblance to the observed
phase, with the simulated phase being offset by up to $\pi/2$ during the
initial forcing period, and the simulated phase remains close to constant
around $\pi$ as long as substantial amplitude remains, missing the gradual
increase in the observed phase.

The addition of the refractive term does little to change the simulated NIW
amplitude in the PM model during most of the event (Fig. 3a,b), indicating
that the phase shifts introduced by refraction are unable to substantially
change the relative alignment between the NIWs and the winds. At later times
the effect of adding refraction is more pronounced in the runs with lower
damping as there is more time for refraction to act on the waves before they
are dissipated. Refraction seems to dominate the phase evolution at all three
damping values, but the phase trends in the opposite direction from what is
observed. The phase tendency is also larger in magnitude than what is
observed.

While the PM model may be able to capture some qualitative features of the NIW
observations, it seems clear that: (i) the physical processes that cause the
NIWs amplitude to decay cannot be accurately captured using a simple linear
drag formulation, and (ii) if refraction is important to NIW evolution, its
effect is not simply to change the phase of the NIWs but must involve
horizontal processes.

The YBJ model captures the observed amplitude and phase evolution much better
than the PM model (Fig. 3c,d). At the central mooring, the YBJ simulation
agrees with observations in terms of peak amplitude and decay timescale.
Similarly, the phase evolution is much closer to observations than the PM
results, in terms of both its value after the forcing and its trend afterward.
Again, we emphasize that, unlike the PM models, there is no tunable parameter
in the YBJ simulations.

The YBJ model also captures observed lateral variations of the NIW signal
across the mooring array (Fig. 3c,d). At the OSW mooring, which at the start
of the event is in a region of anti-cyclonic vorticity (Fig. 2), the YBJ
simulation successfully predicts a substantial enhancement in the NIW
amplitude compared to the central mooring. At the ONE mooring, which at the
start of the event is in a cyclonic region (Fig. 2), the YBJ simulation
successfully shows a reduction in the NIW amplitude compared to the central
mooring. The YBJ model also captures that the NIW phase is much more uniform
across the mooring array than the amplitude.

####  4.1.2 NIW Kinetic-Energy Budget

Figure 4: (a) NIW kinetic energy budget terms. The kinetic energy tendency
(dashed line) is decomposed into the 5 processes in the model which can change
the kinetic energy: advection (blue), horizontal flux divergence (orange),
vertical flux divergence (green), hyperviscosity (brown) and wind forcing
(pink). The budgets are evaluated at the horizontal position of the moorings
and at fixed depth of 25. (b) As in (a) but for the ONE mooring. To better
visualize the terms we only plot the budget for the first 3/4 of the event.
Figure 5: (a) Color map of vorticity (normalized by $f$) at the peak of event
1. Arrows indicate the horizontal energy flux field at $z=-\qty{25}{}$,
showing a flux of NIW kinetic energy out of cyclones and into anticyclones.
(b) Horizontal sections of NIW amplitude at a depth of 25 and three different
times. Sections show the inner half of the simulation domain. There is no
horizontal structure to the forcing and hence the NIWs are initially forced
uniformly throughout the domain. Refraction is the only process which can
impose structure on a uniform NIW field in the YBJ framework. NIWs begin to be
concentrated into anticyclones. Brighter regions in the Oct. 20 plot correlate
with anticyclonic regions in the altimetry. Once refraction creates horizontal
structure other processes can act. Dispersion will eventually counteract
concentration into anticyclones. Advection by the mesoscale eddies will also
stir horizontal structure created by refraction. A signal of advective
stirring is clearly visible in the upper right quadrant of the Oct. 27 plot.

The observations show a higher NIW amplitude at the mooring in the anticyclone
than at the mooring in the cyclone (Fig. 3c,d). This is suggestive of $\zeta$
refraction concentrating NIW kinetic energy in anticyclones. It is hard to
draw this conclusion from the observations alone, however, because other
factors could give rise to the observed amplitude differences. For example,
the current meters were located at slightly different depths (between
$z=-\qty{44}{}$ and $z=-\qty{62}{}$).

With the YBJ simulations capturing the observed differences between the
moorings, we can use the YBJ framework to identify the processes giving rise
to these lateral variations. We construct a point-wise kinetic-energy budget
(4) at the OSW (anticyclonic region) and ONE (cyclonic region) moorings, which
allows us to separate the kinetic-energy tendency into advection, dispersive
flux divergence, dissipation, and forcing (Fig. 4). During the initial forcing
period (up to Oct. 20), the tendency due to wind forcing is similar for both
moorings. At both locations, there is a small positive advective tendency that
turns slightly negative toward the end of the forcing period. The vertical
flux divergence is also similar at the two locations and smaller in magnitude
than the advective tendency. Dissipation is negligible. The most notable
difference between the two mooring locations is in the horizontal flux
divergence term. There is horizontal flux convergence at the OSW mooring and
flux divergence at the OSE mooring. This causes the total tendency to be
larger than the wind forcing at the OSW mooring and smaller than the wind
forcing at the ONE mooring.

At these early times, $\zeta$ refraction is the primary driver of the
horizontal energy flux and causes concentration of NIW kinetic energy into
anticyclonic regions. During the peak of the event, the horizontal energy flux
is directed from cyclonic to anticyclonic regions (Fig. 5a). This arises from
an interplay of the refraction and dispersion terms in the YBJ dynamics:
refraction sets up phase gradients, which cause a dispersive energy flux as
described by (5). As horizontal structure develops, advection can also become
important as it stirs the existing horizontal structure. This sequence of
events was described by Rocha et al. (2018) and captures the early evolution
during this event (Fig. 5b).

We again emphasize that the PM model, even if the refractive term is included,
cannot capture these dynamics because lateral energy transport originating
from the dispersion term in the YBJ equation is crucial. Once dispersion is
included, the YBJ model captures observed lateral variations in the NIW
amplitude between anticyclonic and cyclonic regions and offers a clear
dynamical explanation.

####  4.1.3 NIW Potential-Energy Budget

Figure 6: Vertically integrated and horizontally averaged potential energy
budget showing the potential energy tendency (dashed line) split into
contributions from refractive production (purple), advective production (red)
and hyperviscosity (gray).

The NIW potential-energy budget is of interest because it provides insight
into the energy exchange with mesoscale eddies (Xie and Vanneste, 2015; Rocha
et al., 2018). While we prescribe the mesoscale eddy field using altimetry
rather than evolving a coupled system, we still interpret the sources of NIW
potential energy as estimates of the energy transfer from mesoscale eddies
facilitated by NIW refraction and advection.

The vertically integrated and domain-averaged budget for event 1 shows that
NIW potential energy is generated by both refraction and advection (Fig. 6).
The potential-energy tendency rises sharply from zero to a peak value after
the main forcing period of the event, lagging the peak in the kinetic-energy
tendency by a few days. At these early times, almost all of the NIW potential
energy is created by refractive production. As horizontal structure is
created, advective production ramps up. The refractive production term
decreases throughout the rest of the event and by the end is overwhelmed by
advective production. The initial dominance of refractive production followed
by an increase in advective production is very similar to the succession of
events Rocha et al. (2018) described for stimulated generation in idealized
simulations that included the full coupling with the mesoscale dynamics.

The potential-energy production—and presumed sink of mesoscale energy—peaks at
about 20□. To put this number in context, a global input of 1 into the
mesoscale eddy field (Wunsch and Ferrari, 2004) corresponds to about 3□ on
average. For stimulated generation to be important in the global energy budget
of mesoscale eddies, it must be much stronger elsewhere.

####  4.1.4 NIW Wind Work

Figure 7: Wind-work as estimated from observations at the central mooring
(dashed black), PM (solid green), PM+$\zeta$(dotted green) and YBJ (solid
black). For observations the wind-work is defined as
$\bf{\tau}_{NIW}\cdot\bf{u}_{NIW}$, while for the simulations it is defined as
the kinetic energy production term from the respective kinetic energy budget.
The PM models use $r^{-1}=4$ days which best matches the peak amplitude.

The interaction between NIWs and mesoscale eddies also affects the wind work
in the near-inertial band (Fig. 7). The PM and PM+$\zeta$ models with
$r^{-1}=\qty{4}{}$ best match the peak amplitude in the NIW evolution (Fig.
3a) but overestimate the wind work at the location of the central mooring by a
factor of more than two. We could also tune $r$ to match the wind work
estimated directly from observations integrated over the event. This is
achieved with $r^{-1}=\qty{0.576}{}$, but this means that the peak NIW
amplitude is underestimated by a factor of three and violates $r\ll f$. The
effect of refraction in the PM model is to slightly increase the wind work. In
the YBJ model, in contrast, the wind work matches the observations well. This
is a consequence of the YBJ model’s ability to closely reproduce the observed
NIW evolution. Whether this reveals a substantial bias in wind work estimates
based on the PM model remains unclear, however, because this particular
location and event may not be representative of a time and space average.

###  4.2 Further Events in the Time Series

Figure 8: NIW amplitude (left column) and NIW phase (right column) in
observations at the central mooring (dashed black line) as compared to PM
(solid green) PM+$\zeta$ (dotted green) and YBJ (solid black) for (a) event 2,
(b) event 3, (c) event 4a, (d) event 4b and (e) event 5.

Simulations of the remainder of the NIW events (Fig. 8) confirm that the YBJ
model better captures the observed NIW amplitude and phase evolution than the
PM model does (Fig. 8). Event 4 is different from the others in that it
consists of a double peak in NIW amplitude without a decay to background
levels in between. The mixed layer depth was much more variable than during
other events, especially towards the end of the event. For that reason, we
simulate event 4 in two parts (events 4a and 4b) in order to minimize the
variations in the mixed layer depth over a simulation period.

In general, we see that the YBJ model performs better than the slab models in
reproducing the observed evolution of the NIW amplitude and phase. In event 2,
the YBJ model captures the slow rise of the NIW amplitude as well as the
decrease at late times. The PM models are not able to capture this behavior as
well (Fig. 8a). All of the models, however, have substantial errors in the
phase at early times. This could be due to pre-existing NIWs that we do not
capture because we initialize the models at rest. The notable exception to the
trend that YBJ performs better than the PM models is event 3 (Fig. 8b), where
the YBJ model not only overestimates the peak but also the decay time of the
waves. YBJ does, however, capture the slow evolution of the phase better than
the PM models. We discuss some of the potential reasons for these
disagreements below. For event 4a, all models predict a double peak in the
amplitude that is not seen in observations (Fig. 8c). The YBJ model does
better at later times in both amplitude and phase. In event 4b, the YBJ model
does rather well in predicting the complete evolution of both the NIW
amplitude and phase (Fig. 8d). The YBJ model predicts the timing of the peak
better than the PM models in event 5, but the waves persist for longer than
observed (Fig. 8e).

##  5 Discussion

The YBJ model does well in reproducing the mooring observations of NIWs. There
are still differences between the model and the observations, however, as well
as a couple of events where the YBJ model does less well. Given the
observational inputs to the model, it would be surprising if this were not the
case. One major limitation is the use of altimetry for the mesoscale
streamfunction. As mentioned previously, the altimetry data are a heavily
smoothed version of the real field. We suggested above that some of the
discrepancies between the YBJ model and the observations were due to this
smoothing. We also interpolated the mesoscale vorticity onto a smaller
simulation grid. The resolution of the altimetry product is
$0.25^{\circ}\times 0.25^{\circ}$, which is larger than the size of the
mooring array. Differences across the mooring array come from the
interpolation between neighboring altimetry grid cells, which will have
introduced interpolation errors.

Furthermore, we assumed that the mesoscale eddy field was barotropic. This is
a reasonable assumption if the vertical scale of the waves is much smaller
than the vertical scale of the eddies. There are certainly errors in the YBJ
evolution, however, that arise from neglecting the baroclinicity of mesoscale
eddies. These effects could be investigated in the future by running similar
simulations using an in-situ data set that resolves the vertical structure of
one or more eddies.

In event 1, the vorticity has large variations across the mooring array (Fig.
2a). This is reflected in the NIW observations as large differences in the
amplitude across the mooring array. However, this need not be the case for all
events. Event 3, by contrast, shows weak variations in the vorticity (as
diagnosed from altimetry) across the mooring array. The result is that the YBJ
simulations also show weak variations in the NIW field across the mooring
array region. The altimetry vorticity field is a heavily smoothed
representation of the real vorticity field of the ocean; smaller scale
vorticity features are invisible to measurements from altimetry. If there is
little variation in the larger-scale vorticity across the mooring array, then
these smaller-scale vorticity features may play a more important role in
governing the variations in $\partial_{z}M$. While we do not have a spatial
map of smaller-scale vorticity features, we can estimate their magnitude at
the mooring array by calculating the vorticity by applying Stokes’ theorem to
the area bounded by the outer moorings (Buckingham et al., 2016). We use a
low-pass filtered velocity to extract the balanced signal and give us an
estimate of the mesoscale vorticity. This vorticity agrees in general with the
vorticity calculated from SSH but shows more high-frequency variability (Fig.
9). Event 3 shows the biggest disagreement between the two estimates of all
the events. The mooring estimate shows cyclonic vorticity, whereas the
altimetry estimate shows anticyclonic vorticity. This likely explains why we
see the amplitude decay much quicker in observations compared to simulations
because there would have been a horizontal energy flux out of the region while
the simulations have a flux into the region, which acts to maintain the
amplitude. The other events show better agreement between the two vorticity
estimates, although there are times where the deviation is larger. Events 1
and 2 specifically show two anticyclonic periods that are not captured by
altimetry. These may explain some of the mismatch between simulations and
observations. For example, the anticyclonic excursion in event 1 probably
explains why the YBJ simulation results decay more quickly than observations.
For event 2, the excursion occurs primarily near the start of the event, where
the NIW amplitude is weak, so its effect is tempered. Many of these anomalies
are short in duration, which limits the error in using the altimetric
vorticity.

The vorticity error for event 4a is minimal (Figure 9). We suspect the
mismatch between the YBJ model and the observations in this event is due to
the forcing. The reanalysis product used is not the exact wind-forcing felt at
the OSMOSIS site. Both YBJ and PM show an initial peak in the NIW amplitude
that is not seen in the observations, indicating that the wind-stress may be
wrong at this point. The forcing event seen in observations then has to
destroy these waves before forcing new ones, which causes the lower amplitude
of the NIW peak in the YBJ simulation compared to observations.

Figure 9: Comparison of the vorticity as calculated from Stokes’ theorem
applied to the outer moorings (orange) and from altimetry (blue). The
velocities used in the Stokes’ theorem calculation were obtained by low-pass
filtering with a second order Butterworth filter and a cutoff frequency
corresponding to a period of approximately 5 days. The blue ribbon shows the
spread in the vorticity when interpolated onto the outer mooring positions.

Using a time-averaged stratification profile likely also contributes to
differences between the observations and simulations, despite our attempts to
keep simulation times as short as possible to avoid this effect. The
assumption that the stratification does not vary in time is inherent in the
YBJ scaling assumptions. Related to this is the question of what depth to
force the model over. We used a fixed forcing depth, although this certainly
varied over the simulation periods. This problem could be side-stepped by
representing the wind forcing as a surface stress and adding a parameterized
turbulent vertical momentum flux to the YBJ equation to transfer momentum from
the wind downwards.

The YBJ equation contains no term to represent the breaking of NIWs. The
scaling $H^{*}=u^{*}/\sqrt{N^{*}f}$, where $u^{*}=\sqrt{\tau/\rho_{0}}$ is the
friction velocity and $N^{*}$ is the stratification at the base of the mixed
layer, indicates the depth to which wind-driven turbulence could deepen the
mixed layer (Pollard et al., 1973). If the mixed layer depth is shallower than
$H^{*}$, the effect of NIW dissipation would be important as there is
potential for mixed layer deepening by breaking of NIWs (Plueddemann and
Farrar, 2006). During the OSMOSIS study, Yu et al. (2022) showed that the
mixed layer was generally deeper than $H^{*}$. This is especially true outside
of summer and during the events we consider here. This is due to strong
surface buoyancy forcing deepening the mixed layer by convection to a depth
beyond what the winds alone could do. While one could include a critical
Richardson number criterion to add NIW breaking to the YBJ equation, we
believe that it is justifiable to ignore dissipation of NIWs in the surface
layer in this study.

In event 1, the wind work predicted by YBJ agreed well with that calculated by
observations. The PM model greatly overestimated the wind-work. The YBJ and PM
models also disagree on the wind work for all the other events (not shown). If
mesoscale interactions change the alignment of the NIWs with the winds, then
this results in a very different wind work. The difference in wind work
between the YBJ and PM models depends on the degree to which the mesoscale
changes the alignment of the waves with the wind and the power of the winds at
different frequencies. The differences in wind work can be large for
individual events but the five events we simulated are not enough to estimate
what the average difference might be when averaged over many events.
Nonetheless, this may provide motivation to revisit previous estimates of the
NIW wind work in order to determine whether the mesoscale may modulate wind
work into the NIW band globally.

##  6 Conclusion

The evolution of the NIW field at the OSMOSIS site in the Northeast Atlantic
Ocean is strongly modulated by mesoscale eddies. The observed evolution can be
captured by the YBJ model, which includes NIW refraction, advection, and
dispersion. If these processes are omitted as in the PM model, the
observations cannot be reproduced as well, even if the parameterized damping
rate is tuned.

The YBJ model provides a powerful interpretive framework for understanding
observations of NIWs in the upper ocean. It allows us to attribute the
observed evolution to specific physical processes. Lateral differences in the
NIW amplitude across the mooring array are caused by $\zeta$ refraction, which
causes NIW kinetic energy to be fluxed into anticyclonic regions. While
correlations between NIW amplitudes and mesoscale vorticity can be diagnosed
from the observations alone (e.g., Yu et al., 2022), the YBJ framework
produces a quantitative prediction for the NIW amplitudes that matches
observations and provides a physical interpretation.

The YBJ model also allows us to calculate changes in the NIW potential energy,
which are expected to arise from transfers of balanced mesoscale energy (Xie
and Vanneste, 2015; Rocha et al., 2018). For the strongest NIW event observed
during the OSMOSIS campaign, the NIW potential energy gain is at least two
orders of magnitude smaller than the global average energy input into
mesoscale kinetic energy. Stimulated generation therefore does not appear to
have a major impact on the mesoscale eddy field in this part of the ocean.

###### Acknowledgements.

This material is based upon work supported by the National Science Foundation
under Grant No. OCE-1924354 and by the National Aeronautics and Space
Administration under Grant No. 80NSSC22K1445 issued through the Science
Mission Directorate (Future Investigators in NASA Earth and Space Science and
Technology). The code to run the 3D YBJ model is available at
https://github.com/scott-conn/3DYBJ. All of the OSMOSIS data used as part of
this study is available from the British Oceanographic Data Centre. The
OSMOSIS mooring data is available at
https://www.bodc.ac.uk/data/bodc_database/nodb/data_collection/6093/. The
OSMOSIS glider data is available at
https://doi.org/10.5285/6cf0b33e-a192-549f-e053-6c86abc01204. The SSH data is
available from the E.U.’s Copernicus Marine Service at
https://doi.org/10.48670/moi-00148. The ERA5 reanalysis data is available from
the Copernicus Climate Change Service (C3S) Climate Data Store at
https://doi.org/0.24381/cds.adbb2d47. [A] YBJ Kinetic-Energy Budget The YBJ
kinetic-energy budget may be formed by multiplying (2) by $-M^{*}/2$ and
adding the complex conjugate. Using integration by parts, the resulting
kinetic-energy budget can be written as

$\partial_{t}\mathcal{K}+J(\psi,\mathcal{K})+\nabla\cdot\bm{F}^{\prime}_{\mathcal{K}}+\partial_{z}G^{\prime}_{\mathcal{K}}=d_{\mathcal{K}}+\gamma_{\mathcal{K}}^{F},$
(14)

where

$\displaystyle\mathcal{K}$ $\displaystyle=\frac{1}{2}|\partial_{z}M|^{2},$
(15a) $\displaystyle\bm{F}^{\prime}_{\mathcal{K}}$
$\displaystyle=\frac{iN^{2}}{4f}\left(M\nabla M^{*}-M^{*}\nabla M\right),$
(15b) $\displaystyle G^{\prime}_{\mathcal{K}}$
$\displaystyle=-\frac{M^{*}}{2f}\left(\partial_{zt}M+J(\psi,\partial_{z}M)+\frac{i\zeta}{2}\partial_{z}M+\nu\nabla^{4}\partial_{z}M-\partial_{z}\mathcal{F}\right)+\text{
c.c. },$ (15c) $\displaystyle d_{\mathcal{K}}$
$\displaystyle=\frac{\nu}{2}\left(\partial_{z}M^{*}\nabla^{4}\partial_{z}M+\partial_{z}M\nabla^{4}\partial_{z}M^{*}\right),$
(15d) $\displaystyle\gamma_{\mathcal{K}}^{F}$
$\displaystyle=\frac{1}{2}\left(\partial_{z}M^{*}\partial_{z}\mathcal{F}+\partial_{z}M\partial_{z}\mathcal{F}^{*}\right).$
(15e)

In this form, however, the vertical flux does not vanish on the boundaries at
$z=0$ and $z=-H$. In order to meaningfully separate horizontal and vertical
fluxes, we must first choose a more sensible gauge for the kinetic-energy
fluxes. Let $\mathbf{H}^{\prime}=[\mathbf{F}^{\prime},G^{\prime}]^{T}$ be the
3D kinetic energy flux vector. Consider a gauge transformation
$\mathbf{H}^{\prime}\rightarrow\mathbf{H}=\mathbf{H}^{\prime}+\nabla_{3}\times\bm{\chi}$
with some vector field $\bm{\chi}$. Because the divergence of a curl is zero,
the kinetic-energy budget is invariant under such a transformation. We propose
to pick a $\bm{\chi}$ such that the transformed $G$ is zero on both
boundaries. To determine $\bm{\chi}$, we first integrate (2) from $-H$ to $z$
and then introduce a new field $A$ such that $\partial_{z}A=N^{2}M/f^{2}$,
which allows us to evaluate the integral of the dispersion term. The resulting
equation is

$\partial_{zt}M+J(\psi,\partial_{z}M)+\frac{if}{2}\nabla^{2}A+\frac{i\zeta}{2}\partial_{z}M-\partial_{z}\mathcal{F}+\nu\nabla^{4}\partial_{z}M=\mathcal{C}$
(16)

with

$\mathcal{C}(x,y,t)=\partial_{zt}M+J(\psi,\partial_{z}M)+\frac{if}{2}\nabla^{2}A+\frac{i\zeta}{2}\partial_{z}M-\partial_{z}\mathcal{F}+\nu\nabla^{4}\partial_{z}M\quad\text{at}\quad
z=-H.$ (17)

We are free to choose $\mathcal{C}$ as a boundary condition on the new field
$A$, i.e., we can choose $A$ at $z=-H$ as the solution of (17) for some
specified $\mathcal{C}$. A sensible gauge is found by choosing $\mathcal{C}=0$
and setting

$\bm{\chi}=\frac{if}{4}\begin{pmatrix}M^{*}\partial_{y}A-M\partial_{y}A^{*}\\\
-M^{*}\partial_{x}A+M\partial_{x}A^{*}\\\ 0\end{pmatrix}.$ (18)

Under this transformation, the kinetic-energy budget becomes (4), with the
terms given by

$\displaystyle\mathbf{F}_{\mathcal{K}}$
$\displaystyle=\frac{if}{4}\left(\partial_{z}M^{*}\nabla A-\partial_{z}M\nabla
A^{*}\right),$ (19a) $\displaystyle G_{\mathcal{K}}$
$\displaystyle=\frac{if}{4}\left(\nabla M\cdot\nabla A^{*}-\nabla
M^{*}\cdot\nabla A\right).$ (19b)

It is clear that the vertical flux $G_{\mathcal{K}}$ vanishes at the
boundaries, as desired, because $\nabla M=0$ there.

The divergence of the flux terms vanishes when the budget is integrated over
the entire domain, and the only terms that can alter the domain-integrated
kinetic energy are the wind work and dissipation:

$\partial_{t}\langle\mathcal{K}\rangle=\mathcal{D}_{\mathcal{K}}+\Gamma_{\mathcal{K}}^{F},$
(20)

where $\mathcal{D}_{\mathcal{K}}=\langle d_{\mathcal{K}}\rangle$ and
$\Gamma_{\mathcal{K}}^{F}=\langle\gamma_{\mathcal{K}}^{F}\rangle$.

[B] YBJ Upper Boundary Condition

Beginning from (16) with the choice $\mathcal{C}=0$, we can vertically
integrate from $z=-H$ to $z=0$ and use $M=0$ at $z=-H$ to arrive at

$\left[\partial_{t}M+J(\psi,M)+\frac{i\zeta}{2}M-\mathcal{F}+\nu\nabla^{4}M\right]_{z=0}=-\frac{if}{2}\nabla^{2}\int_{-H}^{0}A\,dz.$
(21)

The no-normal flow condition is imposed by requiring $\nabla M=0$ at $z=0$
(Young and Ben Jelloul, 1997), which eliminates the advection and dissipation
terms. We then horizontally average (denoted by
$\overline{\,\cdot\,\vphantom{a}}$) equation (21). Because $M(x,y,0,t)$ has no
horizontal structure, it is equal to its horizontal average. On a horizontally
periodic domain, all but two terms vanish in the averaged equation:

$\partial_{t}M(x,y,0,t)=\overline{\mathcal{F}}(0,t).$ (22)

Because the forcing is horizontally uniform in all of our simulations, this
reduces to (13). Note that subtracting (22) from (21) yields a condition on
the integral of $A$:

$\frac{if}{2}\nabla^{2}\int_{-H}^{0}A\,dz=\mathcal{F}^{\prime}-\frac{i\zeta}{2}M(x,y,0,t),$
(23)

where $\mathcal{F}^{\prime}=\mathcal{F}-\overline{\mathcal{F}}$.

[C] YBJ Potential Energy Budget

The YBJ potential-energy budget can be formed by multiplying (2) by
$i\partial_{t}M^{*}/2f$ and adding the complex conjugate. However, a more
transparent derivation begins with dotting the gradient of (16) with $-\nabla
A^{*}/4$ and adding the complex conjugate (setting $\mathcal{C}=0$ as before).
The resulting potential-energy budget is:

$\partial_{t}\mathcal{P}+J(\psi,\mathcal{P})+\nabla\cdot\mathbf{F}_{\mathcal{P}}+\partial_{z}G_{\mathcal{P}}=d_{\mathcal{P}}+\gamma_{\mathcal{P}}^{F}+\gamma_{\mathcal{P}}^{R}+\gamma_{\mathcal{P}}^{A},$
(24)

where

$\displaystyle\mathcal{P}$ $\displaystyle=\frac{N^{2}}{4f^{2}}|\nabla M|^{2},$
(25a) $\displaystyle\mathbf{F}_{\mathcal{P}}$
$\displaystyle=\frac{if}{8}\left[(\nabla^{2}A^{*})\nabla A-(\nabla^{2}A)\nabla
A^{*}\right],$ (25b) $\displaystyle G_{\mathcal{P}}$
$\displaystyle=-\frac{1}{4}\left[\nabla A^{*}\cdot\left(\partial_{t}\nabla
M+J(\psi,\nabla M)\right)+\text{c.c.}\right],$ (25c) $\displaystyle
d_{\mathcal{P}}$ $\displaystyle=\frac{\nu}{4}\left[\nabla
A^{*}\cdot\nabla^{4}\nabla\partial_{z}M+\nabla
A\cdot\nabla^{4}\nabla\partial_{z}M^{*}\right],$ (25d)
$\displaystyle\gamma_{\mathcal{P}}^{F}$
$\displaystyle=-\frac{1}{4}\left[\nabla
A^{*}\cdot\nabla\partial_{z}\mathcal{F}+\nabla
A\cdot\nabla\partial_{z}\mathcal{F}^{*}\right],$ (25e)
$\displaystyle\gamma_{\mathcal{P}}^{R}$ $\displaystyle=\frac{i}{8}\left[\nabla
A^{*}\cdot\nabla(\zeta\partial_{z}M)-\nabla
A\cdot\nabla(\zeta\partial_{z}M^{*})\right],$ (25f)
$\displaystyle\gamma_{\mathcal{P}}^{A}$ $\displaystyle=\frac{1}{4}\left[\nabla
A^{*}\cdot J(\nabla\psi,\partial_{z}M)+\nabla A\cdot
J(\nabla\psi,\partial_{z}M^{*})\right].$ (25g)

By making a plane wave ansatz in the vertical and using integration by parts,
this energy budget can be brought into the same form as the potential-energy
budget in Rocha et al. (2018). In this paper we are primarily concerned with
the domain-integrated potential-energy budget. Under domain integration terms
which differ by a divergence are the same and hence the interpretation of
terms also remains the same. The domain-integrated potential energy budget is
(6), where again $\mathcal{D}_{\mathcal{P}}=\langle d_{\mathcal{P}}\rangle$
and $\Gamma_{\mathcal{P}}^{i}=\langle\gamma_{\mathcal{P}}^{i}\rangle$. In the
case of horizontally uniform forcing, there is no generation of potential
energy by the winds and hence $\gamma_{\mathcal{P}}^{F}=0$.

## References

  * Alford (2001) Alford, M. H., 2001: Internal swell generation: The spatial distribution of energy flux from the wind to mixed layer near-inertial motions. Journal of Physical Oceanography, 31 (8), 2359–2368.
  * Alford (2020) Alford, M. H., 2020: Revisiting near-inertial wind work: Slab models, relative stress, and mixed layer deepening. Journal of Physical Oceanography, 50 (11), 3141–3156.
  * Alford et al. (2016) Alford, M. H., J. A. MacKinnon, H. L. Simmons, and J. D. Nash, 2016: Near-inertial internal gravity waves in the ocean. Annual review of marine science, 8, 95–123.
  * Arbic and Flierl (2004) Arbic, B. K., and G. R. Flierl, 2004: Baroclinically Unstable Geostrophic Turbulence in the Limits of Strong and Weak Bottom Ekman Friction: Application to Midocean Eddies. Journal of Physical Oceanography, 34 (10), 2257–2273, 10.1175/1520-0485(2004)034¡2257:BUGTIT¿2.0.CO;2.
  * Asselin and Young (2020) Asselin, O., and W. R. Young, 2020: Penetration of wind-generated near-inertial waves into a turbulent ocean. Journal of Physical Oceanography, 50 (6), 1699–1716.
  * Balmforth et al. (1998) Balmforth, N., S. G. Llewellyn Smith, and W. Young, 1998: Enhanced dispersion of near-inertial waves in an idealized geostrophic flow. Journal of marine research, 56 (1), 1–40.
  * Balmforth and Young (1999) Balmforth, N. J., and W. R. Young, 1999: Radiative damping of near-inertial oscillations in the mixed layer. Journal of marine research, 57 (4), 561–584.
  * Buckingham et al. (2016) Buckingham, C. E., and Coauthors, 2016: Seasonality of submesoscale flows in the ocean surface boundary layer. Geophysical Research Letters, 43 (5), 2118–2126.
  * Burns et al. (2020) Burns, K. J., G. M. Vasil, J. S. Oishi, D. Lecoanet, and B. P. Brown, 2020: Dedalus: A flexible framework for numerical simulations with spectral methods. Phys. Rev. Res., 2, 023 068.
  * Callies et al. (2020) Callies, J., R. Barkan, and A. N. Garabato, 2020: Time scales of submesoscale flow inferred from a mooring array. Journal of Physical Oceanography, 50 (4), 1065–1086.
  * Damerell et al. (2020) Damerell, G. M., K. J. Heywood, D. Calvert, A. L. Grant, M. J. Bell, and S. E. Belcher, 2020: A comparison of five surface mixed layer models with a year of observations in the north atlantic. Progress in Oceanography, 187, 102 316.
  * Damerell et al. (2016) Damerell, G. M., K. J. Heywood, A. F. Thompson, U. Binetti, and J. Kaiser, 2016: The vertical structure of upper ocean variability at the porcupine abyssal plain during 2012–2013. Journal of Geophysical Research: Oceans, 121 (5), 3075–3089.
  * D’Asaro (1985) D’Asaro, E. A., 1985: The energy flux from the wind to near-inertial motions in the surface mixed layer. Journal of Physical Oceanography, 15 (8), 1043–1059.
  * D’Asaro (1995) D’Asaro, E. A., 1995: Upper-ocean inertial currents forced by a strong storm. part iii: Interaction of inertial currents and mesoscale eddies. Journal of physical oceanography, 25 (11), 2953–2958.
  * D’Asaro et al. (1995) D’Asaro, E. A., C. C. Eriksen, M. D. Levine, P. Niiler, P. Van Meurs, and Coauthors, 1995: Upper-ocean inertial currents forced by a strong storm. part i: Data and comparisons with linear theory. Journal of Physical Oceanography, 25 (11), 2909–2936.
  * Dewar and Flierl (1987) Dewar, W. K., and G. R. Flierl, 1987: Some Effects of the Wind on Rings. Journal of Physical Oceanography, 17 (10), 1653–1667, 10.1175/1520-0485(1987)017¡1653:SEOTWO¿2.0.CO;2.
  * Essink et al. (2022) Essink, S., E. Kunze, R.-C. Lien, R. Inoue, and S.-i. Ito, 2022: Near-inertial wave interactions and turbulence production in a kuroshio anticyclonic eddy. Journal of Physical Oceanography, 52 (11), 2687–2704.
  * Ferrari and Wunsch (2009) Ferrari, R., and C. Wunsch, 2009: Ocean circulation kinetic energy: Reservoirs, sources, and sinks. Annual Review of Fluid Mechanics, 41, 253–282.
  * Flexas et al. (2019) Flexas, M. M., A. F. Thompson, H. S. Torres, P. Klein, J. T. Farrar, H. Zhang, and D. Menemenlis, 2019: Global estimates of the energy transfer from the wind to the ocean, with emphasis on near-inertial oscillations. Journal of Geophysical Research: Oceans, 124 (8), 5723–5746.
  * Gill (1984) Gill, A. E., 1984: On the behavior of internal waves in the wakes of storms. Journal of Physical Oceanography, 14 (7), 1129–1151, 10.1175/1520-0485(1984)014¡1129:OTBOIW¿2.0.CO;2.
  * Gnanadesikan et al. (2015) Gnanadesikan, A., M.-A. Pradal, and R. Abernathey, 2015: Isopycnal mixing by mesoscale eddies significantly impacts oceanic anthropogenic carbon uptake. Geophysical Research Letters, 42 (11), 4249–4255, 10.1002/2015GL064100.
  * Guan et al. (2014) Guan, S., W. Zhao, J. Huthnance, J. Tian, and J. Wang, 2014: Observed upper ocean response to typhoon megi (2010) in the northern south china sea. Journal of Geophysical Research: Oceans, 119 (5), 3134–3157.
  * Hersbach et al. (2018) Hersbach, H., and Coauthors, 2018: Era5 hourly data on single levels from 1979 to present. Copernicus climate change service (c3s) climate data store (cds), 10 (10.24381).
  * Jayne and Marotzke (2002) Jayne, S. R., and J. Marotzke, 2002: The Oceanic Eddy Heat Transport. Journal of Physical Oceanography, 32 (12), 3328–3345, 10.1175/1520-0485(2002)032¡3328:TOEHT¿2.0.CO;2.
  * Jochum et al. (2013) Jochum, M., B. P. Briegleb, G. Danabasoglu, W. G. Large, N. J. Norton, S. R. Jayne, M. H. Alford, and F. O. Bryan, 2013: The impact of oceanic near-inertial waves on climate. Journal of Climate, 26 (9), 2833–2844.
  * Kato and Phillips (1969) Kato, H., and O. Phillips, 1969: On the penetration of a turbulent layer into stratified fluid. Journal of Fluid Mechanics, 37 (4), 643–655.
  * Kunze (1985) Kunze, E., 1985: Near-inertial wave propagation in geostrophic shear. Journal of Physical Oceanography, 15 (5), 544–565.
  * Large and Pond (1981) Large, W., and S. Pond, 1981: Open ocean momentum flux measurements in moderate to strong winds. Journal of physical oceanography, 11 (3), 324–336.
  * Müller et al. (2011) Müller, P., J. McWilliams, and J. Molemaker, 2011: Routes to Dissipation in the Ocean: The 2D/3D Turbulence Conundrum. Marine Turbulence: Theories, Observations, and Models, Cambridge University Press.
  * Nikurashin et al. (2013) Nikurashin, M., G. K. Vallis, and A. Adcroft, 2013: Routes to energy dissipation for geostrophic flows in the Southern Ocean. Nature Geoscience, 6 (1), 48–51, 10.1038/ngeo1657.
  * Plueddemann and Farrar (2006) Plueddemann, A., and J. Farrar, 2006: Observations and models of the energy flux from the wind to mixed-layer inertial currents. Deep Sea Research Part II: Topical Studies in Oceanography, 53 (1-2), 5–30.
  * Pollard and Millard (1970) Pollard, R. T., and R. Millard, 1970: Comparison between observed and simulated wind-generated inertial oscillations. Deep Sea Research and Oceanographic Abstracts, 17, 813–821.
  * Pollard et al. (1973) Pollard, R. T., P. B. Rhines, and R. O. Thompson, 1973: The deepening of the wind-mixed layer. Geophysical Fluid Dynamics, 4 (4), 381–404.
  * Renault et al. (2016) Renault, L., M. J. Molemaker, J. C. McWilliams, A. F. Shchepetkin, F. Lemarié, D. Chelton, S. Illig, and A. Hall, 2016: Modulation of wind work by oceanic current interaction with the atmosphere. Journal of Physical Oceanography, 46 (6), 1685–1704.
  * Rocha et al. (2018) Rocha, C. B., G. L. Wagner, and W. R. Young, 2018: Stimulated generation: extraction of energy from balanced flow by near-inertial waves. Journal of Fluid Mechanics, 847, 417–451, 10.1017/jfm.2018.308.
  * Taburet et al. (2019) Taburet, G., A. Sanchez-Roman, M. Ballarotta, M.-I. Pujol, J.-F. Legeais, F. Fournier, Y. Faugere, and G. Dibarboure, 2019: Duacs dt2018: 25 years of reprocessed sea level altimetry products. Ocean Science, 15 (5), 1207–1224.
  * Thomas et al. (2020) Thomas, L. N., L. Rainville, O. Asselin, W. R. Young, J. Girton, C. B. Whalen, L. Centurioni, and V. Hormann, 2020: Direct observations of near-inertial wave $\zeta$-refraction in a dipole vortex. Geophysical Research Letters, 47 (21), e2020GL090 375.
  * Thompson et al. (2016) Thompson, A. F., A. Lazar, C. Buckingham, A. C. N. Garabato, G. M. Damerell, and K. J. Heywood, 2016: Open-ocean submesoscale motions: A full seasonal cycle of mixed layer instabilities from gliders. Journal of Physical Oceanography, 46 (4), 1285–1307.
  * von Storch and Lüschow (2023) von Storch, J.-S., and V. Lüschow, 2023: Wind power input to ocean near-inertial waves diagnosed from a 5-km global coupled atmosphere-ocean general circulation model. Journal of Geophysical Research: Oceans, 128 (2), e2022JC019 111, 10.1029/2022JC019111.
  * Wagner and Young (2016) Wagner, G., and W. Young, 2016: A three-component model for the coupled evolution of near-inertial waves, quasi-geostrophic flow and the near-inertial second harmonic. Journal of Fluid Mechanics, 802, 806–837.
  * Wunsch and Ferrari (2004) Wunsch, C., and R. Ferrari, 2004: Vertical mixing, energy, and the general circulation of the oceans. Annu. Rev. Fluid Mech., 36, 281–314.
  * Xie and Vanneste (2015) Xie, J.-H., and J. Vanneste, 2015: A generalised-lagrangian-mean model of the interactions between near-inertial waves and mean flow. Journal of Fluid Mechanics, 774, 143–169.
  * Young and Ben Jelloul (1997) Young, W., and M. Ben Jelloul, 1997: Propagation of near-inertial oscillations through a geostrophic flow. JMR, 55, 735–766.
  * Yu et al. (2022) Yu, X., A. C. Naveira Garabato, C. Vic, J. Gula, A. C. Savage, J. Wang, A. F. Waterhouse, and J. A. MacKinnon, 2022: Observed equatorward propagation and chimney effect of near-inertial waves in the midlatitude ocean. Geophysical Research Letters, 49 (13), e2022GL098 522.
  * Zhai et al. (2010) Zhai, X., H. L. Johnson, and D. P. Marshall, 2010: Significant sink of ocean-eddy energy near western boundaries. Nature Geoscience, 3 (9), 608–612, 10.1038/ngeo943.