## Abstract

In a series of large-eddy simulations with different forcing, we study the generation of internal gravity waves at the base of the surface mixed layer. If turbulent eddies act as obstacles and undulate the base of the mixed layer, horizontal velocities associated with inertial oscillations and Ekman dynamics can move the obstacles relative to the stratified interior, exciting internal gravity waves similar to lee waves. We find strong evidence that the “obstacle mechanism” is able to excite large parts of the internal wave spectrum, including near inertial waves. The high-frequency part of the excited wave spectrum is filtered by the increased stratification in the transition layer between the mixed layer and lower stratified interior, but a substantial part of the wave spectrum is able to overcome this barrier, hence contributing to interior mixing. The magnitude of the downward-radiated energy below the transition layer depends on the source of turbulence, but we show that the obstacle mechanism, especially under destabilizing heat fluxes, has the potential to contribute considerably to the internal wave energy in the interior ocean.

## 1. Introduction

Breaking of internal waves is a main source of energy for turbulence in the ocean interior. Sources of internal waves include interaction of tidal or balanced flow with bottom topography, loss of balance and wind stress, in particular storms, acting on the surface. Wind stress generated internal waves are often associated with frequencies near the inertial frequency. A prominent mechanism is the so-called “inertial pumping”: temporal fluctuations in the wind stress excite currents in the surface mixed layer which oscillate with the inertial frequency. If these inertial oscillations are horizontally divergent, the vertically fluctuating base of the mixed layer “pumps” near inertial waves in the stratified interior (Price 1983; Gill 1984).

Several studies estimate the flux of near inertial energy into the ocean to lie between 0.34 and 1.4 TW (e.g., Watanabe and Hibiya 2002; Alford 2003). The estimated fraction that leaves the mixed layer and is available for mixing ranges between 5% and 15% (Furuichi et al. 2008; Simmons and Alford 2012; Rimac et al. 2016). The large uncertainties reflect the dependence on details in the model and wind stress products in use. The fate of the remaining ~90% is even less clear as all estimates are based on simple slab or general circulation models, which resolve only a small part of the internal wave spectrum, and where the turbulence in the mixed layer is not resolved but often crudely parameterized.

Besides the inertial pumping, which is considered as the source of internal wave energy in the above cited studies, there are other mechanisms leading to excitement of internal gravity waves. One example is the “mechanical oscillator effect” (Fovell et al. 1992; Ansong and Sutherland 2010): buoyant plumes in the mixed layer overshoot into the adjacent stratified layer where they excite internal gravity waves, a sudden retreat leads to vertical oscillation of the plumes. However, here we focus on the “obstacle mechanism” (Bell 1978; Clark et al. 1986): turbulent eddies undulate the base of the mixed layer and horizontal velocities set the eddies (obstacles) in relative motion to the stratified interior and excite internal gravity waves (Fig. 1). Physically, the mechanism resembles the generation of lee waves with mixed layer eddies instead of topographic obstacles. Regions of high vertical shear in the horizontal velocities where the obstacle mechanism is possibly active includes the Equatorial Undercurrent (Wijesekera and Dillon 1991), tidal velocities in the bottom boundary layer (Gayen et al. 2010), and the inertial oscillations in the surface mixed layer (Bell 1978), which is also the focus of the present study. Polton et al. (2008) find that Langmuir turbulence in the mixed layer generates high-frequency internal gravity waves in the adjacent well stratified transition layer. These waves have frequencies close to the local buoyancy frequency *N* and evolve such that their vector phase velocity matches the depth-averaged mixed layer velocity that rotates as an inertial oscillation, a clear indication of the obstacle mechanism. But as *N* drops below the transition layer, the waves are trapped in the transition layer and might not contribute to interior mixing.

We study the obstacle mechanism using large-eddy simulations (LESs), which are introduced in section 2. In section 3, we present the results beginning with a discussion of the energetics in the mixed layer, while the following subsections show the results for different sources of coherent eddies. We start with eddies generated through shear and convection followed by Langmuir turbulence. Section 4 summarizes and discusses the results.

## 2. The model

The LES we use is based on the nonhydrostatic version of the MIT general circulation model (Marshall et al. 1997). The LES approach demands that the large energetic eddies are resolved by the model grid. The unresolved motions are parameterized by an eddy viscosity which is determined by a subgrid-scale (SGS) model. LES models often do not resolve the stress-carrying small near-wall eddies in wall bounded turbulence. In planetary boundary layers this leads to deviation from Monin–Obukhov similarity theory in the surface layer (Sullivan et al. 1994). To overcome this problem, we apply the two-part eddy viscosity model of Sullivan et al. (1994) where the subgrid momentum fluxes are given by

In Eq. (1) and elsewhere, the Einstein notation is used. The overline denotes a horizontal average. The first term on the right-hand side (rhs) of Eq. (1) refers to the isotropic, fluctuating part of the flow. Parameter $Sij$ is the strain-rate tensor, and $\nu sgs$ is the turbulent eddy viscosity provided by the SGS model as discussed below. The isotropy factor *γ* relates strain of the fluctuations to strain of the mean flow and accounts for the anisotropy in the SGS motions near the wall, that is, *γ* has a value of 1 in the fully isotropic interior and reduces to zero toward the nonisotropic surface boundary. The second term on the rhs of Eq. (1) is the contribution of the mean shear or the nonisotropic part. The mean eddy viscosity $\nu non-isotropic$ is determined so that the vertical shear of the horizontal mean velocities matches Monin–Obukhov similarity theory near the surface (for further details, see Sullivan et al. 1994). Away from the boundary in the mid–mixed layer $\gamma \u21921$ and $\nu non-isotropic\u226a\nu sgs$, following Sullivan et al. (1994), we set $\gamma =1$ and $\nu non-isotropic=0$ for *z* deeper than half the mean mixed layer depth.

