# Energy transfers in surface wave-averaged equations

###### Abstract

Ocean surface gravity waves play an important role for the air-sea momentum
fluxes and the upper ocean mixing, and knowledge of the sea state leads in
general circulation models to improved estimates of the ocean energy budget
and allows to incorporate surface wave impacts, such as Langmuir turbulence.
However, including the Stokes drift ${\mathbf{u}}^{S}$, in phase-averaged
equations for the Eulerian mean motion leads to an Eulerian energy budget
which is physically difficult to interpret. In this note, we show that a
Lagrangian energy budget allows for a closed energy budget, in which all terms
connecting the different energy compartments correspond to well known energy
transfer terms. We show that the so-called Coriolis-Stokes force does not lead
to an energy transfer between surface gravity waves and oceanic mean motions
as previously suggested. Instead, the Coriolis-Stokes force transfers energy
between the Eulerian mean kinetic energy,
$MKE_{E}=\frac{1}{2}\overline{{\mathbf{u}}}\cdot\overline{{\mathbf{u}}}$, and
a mean energy compartment which is the product of the mean Eulerian velocity
$\overline{{\mathbf{u}}}$ and the mean Stokes drift
$\overline{{\mathbf{u}}^{S}}$,
$MKE_{ES}=\overline{{\mathbf{u}}}\cdot\overline{{\mathbf{u}}^{S}}$. Both
energy forms are a result of the unnatural split-up of the Lagrangian velocity
into Eulerian velocity and the Stokes drift. In an energy budget for the
Lagrangian mean kinetic energy, the work done by the Coriolis-Stokes force
does not contribute, and should be used to estimate the kinetic energy balance
in the wave‐affected surface mixed layer. The Lagrangian energy budget is used
to discuss an energetically consistent framework which can be used to couple a
general circulation ocean model to a surface wave model.

##  1 Intoduction

In a non-rotating frame, Stokes (1847) established that surface gravity waves
induce a mean flow in the direction of wave propagation, known as the Stokes
drift ${\mathbf{u}^{S}}$. Craik and Leibovich (1976) were able to incorporate
${\mathbf{u}^{S}}$ in wave averaged Boussinesq momentum equations for the
Eulerian velocity ${\mathbf{u}}$. Including the Coriolis force, the equations
are given by (e.g. Huang (1979), Leibovich (1980)):

$\partial_{t}{\mathbf{u}}+({\mathbf{u}}\cdot\nabla){\mathbf{u}}=-{\mathbf{f}}\times{\mathbf{u}}-\underbrace{{\mathbf{f}}\times{\mathbf{u}^{S}}}_{F_{CS}}+b{\mathbf{z}}-\nabla
p^{*}+\underbrace{{\mathbf{u}^{S}}\times{\bm{\omega}}}_{F_{V}}+{\mathbf{D}_{u}}$
(1)

Here, $b=-g\rho/\rho_{0}$ denotes buoyancy,
${\bm{\omega}}=\nabla\times{\mathbf{u}}$ is the Eulerian vorticity,
${\mathbf{f}}$ is the Coriolis parameter, and ${\mathbf{D}_{u}}$ indicates
dissipation of Eulerian momentum. Surface wave impacts enter the equation
through ${\mathbf{u}^{S}}$ and are given by the Coriolis-Stokes force
$F_{CS}$, the vortex force $F_{V}$, and the modified pressure
$p^{*}=p/\rho_{0}+\frac{1}{2}[({\mathbf{u}}-{\mathbf{u}^{S}})^{2}-{\mathbf{u}}^{2}]$.
The equations are usually referred to as ’Craik-Leibovich equations’ and are
widely used to study the impact of surface waves on the oceanic surface mixed
layer and Ekman-spiral solution (e.g. Skyllingstad and Denbo (1995),
McWilliams et al. (1997), Polton et al. (2005)).

The Coriolis Stokes force $F_{CS}$ originates from the wave-induced Reynolds
stresses, and leads in an inviscid ocean to an Eulerian flow which is exactly
opposite to the Stokes drift, resulting in a vanishing Lagrangian mean flow
Hasselmann (1970). This result is in agreement with previous findings, that in
a rotating ocean, a Lagrangian mass transport cannot arise from a steady wave
field (Ursell and Deacon (1950), Pollard (1970)). Hasselmann (1970) further
established that the Coriolis Stokes force leads to surface wave driven
inertial oscillations. The vortex force $F_{V}$ causes vorticity to tilt in
the direction of the Stokes drift. The result are coherent vortices known as
Langmuir circulation Craik and Leibovich (1976). The associated Langmuir
turbulence is often a dominant source for turbulent motions and mixing in the
oceanic mixed layer Belcher et al. (2012). Note, that the term ’Langmuir
turbulence’ sometimes include the shear driven turbulence of the Eulerian
return flow, called anti-Stokes flow.

Global and regional climate models do not resolve Langmuir turbulence, however
numerous parameterizations exist and have been tested in such models (e.g. Fan
and Griffies (2014), Ali et al. (2019)). The parameterizations usually rely on
the knowledge of the Stokes drift. Under the assumption of a fully developed
sea, the Stokes drift can be approximated by using the local wind. However,
such a sea state in equilibrium seems to be a rather poor assumption Hanley et
al. (2010). Another possibility is to use a third generation surface wave
model like WAM Komen et al. (1996) or WAVEWATCH III (2016). Several studies
have coupled such a surface wave model to regional or global climate models by
using some form of the Craik-Leibovich equation (e.g. Breivik et al. (2015),
Li et al. (2016), Sun et al. (2022)) or phase-averaged equations of even
higher order in vertical shear (Couvelard et al. (2020)).

The Stokes drift obtained from surface wave models has also been utilized to
estimate the energy input into the ocean. If the Craik-Leibovich equation (Eq.
(1)) is used to form an energy equation for the Eulerian mean kinetic energy,
$MKE_{E}=\frac{1}{2}\overline{{\mathbf{u}}}\cdot\overline{{\mathbf{u}}}$, a
transfer term appears of the form
$\partial_{t}MKE_{E}=-\overline{{\mathbf{u}}}\cdot({\mathbf{f}}\times\overline{{\mathbf{u}^{S}}})$.
The overbar denotes an adequate averaging. The Coriolis force, despite being a
fictitious force, contributes then in this energy equation to the energy
budget which is difficult to understand. The above transfer term is considered
to be an energy transfer from the surface waves to the Eulerian kinetic energy
(e.g. Suzuki and Fox-Kemper (2016), hereafter SFK16) and has been used to
calculate the energy input into the mixed layer (Liu et al. (2009), Sayol et
al. (2016), Zhang et al. (2019), among others). On the global scale Liu et al.
(2009) estimated an energy input of $0.29$ Tw in the Ekman layer through the
work done by the Coriolis-Stokes force, which is a significant share of other
important transfer rates in the global ocean energy cycle. However, at least
to our knowledge, no energy equation for surface waves was derived which shows
the corresponding transfer term of opposite sign. Broström et al. (2014) and
Weber et al. (2015) discussed some inconsistencies in the energy budget
related to the Coriolis-Stokes force. They both conclude that the Coriolis-
Stokes force plays no role in the energy budget, if the budget is integrated
vertically to the moving material surface, i.e. in a Lagrangian framework.

If Eq. (1) is used to form an energy equation for the Eulerian turbulent
kinetic energy,
$TKE_{E}=\frac{1}{2}\overline{{\mathbf{u}^{\prime}}\cdot{\mathbf{u}^{\prime}}}$,
a transfer term appears of the form
$\partial_{t}TKE_{E}=-\overline{{\mathbf{u}^{\prime}}w^{\prime}}\cdot\partial_{z}\overline{{\mathbf{u}^{S}}}$
McWilliams et al. (1997). It is interpreted as a transfer from surface wave
energy into $TKE_{E}$ and originates from the Stokes term in the modified
pressure $p^{*}$. The transfer term was also derived from rapid distortion
theory Teixeira and Belcher (2002) and from generalised Lagrangian mean theory
Ardhuin and Jenkins (2006). Apart from breaking waves, this transfer often
dominates the $TKE$ budget in the ocean mixed layer Belcher et al. (2012) and
is associated with Langmuir circulation, and thus may also play an important
role in the ocean energy cycle. It is the aim of this study to integrate all
such energy transfers into a meaningful and consistent Lagrangian framework.

In Section 2 we discuss the mean kinetic energy equations in this Lagrangian
framework which allow for a closed energy budget with well known energy
transfer terms. Large eddy simulations are used to visualise some important
energy transfers in idealised experiments in Section 3. As the model community
starts to become aware of energy consistency (e.g. Eden et al. (2014)), we
discuss in Section 4, how such a consistent framework could be realised in a
general circulation ocean model coupled to a surface wave model.

##  2 Energy budgets

###  2.1 Mean kinetic energy

The surface-wave averaged Boussinesq momentum Eq. (1) can be rewritten in a
mathematical identical form (SFK16):

$\partial_{t}{\mathbf{u}}+({\mathbf{u}^{L}}\cdot\nabla){\mathbf{u}}+{\mathbf{f}}\times{\mathbf{u}^{L}}=b{\mathbf{z}}-\nabla
p-u^{L}_{i}\nabla u^{S}_{i}+{\mathbf{D}_{u}}$ (2)

Here, ${\mathbf{u}}$ is the wave-averaged Eulerian velocity,
${\mathbf{u}^{L}}={\mathbf{u}}+{\mathbf{u}^{S}}$ the Lagrangian velocity, i.e.
the sum of the Eulerian velocity ${\mathbf{u}}$ plus the Stokes drift
${\mathbf{u}^{S}}$ with components $u^{L}_{i}$ and $u^{S}_{i}$. Surface wave
effects enter the momentum equations via the Stokes drift and modify the
advection and the Coriolis term, i.e. the Coriolis-Stokes force adds to the
”traditional” Coriolis force. The third term on the r.h.s. is the Stokes shear
force and is responsible for Langmuir turbulence. Note, that Einstein
summation convention is used here. For now we consider only molecular
dissipation acting on the Lagarangian motion, so that
${\mathbf{D}_{u}}=\mu{\mathbf{\nabla}^{2}}{\mathbf{u}^{L}}$, with $\mu$ being
the molecular viscosity.

A mean kinetic energy equation for the Eulerian velocity,
$MKE_{E}=\frac{1}{2}\overline{u}_{i}\,\overline{u}_{i}$, can be derived by
multiplying equation (2) by $\cdot\overline{{\mathbf{u}}}$, followed by a
suitable averaging denoted by an overbar. The averaging should satisfy
$\overline{\nabla_{\alpha}A}=\nabla_{\alpha}\overline{A}$ for
$\alpha=(t,x,y,z)$, $\overline{\overline{A}}=\overline{A}$, and
$\overline{\overline{A}B}=\overline{A}\;\overline{B}$ for any quantities $A$
and $B$. Suitable methods are therefore ensemble and horizontal averages. Here
we choose horizontal averaging, as we will later show model results in a
horizontally periodic domain. The velocity ${\mathbf{u}}$ can be therefore
split into ${\mathbf{u}}=\overline{{\mathbf{u}}}+{\mathbf{u}}^{\prime}$, where
the turbulent velocity ${\mathbf{u}}^{\prime}$ is the deviation from the mean
velocity $\overline{{\mathbf{u}}}$. In the same way we split
${\mathbf{u}^{L}}$ and ${\mathbf{u}^{S}}$. As all our velocities are phase-
averaged with respect to surface gravity waves, the above horizontal averaging
relies on the non-trivial assumption that turbulent quantities are de-
correlated from the wave phase. This might be especially problematic with
respect to ${\mathbf{u}}^{S\prime}$, however, we here follow SFK16 and keep
this assumption for now.

$\partial_{t}MKE_{E}$ is then given by:

$\begin{split}\frac{\partial}{\partial
t}{MKE_{E}}+\underbrace{\overline{u^{L}_{j}}\frac{\partial}{\partial
x_{j}}MKE_{E}}_{1}+\underbrace{\vphantom{\frac{\partial K_{ME}}{\partial
x_{j}}}\frac{\partial}{\partial
x_{j}}(\overline{u_{j}}\,\overline{p})}_{2}+\underbrace{\frac{\partial}{\partial
x_{j}}\left(\overline{u_{i}}\,\overline{u_{i}^{\prime}u_{j}^{L\prime}}\right)}_{3}-\underbrace{\frac{\partial}{\partial
x_{j}}\left(\mu\,\overline{u_{i}}\frac{\partial\overline{u^{L}_{i}}}{\partial
x_{j}}\right)}_{4}=\\\ \
\underbrace{\overline{u_{i}^{\prime}u_{j}^{L\prime}}\frac{\partial\overline{u_{i}}}{\partial
x_{j}}}_{5}+\underbrace{\vphantom{\frac{\partial\overline{u_{i}}}{\partial
x_{j}}}\delta_{i,3}\overline{b}\overline{u_{i}}}_{6}-\underbrace{\vphantom{\frac{\partial\overline{u_{i}}}{\partial
x_{j}}}\mu\frac{\partial\overline{u_{i}}}{\partial
x_{j}}\frac{\partial\overline{u^{L}_{i}}}{\partial
x_{j}}}_{7}-\underbrace{\vphantom{\frac{\partial\overline{u}_{i}}{\partial
x_{j}}}\epsilon_{ijk}\overline{u_{i}}f_{j}\overline{u^{S}_{k}}}_{8}-\underbrace{\vphantom{\frac{\partial\overline{u}_{i}}{\partial
x_{j}}}\overline{u}_{i}\overline{u^{L}_{j}}\frac{\partial\overline{u^{S}_{j}}}{\partial
x_{i}}}_{9}-\underbrace{\vphantom{\frac{\partial\overline{u}_{i}}{\partial
x_{j}}}\overline{u}_{i}\overline{u^{L\prime}_{j}\frac{\partial
u^{S\prime}_{j}}{\partial x_{i}}}}_{10}\end{split}$ (3)

Eq. (3) distinguishes transport terms on the l.h.s. from exchange terms on the
r.h.s. The terms 1 to 4 are advection of $MKE_{E}$, work done by pressure,
transport by the Reynolds stresses, and transport by viscous stresses,
respectively. These terms redistribute mean kinetic energy. Note, that here
$MKE_{E}$ is advected by the Lagrangian mean velocity
$\overline{{\mathbf{u}^{L}}}$ and that the Reynolds stresses combine the
deviations from the Eulerian and Lagrangian mean, i.e.
$\overline{u_{i}^{\prime}u_{j}^{L\prime}}$. Term 5 gives the exchange with
turbulent kinetic energy, followed by the exchange with mean potential energy
(term 6). The molecular dissipation of $MKE_{E}$ is given by term 7. The
dissipation is not positive definite, which depends whether the dissipation
acts on the Eulerian or the Lagrangian velocity. The decision has consequences
of either having a non-positive definite dissipation in the Eulerian energy or
helocity budget Holm (1996). We assume that the molecular dissipation acts on
the Lagrangian velocity, which leads to a positive definite dissipation in the
Lagrangian energy budget introduced further below.

The term 8 in Eq. (3) is the work done by the Coriolis-Stokes force. Finally,
we have the somewhat unfamiliar terms 9 and 10 which both originate from the
Stokes shear force. Eq. (3) is first given in a form which is similar to
SFK16, who already derived the terms 8 to 10. These terms are absent in an
energy equation without the presence of surface waves and therefore represent
energy exchanges with surface wave energy (SFK16). To our knowledge, the full
$MKE_{E}$ equation in the current form is here presented for the first time.

Eq. (3) can be reformulated in the more familiar form:

$\begin{split}\frac{\partial}{\partial
t}{MKE_{E}}+\underbrace{\overline{u_{j}}\frac{\partial}{\partial
x_{j}}MKE_{L}}_{T1^{E}}+\underbrace{\vphantom{\frac{\partial K_{ME}}{\partial
x_{j}}}\frac{\partial}{\partial
x_{j}}\left(\overline{u_{j}}\,\overline{p}\right)}_{T2^{E}}+\underbrace{\frac{\partial}{\partial
x_{j}}\left(\overline{u_{i}}\,\overline{u_{i}^{L\prime}u_{j}^{L\prime}}\right)}_{T3^{E}}-\underbrace{\frac{\partial}{\partial
x_{j}}\left(\mu\,\overline{u_{i}}\frac{\partial\overline{u^{L}_{i}}}{\partial
x_{j}}\right)}_{T4^{E}}=\\\ \
\underbrace{\overline{u_{i}^{L\prime}u_{j}^{L\prime}}\frac{\partial\overline{u_{i}}}{\partial
x_{j}}}_{E1^{E}}+\underbrace{\vphantom{\frac{\partial\overline{u}_{i}}{\partial
x_{j}}}\delta_{i,3}\overline{b}\overline{u_{i}}}_{E2^{E}}-\underbrace{\vphantom{\frac{\partial\overline{u_{i}}}{\partial
x_{j}}}\mu\frac{\partial\overline{u_{i}}}{\partial
x_{j}}\frac{\partial\overline{u^{L}_{i}}}{\partial
x_{j}}}_{E3^{E}}-\underbrace{\vphantom{\frac{\partial\overline{u}_{i}}{\partial
x_{j}}}\epsilon_{ijk}\overline{u_{i}}f_{j}\overline{u^{S}_{k}}}_{E^{4}}\end{split}$
(4)

The reformulation allows for an easier interpretation of the individual terms
and is given here for the first time. The labelling of the terms in Eq. (4)
distinguishes transport terms $T$ from exchange terms $E$. The advection term
$T1^{E}$ in Eq. (4) is now expressed as the advection of Lagrangian mean
kinetic energy,
$MKE_{L}=\frac{1}{2}\overline{{\mathbf{u}^{L}}}\cdot\overline{{\mathbf{u}^{L}}}$,
advected by the eulerian velocity $\overline{{\mathbf{u}}}$. The Reynolds
stresses in term $T3^{E}$ and $E1^{E}$ are now the deviations from Lagrangian
mean velocities. The unfamiliar terms 9 and 10 of Eq. (3) are absorbed into
the more familiar terms $T1^{E}$, $T3^{E}$ and $E1^{E}$ of Eq. (4), i.e. in
the latter form they can be interpreted as either transport terms or in case
of $E1^{E}$ as an exchange with turbulent kinetic energy.

Adding $\partial_{t}{\mathbf{u}^{S}}$ on both sides of Eq. (2) allows to write
a tendency equation for the Lagrangian velocity Holm (1996):

$\partial_{t}{\mathbf{u}^{L}}+({\mathbf{u}^{L}}\cdot\nabla){\mathbf{u}^{L}}+{\mathbf{f}}\times{\mathbf{u}^{L}}=b{\mathbf{z}}-\nabla
p-{\mathbf{u}^{L}}\times(\nabla\times{\mathbf{u}^{S}})+{\mathbf{D}_{u}}+\partial_{t}{\mathbf{u}^{S}}$
(5)

The momentum equation Eq. (5) will be the cornerstone, on which our Lagrangian
energy budget and suggested model framework relies. In contrast to Eq. (2),
velocities are fully written in terms of ${\mathbf{u}^{L}}$ and
${\mathbf{u}^{S}}$, i.e. the Eulerian velocity is absent. The term
$-{\mathbf{u}^{L}}\times(\nabla\times{\mathbf{u}^{S}})=-u^{L}_{i}\nabla
u^{S}_{i}+({\mathbf{u}^{L}}\cdot\nabla){\mathbf{u}^{S}}$ contains the Stokes
shear term and the advection of ${\mathbf{u}^{S}}$. A temporal change in the
Stokes drift $\partial_{t}{\mathbf{u}^{S}}$ can be interpreted as a forcing
term for the Lagrangian velocity.

Multiplying Eq. (5) by $\cdot\overline{\mathbf{u}^{L}}$ and averaging leads
then to a tendency equation for the Lagrangian mean kinetic energy
$MKE_{L}=\frac{1}{2}\overline{u^{L}_{i}}\,\overline{u^{L}_{i}}$:

$\begin{split}\frac{\partial}{\partial
t}{MKE_{L}}+\underbrace{\overline{u_{j}^{L}}\frac{\partial}{\partial
x_{j}}MKE_{L}}_{T1^{L}}+\underbrace{\vphantom{\frac{\partial K_{ME}}{\partial
x_{j}}}\frac{\partial}{\partial
x_{j}}\left(\overline{u_{j}^{L}}\,\overline{p}\right)}_{T2^{L}}+\underbrace{\frac{\partial}{\partial
x_{j}}\left(\overline{u_{i}^{L}}\,\overline{u_{i}^{L\prime}u_{j}^{L\prime}}\right)}_{T3^{L}}-\underbrace{\frac{\partial}{\partial
x_{j}}\left(\nu\frac{\partial}{\partial x_{j}}{MKE_{L}}\right)}_{T4^{L}}=\\\ \
\underbrace{\overline{u_{i}^{L\prime}u_{j}^{L\prime}}\frac{\partial\overline{u_{i}^{L}}}{\partial
x_{j}}}_{E1^{L}}+\underbrace{\vphantom{\frac{\partial\overline{u}_{i}}{\partial
x_{j}}}\delta_{i,3}\overline{b}\overline{u_{i}^{L}}}_{E2^{L}}-\underbrace{\vphantom{\frac{\partial\overline{u}_{i}}{\partial
x_{j}}}\mu\frac{\partial\overline{u^{L}_{i}}}{\partial
x_{j}}\frac{\partial\overline{u^{L}_{i}}}{\partial
x_{j}}}_{E3^{L}}+\underbrace{\vphantom{\frac{\partial\overline{u}_{i}}{\partial
x_{j}}}\overline{u^{L}_{i}}\frac{\partial\overline{u^{S}_{i}}}{\partial
t}}_{E7^{L}}\end{split}$ (6)

The labelling of Eq. (6) follows Eq. (4), i.e. the number indicates the
physical interpretation of the term. A short version of Eq. (6), without the
terms $T3^{L}$, $E1^{L}$ and $E2^{L}$, is also given in Holm (1996). We first
notice that apart from the new forcing term $E7^{L}$, Eq. (6) represents the
standard textbook form of a mean kinetic energy equation (e.g. Olbers et al.
(2012), their Eq. 11.63), but here purely expressed using the Lagrangian
velocity ${\mathbf{u}^{L}}$. The Coriolis force is absent in the energy budget
for $MKE_{L}$. We further notice that the exchange term with turbulent
Lagrangian energy, $E1^{L}$, includes
$\overline{u_{i}^{\prime}u_{j}^{\prime}}\partial_{x_{j}}\overline{u_{i}^{S}}$,
which is missing in $E1^{E}$ of Eq. (4).