To estimate the SGS viscosity $\nu sgs$ and SGS diffusivity $\kappa sgs$, we use the approach of Deardorff (1980), where we cointegrate a prognostic equation for the SGS turbulent kinetic energy. This type of SGS model was successfully applied to the oceanic boundary layer by several studies (e.g., Denbo and Skyllingstad 1996). The SGS turbulent kinetic energy *e* is given by

where $\Delta $ is the averaged grid spacing $\Delta =\u2061(\Delta x\Delta y\Delta z)1/3$. The diagnostic length scale *l* is given by $l=\Delta $ in unstable conditions and $l=Cle1/2N\u22121/2$ for positive values of the Brunt–Väisälä frequency *N*. The dimensionless coefficient *C* is given by $C=C\epsilon l+C\epsilon s\u2061(l/\Delta )$. For all such coefficients we use the values of Deardorff (1980), that is, $Cdg=2.0$, $C\epsilon l=0.19$, $C\epsilon s=0.51$, $Cl=0.76$, $Ckm=0.1$, $Chl=1$, and $Chs=2$. The various terms on the rhs of Eq. (6) are advection, diffusion, shear production *P*, buoyancy *B*, and dissipation *ε*. We follow the approach of Sullivan et al. (1994) and compute the shear production *P* in Eq. (7) solely from the fluctuating velocities.

The model domain is double periodic and spans an area of 180 m × 512 m in the horizontal. The resolution is 1 m in the horizontal and the vertical resolution increases from 0.25 m in the surface layer to 0.5 m below a depth of *z* = −15 m. The applied Coriolis parameter of $f=2\Omega \u2009sin\phi \u22487.27\xd710\u22125\u2009s\u22121$ leads to inertial periods of ≈1 day in the horizontal mixed layer velocities. We use a linear density equation with an initial stratification of *N*^{2} = 1.2 × 10^{−4} s^{−2} below an unstratified surface mixed layer of 30 m. The model has a uniform depth of 120 m. The deepest 20 m consists of a sponge layer, which relaxes temperature to initial values and vertical velocities to zero to prevent reflection of internal waves.

The forcing in all experiments is constant in time and horizontally homogenous. An exception is the wind stress which is calculated with a constant zonal atmospheric wind $U10$, a constant drag coefficient $cd$, and a constant atmospheric density $\rho a$ but modified by the surface velocity, that is,

We include the surface velocity dependence on the wind stress as it is an important damping mechanism of inertial oscillations (Rath et al. 2013). The wind stress in the case of zero surface velocity is $(\tau x,\tau y)=\u2061(0.08,0)$ N m^{−2}. Depending on the experiment, the forcing is a combination of cooling (C), that is, surface heat loss of 100 W m^{−2}, (zonal) wind stress (W), and Langmuir forcing (L), that is, a Stokes drift $S\u2061(z)$ aligned with the wind stress. The Stokes drift velocity is incorporated by using the wave filtered Craik–Leibovich equations (Craik and Leibovich 1976). The vertical profile of the Stokes drift is computed following Sullivan et al. (2007). The Stokes drift at the surface is *U*_{s} = 0.1561 m s^{−1}, which results in a turbulent Langmuir number of $La=U*/Us=0.24$, that is, Langmuir dominates shear turbulence. The depth where turbulent kinetic energy production (TKE) by shear and buoyancy (W and C) are equally important is the Monin–Obukhov length scale. In our setup *L*_{MO} = −36.5 m, which is close to our initial mixed layer depth of *z* = −30 m. A similar Langmuir stability length scale can be defined for the impact of Stokes drift compared to surface buoyancy (L and C) (Belcher et al. 2012). The Langmuir stability length scale *L*_{L} = −254 m in our setup which implies that Langmuir turbulence also dominates the buoyant TKE production. We neglect forcing through breaking surface waves. According to the applied forcing the experiments are called W, C+W, W+L, and C+W+L.

All experiments start after a spinup of 3 h in which initial turbulence develops. During the spinup the wind stress forcing is reduced to 25%, that is, 0.02 N m^{−2}, and no Coriolis force is present. For the (negative) heat flux during the spinup, we apply ~13 W m^{−2} so that the Monin–Obukhov length scale is kept constant over the whole integration. After the spinup we switch on the Coriolis force and increase the wind stress and heat flux to force comparable inertial oscillations in the mixed layer. Note that we did no sensitivity runs to prove that such a spinup is necessary. All experiments are integrated for four inertial periods.

In all experiments, the Ekman transports are approximately nondivergent. Consequently, the inertial pumping mechanism is absent in our model runs and the horizontal mean of the vertical velocity $w\xaf$ vanishes.

## 3. Results

### a. Mixed layer energetics

In this section we discuss the budget of the mixed layer kinetic energy which provides the energy source for the internal gravity waves in our experiments. We split the velocity $u$ into a horizontal mean and deviations from it so that $u=u\xaf+u\u2032$. The associated kinetic energies are the mean kinetic energy $Mke=\u2061(1/2)uj\xaf\u2009uj\xaf$ and the turbulent kinetic energy $Tke=\u2061(1/2)uj\u2032uj\u2032\xaf$. We further split the mean velocity into $u\xaf=\u2329u\u232a+u^$, where the angle brackets denote a time average over one inertial period and deviations from it are marked with a hat. Figures 2a and 2b show the horizontal mean velocities $\u2329u\u232a$ and $\u2329\upsilon \u232a$ averaged over the second inertial period (day 2) after the spinup. The results are qualitatively similar to the solution of McWilliams et al. (2012), who use a similar model setup. The vertical integrals of the mean momentum balance are given by $\u222b\u2329u\u232a\u2009dz=\u2212\u222bus\u2009dz$, $\u222b\u2329\upsilon \u232a\u2009dz=\u2212\tau x/\rho 0f$ (not shown). Here, $us$ is the constant zonal Stokes drift (see Fig. 2a) and $\tau x$ is the constant zonal wind stress. In experiments without Stokes drift, that is, experiment C+W and W, the vertical integral of the zonal momentum is zero, whereas in experiment W+L and C+W+L the vertical integral of the Eulerian zonal velocity compensates the Stokes drift and is called anti-Stokes flow. The vertically integrated meridional transport is equal in all experiments and is in accordance with the classical Ekman transport.