We propose energy equation Eq. (6) for $MKE_{L}$ to be used for interpretation
and quantification of energy transfers and budget in the surface wave effected
ocean. Eq. (6) will be exploited in section 4, in order to establish an
energetically consistent coupling between a surface wave and ocean model.

$MKE_{L}$ can be split into different energy compartments:

$MKE_{L}=\frac{1}{2}(\overline{{\mathbf{u}}}+\overline{{\mathbf{u}^{S}}})\cdot(\overline{{\mathbf{u}}}+\overline{{\mathbf{u}^{S}}})=MKE_{E}+MKE_{S}+MKE_{ES}$
(7)

$MKE_{S}=\frac{1}{2}\overline{{\mathbf{u}^{S}}}\cdot\overline{{\mathbf{u}^{S}}}$
is the kinetic energy in the Stokes drift, its evolution is solely given by
our prescribed forcing, i.e.
$\partial_{t}MKE_{S}=\overline{{\mathbf{u}^{S}}}\cdot\partial_{t}\overline{{\mathbf{u}^{S}}}$.
$MKE_{ES}=\overline{{\mathbf{u}}}\cdot\overline{{\mathbf{u}^{S}}}$ is some
mixed Eulerian velocity/ Stokes drift energy. By construction $MKE_{ES}$ is
not always positive. A possible negative energy already shows that the split
of $MKE_{L}$ into the different compartments is a purely mathematical
construct which cannot be based on physical arguments. However, we follow this
route for a moment here. A tendency equation for $MKE_{ES}$ can be derived by
multiplying Eq. (2) by $\cdot{\mathbf{u}^{S}}$ followed by averaging, and is
given by

$\begin{split}\frac{\partial}{\partial
t}{MKE_{ES}}+\underbrace{\overline{u_{j}^{S}}\frac{\partial}{\partial
x_{j}}MKE_{L}}_{T1^{ES}}+\underbrace{\vphantom{\frac{\partial K_{ME}}{\partial
x_{j}}}\frac{\partial}{\partial
x_{j}}\left(\overline{u_{j}^{S}}\,\overline{p}\right)}_{T2^{ES}}+\underbrace{\frac{\partial}{\partial
x_{j}}\left(\overline{u_{i}^{S}}\,\overline{u_{i}^{L\prime}u_{j}^{L\prime}}\right)}_{T3^{ES}}-\underbrace{\frac{\partial}{\partial
x_{j}}\left(\mu\,\overline{u^{S}_{i}}\frac{\partial\overline{u^{L}_{i}}}{\partial
x_{j}}\right)}_{T4^{ES}}=\\\ \
\underbrace{\overline{u_{i}^{L\prime}u_{j}^{L\prime}}\frac{\partial\overline{u_{i}^{S}}}{\partial
x_{j}}}_{E1^{ES}}+\underbrace{\vphantom{\frac{\partial\overline{u}_{i}}{\partial
x_{j}}}\delta_{i,3}\overline{b}\overline{u_{i}^{S}}}_{E2^{ES}}-\underbrace{\vphantom{\frac{\partial\overline{u}_{i}}{\partial
x_{j}}}\mu\frac{\partial\overline{u^{S}_{i}}}{\partial
x_{j}}\frac{\partial\overline{u^{L}_{i}}}{\partial
x_{j}}}_{E3^{ES}}-\underbrace{\vphantom{\frac{\partial\overline{u}_{i}}{\partial
x_{j}}}\epsilon_{ijk}\overline{u^{S}_{i}}f_{j}\overline{u_{k}}}_{E4^{ES}}+\underbrace{\vphantom{\frac{\partial\overline{u}_{i}}{\partial
x_{j}}}\overline{u_{i}}\frac{\partial\overline{u^{S}_{i}}}{\partial
t}}_{E7^{ES}}\end{split}$ (8)

By construction, the individual terms of Eq. (8) can be added to the terms of
the tendency equations for $MKE_{E}$ and $MKE_{S}$ to give the tendency
equation for $MKE_{L}$. The triple product of term $E4^{ES}$ in Eq. (8)
corresponds to the triple product of term $E4^{E}$ in Eq. (4) only with
opposite sign, i.e. the work done by the Coriolis-Stokes force exchange energy
between $MKE_{ES}$ and $MKE_{E}$. As both are part of the Lagrangian energy
$MKE_{L}$, the work done by the Coriolis-Stokes force does not appear in the
budget for $MKE_{L}$.

A physical interpretation and/or justification for $MKE_{ES}$ is lacking, and
so is any exchange between $MKE_{ES}$ and $MKE_{E}$. We therefore encourage to
use the the Lagrangian energy $MKE_{L}$ budget, were such an energy transfer
term based on the fictitious Coriolis force is absent.

###  2.2 Turbulent kinetic energy

An equation for the turbulent Lagrangian kinetic energy,
$TKE_{L}=\frac{1}{2}\overline{{\mathbf{u}^{L\prime}}\cdot{\mathbf{u}^{L\prime}}}$
can be derived by multiplying Eq. (5) by $\cdot\mathbf{u}^{L\prime}$ and
averaging:

$\begin{split}\frac{\partial}{\partial
t}{TKE_{L}}+\underbrace{\overline{u_{j}^{L}}\frac{\partial}{\partial
x_{j}}TKE_{L}}_{T1^{L\prime}}+\underbrace{\vphantom{\frac{\partial
K_{ME}}{\partial x_{j}}}\frac{\partial}{\partial
x_{j}}\overline{u_{j}^{L\prime}\,p^{\prime}}}_{T2^{L\prime}}+\underbrace{\frac{\partial}{\partial
x_{j}}\left(\frac{1}{2}\overline{u_{j}^{L\prime}\,u_{i}^{L\prime}u_{i}^{L\prime}}\right)}_{T3^{L\prime}}-\underbrace{\frac{\partial}{\partial
x_{j}}\left(\mu\frac{\partial}{\partial
x_{j}}{TKE_{L}}\right)}_{T4^{L\prime}}=\\\ \
-\underbrace{\overline{u_{i}^{L\prime}u_{j}^{L\prime}}\frac{\partial\overline{u_{i}^{L}}}{\partial
x_{j}}}_{E1^{L\prime}}+\underbrace{\vphantom{\frac{\partial\overline{u}_{i}}{\partial
x_{j}}}\delta_{i,3}\overline{b^{\prime}u_{i}^{\prime
L}}}_{E2^{L\prime}}-\underbrace{\vphantom{\frac{\partial\overline{u}_{i}}{\partial
x_{j}}}\mu\overline{\frac{\partial u^{L\prime}_{i}}{\partial
x_{j}}\frac{\partial u^{L\prime}_{i}}{\partial
x_{j}}}}_{E3^{L\prime}}+\underbrace{\vphantom{\frac{\partial\overline{u}_{i}}{\partial
x_{j}}}\overline{u^{L\prime}_{i}\frac{\partial u^{S\prime}_{i}}{\partial
t}}}_{E7^{L\prime}}\end{split}$ (9)

The tendency equation for $TKE_{L}$ represents again the textbook form of a
turbulent kinetic energy equation (e.g. Olbers et al. (2012), their Eq.
11.62), with the addition of a forcing term $E7^{L\prime}$. The transport
terms $T1^{L\prime}$ to $T4^{L\prime}$ are Lagrangian advection of $TKE_{L}$,
work done by pressure fluctuation, transport by Reynolds-, and by viscous
stresses, respectively. Term $E1^{L\prime}$ represents the exchange term with
$MKE_{L}$ and corresponds to term $E1^{L}$ of Eq. (6). Term $E2^{L\prime}$
gives the exchange with turbulent potential energy, and $E3^{L\prime}$
describes the molecular dissipation of $TKE_{L}$.

Similar to the mean kinetic energy, we can split $TKE_{L}$ into different
compartments , i.e. $TKE_{L}=TKE_{E}+TKE_{S}+TKE_{ES}$.
$TKE_{S}=\frac{1}{2}\overline{{\mathbf{u}^{S\prime}}\cdot{\mathbf{u}^{S\prime}}}$
is the kinetic energy in the Stokes drift fluctuations, its evolution is given
by the forcing term, i.e.
$\partial_{t}TKE_{S}=\overline{{\mathbf{u}^{S\prime}}\cdot\partial_{t}{\mathbf{u}^{S\prime}}}$.
$TKE_{ES}=\overline{{\mathbf{u}^{\prime}}\cdot{\mathbf{u}^{S\prime}}}$ is a
unfamiliar mixed turbulent Eulerian/Stokes energy which is not positive
definite. Finally,
$TKE_{E}=\frac{1}{2}\overline{{\mathbf{u}^{\prime}}\cdot{\mathbf{u}^{\prime}}}$
is the Eulerian turbulent kinetic energy, its evolution is given by:

$\begin{split}\frac{\partial}{\partial
t}{TKE_{E}}+\underbrace{\vphantom{\left(\nu\,\overline{u_{i}^{\prime}\frac{\partial
u^{L\prime}_{i}}{\partial
x_{j}}}\right)}\overline{u_{j}}\frac{\partial}{\partial
x_{j}}TKE_{L}}_{T1^{E\prime}}+\underbrace{\vphantom{\left(\nu\,\overline{u_{i}^{\prime}\frac{\partial
u^{L\prime}_{i}}{\partial x_{j}}}\right)}\frac{\partial}{\partial
x_{j}}\overline{u_{j}^{\prime}\,p^{\prime}}}_{T2^{E\prime}}+\underbrace{\vphantom{\left(\nu\,\overline{u_{i}^{\prime}\frac{\partial
u^{L\prime}_{i}}{\partial x_{j}}}\right)}\frac{\partial}{\partial
x_{j}}\left(\frac{1}{2}\overline{u_{j}^{\prime}\,u_{i}^{L\prime}u_{i}^{L\prime}}\right)}_{T3^{E\prime}}-\underbrace{\frac{\partial}{\partial
x_{j}}\left(\mu\,\overline{u_{i}^{\prime}\frac{\partial
u^{L\prime}_{i}}{\partial x_{j}}}\right)}_{T4^{E\prime}}=\\\ \
-\underbrace{\overline{u_{i}^{\prime}u_{j}^{L\prime}}\frac{\partial\overline{u_{i}^{L}}}{\partial
x_{j}}}_{E1^{E\prime}}+\underbrace{\vphantom{\frac{\partial\overline{u}_{i}}{\partial
x_{j}}}\delta_{i,3}\overline{b^{\prime}u_{i}^{\prime}}}_{E2^{E\prime}}-\underbrace{\vphantom{\frac{\partial\overline{u}_{i}}{\partial
x_{j}}}\mu\overline{\frac{\partial u^{\prime}_{i}}{\partial
x_{j}}\frac{\partial u^{L\prime}_{i}}{\partial
x_{j}}}}_{E3^{E\prime}}-\underbrace{\vphantom{\frac{\partial\overline{u}_{i}}{\partial
x_{j}}}\epsilon_{ijk}\overline{u_{i}^{\prime}f_{j}u^{S\prime}_{k}}}_{E4^{E\prime}}\end{split}$
(10)

The term $E4^{E\prime}$ represents work done by the fluctuating part of the
Coriolis-Stokes force. It does not show up in the $TKE_{L}$ budget, as it
exchanges energy with $TKE_{ES}$ in the same way as its mean part exchanges
energy between $MKE_{E}$ and $MKE_{ES}$. The shear production term
$E1^{E\prime}$ includes the full Lagrangian shear, i.e. it represents only
partly an exchange with $MKE_{E}$ given by the Eulerian shear production, but
additionally exchange energy with $MKE_{ES}$ given by the Stokes shear
production term $E1^{ES}$.

###  2.3 Lagrangian vs. Eulerian energy budget

Figure 1: Schematic of the energy exchanges between different compartments
(boxes) for a) the Lagrangian and b) the Eulerian framework. Blue transfer
terms give transfers within the mechanical energy budget of the ocean. Green
transfer terms give exchanges with surface wave energy through temporal
changes in the Stokes drift, and black transfer terms give exchanges with
internal energy through molecular dissipation (E3 terms). Finally, red
transfer terms give exchanges with the physically difficult to interpret
compartment $MKE_{ES}$. Energy fluxes at the vertical boundaries and lateral
flux divergences are not shown.

In order to simplify the comparison between Lagrangian and Eulerian framework,
we consider the plausible assumption that ${\mathbf{u}^{S\prime}}=0$. The
result is that $TKE_{L}$ and $TKE_{E}$, and their corresponding tendency Eq.
(9) and (10) are identical. Schematics of the vertically integrated energy
budgets which then follow are given in Fig. (1). Our advocated Lagrangian
framework has a closed mechanical energy budget. Compartments and the well
known energy transfer terms are given in black in Fig. (1a). Exchange terms
with other energy compartments are colored in blue. Exchanges with internal
energy through molecular dissipation are given by the terms $E3^{L}$ and
$E3^{L^{\prime}}$ of Eq. (6) and Eq. (9), respectively. The only term, which
exchange energy with ”external” surface wave energy in this framework, is
given by the forcing term
$\int\overline{\mathbf{u^{L}}}\cdot\partial_{t}\overline{\mathbf{u^{S}}}\;dz$.
Note, that the term ”external” surface wave energy is used here, as the
inclusion of the Stokes drift in $MKE_{L}$ represents already some form of
phase-averaged kinetic wave energy. Typically the wind stress, Stokes drift
and near surface Lagrangian velocity are roughly aligned in the same
direction, so that one can expect an overall energy transfer from the waves to
oceanic motions by the forcing term, but exceptions seem possible.

Previous studies have interpreted the Eulerian mechanical energy budget (Fig.
(1b)). The red energy transfer terms were interpreted as an energy exchange
with surface wave energy (e.g. SFK16). However, surface wave energy equations
showing these terms are lacking. A notable exception is the Stokes shear
production term in the wave energy equation of Teixeira and Belcher (2002). By
splitting up $MKE_{L}$ in the different compartments $MKE_{E}$, $MKE_{ES}$,
and $MKE_{S}$, we showed that the red transfer terms can be interpreted as an
exchange of $MKE_{E}$ and $TKE_{E}$ with $MKE_{ES}$. The
$MKE_{ES}=\overline{{\mathbf{u}}}\cdot\overline{{\mathbf{u}^{S}}}$ compartment
is given in red in (Fig. (1b)) because of its dubious physical meaning, for
example, it is not positive definite. As the definition of $MKE_{ES}$ includes
the Eulerian velocity, there seems no reason to consider $MKE_{ES}$ as part of
the surface wave energy. Any physical interpretation of $MKE_{ES}$ or the
associated energy exchanges is very difficult. We therefore recommend to use
the Lagrangian energy budget for quantification and interpretation of the
energy transfers. Using idealized numerical experiments, we will show in the
next section, that a Lagrangian framework can lead to complete different
results in the energy transfers and budgets compared to the Eulerian
framework.

In contrast to the energy transfer associated with the Coriolis-Stokes force,
we do not wish to challenge the well-established interpretation of the Stokes
shear production term
$\overline{u_{i}^{\prime}u_{j}^{\prime}}\partial_{x_{j}}\overline{u_{i}^{S}}$
as an exchange between turbulence and surface wave energy. In the Lagrangian
energy framework this can be interpreted as follows: the Stokes shear
production term removes surface wave energy, which causes a change in the
Stokes drift. The same amount of energy transfer is thus also contained in
${\mathbf{u}^{L}}\cdot\partial_{t}{\mathbf{u}^{S}}$. Note, however, that
changes in the surface wave energy due to wave growth and breaking are
considered to be much larger than due to the Stokes shear production term
Ardhuin et al. (2010).

The energy transfers through wave breaking and surface wave stresses can be
interpreted as vertical boundary conditions for the transport terms labelled
T3 and are not contained in the schematic in Fig. (1). Implications for a
large-scale numerical model framework will be discussed in Section 4.

##  3 Numerical model experiments

Although the results in this study are essentially analytically, we visualise
some of our findings using a numerical model. The model is designed to do
large eddy simulations and is fully three dimensional with lateral cyclic
boundary conditions. At the vertical boundaries, we use no-flux boundary
conditions. A sponge layer damps velocities near the bottom, in order to
prevent reflection of internal waves. The model integrates the surface-wave
averaged Boussinesq momentum given by Eq. (2) and a buoyancy equation of the
form:

$\partial_{t}b+({\mathbf{u}^{L}}\cdot\nabla)b={D_{b}}$ (11)

The sub-grid scale closure follows Ducros et al. (1996), with a turbulent eddy
viscosity $\mu_{t}$ operating on the Eulerian velocity, so that
${\mathbf{D}_{u}}={\mathbf{\nabla}}\mu_{t}{\mathbf{\nabla}}{\mathbf{(}u^{L}-u^{S}})$.
${D_{b}}={\mathbf{\nabla}}\kappa_{t}{\mathbf{\nabla}}{b}$, with a turbulent
eddy diffusivity is $\kappa_{t}=\mu_{t}/Pr$, and a Prandtl number of $Pr=0.7$.

The wave forcing is given by a prescribed Stokes drift which corresponds to
monochromatic, uni-directional deep water wave:

$u^{S}_{eq}=u_{0}^{S}\;\exp\left(\frac{z}{D_{s}}\right)$ (12)

$D_{s}=(2k)^{-1}$ denotes the depth penetration scale, $k$ the wave number,
and $u_{0}^{S}=a^{2}\omega\,k$ the surface Stokes drift, $a$ as the amplitude,
$g$ the gravitational acceleration, and $\omega=\sqrt{gk}$ the frequency. We
choose typical swell conditions with a wave length of $\lambda=100m$, an
amplitude of $a=1.5m$, and a surface Stokes drift in positive x-direction of
$u_{0}^{S}=0.11\frac{m}{s}$, leading to a depth penetration scale of $\sim
8m$. For the growth of swell, we follow the analytical solution of Wagner et
al. (2021):

$u^{S}(z,t)=u^{S}_{eq}(z)\left[1-\exp\left(\frac{-t^{2}}{2T^{2}_{w}}\right)\right]$
(13)

Here, $u^{S}_{eq}(z)$ is the equilibrated Stokes drift as given by Eq. (12),
and $T_{w}$ is a growth time scale, we chose $T_{w}=2h$. The model uses a
rigid-lid and we do not consider a wavy surface. A discussion of the Eulerian
Stokes drift, which lies between the crests and troughs, is given by Broström
et al. (2014).

###  3.1 Laminar flow

The focus of the first experiment is on mean kinetic energies. The model is
initialised with a constant stratification of $N^{2}=5.0\times 10^{-4}s^{-2}$.
Although the Stokes shear term is included, the setup is chosen in a way, that
the model does not generate Langmuir turbulence, i.e. shear and Langmuir
instabilities are not able to overcome the initial stratification. This also
keeps the turbulent viscosity and diffusivity at the lower limit, which
corresponds to molecular friction and diffusion. The results are basically one
dimensional in the vertical, and all velocities can be considered as mean
quantities. The lateral cyclic boundary conditions leads to a vanishing mean
vertical velocity. The Coriolis frequency of
${\mathbf{f}}=f{\mathbf{z}}=7.29\times 10^{-5}s^{-1}$ corresponds to an
inertial period of one day.

In order to see how surface waves drive inertial oscillations theoretically,
we make the following assumptions. We consider a linear, inviscid ocean away
from lateral boundaries, so that no horizontal pressure gradient can be
established on the considered scales. If we choose a horizontal Stokes drift
of the form $u^{S}(z,t)$ the horizontal components of Eq. (5) simplify to

$\displaystyle\partial_{t}{u^{L}}-fv$ $\displaystyle=$
$\displaystyle\partial_{t}{u^{S}}$ (14)
$\displaystyle\partial_{t}{v^{L}}+fu+fu^{S}$ $\displaystyle=$ $\displaystyle
0$ (15)

A time dependent Stokes drift will lead to inertial oscillations as shown by
Hasselmann (1970). In steady state, the mass transport due to surface waves is
exactly zero as the Stokes drift is balanced by the Eulerian anti-Stokes flow,
i.e. $u^{L},v^{L}=0$ and $u^{S}=-u$.

Figure 2: Time evolution of surface velocities: Eulerian zonal velocity
(orange), Lagrangian zonal velocity (blue), Eulerian/Lagrangian meridional
velocity (green), and zonal Stokes drift (red)

Although the non-linear and inviscid assumptions are not made in our model
experiment, the results corresponds to Eqs. (14) and (15). Figure (2) shows
the different velocity compartments at the sea surface. The model is forced by
the arrival of swell, i.e. $\partial_{t}u^{S}$ in the first few hours
according to Eq. (13) with a maximum forcing at $T_{W}=2h$. The forcing leads
to nearly ideal inertial oscillations with vanishing time mean in the
Lagrangian velocities. As the Stokes drift is only in zonal direction, the
meridional component of the Eulerian and Lagrangian velocities are identical.
The amplitude of the inertial oscillations reaches close to the magnitude of
the Stokes drift and depends on the ratio between $T_{w}$ and the inertial
period. If $T_{w}$ goes to zero, the amplitudes of the Stokes drift and
inertial oscillations will be equal. The Eulerian velocity is therefore always
negative for a positive Stokes drift. If we average over one inertial period,
the Eulerian velocity will oppose the Stokes drift $u=-u^{S}$, known as the
anti-Stokes flow. The model therefore reproduces the findings of Hasselmann
(1970), and averaged over an inertial oscillation, the no net mass flux of
Ursell and Deacon (1950).

In our simplified setting the Eulerian and Lagrangian mean kinetic energy Eqs.
(4) and (6) reduce to:

$\displaystyle\partial_{t}{MKE_{E}}$ $\displaystyle=$
$\displaystyle-\overline{{\mathbf{u}}}\cdot
f{\mathbf{z}}\times\overline{{\mathbf{u}^{S}}}$ (16)
$\displaystyle\partial_{t}{MKE_{L}}$ $\displaystyle=$
$\displaystyle\overline{{\mathbf{u}^{L}}}\cdot\partial_{t}\overline{{\mathbf{u}^{S}}}$
(17)