The vertical mixing of horizontal momentum agrees with the scaling arguments of the previous section, that is, a small Langmuir number of La = 0.24 and a Langmuir stability length of *L*_{L} = −245 m. The higher Tke production in the Stokes-driven experiments leads to weaker vertical gradients in the mixed layer. The strongest gradients are found in experiment W, where the Tke production is based exclusively on wind driven shear. At the base of the mixed layer, on the other hand, the vertical shear is rather small in experiment W, which is most evident in the meridional component.

There exist several methods to estimate the depth of the mixed layer *D*. Throughout this study, we define the mixed layer depth as the depth where the deviation from the sea surface temperature exceeds 0.1 K. The difference to the more physical based definition using the minimum turbulent buoyancy flux $b\u2032w\u2032\xaf$ is usually small in our setup, and on the order of the vertical grid resolution. We here use the temperature criterion based definition, when we later use the mixed layer depth *D* as a function of (*x*, *y*, *t*).

The depth of the mixed layer after ~1 day is already quite different among the experiments, as seen by the depths of the sudden momentum decline at the base of the mixed layer (Figs. 2a,b). Averaged over day 2, the mixed layer depth $D=\u2061(\u221234.7,\u221235.67,\u221237.2,\u221237.9)$ for experiments W, W+L, C+W and C+W+L, respectively. The entrainment rate, that is, $\u2212\u2202D/\u2202t$, is larger in experiments including buoyancy forcing (cooling). Langmuir turbulence increases the entrainment and mixed layer deepening far beyond the penetrating depth scale of the Stokes drift (e.g., Polton and Belcher 2007). Consequently the entrainment and mixed layer deepening is strongest in experiment C+W+L followed by experiments C+W, W+L, and W.

The evolution of mean kinetic energy $Mke=\u2061(1/2)uj\xaf\u2009uj\xaf$ assuming horizontal homogeneity is given by

The first term on the rhs is the shear production term and represents the exchange with TKE, which is usually a sink of MKE. The second term appears only in experiments with Stokes forcing and is called “Coriolis–Stokes work.” It represents an exchange between Mke and wave energy (Hasselmann 1970; Suzuki and Fox-Kemper 2016). The third term is the transport of Mke through Reynolds stresses. Finally, $\Gamma $ represents the dissipation by the subgrid scale closure. Inertial oscillations (IOs) are forced either by a change in the wind stress or a change in the Stokes drift (Hasselmann 1970). The dominant balance in the Mke budget of the mixed layer is between work done by the wind and the shear production term for experiment C+W and W. In experiments with Stokes forcing the shear production term is mostly balanced by the Coriolis–Stokes work. We do not show the Mke balances here, but they are in agreement with previous estimates, (e.g., McWilliams et al. 2012). Note that the high Tke production in experiments with Stokes drift is caused by the Stokes shear production, which is an energy conversion from wave energy directly to Tke and does not show up in the Mke equation. An exchange with potential energy is missing in our experiments as we have no mean vertical velocity.

The mixed layer Mke averaged over an inertial period can be split into

The first term on the rhs corresponds to the evolving steady Ekman spiral and is characterized by strong vertical shear within the mixed layer. Under the assumption that the damping scale of inertial oscillations (IOs) is much larger than the inertial period, the second term represents the energy in the IOs, $Mkeio=\u2061(1/2)\u2329u^ju^j\u232a$, which is shown in Fig. 2c. Common to all experiments is that the inertial velocities $u^$ show very weak vertical shear within the mixed layer but strong shear at the base. A plausible explanation for the weak vertical shear is the large time scale separation between turbulence and inertial period. The eddy turnover scale $t*$ can be estimated from the velocity scales involved, that is, the friction velocity $u*=|\tau |/\rho $, the surface Stokes drift $uz=0s$, or Deardorff’s vertical convection scale $w*=\u2061(\u2212B0D)1/3$ (Deardorff 1970), where *D* is the mixed layer depth and $B0$ the surface buoyancy flux. The minimum of each possible combination $t*=D/\u2061[min\u2061(u*,uz=0s,w*)]=3355\u2009s\u22481\u2009h$ is much faster than an inertial period. The energy in the inertial oscillations is forced at the surface with $u^j\tau ^j$. Assuming a constant or slowly varying *τ*, this energy input is oscillating with the inertial period, and $\u2329u^j\tau j\u232a\u22480$. The fast vertical turbulent momentum transfer is able to redistribute the inertial energy homogeneously over the mixed layer. In contrast, the mean energy input, associated with the Ekman dynamic, $\u2329uj\u232a\u2329\tau j\u232a$, is constant in time, and energy loss by shear turbulence is needed to adjust for the constant energy input.

Figure 2c shows that the inertial energy is quite different in the different experiments. Experiment W+L has the most inertial energy, experiment W the least. The reasons for this difference remain unclear, but it seems that the wave forcing is more efficient in driving IOs probably because the “Coriolis–Stokes force” acts as a body force rather than a surface stress. Chaigneau et al. (2008) estimate the energy density of the IOs integrated over the mixed layer from surface drifter, that is, $HKE=\rho \u222b\u2212D0Mkeio\u2009dz$. Their global average of the mean mixed layer inertial velocities $|u^|=9.9\u2009cm\u2009s\u22121$ and the inertial energy density HKE = 286 J m^{−2}. The corresponding values for our experiments range from $|u^|=4.9\u2009cm\u2009s\u22121$ to 5.9 cm s^{−1} and HKE = 43.57–63.95 J m^{−2}. Our HKE, corresponds therefore to values in rather calm regions like the subtropics away from intense storm tracks [see Chaigneau et al. (2008) for global maps].