Changes in $MKE_{E}$ are induced by the work done by the Coriolis-Stokes
force. As noticed by Polton (2009), the term ${\mathbf{u}}\cdot
f{\mathbf{z}}\times{\mathbf{u}^{S}}$ is a scalar product between a phase-
averaged velocity and a phase-averaged non-linear momentum term (wave Reynolds
stress), where the latter only gives $f{\mathbf{z}}\times{\mathbf{u}^{S}}$
after phase-averaging. Therefore, the energy equation should be derived before
phase-averaging, as wave correlated terms could give rise to an additional
term in the energy budget. We checked that for our monochromatic wave, to find
out that these contribution can be safely neglected here. In the tendency
equation for the $MKE_{L}$ no Coriolis-Stokes force appears and the evolution
is dependent on the forcing, i.e.
$\overline{{\mathbf{u}^{L}}}\cdot\partial_{t}\overline{{\mathbf{u}}^{S}}$. As
outlined in the previous section $MKE_{L}=MKE_{E}+MKE_{ES}+MKE_{S}$. The two
remaining compartments are given by

$\displaystyle\partial_{t}MKE_{ES}$ $\displaystyle=$
$\displaystyle\overline{{\mathbf{u}}}\cdot
f{\mathbf{z}}\times\overline{{\mathbf{u}^{S}}}+\overline{{\mathbf{u}}}\frac{\partial\overline{{\mathbf{u}^{S}}}}{\partial
t}$ (18) $\displaystyle\partial_{t}MKE_{S}$ $\displaystyle=$
$\displaystyle\overline{{\mathbf{u}^{S}}}\frac{\partial\overline{{\mathbf{u}^{S}}}}{\partial
t}$ (19)

Figure 3: Time evolution of vertically integrated mean kinetic energy
compartments: Eulerian mean kinetic energy $MKE_{E}$ (orange), Lagrangian mean
kinetic energy $MKE_{L}$ (blue), the mixed kinetic energy $MKE_{ES}$ (green),
and kinetic energy of the Stokes drift $MKE_{S}$ (red). See text for details
and definitions.

The vertically integrated energy budgets for the compartments are shown in
Fig. (3) . The budget for $MKE_{E}$ shows strong undulations during an
inertial cycle. Averaged over an inertial period the budget is 0.46
$m^{3}s^{-2}$ and much higher than for the Lagrangian energy $MKE_{L}$, which
is 0.22 $m^{3}s^{-2}$. The exchange between $MKE_{E}$ and $MKE_{ES}$ through
the Coriolis-Stokes term is the dominant signal. The sum of the $MKE_{E}$ and
the $MKE_{ES}$ budget is still slightly negative as it is determined by the
negative contribution of the r.h.s. of Eq. (18). As no further energy exchange
terms are given, the difference in the budget between $MKE_{E}$ and $MKE_{ES}$
is identical to the difference between the two other compartments $MKE_{L}$
and $MKE_{S}$. $MKE_{L}$ gives the energy budget of the inertial oscillations
(see also Fig. (2)) in this experiment, which in a more complex setting will
decrease over time mainly through form stress at the base of the mixed layer
and dissipation (see Czeschel and Eden (2019) and references therein). If all
energy in the inertial oscillations is dissipated, $MKE_{L}$ would be zero and
any movement of fluid particles would stop. However, $MKE_{E}$ would be still
positive in such a steady state, as the Eulerian velocity would exactly oppose
the Stokes drift. The physical interpretation of such an $MKE_{E}$ budget is
difficult, as the energy cannot be transferred to, for example, $TKE$.

Multiplying the integrated mean energies by reference density, e.g.
$\rho_{0}=1000\,kg/m^{3}$, allows a comparison with the total kinetic energy
of the surface waves of $E_{kin}=\frac{1}{4}g\rho_{0}a^{2}=5518\;J/m^{2}$,
with $a=1.5m$ being the amplitude of our prescribed swell. After the initial
forcing period,
$\int\rho_{0}\;MKE_{S}\;dz\approx\int\rho_{0}\;MKE_{L}\;dz\approx
25\;J/m^{2}$, i.e. it is just a small fraction of the total kinetic energy of
the surface waves. $E_{kin}$ accounts for the full orbital motion, whereas
$MKE_{S}$ accounts only for the Stokes drift, i.e. the deviation from a closed
orbital loop. The $MKE_{L}$ budget of $\approx 25\;J/m^{2}$ corresponds to the
energy loss of $E_{kin}$ to inertial oscillation within the first $\approx 4$
hours as given by the r.h.s of Eq. (17). In contrast to the Eulerian energy
budget, the Lagrangian energy budget therefore allows for a clear physical
interpretation of the exchange terms.

###  3.2 Turbulent flow

The impact of turbulence on surface waves and the associated Stokes drift is
largely unknown. Applying the same eddy viscosity on the Stokes drift as used
in mixing parameterizations for the upper ocean Eulerian currents seems not
appropriate, as the energy loss would be much too strong Ardhuin and Jenkins
(2006). A physical explanation for the different impact of turbulence on the
Stokes drift and the Eulerian current might be given by the overlapping time
and spatial scales in wave dynamics and turbulence. For example, the time
scale associated with Langmuir turbulence is often larger than the wave
periods of typical wind waves, i.e. Langmuir turbulence might have no impact
on such waves. The consequence for the energy budget would be that the
exchange terms between $MKE_{L}$ and $TKE_{L}$ in Eqs. (6) and (9), i.e.
$\overline{u_{i}^{L\prime}u_{j}^{L\prime}}\frac{\partial\overline{u_{i}^{L}}}{\partial
x_{j}}$, do not share the same Reynolds stresses as they act on the Eulerian
or the Stokes drift shear. The Lagrangian shear production term should be then
reformulated as
$\overline{u_{i}^{L\prime}u_{j}^{L\prime}}^{E}\frac{\partial\overline{u_{i}}}{\partial
x_{j}}+\overline{u_{i}^{L\prime}u_{j}^{L\prime}}^{S}\frac{\partial\overline{u_{i}^{S}}}{\partial
x_{j}}$ with
$\overline{u_{i}^{L\prime}u_{j}^{L\prime}}^{E}\neq\overline{u_{i}^{L\prime}u_{j}^{L\prime}}^{S}$,
here the different overbars denote different averaging scales. Parameterizing
these Reynolds stresses would then demand different scale-dependent eddy
viscosities.

However, the different Reynolds stresses are difficult to realise in models
using phase-averaged equations, like the Craik-Leibovich equations. Such
models typically prescribe the Stokes drift, and possible impacts of
turbulence on the Stokes drift are neglected. This is achieved by an eddy
viscosity that acts only on the Eulerian velocity and by neglecting the
advection of Stokes drift.

We follow this approach here, and repeat the experiment from the laminar case
but with a uniform mixed layer of 50m ($N^{2}=0\,s^{-2}$) on top of the
stratified interior ($N^{2}=5.0\times 10^{-4}s^{-2}$). We additionally cool
the ocean for 6 hours with 10 $Wm^{-2}$, in order to generate some initial
turbulence. The surface wave forcing follows again Eq. (13) starting at t=0
after the cooling period, i.e. the model starts without mean kinetic energy,
but has a vertically integrated turbulent kinetic energy of
$TKE_{L}=1.73\times 10^{-4}m^{3}s^{-2}$. $TKE_{L}$ equals $TKE_{E}$ in our
set-up, as the Stokes drift is horizontally constant and
${\mathbf{u}^{S\prime}}=0$. The large time and spatial scales of our swell
forcing suggest that all generated turbulence act on the Stokes drift shear.
As we apply lateral cyclic boundary conditions, all vertical mean velocities
are zero, and the evolution of the mean horizontal velocities are governed by

$\displaystyle\partial_{t}\overline{{u^{L}}}-f\overline{v}$ $\displaystyle=$
$\displaystyle\partial_{z}\left(\overline{\mu_{t}\frac{\partial u}{\partial
z}}-\overline{u^{\prime}w^{\prime}}\right)+\partial_{t}\overline{{u^{S}}}$
(20)
$\displaystyle\partial_{t}\overline{{v^{L}}}+f\overline{u}+f\overline{u^{S}}$
$\displaystyle=$
$\displaystyle\partial_{z}\left(\overline{\mu_{t}\frac{\partial v}{\partial
z}}-\overline{v^{\prime}w^{\prime}}\right)$ (21)

here $\mu_{t}$ is the turbulent eddy viscosity acting on the Eulerian
velocity. Note again, that $v^{S}=0$ in our experiment, which also gives
$v^{L}=v$. In steady state, and assuming no-flux boundary conditions at the
surface and somewhere in the stratified interior, only the vertical integrals
of the Stokes drift and the anti-Stokes flow balances, so that
$\int\overline{u^{S}}dz=-\int\overline{u}dz$ and $\int\overline{v}dz=0$. This
should be compared to the laminar case (Eq. 14 and 15), where the steady state
solution was $\overline{u^{S}}={-u}$. Resolved and unresolved turbulence are
therefore shaping the vertical profiles of the Eulerian velocities. The
assumption that turbulence acts solely on the Eulerian velocities has
therefore strong consequences on the vertical gradients of
$\overline{u^{L}},\overline{u}$, and $\overline{u^{S}}$, and therefore also
for the exchange between $MKE_{L}$ and $TKE_{L}$ as given by the Lagrangian
shear production term.

Figure 4:  Vertical profiles of horizontal averaged velocities. The velocities
are averaged over one inertial period, or from model hour 12 to 36.

In our simulation, the vertical shear in the Stokes drift increases within the
first $\sim 6$h due to the growing swell. The result is a burst of turbulence
driven by the developing Langmuir circulation. The $TKE$ budget is very
similar to the findings of Wagner et al. (2021) and is not repeated here.
Although the turbulence is not in equilibrium in such a setup, the ”quasi
steady-state” velocity profiles from our turbulent experiment are given in
Fig. 4. The shaping of the anti-Stokes flow $\overline{u}$ through turbulence
is clearly visible. As predicted from Eq. (20) and Eq. (21) the vertical
integrals of $\overline{u}$ and $\overline{u^{S}}$ cancel each other, and
$\overline{u^{L}}$ and $\overline{v}$ integrate to zero. The differences
between the profiles of $\overline{u^{S}}$ and $(-)\,\overline{u}$ depend on
the amount of turbulence, which is rather weak in our experiment only driven
by swell. For example, Langmuir turbulence driven by wind waves is usually
stronger, and so are the differences between $\overline{u}$ and
$\overline{u^{S}}$. The differences become even larger, if other sources of
turbulence also contribute.

Figure 5: Time evolution of vertically integrated mean kinetic energy tendency
terms in $m^{3}/s^{3}$. Shown are temporal changes in $MKE_{L}$ (blue), the
forcing due to a change in the Stokes drift (orange), the Lagrangian shear
production term (purple) , which consists of the Eulerian production term
(shown in green) and the Stokes production term (red), and the The vector
$\mathbf{u}_{h}=(u,v)$ denotes the horizontal components.