Applying a temporal average over one inertial period to Eq. (13), we can derive an equation for the inertial energy $Mkeio$:

The first term of the rhs is the shear production term, that is, the exchange between Tke and $Mkeio$. The second term is the Coriolis–Stokes work representing the exchange between $Mkeio$ and the wave energy followed by the transport through Reynold stresses. The last term, that is, $\Gamma io$, represents the dissipation of $MKEio$ by the subgrid closure.

In our setup, the IOs are forced by a sudden increase in $\tau x$ and $uS$ after the spinup. Under the assumption that the inertial damping time scale is much larger than an inertial period, the inertial velocities $u^$ then follow a sine/cosine function. The approximate validity of this assumption can be seen in Fig. 2d after removing the time mean $\u2329u\u232a$ from the full velocity $u\xaf=\u2329u\u232a+u^$. Consequently, after the initial push, the Coriolis–Stokes work and the work done by the wind, that is, $\u2329\u2061(1/\rho 0)\tau ju^j\u232a$ vanish in thecase of constant $\tau x$ and $uS$. In our setup, this only holds for $uS$, as $\tau x$ is dependent on the surface velocity, which leads to damping of IOs (Rath et al. 2013).

As $\u2202u^j/\u2202z$ is small within the mixed layer (Fig. 2c), the damping of inertial energy is concentrated at the boundaries of the mixed layer. The exchange with Tke within the mixed layer is therefore much smaller than the corresponding term in Eq. (13) for the Mke. Similar arguments hold for the transport of $Mkeio$ through Reynolds stresses and the dissipation by the subgrid closure $Dissio$, which are also small within the mixed layer.

Under constant forcing $Mkeio$ is expected to decrease and the solution is heading toward a steady-state Ekman solution. Typical damping time scales of IOs are 4–20 days (e.g., D’Asaro et al. 1995). Bell (1978) suggested a damping mechanism based on the work done by mixed layer motions against the base of the mixed layer, that is, the obstacle mechanism. An equation for the total kinetic energy $K=\u2061(1/2)ujuj$ integrated over the mixed layer is derived in the appendix and is given by

Here, Greek indices denote horizontal indices (1, 2) and $p*=p/\rho o+\u2061(1/2)\u2061(|u\u2212uS|2\u2212u2)$ is the modified pressure. The first terms on the rhs describe the work done by mean horizontal mixed layer motion against the inclined base of the mixed layer. The next two terms are associated with turbulent motions. The fourth term is a turbulent exchange with potential energy, followed by two exchange terms with the wave energy. Parameter $\Gamma K$ is dissipation by the subgrid scale model and *E* is the flux of *K* through the base of the mixed layer due to entrainment. Often the base of the mixed layer is considered as a material surface, that is, the entrainment of energy into the mixed layer is disregarded (e.g., Bell 1978). This is a valid assumption in our experiments as *K* is small below the mixed layer (not shown). Horizontal fluxes of *K* across a sloping mixed layer base are locally possible but should be small in the horizontal mean, as horizontal mass fluxes vanish in the horizontal mean using a rigid-lid surface. However, our results will not depend on the assumption $E\xaf=0$.

The first three terms on the rhs of Eq. (16) describe the loss of energy through the internal wave field, and the first and second can be associated with the obstacle mechanism. Note that these terms do not show up if we vertically integrate an equation for Tke or Mke [Eq. (13), as $\u222b\u2212D0K\u2009dz\xaf\u2260\u222b\u2212D\xaf0K\xaf\u2009dz$]. Nonetheless, we can split the work done by the mean horizontal mixed layer motion at the base of the mixed layer from Eq. (16) into work done by time mean motions and inertial oscillations:

### b. Sheared convection

It is well known that planetary boundary layers under convective forcing form convection cells with narrow downdrafts in the ocean or updrafts in case of the atmosphere (e.g., Schmidt and Schumann 1989). An additionally applied wind stress leads to sheared convection and the downdrafts tend to organize themselves in alignment with the wind stress direction. The pattern changes from a cellular structure in the “free convection” case into convective rolls for “sheared convection” in a wide range of combined shear and buoyancy forcing (Deardorff 1972; Sykes and Henn 1989; Moeng and Sullivan 1994; Heitmann and Backhaus 2005). These rolls are shown in Fig. 3a for experiment C+W in terms of the vertical velocity at *z* = −15 m. Four downdrafts resulting in eight rolls can be clearly identified. The downdrafts are ~25 m thick and ~125 m apart. As the convective rolls are zonally homogenous the problem in experiment C+W becomes nearly two dimensional. A zonally averaged view is given in Fig. 3b. In most cases the downdrafts extend over the whole mixed layer, which is separated from the interior by a transition layer of increased stratification ranging from *z* ≈ −36 to −41 m. The downdrafts and the broader region of upward motion cause a regular undulation of the base of the mixed layer although of varying magnitude.

The increase in wind stress after the spinup leads to IOs, which set the up-/downdrafts in motion relative to the interior. Additionally, the mean meridional velocity increases in time in agreement with Ekman dynamics. Such a vertical shear in the horizontal motions across the base of the mixed layer is a requirement for the obstacle mechanism (see Fig. 1). The IOs start with an amplitude of ~5 cm s^{−1} in experiment C+W. The zonal component $u\xaf$ of the averaged mixed layer velocities follows a sine function and $\upsilon \xaf$ follows a cosine (Fig. 2d). The nearly two dimensional character of the roll vortices allows us to concentrate on the meridional component. Because of the mean southward Ekman component, $\upsilon \xaf$ is always negative in our setup and returns to a value close to zero at $t=2\pi $.

The pressure anomalies $p\u2032$ for two different stages during an inertial period are shown in Fig. 4. In the mixed layer, that is, above *z* = −40 m, the downdraft are located at the strong negative meridional gradients of $p\u2032$, for example, at *x* = 40, 220, 400, and 480 m in Fig. 4a. Below the mixed layer, pressure anomalies are inclined toward the horizontal. The $p\u2032$ contours represent here lines of constant phases of internal gravity waves, propagating in meridional direction. The vertical phase (energy) propagation is upward (downward) as seen by the Hovmöller diagrams of $p\u2032$ for $t=2\pi $ and *π* at a fixed location (Fig. 5).

Under the assumption that the obstacle mechanism is generating the internal waves seen in Fig. 4, we can use linear theory to anticipate the wave characteristics analytically. The frequency *ω* is given by

Here, only the meridional component of $u\xaf|\u2212D$ is relevant, which evolves like the mean mixed layer velocity $\upsilon \xaf$ (Fig. 2d). Parameter $\upsilon \xaf|\u2212D$ is related to the horizontal phase/group speed of the waves. Parameter $kh$ denotes the horizontal wavenumber, which can be obtained by an approximate wavelength of $\lambda =Ly/2No$, where $No$ is the number of roll vortices over a distance $Ly$. The dispersion relation of internal waves is given by

where *ϕ* is the angle between the wavenumber $k$ and the horizontal plane. The theoretically predicted wavenumber vectors using Eqs. (18) and (19) are also shown in Fig. 4 and agree well with the simulated phase propagation. The frequency *ω* and angle *ϕ* change during an inertial period from $t=0$ to $2\pi $. The angle *ϕ* has a minimum value of $\varphi =76.3\xb0$ at $t=\pi $ (Fig. 4a), which corresponds to a maximum $\upsilon \xaf|\u2212D$. A minimum in *ϕ* is equivalent to a maximum in vertical group velocity which is $cgz\u224840\u2009m h\u22121$ (Fig. 5a). If $\upsilon \xaf|\u2212D$ slows down for $t>\pi $, *ϕ* increases. Figure 4b shows the solution at $t=\pi /2$ when the $\upsilon \xaf|\u2212D$ has a medium strength. Here, linear theory estimates an angle of $\varphi =84.3\xb0$. Internal wave generation is bounded by the frequency limits *N* and *f*, that is, in this setup by

Note that the moderate $\upsilon \xaf|\u2212D$ in our experiment does not reach the upper limit related to *N* but does approach the lower limit related to f during $t=2\pi $, where the troughs and crests of the internal waves are nearly horizontal (Fig. 5b).

The horizontal wavelength $\lambda h$ of the internal waves seem to be set by the number of roll vortices $No$ over $Ly$, which would be in agreement with a generation through the obstacle mechanism. Previous studies show that the number of roll vortices and the spacing between them show some correlation with the mixed layer depth and the shear close to the boundaries (Sykes and Henn 1989; Polton et al. 2008). In our simulations the spacing between the downdrafts and the intensity within the downdrafts varies in time, they also sometimes merge or split but no long-term tendency is detectable. Horizontal wavenumber $kh$ seems to be additionally related to the vertical shear of the meridional velocity in the mixed layer, that is, the tilting of the roll vortices. This might also explain why $No$ switches between six and eight during an inertial period in experiment C+W. In Fig. 4b, $No$ is only six because of the stronger vertical shear of $\upsilon \xaf$ compared to Fig. 4a where $No$ = 8. Nonetheless, the number of clearly detectable downdrafts goes not below three, that is, below six roll vortices, during the integration of experiment C+W.

In our model setup, all motions in the interior below the transition layer can be attributed to internal waves. Their energy source is the kinetic energy in the mixed layer as given by, for example, Eq. (16). The vertical flux of internal wave energy is given by

Note that $p*\u2032=p\u2032/\rho 0$ below the penetration depth scale of the Stokes drift, which is ~15 m. The transfer of energy from the mixed layer motions to the internal wave field is given by the first three terms on rhs of Eq. (16). Quantitative results using these transfer terms depend on the definition of the mixed layer depth. However, in the model and in the real ocean there might be no sharp interface between the surface layer and the strongly stratified transition layer and the definition of the mixed layer depth *D* becomes ambiguous.

In our analysis we define the mixed layer depth as the depth where the deviation from the sea surface temperature exceeds 0.1 K. Figure 6 shows the transfer of mixed layer *K* into the internal wave field at the base of the mixed layer from the first term on the rhs of Eq. (16) compared to the vertical flux as in Eq. (21) at *z* = −50 m, which is well below the mixed and also below the transition layer. The general good agreement suggests that the internal waves in the interior are driven by the mean mixed layer motions. The averaged energy flux is downward as the internal waves will be damped by the sponge layer in the lowermost grid cells thereby preventing reflection. The energy flux is strongest around $t=\pi $ when $\upsilon \xaf|\u2212D$ has its maximum, *ϕ* its minimum, and $cgz$ its maximum during an inertial period. Note that $Eflux=cgzE$, with *E* being the total internal wave energy. At $t=2\pi $, the $Eflux$ is close to zero as $\upsilon \xaf|\u2212D\u22480$ and the generated waves approaching frequencies close to *f* with vanishing vertical group velocity (see also Fig. 5b). The internal or interfacial waves at the base of the mixed layer give some of the energy back to the mixed layer motions around $t=2\pi $.