A detailed discussion of Langmuir turbulence driven by swell can be found in
McWilliams et al. (2014) and Wagner et al. (2021) and is not the scope of the
present study. We here concentrate on the impact on $MKE_{L}$. The evolution
in the vertically integrated tendency terms of $MKE_{L}$ are given in Fig. 5.
Similar to the laminar case, the temporal change in $MKE_{L}$ (blue) is
initially given by the forcing term due to temporal changes in the Stokes
drift (orange). As in the laminar case, the term forces surface wave driven
inertial oscillations. Around model hour three, Langmuir turbulence start to
transfer energy from $MKE_{L}$ to $TKE_{L}$ as given by the Lagrangian shear
production term (purple). After model hour 7 the change in $MKE_{L}$ is solely
given by Lagrangian shear production term. The Lagrangian shear production
consists of the Eulerian (green) and Stokes shear production (red). The Stokes
shear production term changes sign after $\sim$ 9 hours. This is because the
Reynolds stresses also changes sign, as they rotate with the Lagrangian mean
flow, which is effected by the inertial oscillation (see McWilliams et al.
(2014) for details). The Eulerian and Stokes shear production show some high
frequent oscillations, which are largely compensated, so that the evolution of
the Lagrangian shear production is much smoother. The compensation points to
the somewhat artificial split-up of the Lagrangian shear production into
Eulerian- and Stokes shear production in such models. The model is not able to
differentiate between the different energy sources. Remember, that the Stokes
shear production is interpreted as a direct energy exchange between surface
wave energy and $TKE$ Teixeira and Belcher (2002).

Figure 6: Time evolution of vertically integrated mean kinetic energy
compartments similar to figure 2 :Eulerian mean kinetic energy $MKE_{E}$
(orange), Lagrangian mean kinetic energy $MKE_{L}$ (blue), the mixed kinetic
energy $MKE_{ES}$ (green), and kinetic energy of the Stokes drift $MKE_{S}$
(red). Note that $MKE_{S}$ is solely based on the Stokes drift forcing which
is identical to the laminar experiment.

Similar to the energy loss to inertial oscillations, the energy loss to $TKE$
is only of minor importance for the energy budget of the surface waves on the
here considered scales, but has strong consequences for the oceanic motions.
On longer time and spatial scales, however, such an energy loss might be an
important contribution to swell dissipation Ardhuin and Jenkins (2006).

The vertically integrated mean kinetic energy compartments are given in Fig.
6. In contrast to the laminar case (Fig. 3), $MKE_{L}$ drops to a value much
lower than the kinetic energy in the Stokes drift ($MKE_{S}$). The reduction
in $MKE_{L}$ is caused by the energy loss to $TKE$. $MKE_{E}$ and $MKE_{ES}$
again exchange energy through the work done by the Coriolis-Stokes force. This
exchange dominates the evolution of the $MKE_{E}$ budget. The loss to $TKE$
through the Eulerian shear production, i.e. term E1 in Eq. (4), is not visible
from this budget. $MKE_{E}$ is again difficult to interpret, as it contains
the largest energy budget, but most of this energy can not be exchanged with
other physically meaningful energy compartments.

##  4 A framework for an energetically consistent coupling of a wave model to
an ocean model

In this chapter we discuss several issues related to a large-scale general
circulation ocean model coupled to a surface wave model. Special focus is
given on consistent energy transfers. We concentrate on deep water waves and
open ocean dynamics.

###  4.1 Model equations

In the suggested framework, the ocean model integrates the Lagrangian velocity
${\mathbf{u}^{L}}={\mathbf{u}}+{\mathbf{u}^{S}}$. The Stokes drift
${\mathbf{u}^{S}}$, and its evolution, will be provided by the wave model. In
the open ocean, the horizontal gradients in the Stokes drift, as given by wave
models, are governed by the atmospheric synoptic scales Haney et al. (2015),
so that
$\partial_{x},\partial_{y}{\mathbf{u}^{L}}>>\partial_{x},\partial_{y}{\mathbf{u}^{S}}$.
The divergence in the Stokes drift is therefore expected to be small, and
allows us to assume $w^{S}=0$, so that any horizontal divergence of the Stokes
drift is compensated by $w$.

Exploiting the assumptions for the open ocean in momentum equation (5) and
buoyancy equation (11) , leads to the following equations to be used in large
scale primitive equation ocean models:

$\partial_{t}{\mathbf{u}_{h}^{L}}+({\mathbf{u}_{h}^{L}}\cdot\nabla){\mathbf{u}^{L}}+f{\mathbf{z}}\times{\mathbf{u}^{L}}=-\nabla_{h}p+\nabla\mu_{t}{\mathbf{\nabla}}({\mathbf{u}^{L}}-{\mathbf{u}^{S}})+\partial_{t}{\mathbf{u}_{h}^{S}}$
(22)

$\partial_{t}b+({\mathbf{u}_{h}^{L}}\cdot\nabla)b=\nabla\kappa_{t}{\mathbf{\nabla}}b$
(23)

Here, the subscript h denotes horizontal vector components, and $\mu_{t}$ and
$\kappa_{t}$ give turbulent viscosity and diffusivity, respectively. The
Lagrangian velocity ${\mathbf{u}^{L}}$ is the only prognostic velocity in the
ocean model. Existing numerical codes for advection, the Coriolis force,
dissipation, and the buoyancy equation can be used. The prognostic model
velocity is then re-interpreted as the Lagrangian velocity $u^{L}$, similar to
the re-interpretation of the model velocity as residual mean velocity
including the Quasi-Stokes velocity in the transformed residual mean theory of
McDougall and McIntosh (2001). Note, however, that dissipation should operate
on ${\mathbf{u}}={\mathbf{u}^{L}}-{\mathbf{u}^{S}}$.

Although it is well known that $\nabla\cdot{\mathbf{u}^{S}}\neq 0$ McIntyre
(1988), we ignore this divergence effect as it is usually quite small Vanneste
and Young (2022). The continuity equation is therefore given by
$\nabla\cdot{\mathbf{u}^{L}}=0$.

As a consequence of the open ocean assumption, the term
$-{\mathbf{u}}^{L}\times(\nabla\times\mathbf{u}^{S})$ in Eq. (5) can be
neglected in the horizontal momentum equations (Eq. (22)). Using scaling
arguments, SFK16 show the possible importance of the term for the vertical
momentum equation. SFK16 suggest to modify the hydrostatic balance in
primitive equation models to a ”wavy hydrostatic balance” of the form,
$\partial_{z}p-b=-u^{L}\partial_{z}u^{S}-v^{L}\partial_{z}v^{S}$. The
necessity of including wave effects in the hydrostatic balance also depends on
the resolved oceanic scales, and needs to be tested upon realization. Note,
that using such a ”wavy hydrostatic balance” would change the budget for
$MKE_{L}$, and we recommend to use the ”standard” hydrostatic balance, i.e.
$\partial_{z}p=b$, for now.

###  4.2 Energy and momentum fluxes

A wave model typically integrates a version of the wave energy balance
equation. For deep water waves in the open ocean, it reads

$\partial_{t}F+\nabla\cdot({\mathbf{v}_{g}}F)=S_{in}+S_{nl}+S_{diss}$ (24)

where $F(\omega,\theta)$ is the two-dimensional wave energy spectrum, which
gives the energy distribution over angular frequency $\omega$ and propagation
direction $\theta$. $\mathbf{v_{g}}$ is the group velocity. The r.h.s. of Eq.
(24) gives the local source terms, which consists of wind input $S_{in}$, non-
linear transfer $S_{nl}$, and dissipation due to wave breaking $S_{diss}$.

The source terms can be utilised to determine the energy and momentum fluxes
between wind, waves and ocean. The momentum flux ${\boldsymbol{\tau}}_{in}$
and energy flux ${\Phi_{in}}$ from wind to the waves are given by (Janssen,
2012):

$\displaystyle{\boldsymbol{\tau}}_{in}$ $\displaystyle=$
$\displaystyle\rho_{w}g\int_{0}^{2\pi}\int_{0}^{\infty}\frac{\mathbf{k}}{\omega}\;S_{in}\;d\omega
d\theta$ (25) $\displaystyle\Phi_{in}$ $\displaystyle=$
$\displaystyle\rho_{w}g\int_{0}^{2\pi}\int_{0}^{\infty}\;S_{in}\;d\omega
d\theta$ (26)

and the fluxes from the waves to the ocean column by:

$\displaystyle{\boldsymbol{\tau}}_{diss}$ $\displaystyle=$
$\displaystyle\rho_{w}g\int_{0}^{2\pi}\int_{0}^{\infty}\frac{\mathbf{k}}{\omega}\;S_{diss}\;d\omega
d\theta$ (27) $\displaystyle\Phi_{diss}$ $\displaystyle=$
$\displaystyle\rho_{w}g\int_{0}^{2\pi}\int_{0}^{\infty}\;S_{diss}\;d\omega
d\theta$ (28)

Here, $\rho_{w}$ is the water density and ${\mathbf{k}}=(k_{x},k_{y})$ the
wave number. Note, that the momentum fluxes are mostly determined through the
high frequency part of the wave spectrum, as they scale with the inverse of
the phase velocity $g=\omega/k$.

The atmospheric or air-side stress is given by
${\boldsymbol{\tau}}_{a}=\rho_{a}\mathbf{u_{*}^{2}}$, here $\rho_{a}$ is the
air density and $\mathbf{u_{*}}$ the air friction velocity. The air-side
stress is usually determined by a drag coefficient $c_{d}$ and the wind speed
in 10m height. As the drag coefficient $c_{d}$ is dependent on the surface
roughness, it should be modified by the sea state, and the wave model can be
used to determine the surface roughness (see e.g. Breivik et al. (2015) for
details).

The ocean side stress $\mathbf{\tau}_{oc}$ can be then considered as the
atmospheric stress minus the residual momentum flux absorbed or released by
the wave field

${\boldsymbol{\tau}}_{oc}={\boldsymbol{\tau}}_{a}-{\boldsymbol{\tau}}_{in}-{\boldsymbol{\tau}}_{diss}$
(29)

If wind increases over a calm ocean, the waves respond first by growing and
${\boldsymbol{\tau}}_{in}>-{\boldsymbol{\tau}}_{diss}$. As the waves mature,
breaking intensifies, and so does the momentum transfer from the waves to the
ocean ${\boldsymbol{\tau}}_{diss}$. At some point during wave growth,
${\boldsymbol{\tau}}_{diss}$ catch up with the momentum transfer from
atmosphere to waves ${\boldsymbol{\tau}}_{in}$. The wave field then is in
equilibrium and ${\boldsymbol{\tau}}_{in}=-{\boldsymbol{\tau}}_{diss}$ Breivik
et al. (2015). At this point the ocean side stress equals the air-side stress
and Eq. (29) reduces to ${\boldsymbol{\tau}}_{oc}={\boldsymbol{\tau}}_{a}$,
which is the assumption made in classical bulk formulas, where
${\boldsymbol{\tau}}_{oc}$ is a function of the wind speed in 10m height.
However, most of the time such an equilibrium is a poor assumption Hanley et
al. (2010). At some point the wind decreases, and the waves have a net
momentum transfer into the ocean
${\boldsymbol{\tau}}_{in}<-{\boldsymbol{\tau}}_{diss}$.