Averaged over the first 4 days $u\alpha \xafp\u2032\u2061(\u2202D\u2032/\u2202x\alpha )\xaf=\u22121.84\xd710\u22125\u2009N\u2009m\u22121\u2009s\u22121$. Splitting the result according to Eq. (17) into the transfer associated with the Ekman component $\u2329u\alpha \xaf\u232a\u2329p\u2032\u2061(\u2202D\u2032/\u2202x\alpha )\xaf\u232a=\u22121.08\xd710\u22125\u2009N\u2009m\u22121\u2009s\u22121$ and a part associated with the inertial oscillations $\u2329u\alpha \xaf^p\u2032\u2061(\u2202D\u2032/\u2202x\alpha )\xaf^\u232a=\u22120.76\xd710\u22125\u2009N\u2009m\u22121\u2009s\u22121$ shows that the contribution of the Ekman component exceeds the contribution of the IOs in our experiment.

The difference between the $w\u2032p\u2032\xaf$ at *z* = −50 m and the transfer of *K* into the internal wave field from the first on the rhs of Eq. (16) is the result of turbulent motions at the base of the mixed layer. The turbulent motions are highly intermittent and sometimes they reach far into the transition layer. The contribution of the turbulent horizontal motions, that is, $\u2329u\alpha \u2032p\u2032\u2061(\u2202D\u2032/\u2202x\alpha )\xaf\u232a=\u22127.2\xd710\u22127\u2009N\u2009m\u22121\u2009s\u22121$ averaged over the first 4 days, is more than an order of magnitude smaller.

The third term on the rhs of Eq. (22), $\u2329w\u2032p\u2032\xaf\u232a=\u22121.64\xd710\u22125\u2009N\u2009m\u22121\u2009s\u22121$, averaged over the first 4 days contributes considerably to the energy transfer. Possible mechanisms contributing to the turbulent transfer into the internal wave field include Kelvin–Helmholtz instabilities and the mechanical oscillator effect (Fovell et al. 1992; Ansong and Sutherland 2010). Common to all these driving mechanism is that the resulting internal waves have high frequencies, that is, close to the adjacent stratification frequency *N*. Entrainment at the base of the mixed layer leads to an increasing *N* and the formation of a higher stratified transition layer. As *N* gives the upper limit of the internal wave branch, the transition layer acts as a low pass filter and large parts of the high-frequency internal wave energy might not reach the interior. Note that *N*^{2} increases from the initial *N*^{2} = 1.2 × 10^{−4} s^{−2} to maximum values of *N*^{2} = 4.0 × 10^{−4} s^{−2} after 3 days. The lack of a transition layer at the beginning of our integration might also to explain the initial high fluxes of internal wave energy at *z* = −50 m on the first half of day 1.

Nonetheless, some high-frequency internal wave energy reaches the interior also with an existing transition layer. In fact, the internal waves in the interior consists of two distinct frequency bands. Figure 7 shows a section at *x* = 80 m of the pressure anomaly and vertical velocity field. The lower-frequency wave field manifests itself in the pressure anomalies. The horizontal wavelength ranges from *λ*_{h} ≈ 120 to 180 m and the angle between the wavenumber $k$ and the horizontal plane ranges from *ϕ* ≈ 72° to 90° during an inertial oscillation. Much higher frequencies show up in the vertical velocity with *λ*_{h} ≈ 35 m and *ϕ* ≈ 35°–40°. From the polarization equations of internal waves (e.g., Sutherland 2010) the ratio of the amplitudes of pressure anomalies $Ap$ to vertical velocities $Aw$ is given by

High-frequency waves with small *ϕ* are therefore dominating the vertical velocities and lower-frequency waves the pressure (and also horizontal velocity) anomalies. The relatively small high-frequency wave band is found in many modeling and laboratory studies were a stratified interior is adjacent to a turbulent boundary layer (e.g., Sutherland and Linden 1998; Taylor and Sarkar 2007; Dohan and Sutherland 2003; Polton et al. 2008; Gayen et al. 2010). The reason for the narrow band might be that waves with an angle of *ϕ* = 35° (*ϕ* = 45°) are most efficient in transporting energy (momentum) away from the turbulent region (Dohan and Sutherland 2003).

The power spectra of the vertical velocities and spectra of the vertical energy flux at *z* = −50 m are shown in Fig. 8. Figures 8a and 8c give the variance conserving form showing that most of the energy is contained in the meridional propagating waves, that is, normal to the roll vortices. The differences between the thick and thin contours give the sudden drop of energy especially in the vertical velocities after day 1 caused at least partly by the emerging transition layer. The energy-rich wavenumbers can be identified in the corresponding spectra which are not scaled by the wavenumber (Figs. 8c,d). On the first day the energy in the vertical velocities peak at a wavelength of 42 m in agreement with Fig. 7b. The average over the following 2 days show the strong reduction in energy but beside the peak at the short waves another maximum with wavelength from 128 to 170 m show up. These large wavelength dominate the vertical energy flux (Figs. 8d,7a). A wavelength of 170 m, for example, is the result of 6 roll vortices or 3 downdrafts as shown in Fig. 7a (also compare schematic in Fig. 1). As discussed above the number of roll vortices switches between 6 and 8, but we also find peaks in the lower harmonics as often one roll vortex is more dominant than the other. For example the peak at the first day in Fig. 8d gives a wavelength of 256 m.

### c. Langmuir turbulence

Another type of coherent eddies in the mixed layer is associated with Langmuir turbulence. Langmuir turbulence is the result of an interaction between surface waves and wind driven currents. The resulting Langmuir cells are aligned between the direction of the wind and the surface waves (Skyllingstad and Denbo 1995; McWilliams et al. 1997). These “classical” Langmuir rolls have depth scales associated with the surface wave-driven Stokes drift, which is ~12 m and well above the mixed layer depth in our setup. Figures 9a and 9b show the Langmuir rolls by means of the vertical velocity at *z* = −5 m for experiment W+L and C+W+L. The downdrafts are ~10 m thick and the rolls have a horizontal scale of ~30 m. At these depths the turbulence is mainly driven by the Stokes shear production, that is, $uh\u2032w\u2032\xaf\u2061(\u2202us\xaf/\u2202z)$, with $us$ being the Stokes velocity. The additional cooling in experiment C+W+L leads to only an marginal increase of 2% in the vertical velocity variance $w\u20322\xaf$ compared to experiment W+L at *z* = −5 m. Note that the increase compared to experiment C+W is more than 300% in both experiment driven by the Stokes forcing. The strong turbulent kinetic energy production cannot be dissipated locally and has to be transported to greater depths. The resulting deeper jets impact the mixed layer evolution below the direct influence of the Stokes drift (Polton and Belcher 2007) and are shown in Figs. 9c and 9d at *z* = −15 m. These jets are also aligned in the direction of the wind stress and the surface waves. The separation scale of ~60 m is greater compared to the classical Langmuir cells in Figs. 9a and 9b. In experiment C+W+L the deep jets are more coherent and more intensified with an increase of 31% in $w\u20322\xaf$ compared to experiment W+L.

The intermittent downwelling deep jets are reminiscent of the downwelling plumes in case of shear convection. The difference is described by Li and Fox-Kemper (2017): though unlike buoyancy force, which acts on a less buoyant water parcel and accelerates the downwelling motion all the way down to the base of the mixed layer, the Stokes shear force is confined in the Stokes shear layer, acting as an initial pump. Li and Fox-Kemper (2017) also show that Langmuir turbulence effects shear turbulence but has less influence in convective turbulence, so that Langmuir turbulence and convective turbulence are approximately additive. However, the deep jets in experiment W+L and C+W+L are less coherent, than the downwelling plumes in experiment C+W, probably because of the high amount of turbulence in the Stokes shear layer.

Nonetheless, the deep jets in experiment W+L and C+W+ lead to an undulation of the mixed layer base and possible internal wave generation through the obstacle mechanism. Similar to the shear convection case the internal wave field is dominated by two distinct frequency bands. Figure 10 shows the spectra of the vertical energy flux for the experiment including Langmuir turbulence at *z* = −50 m below the transition layer averaged over 3 days. In agreement with the dominant zonal orientation of the mixed layer undulations, the energy in the meridional exceeds the zonal propagating waves. The mixed layer undulations are more pronounced in experiment C+W+L which leads to more energetic internal waves compared to experiment W+L. The meridional wavelength peak around 100 m for experiment C+W+L and 50–70 m for experiment W+L (see also Figs. 8c,d).

In cases where Stokes drift and the wind stress are not aligned, the Langmuir rolls are orientated between wind and surface wave direction (Van Roekel et al. 2012). We have not performed any experiments with misaligned forcing, but the obstacle mechanism should be present also here.

The evolution of the vertical energy flux at *z* = −50 m for all experiments is given in Fig. 11. All experiments show maximum values in first couple of hours. One reason is the lack of a transition layer in the beginning of each integration but the initial transition to turbulent motions at the base of the mixed layer after a change in the forcing might also contribute. After the initial response all experiments follow the inertial cycle with maximum values at $t=\pi $ which corresponds to a maximum $\upsilon \xafz=\u2212D$. An exception is the pure shear case experiment W, where the maximum in the inertial cycle, also much lower in amplitude, starts somewhat earlier and corresponds to the maximum shear of zonal momentum at the base of the mixed layer. In experiment C+W, the energy flux is an order of magnitude stronger compared to experiment W. In experiment W, the eddies in the mixed layer are smaller and more confined to the surface as the shear stress follows the classical law of the wall. The energy flux is much larger in experiments with destabilizing buoyancy forcing. Maximum energy fluxes are found in experiment C+W, which also exceeds experiment C+W+L. One reason might be the less coherent roll vortices in the experiment with Stokes forcing, but the undulations in base of the mixed layer are also higher in experiment C+W. The time mean standard deviation of anomalous mixed layer depths is 0.24 m in experiment W, 0.35 m in experiment C+W, 0.31 m in experiment C+W+L, and 0.26 m in experiment W+L. The standard deviation is therefore a good indicator for the magnitude of the internal wave energy flux.

## 4. Summary and discussion

Using an LES model, we study the internal wave generation caused by turbulent motions in the mixed layer driven by different forcing functions. Common in all experiments is the existence of coherent turbulent eddies in the mixed layer, which are advected by horizontal mean motions. The horizontal mean motion consists of a rather steady but vertically sheared Ekman component and vertically homogenous inertial oscillations. All experiments show internal wave activity in the interior of the ocean, that is, below the transition layer, which acts as a low-pass filter for waves generated at the base of the mixed layer. The internal waves in the interior consist of two distinct wave bands. A lower-frequency wave band with angles between the wavenumber $k$ and the horizontal plane of *ϕ* ≈ 72°–90°, that is, including near inertial waves, and higher frequency wave band with angles of *ϕ* ≈ 35°–40°. The vertical energy flux is dominated by the low-frequency wave band.

There is convincing evidence that the low-frequency band is driven by the horizontal mean mixed layer motions through the obstacle mechanism:

The horizontal wavelength of the internal waves corresponds to the eddy/obstacle size in the mixed layer.

The internal waves propagate predominantly normal to the roll vortices in the direction of $\upsilon \xaf|\u2212D$, which we would not expect in case of other driving mechanisms.

The frequency of the waves follow a cycle during one inertial period, according to $\omega =u\xaf|\u2212D\u22c5kh$.