The prognostic frequency range in wave models has an upper limit $\omega_{c}$,
above $\omega_{c}$ the wave spectrum is given by a widely accepted
$\omega^{-5}$ power law. In the high-frequent diagnostic range, we assume that
the wave field is always in equilibrium, so that
${\boldsymbol{\tau}}_{oc}={\boldsymbol{\tau}}_{a}$ for $\omega>\omega_{c}$.
The practical consequence is that Eq. (29) still holds, if we integrate
${\boldsymbol{\tau}}_{in}$ and ${\boldsymbol{\tau}}_{diss}$ over the
prognostic range. A possible physical justification is provided by Chalikov
and Belevich (1993), they argue that the small waves are sheltered by the
large waves and by-pass wave growth, thereby directly driving mean motions.
Such an interpretation is of course an over-simplification of the problem. How
exactly momentum fluxes enter the ocean is subject of active research, but
beyond the scope of the present study.

Similar to the momentum fluxes, the energy flux into the ocean is given by
(Janssen, 2012):

${\Phi}_{oc}={\Phi}_{in}-\rho_{w}g\int_{0}^{2\pi}\int_{0}^{\omega_{c}}\;(S_{in}+S_{diss})\;d\omega
d\theta$ (30)

which can be also written as

${\Phi}_{oc}=\rho_{w}g\int_{0}^{2\pi}\int_{\omega_{c}}^{\infty}\;S_{in}\;d\omega
d\theta\,-\,{\Phi}_{diss}$ (31)

i.e. the energy input consists of the direct energy gain from air in the
diagnostic range and the dissipation of wave energy in the prognostic range,
the latter is mainly the result of white capping. Note, the change in the
limits of integration from Eq. (30) to Eq. (31).

${\Phi}_{oc}$ is the energy transfer from the surface waves available to drive
oceanic mean motions and turbulence. The kinetic energy of the ocean model in
the resolved velocities is governed by the tendency Eq. (6) for $MKE_{L}$. To
allow for a consistent energy transfers to sub-grid scales, we suggest using a
second moment closure, for example a $k-\epsilon$ model. Such a closure
integrates a TKE equation similar to Eq. (9) for $TKE_{L}$, however the
transfer terms are parameterized. For simplicity we assume here a rigid-lid
surface boundary condition, with $w^{L}=0$ and
${\boldsymbol{\tau}}_{oc}=\overline{{\mathbf{u}}^{L\prime}w^{\prime}}$ at
$z=0$. The energy gain of the ocean column is then governed by the transport
terms T3 and the forcing terms E7 in Eqs. (6) and (9). ${\Phi}_{oc}$ should
then match:

${\Phi}_{oc}=\overline{{\mathbf{u}}^{L}}\cdot{\boldsymbol{\tau}}_{oc}+\int\overline{{\mathbf{u}^{L}}}\cdot\partial_{t}\overline{{\mathbf{u}^{S}}}\;dz+\overline{{\mathbf{u}}^{L\prime}\cdot{\boldsymbol{\tau}}^{\prime}_{oc}}+\int\overline{{\mathbf{u}^{L\prime}}\cdot\partial_{t}{\mathbf{u}^{S\prime}}}\;dz+\epsilon_{break}$
(32)

where $\epsilon_{break}$ denotes injection of $TKE$ by breaking waves.

It is unclear, if either
$\overline{{\mathbf{u}}^{L\prime}\cdot{\boldsymbol{\tau}}^{\prime}_{oc}}$ or
$\overline{{\mathbf{u}^{L\prime}}\cdot\partial_{t}{\mathbf{u}^{S\prime}}}$
show any significant correlation in the primed terms, so that they could be
possibly neglected. However, we combine them together with the also unknown
dissipation due to breaking of waves, so that

$\Gamma_{break}=\overline{{\mathbf{u}}^{L\prime}\cdot{\boldsymbol{\tau}}^{\prime}_{oc}}+\int\overline{{\mathbf{u}^{L\prime}}\cdot\partial_{t}{\mathbf{u}^{S\prime}}}\;dz+\epsilon_{break}$
(33)

$\Gamma_{break}$ can be determined in our coupled framework from Eq. (32), as
the remaining terms can be obtained directly from the ocean or the surface
wave model.

Figure 7: Schematic of the energy exchanges between different compartments
(coloured boxes). PE denotes potential energy and $\epsilon$ gives the
integrated dissipation of subgrid $TKE$, i.e. the exchange with internal
energy. Note, that the sign of the terms are associated with $MKE_{L}$ if
possible.

Figure (7) gives an overview of the different energy compartments and the
involved energy transfers. The energy transfer from the wave model to the
ocean is governed by Eq. (32). $\Gamma_{break}$ goes directly into the subgrid
$TKE$, whereas
$\overline{{\mathbf{u}}^{L}}\cdot\overline{{\boldsymbol{\tau}}_{oc}}$ and
$\int\overline{{\mathbf{u}^{L}}}\cdot\partial_{t}\overline{{\mathbf{u}^{S}}}\;dz$
drive $MKE_{L}$, i.e. resolved oceanic motion.

The Stokes production term
$-\int\overline{{\mathbf{u}^{L\prime}}w^{\prime}}\;\frac{\partial\overline{{\mathbf{u}^{S}}}}{\partial
z}\;dz$ could be interpreted as a direct transfer from wave energy to subgrid
$TKE$. In our framework it is part of the Lagrangian shear production term
$\int\overline{{\mathbf{u}^{L\prime}}w^{\prime}}\;\frac{\partial\overline{{\mathbf{u}^{L}}}}{\partial
z}\;dz$, which transfers energy between $MKE_{L}$ and subgrid $TKE$. To be
fully energetically consistent, the shear production term should remove energy
from the surface wave model and should be included in the wave dissipation
term $\Phi_{diss}$ and therefore $\Phi_{oc}$. Ardhuin et al. (2010) in their
recent update on wave dissipation parameterization, discussed the inclusion of
the shear production term. It was, however, neglected, as its contribution is
considered to be very small.

###  4.3 Parameterizations

As mentioned above, we suggest to use a second moment closure for the subgrid
$TKE$. Harcourt (2013) and Harcourt (2015) provide such a closure using the
Craik-Leibovich equations, i.e. it involves a parameterization for the Stokes
production term. The unknown Reynolds stresses are parameterized with two
distinct eddy viscosities $K_{M}$ and $K^{S}_{M}$, so that, e.g.
$\overline{u^{\prime}w^{\prime}}=-K_{M}\partial_{z}u-K^{S}_{M}\partial_{z}u^{S}$.
The eddy viscosities are derived from stability functions (see Harcourt (2013)
for details). It should be however mentioned, that the distinct eddy
viscosities lead to the same Reynolds stresses in the Eulerian and in the
Stokes production term, i.e. it does not resolve possible issues with
overlapping wave and turbulence scales, as discussed in the previous section.

$\Gamma_{break}$ mainly consists of turbulence injected by breaking waves, the
amount of energy is determined by eq. (32). The energy can be injected as a
flux boundary condition in the equations for the second moment closure
following Burchard (2001).

##  5 Summary and discussion

Based on the wave-averaged Craik-Leibovich equations including a prescribed
Stokes drift, we present a closed Lagrangian energy framework. The only energy
exchange of the Lagrangian ocean flow with surface waves is due to changes in
the Stokes drift forcing. Advantages compared to an Eulerian kinetic energy
budget are that all energy transfer terms are well known and easy to
interpret. In particular the work done by the ”fictitious” Coriolis-Stokes
force is absent in the Lagrangian energy budget. Previous studies have
suggested that the work done by the Coriolis-Stokes force is associated with
an energy transfer from the surface waves to the Eulerian kinetic energy. We
argue that this energy gain is an artefact of the split-up of the Lagrangian
kinetic energy into different compartments. The work done by the Coriolis-
Stokes force is an exchange between the Eulerian kinetic energy and an energy
compartment defined as
$MKE_{ES}=\overline{{\mathbf{u}}}\cdot\overline{{\mathbf{u}}^{S}}$. The
individual compartments of the Lagrangian kinetic energy are physically
difficult to interpret, but the compensation between $MKE_{E}$ and $MKE_{ES}$
suggests that large parts of the Eulerian mean kinetic energy are not
available for a transfer to turbulent kinetic energy and finally for mixing.

The ambiguity of the $MKE_{E}$ budget suggests that previous estimates of the
energy input into the ocean by wind and waves should be interpreted with care.
Our Lagrangian framework suggests to ignore the work done by the Stokes-
Coriolis force in such an estimate. The energy input into mean motions by the
wind stress is ${\mathbf{u}}^{L}\cdot\tau$ in our framework. At first glance
this might be very different from previous estimates which used the Eulerian
velocity or even the surface geostrophic velocity. However, a Stokes drift is
usually accompanied by an Eulerian anti-Stokes flow of similar order (at least
vertically integrated), reducing the impact of the Stokes drift.

The compensation of Stokes drift by an anti-Stokes flow is very dependent on
the impact of turbulence on both. Unfortunately, such an impact on surface
waves and the associated Stokes drift is largely unknown, hence ignored in
most phase-averaged models . Nonetheless, the difference between Stokes drift
shear and anti-Stokes flow shear plays an important role in the shear driven
turbulence. The amount of compensation also modifies the surface Lagrangian
velocity, and thus, the energy transfer from surface waves to mean motions.

In the current phase-averaged ocean models, the same Reynolds stresses act on
the Eulerian and the Stokes shear. Overlapping temporal and spatial scales of
turbulence and surface waves suggest that this assumption might not be very
realistic. It might be valid for Swell though. At present, there is no
solution in phase-averaged models, and how turbulence effect surface waves is
a future task for the phase resolving modelling community.

Using our Lagrangian framework, we suggest an energetically consistent
coupling between a surface wave model and a large-scale ocean model. We
recommend to use the Lagrangian velocity as prognostic velocity in the model
as given by Eq. (22). Other forms of the Craik-Leibovich equation are equally
valid, but might demand more changes in an existing numerical code. The energy
transfer from the waves to the ocean can then be split-up into energy which
goes into mean motions and energy which goes into sub-grid turbulence. The
transfer to mean motions consists of the work done by the surface stress, i.e.
$\overline{{\mathbf{u}}^{L}}\cdot\overline{{\boldsymbol{\tau}}_{oc}}$, and
through temporal changes in the Stokes drift,
$\int\overline{{\mathbf{u}^{L}}}\cdot\partial_{t}\overline{{\mathbf{u}^{S}}}\;dz$.
The latter term is expected to be much weaker than the work done by the
surface stress, but is able to force, for example, strong surface wave driven
inertial oscillations. As the overall energy transfer from the surface waves
model to the ocean $\Phi_{oc}$ can be estimated using Eq. (31), the remainder
energy transfer is related to wave breaking, which goes directly into subgrid
turbulence.

The wave models usually allow to incorporate the sea surface velocity to be
included in the second term of Eq. 24, which accounts for wave refraction by
horizontal shear in $\mathbf{u}^{L}|_{z=0}$. The surface velocity could be
also used in the wave model to compute the relative wind
$\mathbf{u}^{atm}_{10}-\mathbf{u}^{L}|_{z=0}$, with $\mathbf{u}^{atm}_{10}$
being the atmospheric absolute wind at 10m. The relative wind rather than the
absolute wind is often used in bulk formulations for the atmospheric stress
$\tau_{a}$. Note, however, that wave models seem very sensitive to the bulk
formulation (see Couvelard et al. (2020) for a detailed discussion).

We here describe a rather simple method to utilize the Lagrangian velocity in
ocean models. More sophisticated methods are suggested, for example using the
generalized Lagrangian mean Ardhuin et al. (2008) or vertically Lagrangian
coordinates Aiki and Greatbatch (2012). These phase-averaging methods follow
the wave motions and allow for a concise treatment of the wavy surface.
However, they are much more difficult to realize.

###### Acknowledgements.

This paper is a contribution to the Collaborative Research Centre TRR 181
Energy Transfer in Atmosphere and Ocean funded by the Deutsche
Forschungsgemeinschaft (DFG, German Research Foundation) Projektnummer
274762653. The numerical code, relevant model data, and scripts are available
on Zenodo with the identifier https://doi.org/10.5281/zenodo.10043904.

## References

  * Aiki and Greatbatch (2012) Aiki, H., and R. J. Greatbatch, 2012: Thickness-weighted mean theory for the effect of surface gravity waves on mean flows in the upper ocean. J. Phys. Oceanogr., 42 (5), 725–747.
  * Ali et al. (2019) Ali, A., K. H. Christensen, Ø. Breivik, M. Malila, R. P. Raj, L. B. an E. P. Chassignet, and M. Bakhoday-Paskyabi, 2019: A comparison of Langmuir turbulence parameterizations and key wave effects in a numerical model of the North Atlantic and Arctic Oceans. Ocean Modelling, 137, 76–97.
  * Ardhuin and Jenkins (2006) Ardhuin, F., and A. D. Jenkins, 2006: On the interaction of surface waves and upper ocean turbulence. J. Phys. Oceanogr., 36 (3), 551–557.
  * Ardhuin et al. (2008) Ardhuin, F., N. Rascle, and K. A. Belibassakis, 2008: Explicit wave-averaged primitive equations using a generalized lagrangian mean. Ocean Modelling, 20 (1), 35–60.
  * Ardhuin et al. (2010) Ardhuin, F., and Coauthors, 2010: Semiempirical dissipation source functions for ocean waves. Part I: Definition, calibration, and validation. J. Phys. Oceanogr., 40 (9), 1917–1941.
  * Belcher et al. (2012) Belcher, S. E., and Coauthors, 2012: A global perspective on Langmuir turbulence in the ocean surface boundary layer. Geophys. Res. Lett., 39 (18).
  * Breivik et al. (2015) Breivik, Ø., K. Mogensen, J.-R. Bidlot, M. A. Balmaseda, and P. A. E. M. Janssen, 2015: Surface wave effects in the NEMO ocean model: Forced and coupled experiments. J. Geophys. Res., Oceans, 120 (4), 2973–2992.
  * Broström et al. (2014) Broström, G., K. H. Christensen, M. Drivdal, and J. E. H. Weber, 2014: Note on Coriolis-Stokes force and energy. Ocean Dynamics, 64, 1039–1045.
  * Burchard (2001) Burchard, H., 2001: Simulating the wave-enhanced layer under breaking surface waves with two-equation turbulence models. J. Phys. Oceanogr., 31 (11), 3133–3145.
  * Chalikov and Belevich (1993) Chalikov, D., and M. Y. Belevich, 1993: One-dimensional theory of the wave boundary layer. Boundary-Layer Meteorology, 63, 65–96.
  * Couvelard et al. (2020) Couvelard, X., F. Lemarié, G. Samson, J.-L. Redelsperger, F. Ardhuin, R. Benshila, and G. Madec, 2020: Development of a two-way-coupled ocean–wave model: assessment on a global nemo (v3. 6)–ww3 (v6. 02) coupled configuration. Geosci. Model Dev., 13 (7), 3067–3090.
  * Craik and Leibovich (1976) Craik, A. D. D., and S. Leibovich, 1976: A rational model for Langmuir circulations. J. Fluid Mech., 73 (03), 401–426.
  * Czeschel and Eden (2019) Czeschel, L., and C. Eden, 2019: Internal wave radiation through surface mixed layer turbulence. J. Phys. Oceanogr., 49 (7), 1827–1844.
  * Ducros et al. (1996) Ducros, F., P. Comte, and M. Lesieur, 1996: Large-eddy simulation of transition to turbulence in a boundary layer developing spatially over a flat plate. J. Fluid Mech., 326, 1–36.
  * Eden et al. (2014) Eden, C., L. Czeschel, and D. Olbers, 2014: Toward energetically consistent ocean models. J. Phys. Oceanogr., 44 (12), 3160–3184.
  * Fan and Griffies (2014) Fan, Y., and S. M. Griffies, 2014: Impacts of parameterized Langmuir turbulence and nonbreaking wave mixing in global climate simulations. J. Climate, 27 (12), 4752–4775.
  * Haney et al. (2015) Haney, S., B. Fox-Kemper, K. Julien, and A. Webb, 2015: Symmetric and geostrophic instabilities in the wave-forced ocean mixed layer. J. Phys. Oceanogr., 45 (12), 3033–3056.
  * Hanley et al. (2010) Hanley, K. E., S. E. Belcher, and P. P. Sullivan, 2010: A global climatology of wind–wave interaction. J. Phys. Oceanogr., 40 (6), 1263–1282.
  * Harcourt (2013) Harcourt, R. R., 2013: A second-moment closure model of Langmuir turbulence. J. Phys. Oceanogr., 43 (4), 673–697.
  * Harcourt (2015) Harcourt, R. R., 2015: An improved second-moment closure model of Langmuir turbulence. J. Phys. Oceanogr., 45 (1), 84–103.
  * Hasselmann (1970) Hasselmann, K., 1970: Wave-driven inertial oscillations. Geophy. and Astrophys. Fluid Dyn., 1 (3-4), 463–502.
  * Holm (1996) Holm, D. D., 1996: The ideal Craik-Leibovich equations. Physica D: Nonlinear Phenomena, 98 (2-4), 415–441.
  * Huang (1979) Huang, N. E., 1979: On surface drift currents in the ocean. J. Fluid Mech., 91 (1), 191–208.
  * Janssen (2012) Janssen, P. A., 2012: Ocean wave effects on the daily cycle in sst. J. Geophys. Res., Oceans, 117 (C11).
  * Komen et al. (1996) Komen, G. J., L. Cavaleri, M. Donelan, K. Hasselmann, S. Hasselmann, and P. Janssen, 1996: Dynamics and modelling of ocean waves. Cambridge University Press.
  * Leibovich (1980) Leibovich, S., 1980: On wave-current interaction theories of Langmuir circulations. J. Fluid Mech., 99 (4), 715–724.
  * Li et al. (2016) Li, Q., A. Webb, B. Fox-Kemper, A. Craig, G. Danabasoglu, W. G. Large, and M. Vertenstein, 2016: Langmuir mixing effects on global climate: WAVEWATCH III in CESM. Ocean Modelling, 103, 145–160.
  * Liu et al. (2009) Liu, B., K. Wu, and C. Guan, 2009: Wind energy input to the ekman-stokes layer: reply to comment by jeff a. polton. Journal of Oceanography, 65, 669–673.
  * McDougall and McIntosh (2001) McDougall, T. J., and P. C. McIntosh, 2001: The temporal-residual-mean velocity. part ii: Isopycnal interpretation and the tracer and momentum equations. J. Phys. Oceanogr., 31 (5), 1222–1246.
  * McIntyre (1988) McIntyre, M., 1988: A note on the divergence effect and the Lagrangian-mean surface elevation in periodic water waves. J. Fluid Mech., 189, 235–242.
  * McWilliams et al. (2014) McWilliams, J. C., E. Huckle, J. Liang, and P. P. Sullivan, 2014: Langmuir turbulence in swell. J. Phys. Oceanogr., 44 (3), 870–890.
  * McWilliams et al. (1997) McWilliams, J. C., P. P. Sullivan, and C.-H. Moeng, 1997: Langmuir turbulence in the ocean. J. Fluid Mech., 334, 1–30.
  * Olbers et al. (2012) Olbers, D., J. Willebrand, and C. Eden, 2012: Ocean dynamics. Springer Science & Business Media.
  * Pollard (1970) Pollard, R. T., 1970: Surface waves with rotation: An exact solution. J. Geophys. Res., 75 (30), 5895–5898.
  * Polton (2009) Polton, J. A., 2009: A wave averaged energy equation: Comment on “global estimates of wind energy input to subinertial motions in the Ekman-Stokes layer” by bin liu, kejian wu and changlong guan. Journal of Oceanography, 65, 665–668.
  * Polton et al. (2005) Polton, J. A., D. M. Lewis, and S. E. Belcher, 2005: The role of wave-induced coriolis–stokes forcing on the wind-driven mixed layer. J. Phys. Oceanogr., 35 (4), 444–457.
  * Sayol et al. (2016) Sayol, J. M., A. Orfila, and L.-Y. Oey, 2016: Wind induced energy–momentum distribution along the Ekman–Stokes layer. application to the Western Mediterranean Sea climate. Deep Sea Research Part I: Oceanographic Research Papers, 111, 34–49.
  * Skyllingstad and Denbo (1995) Skyllingstad, E. D., and D. W. Denbo, 1995: An ocean large-eddy simulation of Langmuir circulations and convection in the surface mixed layer. J. Geophys. Res., 100 (C5), 8501–8522.
  * Stokes (1847) Stokes, G. G., 1847: On the theory of oscillatory waves. Trans. Cam. Philos. Soc., 8, 441–455.
  * Sun et al. (2022) Sun, R. A., and Coauthors, 2022: Waves in SKRIPS: WaveWatch III coupling implementation and a case study of cyclone Mekunu. EGUsphere.
  * Suzuki and Fox-Kemper (2016) Suzuki, N., and B. Fox-Kemper, 2016: Understanding Stokes forces in the wave-averaged equations. J. Geophys. Res., 121 (5), 3579–3596.
  * Teixeira and Belcher (2002) Teixeira, M. A. C., and S. E. Belcher, 2002: On the distortion of turbulence by a progressive surface wave. J. Fluid Mech., 458, 229–267.
  * Ursell and Deacon (1950) Ursell, F., and G. E. R. Deacon, 1950: On the theoretical form of ocean swell. On a rotating earth. Geophysical Journal International, 6, 1–8.
  * Vanneste and Young (2022) Vanneste, J., and W. R. Young, 2022: Stokes drift and its discontents. Philos. Trans. Roy. Soc. London, 380 (2225), 20210 032.
  * Wagner et al. (2021) Wagner, G. L., G. P. Chini, A. Ramadhan, B. Gallet, and R. Ferrari, 2021: Near-inertial waves and turbulence driven by the growth of swell. J. Phys. Oceanogr., 51 (5), 1337–1351.
  * WAVEWATCH III (2016) WAVEWATCH III, R., 2016: Development group (WW3DG): User manual and system documentation of WAVEWATCH III R version 5.16. Tech. Note 329, NOAA/NWS/NCEP/MMAB, College Park, MD, USA.
  * Weber et al. (2015) Weber, J. E. H., M. Drivdal, K. H. Christensen, and G. Broström, 2015: Some aspects of the Coriolis-Stokes forcing in the oceanic momentum and energy budgets. Journal of Geophysical Research: Oceans, 120 (8), 5589–5596.
  * Zhang et al. (2019) Zhang, Y., Z. Song, K. Wu, and Y. Shi, 2019: Influences of random surface waves on the estimates of wind energy input to the Ekman layer in the Antarctic circumpolar current region. Journal of Geophysical Research: Oceans, 124 (5), 3393–3410.