The energy flux in the interior is in very good agreement with the work done by the mean mixed layer motions against an inclined mixed layer base.

The amplitude of the low-frequency waves is linked to the strength of the undulations at the base of the mixed layer. These undulations depend on the source of turbulence with convective eddies being most efficient.

The energy sources of the low-frequency wave band can be separated into an Ekman component which accounts for 58% and inertial oscillations which accounts for 42% of the energy in experiment C+W. The inertial energy in the mixed layer in our experiments corresponds to values in rather calm regions like the subtropics away from intense storm tracks [see Chaigneau et al. (2008) for global maps]. Stronger wind stress forcing or a turning of the wind stress, which amplifies the inertial oscillations, would increase $\u2329u\u232a$ and $u^$. Note that the estimated energy flux by the inertial pumping mechanism is mostly caused by storms (D’Asaro 1985; Alford et al. 2016). Stronger $\u2329u\u232a$ and $u^$ would increase the frequency, vertical group velocity and the vertical energy flux of internal waves generated by the obstacle mechanism. But it is likely that during a storm, the frequency of the wave generation by the obstacle mechanism would exceed *N*, the natural limit of internal waves.

As outlined in the introduction, global estimates of internal wave energy radiated from the surface into the interior is based on the inertial pumping mechanism. The estimates range from 0.34 to 1.4 TW through the surface (Watanabe and Hibiya 2002; Alford 2003, among other) from which ~5%–15% reach the interior (Furuichi et al. 2008; Simmons and Alford 2012; Rimac et al. 2016) resulting in 17–210 GW. To assess the potential role of the obstacle mechanism for interior mixing, we added a dashed line in Fig. 11 corresponding to a globally integrated vertical energy flux of 10 GW for a surface of 3.61 × 10^{8} km^{2} of the world’s ocean. Although our results are based on idealized experiments under constant forcing and for one given initial mixed layer, the comparison shows that at least for destabilizing heat fluxes, the turbulence generated energy flux in the internal wave field should be quantified more precisely.

Polton et al. (2008) run a LES simulation similar to our experiment W+L. They also found high-frequency internal waves closer to *N* driven by Langmuir turbulence. They report that these waves evolve such that their vector phase velocity matches the depth-averaged mixed layer velocity that rotates as an inertial oscillation. The difference to our study is that their analysis is based on the vertical velocity of the first ~12 h of the simulation, that is, the time period where the transition layer is not fully established. Our focus is on the energy flux below the transition layer, so that we do not analyze the first ~12 h in detail, which shows strong differences to the energy flux with established transition layer (Fig. 11).

Undulation of the mixed layer base are also caused by other dynamics like submesoscale eddies with length scales on the order of the mixed layer Rossby radius of deformation, Ro = *O*(1) km. Internal waves generated through the obstacle mechanism would have frequencies of $\omega =u\xaf|\u2212D\u22c5kh=|u\xaf|\u2212D|Ro\u22121$. To reach the limit of $\omega >f$ the velocity has to be larger than *O*(0.1) m s^{−1} in the midlatitudes. The obstacle mechanism might therefore not applicable for undulations with length scales larger than 10 km.

Recent observations of the middepth internal wave spectrum at 16°N in the North Atlantic show increased wave energy in the frequencies above 6 cpd during the winter season (Köhler et al. 2018). The authors show that the internal wave energy in this frequency band correlates with increased surface swell induced by winds north of the trade zone. The present study suggests that the obstacle mechanism is a plausible mechanism contributing to the observed seasonal cycle in the internal wave energy. Swell-generated Langmuir turbulence and, most notably, convective eddies caused by wintertime cooling, are capable of generating the observed increase in internal wave energy.

## Acknowledgments

The model integrations have been performed on the “High Performance Computing system for Earth system research (HLRE-3)” at the Deutsches Klimarechenzentrum (DKRZ), Hamburg, Germany. We thank two anonymous reviewers who helped to improve the manuscript. This paper is a contribution to the project T2 (Energy budget of the surface mixed layer) of the Collaborative Research Centre TRR 181 “Energy Transfer in Atmosphere and Ocean” funded by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) - Projektnummer 274762653.

### APPENDIX

#### Kinetic Energy Equation

The wave-averaged Boussinesq momentum equation is given by (Craik and Leibovich 1976; McWilliams et al. 1997)

Here, $u$ is the wave-filtered Eulerian velocity, $uS$ is the Stokes drift, $g=\u2061(0,0,g)$ is the gravity, $p*=p/\rho o+\u2061(1/2)\u2061(|u\u2212uS|2\u2212u2)$ is the modified pressure, and $\nu sgs$ is the viscosity provided by the subgrid-scale model.

We now form an energy equation for $K=\u2061(1/2)ujuj$, which results in

with $\Gamma K=\u2212uj\u2061(\u2202/\u2202xk)\nu sgs\u2061(\u2202/\u2202xk)uj$ being the dissipation by the subgrid-scale model. Latin indices denote (1, 2, 3) and Greek indices denote horizontal indices (1, 2). We now integrate vertically over the mixed layer. For the boundary conditions we choose a rigid lid surface at $z=0$ with

and a kinematic condition with vanishing Stokes drift at $z=\u2212D$, that is,

where $M=Mflux/\rho $ and $Mflux$ being the mass flux through the lower boundary, that is, the base of the mixed layer.

Vertical integration of Eq. (A2) gives

where $E=MK|\u2212D$ is the flux of *K* through the base of the mixed layer.

## REFERENCES

*Internal Gravity Waves*. Cambridge University Press, 377 pp.

## Footnotes

© 2019 American Meteorological Society. For information regarding reuse of this content and general copyright information, consult the AMS Copyright Policy (www.ametsoc.org/PUBSReuseLicenses).