A wind-driven model of the ocean surface layer with wave radiation physics

Abstract

Surface windstress transfers energy to the surface mixed layer of the ocean, and this energy partly radiates as internal gravity waves with near-inertial frequencies into the stratified ocean below the mixed layer where it is available for mixing. Numerical and analytical models provide estimates of the energy transfer into the mixed layer and the fraction radiated into the interior, but with large uncertainties, which we aim to reduce in the present study. An analytical slab model of the mixed layer used before in several studies is extended by consistent physics of wave radiation into the interior. Rayleigh damping, controlling the physics of the original slab model, is absent in the extended model and the wave-induced pressure gradient is resolved. The extended model predicts the energy transfer rates, both in physical and wavenumber-frequency space, associated with the wind forcing, dissipation in the mixed layer, and wave radiation at the base as function of a few parameters: mixed layer depth, Coriolis frequency and Brunt-Väisälä frequency below the mixed layer, and parameters of the applied windstress spectrum. The results of the model are satisfactorily validated with a realistic numerical model of the North Atlantic Ocean.

1 Introduction

The generation of near-inertial waves by windstress in the upper ocean has extensively been studied using observations and analytical and numerical models. Notable observational results were obtained in the 1970s and 1980s (e.g., Leaman and Sanford 1975; Leaman 1976; Kundu 1976; Fu 1981; Price 1983; D’Asaro 1985) and later by the pioneering work of D’Asaro et al. (1995). These studies demonstrate that velocity fluctuations in the upper ocean below the mixed layer base are dominated by near-inertial frequency oscillations and are qualitatively consistent with internal wave kinematics. Models that have been used to understand the energy budget in the mixed layer and the wave radiation in the stratified ocean below show a broad variety: They range from extensions of the classical Ekman spiral (e.g., Kroll 1975; Kim et al. 2014; Wenegrat and McPhaden 2016) in the attempt to pin down the transfer function between the observed windstress and near-surface currents, via the applications of slab models of varying complexity (Pollard and Millard 1970; Whitt and Thomas 2015; Jing et al. 2017), to diagnostics of high-resolution ocean models (e.g., Simmons and Alford 2012; Rimac et al. 2016). The physics of near-inertial waves in the ocean and the observations, theory, and models have recently been reviewed by Alford et al. (2016).

In the center of the renewed interest in near-inertial waves in the recent decade is the power that near-inertial internal waves (NIW) supply to the ocean interior. Most of this energy supply is expected to be dissipated in the ocean surface layer and used there for mixing and entrainment. The fraction of the radiated flux to the surface flux of energy is thus of particular interest to ocean modelers. The radiative energy flux at the mixed layer base determines how much energy can be converted to turbulence in the interior of the ocean and made available for mixing the stratification. Jochum et al. (2013) and others used a constant value for this fraction to parameterize wave-induced mixing in circulation models. In the wave energy balance and mixing model IDEMIX (Olbers and Eden 2013; Eden and Olbers 2014), a constant fraction is implemented as well. Rimac et al. (2013) and Rimac et al. (2016) give a first estimate of the parametric dependencies (on the local mixed layer depth and windstress amplitude) of the fraction using results from a global numerical ocean model. The present study has the aim to reveal and understand such dependencies with the help of an analytical model.

In extension of the slab model of Pollard and Millard (1970), we develop a simple model of the energetics of wind-driven near-inertial oscillations (NIO) in the mixed layer and the subsequent radiation of near-inertial internal waves (NIW) from the mixed layer base into the ocean. Unlike models of NIO and NIW which are formulated in the time domain as initial value problems, the model presented here evaluates the energy balance of the mixed layer in wavenumber and frequency space, aiming at the spectral flux induced by the windstress at the surface and the flux associated with wave radiation at the mixed layer base. The physics of the model, however, mostly follows other analytical linear models investigating the energetics of NIO in the mixed layer and radiation of NIW from the mixed layer base.

Various types of analytical models have been proposed and most simulate the mixed layer response of particular storm events or a historical windstress climatology. The models by Pollard and Millard (1970) and Gill (1984) have become standards in research on NIO and NIW. Pollard and Millard’s model is limited to the inertial response of the mixed layer currents, and because horizontal homogeneity in the dynamics (no pressure gradients) is assumed, wave radiation into the interior is absent. Gill’s model resolves the pressure forces and, as it is formulated in terms of vertical modes, the slow de-phasing of the modal superposition of an inertial current left in the mixed layer behind a passing storm simulates a vertically propagating signal in the initially quiet and stratified layer below. The radiation is associated with a vertical displacement of the mixed layer base with a near-inertial frequency, a mechanism called by Gill “inertial pumping.” A simple extension of the model by Pollard and Millard (1970) to include radiation physics has been proposed by Olbers et al. (2012) where the inertial pumping is calculated from the divergence and curl of the applied windstress.

Käse and Olbers (1980) base their study on the time-dependent Ekman model to simulate observed wind-driven inertial currents in the subtropical Atlantic. Friction is parameterized by a viscous stress and the model resolves the current profile. The model is then coupled to the wave guide below the mixed layer, anticipating Gill’s inertial pumping. Physically plausible time-dependent Ekman models have in fact been proposed very rarely (notable exceptions are Kroll 1975; Lewis and Belcher 2004; Kim et al. 2014; Wenegrat and McPhaden 2016). The model by Pollard and Millard (1970) with Rayleigh friction as a non-physical parameterization of the flux divergence of momentum has gained more attention (see e.g. Alford 2001, 2003); the likely reason is the lack of knowledge about the profile of viscosity (or stress) which drops in a vertically integrated slab model.

Spectral transfer models have been presented by Käse (1979) and Rubenstein (1994), among others. The windstress is modelled as a body force residing in the mixed layer and the waves are represented by the vertical modes corresponding to the prescribed profile of the Brunt-Väisälä frequency, similar to Gill’s work. In contrast to the initial value problem, however, the wind forcing is assumed to be a random process with a given spectrum in wavenumber-frequency space and the resulting rate of change of the energy spectrum of the forced wave field is calculated (if friction is included as in Rubenstein’s work, the spectrum of energy can be computed as well). Phase averaging is an essential ingredient of the solution technique. Individual wind events can thus not be followed and the evolving stage of vertical propagation from the mixed layer cannot be studied. The time-mean energy transfer from the wind to the wave field is the main outcome of these models.

As these models, our analytical model is linear and works on the f -plane. The wavelengths of waves radiating from the mixed layer base are thus entirely prescribed by the spatial variations of the wind forcing. If these scales are large, the implied vertical propagation speeds may be too small for the waves to leave the upper ocean as observed within a few tens of days (Gill 1984). Thus, if these long scales are imprinted on the mixed layer, there must be a scale transformation. A responsible process is the β-effect (Gill 1984; D’Asaro 1989), which leads to a quick shrinking of the horizontal coherence scale to less than 100 km (Garrett 2001). It is worth mentioning that Kroll (1975) has anticipated much of this work. The β-induced transformation is also the main agent in the NIO equation of Young and Ben Jelloul (1997). The β-effect in combination with balanced shear flow in the mixed layer is shown to be responsible for scale reduction in models like Balmforth et al. (1998) and Balmforth and Young (1999). Lateral propagation and scale reduction are, however, not part of our analytical model.

The new model is developed in steps. We start in Section 2 with the Pollard-Millard model in Fourier space and concentrate on the role of the Rayleigh friction. Wave radiation is then added in Section 3, arising from non-homogeneous forcing, still in a rudimentary approach with inconsistent energetics. The physical ingredients of the new analytical model are developed in Section 4. We solve the linearized equations of motion in the mixed layer in spectral space, including pressure forces and without artificial friction. The energy balance is derived and a complete solution for a given stress profile in the mixed layer is presented. Further analysis is then developed in a long-wave approximation. This enables us to compute the spectrum of the radiation flux in terms of the windstress cross-spectrum (Section 5). The surface energy flux, on the contrary, is shown to depend on interior properties of the mixed layer stress and a less accurate treatment must be applied (Section 4). A comparison with corresponding quantities extracted from a realistic numerical model of the North Atlantic, described in Section 6, is made at various stages of the study (Sections 2, 7, and 8). A summary and discussion follows in Section 9.

2 The role of Rayleigh damping in the Pollard-Millard slab model

We consider a mixed layer (ML) of depth d on top of a stratified ocean. The fluid is driven by a windstress τ₀ applied at z = 0, the ocean surface. The model of Pollard and Millard (1970) (henceforth PM) is given by¹:

$$ \begin{array}{@{}rcl@{}} \partial_t{\textbf{ U}} + f\underset{\neg}{\textbf{U}} = \boldsymbol{\tau}_0 - r \textbf{ U} \end{array} $$

(1)

In this equation, $\textbf { U }=(U,V)={\int \limits }_{-d}^{0} \textbf { u} dz$ is the ML transport, i.e., depth integrated velocity, and τ₀ is the windstress (divided by a reference density). Some assumptions must be made to derive the PM model from Boussinesq equations: (1) a small Rossby number to abandon momentum advection, (2) horizontal homogeneity in the dynamics to ignore the pressure gradient, (3) constant (in time) mixed layer depth (MLD) d to shift the vertical integral under the time derivative, and (4) vanishing stress at the mixed layer base (MLB), τ(z = −d) = 0, so that the ML absorbs the entire turbulent stress in the water column. In the following, we discuss the particular role of the Rayleigh friction term rU in the PM model.

The PM model is readily solved in Fourier space where $\textbf { U}(t)={\int \limits }_{-\infty }^{\infty } d\omega \hat {\textbf {U}}(\omega ) \exp {(-i\omega t)}$, hence:

$$ \begin{array}{@{}rcl@{}} \hat{\mathbf{U}}(\omega) = \frac{(i\omega - r) \boldsymbol{\hat \tau}_0(\omega) + f \underset{\neg}{\boldsymbol{\hat\tau}_0} (\omega) }{\omega^2 - f^2 - r^2 + 2i r\omega} \end{array} $$

(2)

Note that $\hat {\mathbf {U}}^{*}(\omega ) = \hat {\mathbf {U}}(-\omega )$ because of reality of U, correspondingly for the stress. The energetics of the PM model is the balance between the energy flux due to transfer of momentum through the surface by the windstress and the frictional dissipation from the damping term. The surface flux F(ω) is evaluated as:

$$ \begin{array}{@{}rcl@{}} F(\omega) &=& \frac{1}{2} \left( \langle \hat {\textbf{ U}}/d\cdot {\boldsymbol{\hat \tau}_0}^{*} \rangle + c.c \right)\\& =& \frac{r}{d} \frac{(\omega^2+f^2+r^2)\langle\boldsymbol{\hat \tau}_0\cdot \boldsymbol{\hat \tau}_0^{*}\rangle +4 \omega f \Im \langle \hat \tau_{0}^{x} (\hat \tau_{0}^{y})^{*}\rangle}{(\omega^2-f^2+r^2)^2 + 4r^2f^2} \end{array} $$

(3)

where I denotes the imaginary part. The same expression is found for the dissipation:

$$ \begin{array}{@{}rcl@{}} 2 r E(\omega)= \frac{r}{d} \langle \hat {\textbf{ U}} \cdot\hat{\mathbf{U}}^{*} \rangle = F(\omega) \end{array} $$

(4)

closing the balance of energy for PM. The equation can be read with $\langle \boldsymbol {\hat \tau }_{0}\cdot \boldsymbol {\hat \tau }_{0}^{*}\rangle$ as being the periodogram of a particular time series of τ₀(t) or the spectral estimate obtained by ensemble averaging (the other angled expressions likewise). In the latter case, the model becomes stochastic.

In the original work by PM, the Rayleigh coefficient r is used to fit the simulated U(t) for a history of observed τ₀(t) to observations of ML currents at a particular mooring, and values of the decay time 1/r of some days (2 to 20) were obtained. Trivially, the stochastic PM model becomes meaningless if friction is absent: both terms in the PM energy balance vanish identically for r = 0. As can be inferred from Eq. 2, friction in the PM model is important to detune the phase of the horizontal velocity from being in quadrature with the applied stress, i.e., only for non-zero r a non-zero surface energy flux becomes possible. It is often mentioned (e.g., D’Asaro 1985; Alford 2001) that the flux is largely insensitive to the value of r. In view of the above analytical form, this statement is surprising and not supported by the following analysis.

The behavior of Eq. 3 for small r is of particular interest. The expression of the transfer coefficient develops a δ-function at ω² = f² in the limits of r approaching zero. Using

$$ \begin{array}{@{}rcl@{}} \delta(x-t) = \underset{\epsilon \rightarrow 0}{\lim}\frac{\epsilon/\pi}{\epsilon^2 + (x-t)^2} \end{array} $$

(5)

the integrated flux and dissipation of the PM model get a finite value. With 𝜖 = 2rf, we arrive at:

$$ \begin{array}{@{}rcl@{}} \underset{r \rightarrow 0}{\lim}{\int}_{-\infty}^\infty d\omega F(\omega) &\! = & \underset{r \rightarrow 0}{\lim}{\int}_{-\infty}^\infty d\omega 2 r E(\omega) \!= \frac{\pi}{2} \left( T(f) + T(-f) \right) \end{array} $$

with $T(\omega ) = \langle \boldsymbol {\hat \tau }_{0}\cdot \boldsymbol {\hat \tau }_{0}^{*}\rangle (\omega )$. The approach of Eq. 3 toward this δ-function limit with decreasing r is demonstrated in Fig. 1. Here, we integrate the transfer term for an artificial spectrum $\langle \boldsymbol {\hat \tau }_{0}\cdot \boldsymbol {\hat \tau }_{0}^{*}\rangle (\omega ) =\sqrt {a}(f/\pi )/(\omega ^{2} + a f^{2})$ (normalized to unity) which can be done analytically. The parameter a controls the level of the spectrum at low frequencies (see Fig. 1 left panel). As demonstrated, the integrated flux (right panel) for a = 1 attains a value close to the δ-function limit quite fast (already r/f = 0.1 is sufficient). But for smaller values of a, shifting the wind variance to smaller frequencies, the approach can be significantly delayed and much smaller r/f are required. Hence, the figure shows that the dependence of the PM model on the Rayleigh friction coefficient might be substantial if the wind has strong variance at frequencies below the inertial frequency.

Fig. 1

a Normalized spectrum $\langle \boldsymbol {\hat \tau }_{0}\cdot \boldsymbol {\hat \tau }_{0}^{*}\rangle (\omega )=\sqrt {a} (f/\pi )/(\omega ^{2}+a f^{2})$ for four values of a = 1,0.1,0.01,0.001 (black to magenta). b The flux for the normalized spectrum, integrated over ω, as function of r/f. The respective limiting δ-function value is shown as black dashed line

The sore point of the Rayleigh friction in the PM model becomes obvious: the limit of vanishing Rayleigh coefficient $r\rightarrow 0$ leads to a finite net flux; the energy in the ML, however, explodes because it is 1/2r times the flux (see Eq. 3). The model in fact becomes meaningless: with decreasing r, the flux and the dissipation always balance with finite values in a system with ever-increasing energy.

What is a reasonable value for the Rayleigh coefficient? We estimate r using spectral data obtained from a high-resolution ocean circulation model, details of which are given in Section 6 where also the spectral estimation method is described. For the analysis presented here, it is important to note that the model resolves the upper ocean layer with up-to-date turbulence parameterization and adequate vertical resolution. There is, however, no Rayleigh damping in the model but rather harmonic diffusion of momentum. The driving windstress product has a temporal resolution of 1 h. The MLD d was identified as described in Section 6 and spectra of the ML transport U and windstress have been calculated from the model output of 6 winter months of 2 consecutive years in 6^∘× 2^∘ boxes along 30^∘W in the North Atlantic (see Section 6). Summer data yield similar results. Figure 2a shows r/f = F(ω)/(2fE(ω)), estimated from the PM flux F(ω) and energy E(ω), both computed from the numerical model data. There is quite a big scatter of two orders of magnitude for the resulting r/f both in latitude and frequency, however, if a single value has to be picked out $r/f \sim 0.01$ seems reasonable. Figure 2b shows the numerical data of the flux F(ω) and dissipation 2rE(ω) together with the expression of Eq. 3, the theoretical PM flux, the latter two for r/f = 0.02. For this value of r, there is good agreement between the numerical F(ω) and the dissipation 2rE(ω): the curves more or less coincide over most of the resolved frequency range, notably better for super-inertial frequencies. Note that the theoretical PM flux expression follows the general decay in frequency but the level is off up to a order of magnitude; a higher r/f = 0.05 could resolve the mismatch to the numerical F(ω) but increases the dissipation in an unacceptable way. There is a quite a narrow range for r to fit the numerical model results.

Fig. 2

ar/f estimated as F(ω)/(2fE(ω)) from the PM energy balance (3). The data of flux and energy are from the numerical model described in Section 6; the spectra are from 4 latitudes 50^∘N (blue) to 20^∘N (purple). b Numerical data of flux F(ω) (black) and dissipation 2rE(ω) (color as in left panel), and the PM model expression in Eq. 3 for r/f = 0.02 (magenta). In b, the spectra from the lower 3 latitudes are displaced by respective factors 10⁵

3 Inertial pumping in the Pollard-Millard slab model

The role of the Rayleigh term has become increasingly overloaded in recent studies using PM. The PM authors attribute dispersion effects to the r-term. In some studies, it is interpreted as parameterizing downward propagating waves, e.g., in Alford (2003), or even called “radiation stress”², e.g., in Plueddemann and Farrar (2006). The physics of radiation, in fact, is absent from PM. Moreover, as shown below, the radiation of internal waves from the ML does not show up as stress divergence in the momentum balance but rather as a flux in the energy balance, arising from pressure work which is not resolved in PM.

Our aim is to install correct physics of wave radiation into the PM model. The next step is to realize that the fluctuating ocean windstress is non-homogeneous in space, τ₀ = τ₀(x, t). Global applications as Alford’s studies (Alford 2001, 2003) use such windstress. The consequence is that the windstress scale is imprinted on the ML currents which become divergent such that pumping by a vertical velocity W at z = −d (MLB) is induced and internal gravity waves (IWs) are excited below the ML. The process has been called inertial pumping by Gill (1984). Using for simplicity rigid lid conditions at the surface, we obtain $\hat W = \nabla \cdot \hat {\textbf { U}}$ and the pumping velocity is then determined by the divergence and curl of the windstress,

$$ \begin{array}{@{}rcl@{}} \hat W(\textbf{ k}, \omega)= - \frac{(\omega + i r) \textbf{ k} +i f \underset{\neg}{\textbf{k}}}{\omega^2 - f^2 - r^2 + 2i r\omega}\cdot\boldsymbol{\hat\tau}_0 (\textbf{ k}, \omega) \end{array} $$

(6)

derived from Eq. 2. Note that we have switched with the hat notation to Fourier transforms of time and space.

The amplitude and energy of a perturbation with wave vector k and frequency ω, induced by $\hat w_{iw} = \hat W(\textbf { k}, \omega )$ at z = −d, are readily evaluated. The vertical velocity $\hat w_{iw}=\hat w_{iw}(\textbf { k}, \omega , z)$ of the wave field is governed by:

$$ \begin{array}{@{}rcl@{}} \partial_{zz}{\hat w_{iw}} + m^2(z) {\hat w_{iw}} = 0 \end{array} $$

(7)

with the squared vertical wavenumber

$$ \begin{array}{@{}rcl@{}} m^2(z)=k^2 \frac{N^2(z)-\omega^2}{\omega^2 - f^2} \end{array} $$

(8)

where N(z) is the Brunt-Väisälä frequency. We may solve (7) with boundary conditions $\hat w_{iw}= \hat W$ at the MLB z = −d and zero at the bottom by which a modal solution is gained. Alternatively, without the bottom condition but a radiation condition at the MLB (see below), we obtain a radiating wave solution in the internal wave range f² < ω² < N²(z). In each case, the remaining field variables, horizontal velocity, buoyancy, and pressure, follow from the polarization equations (Olbers et al. 2012). For pressure, we find:

$$ \begin{array}{@{}rcl@{}} \hat p_{iw}=i \frac{\omega^2-f^2}{\omega k^2} \partial_z \hat w_{iw} \end{array} $$

(9)

For the radiative type of solution (which we prefer in the following), the wave has the vertical dependence $\sim \exp (i {\int \limits }^{z} m(z') dz')$, propagating in the Brunt-Väisälä frequency profile N(z). We require a positive vertical wavenumber m to have downward group velocity (see e.g. Olbers et al. 2012). This is the radiation condition. The wave field supports a vertical flux of energy given by³$(w_{iw} p^{*}_{iw} + c.c.)/2$ which results in the radiative energy flux:

$$ \begin{array}{@{}rcl@{}} {\Phi}(\textbf{ k}, \omega) &=& \frac{1}{2} \langle{\hat W \hat P^{*}}\rangle + c.c.\\&=& - \frac{ (N^2-\omega^2)^{1/2}(\omega^2-f^2)^{1/2}}{\omega k} \langle \hat W(\textbf{ k}, \omega) \hat W^{*}(\textbf{ k}, \omega) \rangle\\ \end{array} $$

(10)

at the MLB where $\hat P=\hat p_{iw}(-d)$ is the pressure at that depth. The spectrum $\langle \hat W(\textbf { k}, \omega ) \hat W^{*}(\textbf { k}, \omega )\rangle$, appearing in this equation, is that of vertical kinetic energy of the wave field at the MLB. It may be expressed in terms of the stress covariance spectrum $\langle \hat {\boldsymbol \tau }_{0}(\textbf { k}, \omega ) \hat {\boldsymbol \tau }_{0}^{*}(\textbf { k}, \omega )\rangle$ using Eq. 6. We continue with the mathematics of this replacement later in Section 5.

The essential deficit of the PM model becomes obvious: the energy flux Φ(k, ω) must be continuous across the MLB as continuity of $\hat w_{iw}$ and $\hat p_{iw} \sim \partial _{z} \hat w_{iw}$ is required. The PM model, however, does not resolve the pressure, and the energy balance of the ML is still the one between wind-driving and Rayleigh dissipation as in Eq. 3. The energy flux into the wave field is not accounted for and the model thus becomes energetically inconsistent. On the other hand, one can show from Eqs. 6 and 10 that roughly:

$$ \begin{array}{@{}rcl@{}} 2 r E /|{\Phi}| \simeq r /(N kd) \end{array} $$

(11)

for ω² ≪ N² and thus one could argue that the radiation flux is contained in the Rayleigh dissipation of energy. Hence, it could be meaningful to consider 2(r − Nkd)E as dissipation in the ML and 2NkdE as a parameterization of wave radiation. For this to hold r > Nkd, leading to r/f > Nkd/f = 0.3 (take d = 100m,N = 5 × 10^− 3s^− 1, and a wavelength of 100 km) which is a heavy constraint for the Rayleigh coefficient, rendering this redefinition of PM as physically unsatisfying. In the following, we equip the PM model with resolution of the pressure gradient.

4 A model with wave radiation

This section is devoted to the development of a new model with radiation physics. Unlike PM, the new model will resolve the ML momentum balance in the vertical direction and thus becomes mathematically much more involved. In the interior of the ML, a turbulent stress τ(z) is present and—as in PM—assumed to be fully embedded, i.e., τ(z ≤−d) = 0, such that the stress sets an inertial current into motion which is restricted to the ML. As outlined in the preceding section, the ocean part below the ML would remain motionless unless we invoke horizontal variations of the windstress or of the ML depth (the latter, however, are ignored in the present study). As explained in the previous section, a divergence of the inertial current then results and generates a vertical velocity at the MLB which is induced by the divergence of the inertial current much like Ekman pumping is induced by a divergent wind-driven Ekman transport (Gill 1984). If the pumping at the MLB has frequencies above f, it generates internal waves which radiate downward into the stratified ocean part. These waves are called near-inertial waves (NIW) because due to the resonance of the system at ω = f the waves have a frequency very close to f. It is our aim to formulate a model for this process in the most simple form and give a solution for the radiation of the NIWs in spectral space.

4.1 Equations of motion and energetics

The linearized equations of motion for the surface ML − d < z < 0 are:

$$ \begin{array}{@{}rcl@{}} \partial_t {\textbf{ u}} + f \underset{\neg}{\mathbf{u}} + \nabla p = \partial_z{\boldsymbol\tau}\\ \partial_t w + \partial_z p = 0\\ \nabla\cdot \textbf{ u }+ \partial_z w = 0 \end{array} $$

(12)

where u is the horizontal velocity, p is the pressure, and w is the vertical velocity. The turbulent interior stress τ is not parameterized as e.g. in the Ekman model but regarded as given with the dynamic surface condition τ = τ₀ at z = 0. As in the PM model (1), there is a singularity in Eq. 12 at the resonance frequency f. This singularity can of course be avoided by adding e.g. horizontal viscous friction or a simple Rayleigh damping term as in PM. However, to keep the analysis and presentation of our model simple, we proceed without friction terms in Eq. 12 and only introduce them later when needed (in Section 8).

As before for PM, we are looking for a solution of the system (12) in Fourier space where the vertical velocity is represented as:

$$ \begin{array}{@{}rcl@{}} w(\textbf{ x}, t,z)\! =\! \int \mathrm{d}^2 k\! {\int}_{-\infty}^\infty\! \mathrm{d}\omega \hat w(\textbf{ k}, \omega, z) e^{\displaystyle{i \left( \textbf{ k}\cdot \textbf{ x }\! -\! \omega t \right)}} \end{array} $$

(13)

and the other fields accordingly. Note that $\hat w^{\ast }(\textbf { k}, \omega , z) = \hat w(-\textbf { k}, -\omega , z)$.

The aim is to investigate the terms of the energy balance derived from Eq. 12, in particular their relation to the windstress spectrum. The balance of kinetic energy results in the vertical integral over the ML in the form:

$$ \begin{array}{@{}rcl@{}} \partial_t{} \frac{1}{2} {\int}_{-d}^0 \mathrm{d} z \langle\textbf{ u}\cdot \textbf{ u} + w^2 \rangle = - \nabla \cdot{\int}_{-d}^0 \mathrm{d} z \langle \textbf{ u }p \rangle + \langle w p \rangle\Big|_{z=-d} + \langle {\textbf{ u}} \cdot {\boldsymbol \tau} \rangle\Big|_{z=0}- {\int}_{-d}^0 \mathrm{d} z \langle\partial_z \textbf{ u}\cdot \boldsymbol \tau \rangle \end{array} $$

(14)

derived from Eq. 12 with rigid lid condition at the surface and vanishing stress at the MLB.⁴ A statistical steady state is assumed (left-hand side equal to zero) as well as horizontal homogeneity of the statistical quantities. Then, there is no contribution from the divergence of the horizontal energy flux. Taking the Fourier transform, we arrive at the balance of integrated kinetic energy in the ML:

$$ \begin{array}{@{}rcl@{}} \mathcal{P}(\textbf{ k,} \omega) + \mathcal{F}(\textbf{ k}, \omega) - \mathcal{D}(\textbf{ k}, \omega) = 0 \end{array} $$

(15)

where

$$ \begin{array}{@{}rcl@{}} \mathcal{P} = \frac{1}{2} \langle\hat w \hat p^\ast \rangle\Big|_{z=-d} + c.c. ~,~ \mathcal{F} = \frac{1}{2} \langle\hat {\textbf{ u}} \cdot \hat {\boldsymbol \tau}^\ast \rangle\Big|_{z=0} + c.c. ~,~ \mathcal{D} = \frac{1}{2} {\int}_{-d}^0 \mathrm{d} z \langle\partial_z{\hat {\textbf{ u}} } \cdot \hat{\boldsymbol \tau}^\ast \rangle +c.c. \end{array} $$

(16)

The above energy balance can as well be derived directly from the Fourier transform of Eq. 12. The pressure flux ${\mathcal{P}}$ includes the wave radiation flux defined in Eq. 10 but, as shown below, it also has contributions from non-wave pressure components. It pumps energy at all frequencies and wavenumbers from the ML into the stratified part below. Note, however, that only perturbations in the internal wave frequency range can propagate freely. The second term ${\mathcal{F}}$ in Eq. 15 denotes the energy flux due to transfer of horizontal momentum through the surface, induced by the windstress, and the last term ${\mathcal{D}}$ is dissipation in the ML due to the vertical stress. The dissipation is a source of turbulent energy in the mixed layer.

4.2 Solution

From the system of Eq. 12 a single equation for the vertical velocity is derived in the standard way by cross-ward differentiation and elimination of the horizontal divergence and vorticity (see e.g. Olbers and Herterich 1979; Olbers 1986). The governing equation for the Fourier transform $\hat w(\textbf { k}, \omega , z)$ of the vertical velocity becomes:

$$ \begin{array}{@{}rcl@{}} \partial_{zz}{\hat w} - q^2 \hat w = \frac{\boldsymbol{\Omega}\cdot \partial_{zz} {\hat { \boldsymbol\tau}} }{\omega^2-f^2} \end{array} $$

(17)

with q² = k²ω²/(ω² − f²) and $\boldsymbol {\Omega } = \omega \textbf { k }+ if \underset {\neg }{\mathbf {k}}$. From $\hat w$, the complete solution of the equations of motion (12) can be computed:

$$ \begin{array}{@{}rcl@{}} \hat{\textbf{ u}}\! =\! \frac{ i}{\omega k^2}\boldsymbol{\Omega}^\ast \partial_z{\hat w} \! +\! \hat{\boldsymbol{U}}~~\text{and}~~ \hat p \! =\! \frac{ i}{\omega k^2} \left( \! \omega^2\! -\! f^2\! \right) \partial_z{\hat w} \! +\! \hat{\Pi} \end{array} $$

(18)

with

$$ \begin{array}{@{}rcl@{}} \hat{\boldsymbol{U}} = \frac{ i}{\omega k^2} (\raisebox{-1.5ex}{{$\stackrel{\textstyle \textbf{\textbf{ k}}}{\neg}$}} \cdot \partial_z{\hat { \boldsymbol\tau}}) \raisebox{-1.5ex}{{$\stackrel{\textstyle \textbf{\textbf{ k}}}{\neg}$}} \quad\quad\text{and}\quad\quad {\hat {\Pi}} = -\frac{ i}{\omega k^2} \boldsymbol{\Omega} \cdot \partial_z{\hat { \boldsymbol\tau}} \end{array} $$

(19)

This form of solution may be constructed by inserting the ansatz $\hat {\textbf { u}} = \alpha \partial _{z} \hat w + \hat {\boldsymbol {U}}, \hat p =\upbeta \partial _{z} \hat w + \hat {\Pi }$ into the Fourier transform of the equations of motion (12) and determining the coefficients α and β. The respective first term in this ansatz is the divergent part; the second is the non-divergent part. Hence, the directly stress-driven velocity $\hat {\boldsymbol {U}}$, associated with the pressure field $\hat {\Pi }$, is non-divergent and has zero vertical velocity.

In addition to the surface rigid lid condition $\hat w = 0$ at z = 0, the system (12), and correspondingly also Eq. 17, requires a further vertical boundary condition, and we demand:

$$ \begin{array}{@{}rcl@{}} \hat w = \hat W ~~ \text{at}~~ z=-d \end{array} $$

(20)

with the yet unknown pumping velocity $\hat W=\hat W(\textbf { k}, \omega )$. To match the ML field to the stratified interior below the MLB, we assume continuity of $\hat w$ and $\partial \hat w/\partial z$ at the MLB z = −d (see previous section). The unknown $\hat W$ is then determined by these matching conditions which become:

$$ \begin{array}{@{}rcl@{}} \hat w = \hat w_{iw} = \hat W ~~\text{and}~~ \partial_z{\hat w} = \partial_z{\hat w_{iw}} ~~\text{at}~~ z=-d ~~\to~~\partial_z{\hat w}\Big|_{z=-d} = i m \hat W \end{array} $$

(21)

where m is the vertical wavenumber just below the MLB. We can now evaluate the energy fluxes ${\mathcal{P}}$ and ${\mathcal{F}}$, appearing in the energy balance (15) of the ML, in terms of $\hat {w}$. From Eq. 18, we find:

$$ \begin{array}{@{}rcl@{}} \mathcal{P} \!\!&=&\!\! \frac{1}{2}\! \left( \! \frac{ i}{\omega k^2}\! \left( \! \omega^2 - f^2 \!\right) \!\langle \hat w^\ast \partial_z{\hat w}\rangle\Big|_{z=-d} + \langle \hat w^\ast \hat {\Pi}\rangle\!\Big|_{z=-d}\right) + c.c.\\ \end{array} $$

(22)

$$ \begin{array}{@{}rcl@{}} \mathcal{F}\!\! &=&\! \frac{1}{2} \left( \frac{ i}{\omega k^2}\boldsymbol{\Omega}^\ast\!\cdot\! \langle\hat{\boldsymbol\tau}^\ast \partial_z{\hat w}\rangle\Big|_{z=0} +\! \langle\hat{\boldsymbol{U}}\cdot \hat {\boldsymbol \tau}^\ast \rangle\Big|_{z=0} \right) \! + c.c. \end{array} $$

(23)

The dissipation ${\mathcal{D}}$ can be expressed accordingly. Using Eqs. 20 and 21, the first term in Eq. 22 then gives for IW frequencies the wave radiation flux Φ(k, ω) as in Eq. 10. The second term in Eq. 22, involving $\hat {\Pi }$, vanishes if ∂_zτ = 0 at the base z = −d. We assume that this condition holds so that the turbulent stress deposits the momentum primarily in the water column above the MLB. This assumption is confirmed by the already mentioned high-resolution numerical model of the North Atlantic circulation, introduced later in Section 6. Figure 3 shows the variance of the stress profiles from the closure by Gaspar et al. (1990) which we use in the numerical model, decaying rapidly within the ML. Note that only the variance of the stress is important for our analytical model.

Fig. 3

Profiles of standard deviation of stress components in the North Atlantic model described in Section 6. The upper ocean is shown, including windstress at the surface. Left: zonal component τ^x, right: meridional component τ^y are shown for January conditions. Profiles in each subplot are from different locations from 20^∘N (blue) to 50^∘N (yellow) in steps of 10^∘, all along 30^∘W; averages over boxes of approximately 2^∘× 2^∘ in the horizontal. Units are m²s^− 2. Frequency integration is done over all ω. Horizontal lines denote the surface mixed layer depth based on a density criterion: temporal maximum in the considered month horizontally averaged inside the boxes. Stress profiles clearly vanish when the surface mixed layer depth is reached

It remains to compute the vertical velocity field $\hat w$ from Eq. 17 and the aforementioned boundary conditions. An exact solution by a Green’s function approach similar to Olbers and Eden (2016) is described in Appendix A. In the following, however, we restrict the analysis to the long-wave approximation. For (qd)² ≪ 1, Eq. 17 can be solved by simple integration. On the other hand, the long-wave approximation is the limit (kd)² ≪ 1, i.e., the wavelength is large compared to the MLD, which is satisfied for any reasonable windstress forcing. Writing the long-wave condition in terms of q, we find (kd)² = (qd)²(1 − f²/ω²). For ω² ≫ f², the condition (kd)² ≪ 1 is obviously equivalent to (qd)² ≪ 1. For the near-inertial frequencies with ω = (1 + 𝜖/2)f and 𝜖 ≪ 1 the factor becomes 1 − f²/ω² ≈ 𝜖. Therefore, we have to assume the stronger condition (kd)² ≪ 𝜖 ≪ 1 in order to gain (qd)² ≪ 1. Taking a lower limit of wavelengths of 10km, which is satisfied even for realistic atmospheric meso-scale systems, and $d\sim 100 \text {m}$ we obtain (kd)² = 0.004, which demonstrates that even for the near-inertial frequencies with 𝜖 ≈ 0.01 the condition (qd)² ≪ 1 indeed holds.

We now assume (qd)² ≪ 1 and integrate (17) twice. Using the rigid lid boundary condition, the condition (20) at the mixed layer base, and the matching condition (21), we find:

$$ \begin{array}{@{}rcl@{}} \hat w(z) = \hat W \left( 1 + im(z+d) \right) + \frac{\boldsymbol {\Omega}\cdot \boldsymbol {\hat \tau}(z) }{\omega^2 - f^2} \end{array} $$

(24)

such that the pumping velocity becomes

$$ \begin{array}{@{}rcl@{}} \hat W &=& - K \frac{ \boldsymbol{\Omega}\cdot \hat {\boldsymbol \tau}_0}{\omega^2-f^2} ~~\text{with}~~ K = \frac{1}{1+imd} \end{array} $$

(25)

For K ≡ 1, it is identical to the one from the divergent PM model, described in Section 2 by Eq. 6, taking there r = 0. The factor K in Eq. 25 is thus a manifestation of the pressure resolution in our new model. It enters the radiation flux

$$ \begin{array}{@{}rcl@{}} {\Phi}(\textbf{k}, \omega)\! \! &=& \!\! - \frac{ (N^2\! -\!\omega^2)^{1/2}(\omega^2\! -\! f^2)^{1/2}}{\omega k} \langle \hat W(\textbf{k}, \omega) \hat W^{*}(\textbf{k}, \omega) \rangle \\ &=& - \frac{(N^2 - \omega^2)^{1/2}}{{\omega k}(\omega^2 - f^2)^{3/2}} |K|^2 \langle|\boldsymbol{\Omega}\cdot \hat {\boldsymbol \tau}_0|^2\rangle \end{array} $$

(26)

derived from Eq. 22 as factor |K|² which is:

$$ \begin{array}{@{}rcl@{}} |K|^2 = \frac{1}{1+(md)^2} = \frac{\omega^2 - f^2}{\omega^2 -f^2 + (kd)^2 (N^2 - \omega^2)} \end{array} $$

(27)

We note that deviations from K = 1 occur predominantly for NIWs (see Fig. 4). Most important, however, is that the radiation flux is completely determined by the surface forcing of the model, the windstress. Other properties of the stress profile do not occur. This is a consequence of the long-wave approximation: without it, the pumping velocity is an integral of the complete profile of stress in the ML (see Eq. 53 in Appendix A).

Fig. 4

a scaled vertical wavenumber m(k, ω)d (logarithm), b factor |K|²(k, ω), and c ratio of radiation flux to surface energy flux ${\Phi }(k, \omega )/{\mathcal{F}}(k, \omega )$ (logarithm) as function of ω/f and k/k_a. For f = 7.29 × 10^− 5s^− 1, N/f = 100,d = 100 m and 2π/k_a = 100 km. For the ratio α = 0.1 was used

The other components of the energy balance (15), namely the energy flux ${\mathcal{F}}$ through the surface and the dissipation ${\mathcal{D}}$, do not share such simplification. The surface flux depends on the derivative of $\hat w$ at the surface. It follows from Eq. 24,

$$ \begin{array}{@{}rcl@{}} \partial_z{\hat w}\Big|_{z=0} &=& i m \hat W + \frac{ \boldsymbol{\Omega} \cdot \partial_z{\hat {\boldsymbol\tau}}|_{z=0}}{\omega^2-f^2} \end{array} $$

(28)

As a surprise, the first term $i m \hat W$ leads in Eq. 23 (first term) to the same expression as the radiation flux (26) with reversed sign so that:

$$ \begin{array}{@{}rcl@{}} \mathcal{F} = - {\Phi} + \mathcal{F}_{diss} \end{array} $$

(29)

where ${\mathcal{F}}_{diss}$ derives from the second term on the right-hand side of Eq. 28 and the $\hat {\boldsymbol {U}}$-term in Eq. 23,

$$ \begin{array}{@{}rcl@{}} \mathcal{F}_{diss}(\textbf{ k}, \omega) & = & \frac{1}{2} \frac{ i}{\omega k^2} \left( \frac{\langle \hat{\boldsymbol\tau}^\ast \cdot\boldsymbol{\Omega}^\ast \boldsymbol{\Omega}\cdot \partial_z\boldsymbol {\hat \tau} \rangle}{\omega^2-f^2} + \langle \hat {\boldsymbol \tau}^\ast \cdot \underset{\neg}{\mathbf{k}} \underset{\neg}{\mathbf{k}} \cdot \partial_z{\hat { \boldsymbol\tau}}\rangle \right)\Big|_{z=0} + c.c. \end{array} $$

(30)

Because of Eq. 29, the energy balance reduces to ${\mathcal{D}} = {\mathcal{F}}_{diss}$. The part ${\mathcal{F}}_{diss}$ of the surface flux ${\mathcal{F}}$ must thus be dissipated in the ML. We may refer to this flux component as “dissipated surface flux” and to the part −Φ as “transmitted surface flux.”

Both terms in Eq. 30 involve the divergence $\partial _{z} \boldsymbol {\hat \tau }$ of the stress at the surface which remains unknown in the present study. Further investigation of the dissipated surface flux requires an additional parameterization of the stress profile. Note that a simple unidirectional stress τ(z, t) = τ₀(t)g(z) in the water column leads to a vanishing ${\mathcal{F}}_{diss}$ which therefore is an unrealistic setting: the stress vector must turn with depth. Hence, to give an order of magnitude estimate, we insert $\partial _{z}\boldsymbol {\hat \tau }_{0} \simeq - i \boldsymbol {\hat \tau }_{0}/(\alpha d)$ into the first (larger) term in Eq. 30. The parameter α controls the divergence of the stress profile. It must be positive such that ${\mathcal{F}}_{diss}$ becomes positive. We find:

$$ \begin{array}{@{}rcl@{}} |{\Phi}\big/\mathcal{F}_{diss}| \simeq \alpha kd\left( \frac{N^2-\omega^2}{\omega^2-f^2}\right)^{1/2} |K|^2 = \frac{\alpha md}{1+ (md)^2} \end{array} $$

(31)

which becomes α/2 at maximum. Of interest is the ratio of the radiative to the total surface flux:

$$ \begin{array}{@{}rcl@{}} |{\Phi}\big/\mathcal{F}| \simeq\frac{\alpha md}{1 + (md)^2 + \alpha md} \end{array} $$

(32)

resulting from the previous estimate. The maximum is α/(2 + α) which is small since α should be less than 1 (see Fig. 3). We conclude that most of the surface flux is dissipated in the ML, as shown in the example in Fig. 4. In that example, the ratio $|{\Phi }\big /{\mathcal{F}}|$ is found to be 0.03 maximum, i.e., only 3% of the wind-driven energy input at the surface leaves the ML by radiation. A larger α increases the ratio, e.g., for α = 0.3 we find 10%. In this example, the vertical wavenumber m is computed for given k and ω from the dispersion relation (8).

In summary, the long-wave approximation (qd)² ≪ 1 of the ML equations of motion allows expressing the radiation flux Φ(k, ω) at the MLB in terms of the windstress cross-spectrum which we regard in this study as given. The other terms in the energy balance of the ML, namely the surface flux and the dissipation, depend on properties of the interior ML stress; in this approximation, it is the divergence of stress at the surface which remains unknown in the present model. A coarse parameterization of the stress divergence at the surface, however, allows investigation of the ratio of radiative to surface fluxes which is then found to be small of order 10%.

In the following, we evaluate the radiation flux for a specific form of the windstress cross-spectrum, derive frequency spectra of the respective quantities, and compare them to data from a numerical ocean circulation model.

5 Windstress model and derivation of frequency transfer functions

To further evaluate the wave radiation flux Φ(k, ω), derived in the preceding section, we express $\langle |\boldsymbol {\Omega } \cdot \hat {\boldsymbol \tau }_{0}|^{2}\rangle$ in terms of the windstress covariance spectrum:

$$ \begin{array}{@{}rcl@{}} F^\tau_{\alpha\upbeta}(\textbf{ k}, \omega)=\langle \hat\tau_{0\alpha}(\textbf{ k}, \omega) \hat\tau_{0\upbeta}^\ast(\textbf{ k}, \omega) \rangle ~, ~ (\alpha, \upbeta)=(1,2) \end{array} $$

(33)

which leads to

$$ \begin{array}{@{}rcl@{}} \langle |\boldsymbol{\Omega}\cdot \hat {\boldsymbol \tau}_0|^2\rangle &=& \omega^2 k^2 \left( (\cos^2 \phi + \frac{f^2}{\omega^2} \sin^2 \phi ) F^\tau_{11} (\textbf{ k}, \omega) + (\sin^2 \phi + \frac{f^2}{\omega^2} \cos^2 \phi ) F^\tau_{22} (\textbf{ k}, \omega) \right. \\ && \left. + 2 \frac{\omega^2 -f^2}{\omega^2} \sin \phi \cos \phi \Re \left( F^\tau_{12}(\textbf{ k}, \omega) \right) + 2 \frac{f}{\omega} \Im \left( F^\tau_{12}(\textbf{ k}, \omega) \right) \right) \end{array} $$

(34)

where $\textbf { k} =k (\cos \limits \phi , \sin \limits \phi )$ and R and I denote the real and imaginary parts. We need a model for the windstress covariance spectrum. Concerning the wavenumber dependence, we assume that the windstress is a random isotropic vector process (see e.g. Batchelor 1990; Olbers et al. 2012), hence

$$ \begin{array}{@{}rcl@{}} F^\tau_{\alpha\upbeta}(\textbf{ k}, \omega) = A(k,\omega) \frac{k_\alpha k_\upbeta}{k^2} + B(k, \omega) \delta_{\alpha\upbeta} \end{array} $$

(35)

with yet unspecified functions A and B of wavenumber modulus k and frequency ω. Note that A + B carries the spectral properties of the non-divergent part of the stress vector while B carries the rotational part. Transforming the spectral densities from (k, ω) to (k, ϕ, ω), we obtain $F^{\tau }_{\alpha \upbeta }(k, \phi , \omega ) = F^{\tau }_{\alpha \upbeta }(\textbf { k}, \omega ) k$. Integration over horizontal angle ϕ, the isotropic model leads to:

$$ \begin{array}{@{}rcl@{}} F^\tau_{\alpha\upbeta}(k, \omega) = \delta_{\alpha\upbeta}\pi k (A(k, \omega) + 2 B(k, \omega)) \end{array} $$

(36)

Further integration over k and ω yields the covariance of the windstress components τ_0α and τ_0β. Figure 7 below demonstrates that the assumption $F^{\tau }_{11}\simeq F^{\tau }_{22} > |F^{\tau }_{12}|$ of the isotropic windstress model is well satisfied for the observed frequency spectra from the windstress that we utilize later in Section 6.

Implementing the isotropic model into Eq. 26 and integrating over the angle ϕ, we find the radiation flux in the form:

$$ \begin{array}{@{}rcl@{}} {\Phi}(k,\omega) = - |K(k, \omega)|^2 \frac{\omega(N^2-\omega^2)^{1/2}}{(\omega^2-f^2)^{3/2}} 2\pi k^2 \left[ A(k,\omega) + \left( 1 + \frac{f^2}{\omega^2}\right) B(k,\omega) \right] \end{array} $$

(37)

If the pressure force is not resolved, as in the divergent PM model (see Section 2), and |K|² ≡ 1 results, friction is required to achieve a finite total radiation flux. In fact, with |K|² ≡ 1 the flux behaves as ${\Phi }(k, \omega ) \sim (\omega ^{2}-f^{2})^{-3/2}$ and the integral over wavenumber and frequency diverges. The behavior of |K|² at low frequencies, as given by Eq. 27, yields ${\Phi }(k, \omega ) \sim (\omega ^{2}-f^{2})^{-1/2}$ with a convergent integral.

We consider a factorized windstress spectrum for both the non-divergent and the rotational parts, for simplicity in identical form, requiring thus:

$$ \begin{array}{@{}rcl@{}} \pi k(A(k, \omega)\! +\! B(k, \omega))\! \approx\! \upbeta \pi k B(k, \omega) \! =\! F^\tau(\omega) G^\tau(k) \end{array} $$

(38)

where G^τ(k), normalized to unity by integration over k, specifies the wavenumber dependence. We will assume $\upbeta \simeq 1$ in the following. Note that:

$$ \begin{array}{@{}rcl@{}} F^\tau(\omega) \simeq \frac{1}{2} F^\tau_{11}(\omega)\simeq \frac{1}{2} F^\tau_{22}(\omega) \simeq \frac{1}{4} \left( F^\tau_{11}(\omega) + F^\tau_{22}(\omega) \right)\\ \end{array} $$

(39)

For the wavenumber part, we adapt the model of Overland and Wilson (1984) to an isotropic form:

$$ \begin{array}{@{}rcl@{}} G^\tau(k) = \frac{2 k/k_a^2}{((k/k_a)^2+1)^2} \end{array} $$

(40)

The wavenumber k_a is a roll-off wavenumber where the k-dependence of the stress spectrum changes from a $\sim k$ to a $\sim k^{-3}$ power law. Other power laws with a well-defined peak wavenumber yield results similar to those presented below.

Figure 5 shows the radiation flux as function of (k, ω) for a particular set of parameters. As expected, the flux spectrum is heavily peaked at the local inertial frequency at wavenumbers around k_a. Our intention is to integrate the flux spectrum over wavenumbers and discuss the frequency spectrum and compare the results with the numerical model later introduced in Section 6. The integration of the quantities over k can be performed by analytical means and yields the frequency spectrum of the flux:

$$ \begin{array}{@{}rcl@{}} {\Phi}(\omega) &=& -T_{\Phi}(\omega) F^\tau(\omega) \\ T_{\Phi}(\omega) &=& 2 \pi k_a \left( \frac{N^2 - \omega^2}{\omega^2 -f^2}\right)^{1/2} \left( 1 + \frac{f^2}{\omega^2}\right) \frac{\omega}{((N^2 - \omega^2)^{1/2} k_ad + (\omega^2 -f^2)^{1/2} )^2} \end{array} $$

(41)

Note that the frequency transfer function T_Φ(ω) of the flux depends on the parameters N, f, d and the peak wavenumber k_a. Further integration of the flux over frequency yields the total flux. It is given in Appendix B in analytical form and discussed in Section 7.

Fig. 5

Radiation flux ${\Phi }(k, \omega )/F_{0}^{\tau }$ (logarithm) as function of ω/f and k/k_a. For f = 7.29 × 10^− 5s^− 1, N/f = 100,d = 100 m and 2π/k_a = 100 km, $F^{\tau }(\omega )=F_{0}^{\tau } (f/\omega )^{2}$ and the Overland spectrum (40). For the right panel, |K|² was set to 1

In the same manner, we evaluate the dissipated surface flux ${\mathcal{F}}_{diss}$, derived in Eq. 30. With the parameterization of the stress divergence at the ocean surface, discussed there, we find:

$$ \begin{array}{@{}rcl@{}} {\int}_0^{2\pi}d\phi \mathcal{F}_{diss}(k, \phi, \omega) &=& T_{diss}(\omega) F^\tau(\omega) G^\tau(k) \\ T_{diss}(\omega) &=& \frac{4}{\alpha d} \frac{\omega}{\omega^2-f^2} \end{array} $$

(42)

Integration over wavenumbers is readily done and the transfer function T_diss(ω) for ${\mathcal{F}}_{diss}(\omega )$ is obtained. It is independent of the parameters of the windstress spectrum. Note that the total energy flux through the surface, ${\mathcal{F}}= -{\Phi } + {\mathcal{F}}_{diss}$ is characterized by the transfer function $T_{{\mathcal{F}}}=T_{\Phi } + T_{diss}$. However, while T_ΦF^τ is generally integrable and thus yields a finite total radiation flux, the dissipated flux T_dissF^τ is non-integrable. Dissipation needs addition of friction to the momentum balance.

The transfer function T_Φ(ω) is displayed in Fig. 6. Note that fT_Φ/k_a is a dimensionless function of the dimensionless parameters ω/f, N/f and k_ad. We display fT_Φ/k_a as function of ω/f for various values of k_ad and N/f. This function rapidly decays with frequency from a peak at the inertial frequency. For the small k_ad, considered here, and large N/f the decay is $\sim (\omega /f)^{-2}$ away from the peak. Close to the peak and for small N/f, the decay is much stronger. The transfer function T_diss of the dissipated flux is simple and not displayed. It decays as ω^− 1 away from the singular behavior at the inertial frequency.

Fig. 6

Dimensionless response function fT_Φ/k_a as function of ω/f and various values for k_ad and N/f. Left: N/f = 100,k_ad = [0.0006,0.0019,0.0063,0.0188,0.0628], realized for instance for a peak wavelength 2π/k_a = 100km and different depths d = [10,30,100,300,1000]m (blue to green). Right: N/f = [10,30,100],k_ad = 0.0063 (blue, red, yellow). The black dashed line indicates the power law (ω/f)^− 2

6 Spectral data from a numerical ocean circulation model

In this section, we report on wave radiation in a wind-driven numerical model with the aim to compare with the analytical radiation model. The numerical model is designed to resolve low-frequency IWs, roughly up to one decade of frequencies above the respective local f. It is fully non-linear, develops a ML from an up-to-date turbulence parameterization, and is forced by a realistic windstress field the statistics of which are non-homogeneous. The physics of the numerical model goes far beyond the linear, homogeneous setting of the analytical model developed in Section 4 and systematic discrepancies between the two models are therefore to be expected.

6.1 The model

We use an eddy-permitting realistic model of the North Atlantic Ocean using the MITgcm code (Marshall et al. 1997). The configuration of the model is similar to the model used in e.g. Eden and Böning (2002) and identical to the one in Thomsen et al. (2014). It is described in detail in Jurgenowski-Wiegandt (2018). The model domain covers a large part of the north and central Atlantic; it reaches approximately from 18^∘S to 70^∘N. The zonal resolution is 1/12^∘ everywhere while the meridional resolution changes in order to have quadratic grid cells with respect to physical length. The model has 45 vertical levels with thickness ranging from 10m at the surface to 250m at the bottom. To parameterize small-scale turbulence in the ocean surface mixed layer, the scheme of Gaspar et al. (1990) is used, which centers around of a prognostic equation for small-scale turbulent energy (TKE). The closure predicts vertical diffusivity and viscosity in the mixed layer from a mixing length assumption, where the characteristic velocity results from the TKE and the mixing length from a diagnostic relation. In the stratified interior, the mixing is essentially given by small prescribed values. It was shown in Eden and Böning (2002) and later studies that this model configuration allows for a realistic simulation of the large-scale circulation of the North Atlantic with a realistic level of meso-scale eddy energy.

In agreement with our analytical considerations, the model has a rigid-lid surface. In a spin-up phase of 16 years, the model is forced by windstress fields with monthly resolution. The model then settles in a quasi-stationary equilibrium. In a 2-year integration following the spin-up phase, the windstress components are changed to a windstress product with time resolution of 1 h while the time step of the model is 5 min in all experiments. These forcing data are taken from the Climate Forecast System Reanalysis (CFSR, Saha et al. (2010)) product for the years 2003 to 2005 which are used for the following analysis. We have deliberately taken the hourly averaged windstress and not the wind data to avoid a bias in the resulting power supply for the ocean, described in Zhai et al. (2012). In the horizontal, the CFSR windstress is linearly interpolated to our model grid. Model output is saved every 30 min in order to resolve the relevant frequency range between Coriolis frequency f and the adjacent one or two decades of frequency.

In the following, we use quantities near the surface as well as at and below the MLD. Surface quantities, e.g., horizontal surface velocities, are simply taken at the model’s first vertical level. MLD d is determined as the depth at which potential density, when compared with its surface value, has increased by an amount equivalent to a temperature drop of 0.8 ^∘C, following Kara et al. (2000). Since MLD is a Lagrangian quantity, we use the monthly maximum of d, and the annual mean of these values is computed. As detailed below, we evaluate the quantities which we need at MLD d, i.e., vertical velocity w and pressure p, at one grid cell below the diagnosed d.

6.2 The spectra

We discuss results of the numerical model for the years 2003 to 2005 in a spectral view. We preferentially present spectra computed from the 6 winter months of 2003/2004 and 2004/2005 but also of the 6 corresponding summer months of 2004 and 2005. Figure 7 shows frequency spectra of the zonal (${\tau _{0}^{x}}$) and meridional (${\tau _{0}^{y}}$) components of the windstress forcing calculated in boxes around 30^∘W and at two latitudes as indicated in the figure. Every box has a size of 2^∘ in longitude and 2^∘ in latitude and thus includes a varying number of horizontal grid points because of the varying model’s resolution. The time series at each grid point are transformed to spectral space before averaging over the grid points in each box. Every month of the time series data is divided into a sequence of half-overlapping windows with durations of approximately 12 days. This duration has been found best in terms of the minimum resolved frequency and statistical significance. In each of these windows, the discrete Fourier transform is performed separately where a Hanning-type weighting function is applied, and the results are averaged for each month. For the present study, 6-month averages for either winter or summer of the monthly spectra are used. The frequency axis is normalized by the respective local inertial frequency f of the box center. The same analysis has been made for other latitudes along 30^∘ W and the data from two boxes exploited here and in the following were found to be typical. Also included in the windstress spectra plots are curves with power laws ω^− 2 and ω^− 3 for comparison. The windstress spectra mostly lie between these power laws, around ω = f and roughly for half a decade above f the spectral slope is close to − 2, dropping then to − 3 and less at higher frequencies. Summer spectra look similar to the winter spectra but the spectral level is significantly lower by almost an order of magnitude. The overall level of power also decreases from north to south (by a factor 10 from 50^∘ to 20^∘). The vertical lines in Fig. 7 (and the following figures) indicate the expected internal wave frequency range allowed by the model grid. We have approximated the upper bound of this range using the long-wave dispersion relation $\omega ^{2} = f^{2}(1 + k^{2} {R_{1}^{2}})$, with R₁ = NH/(fπ) of the first baroclinic mode where N is local stability frequency at the MLB; furthermore, H = 2000m, and k = 2π/4Δx, where Δx is the grid spacing. Note that this upper bound is a qualitative measure since the parameters are set rather arbitrarily. While at low latitudes, the model grid allows waves also at frequencies much larger than f; at higher latitudes, the frequency range gets much smaller.

Fig. 7

Frequency spectra of windstress components for winter conditions. Blue lines denote the variance spectrum of ${\tau _{0}^{x}}$; purple, the variance spectrum of ${\tau _{0}^{y}}$; and yellow, the absolute value of the covariance spectrum of ${\tau _{0}^{x}}$ and ${\tau _{0}^{y}}$. Data from horizontal boxes of approximately 2^∘× 2^∘ extent are considered whose centers lie along 30^∘W. Left: at 30 ^∘N, right: at 50^∘N. Monthly determined spectra are averaged for winter value. Horizontal averages over the boxes are performed after spectral calculations. The black vertical line denotes f; the red one the frequency up to which internal gravity waves are expected to be resolved (see text) which depends on N at d. Data are taken from years 2004 to 2005. Both axes are logarithmic; the frequency axis is normalized by the respective local Coriolis frequency f. Power laws ω^− 2 and ω^− 3 are indicated by black dashed lines

Figure 8 shows frequency spectra of the different contributions to kinetic energy: horizontal kinetic energy $E_{ek}^{h}(\omega ) = \langle \hat u \hat { u}^{\ast } + \hat v\hat { v}^{\ast } \rangle / 2$ slightly below MLD and its vertical analogue $E_{ek}^{v}(\omega ) = \langle \hat w \hat { w}^{\ast } \rangle / 2$, as well as total kinetic energy at the ocean surface. Here, $\hat u$, $\hat v$, and $\hat w$ denote the complex amplitudes of the Fourier transforms of horizontal velocity u, v, and vertical velocity w, as before, and the angular brackets denote the statistical expectation. $E_{ek}^{v}$ is many orders of magnitude lower than $E_{ek}^{h}$, appropriate to the long gravity waves resolved by the model grid. Both $E_{ek}^{h}$ and $E_{ek}^{v}$ at the mixed layer base feature peaks at f and around 2f and there is also some indication of a peak around 3f. This is not seen in the total surface kinetic energy where a peak appears only at f. The spectral peaks near the inertial frequency increase and widen toward the south. The reason is that the inertial peak at a certain latitude comes about by direct local forcing and due to waves propagating to the specific point from farther north, having a larger frequency imprinted there as near-inertial. Remember that NIW cannot move much northward because they encounter turning points slightly north of their generation region. Summer spectral levels are significantly lower than the winter levels shown here. This is particularly valid for the harmonic peaks.

Fig. 8

Same as Fig. 7 but for total kinetic energy at the surface (purple), vertical (red), and horizontal (blue) kinetic energy taken just below MLD (see text for details). Power laws ω^− 4 and ω^− 6 are indicated by black dashed lines. The black full vertical line denotes f; the black dotted one 2f. The yellow vertical line indicates the expected internal wave frequency range allowed by the model grid as detailed in the text

A closer inspection of the peaks near 2f and 3f reveals that they occur precisely at these frequencies at high latitudes (see 50^∘N in Fig. 8), whereas they are shifted to higher frequencies than 2f and 3f at lower latitudes (see 30^∘N in Fig. 8). The shift is generally in good agreement with propagation from higher latitudes as well, e.g., at 30^∘ the “2f” peak is at a frequency roughly equal to 2f at 50^∘. We conclude that the NIW spectral part of the numerical model data is quite severely modified by propagation effects, a process that is completely absent in the analytical model.

Due to the presence of peaks in the energy spectra it appears difficult to point down a dominating power law behavior. However, a rough agreement with (ω/f)^− 6 is reasonable for the spectra below MLD, in particular for the vertical kinetic energy. The surface kinetic energy rather follows (ω/f)^− 4. Such behavior is very different from what is found in the deep ocean where frequency spectra show slopes close to − 2 (see e.g. Fu 1981). Such spectral slopes close to − 2 also apply to upper ocean data (Roth et al. 1981; Levine et al. 1983). Note, however, that these reported slopes refer to the behavior at frequencies much higher than f, which are not resolved by the numerical model, while slopes in the near-inertial range in the reported observations are considerably steeper.

Figure 9 shows energy flux terms important for the energy balance which has been derived in Section 4, namely the cross-spectra $F_{u\tau }(\omega ) = (\langle \hat {\textbf { u}}_{0} \cdot \hat {\boldsymbol \tau }_{0}^{\ast } \rangle + c.c.)/2$ (evaluated at z = 0) as well as $F_{wp}(\omega ) = (\langle \hat {w} \hat {p}^{\ast } \rangle +c.c.)/2$ (evaluated slightly below z = −d) where $\hat {\textbf { u}}, \hat w, \hat p$, and $\hat {\boldsymbol \tau }_{0}$ are the Fourier transforms of horizontal and vertical velocity, pressure, and windstress vector. F_uτ and F_wp are the fluxes of total kinetic energy through the ocean surface and through the mixed layer base, respectively. The power law behavior of the fluxes is similar to the energies: we find roughly (ω/f)^− 6 for the radiation flux and (ω/f)^− 4 for the surface flux of energy. The input of energy from the wind into the mixed layer (F_uτ) is directed into the ocean at all frequencies and exceeds the pumping flux (F_wp) out of the mixed layer through its base at all frequencies. Evidently, in frequency space, there are regions where the energy flux F_wp is directed not down- but upwards (down in Fig. 9 is red, up is blue). Note that in the distinct peak around inertial frequency f, energy is generally transported downwards as expected from the wind-driving scenario, whereas in some peaks, e.g., the one around 2f in winter at 50^∘ N, the direction of the energy flux is upward toward the MLB. Over large parts of the frequency axis away from the peaks, energy fluxes are varying between up- and downward propagation. The source of this upward-propagating energy cannot be directly forced by the wind and we conclude that this is a sign for non-linear interactions of the directly wind-forced waves by which energy is redistributed from the near-inertial frequency to its higher harmonics. The conversion is likely occurring deeper in the water column among the incident and reflected wave components to explain the upward direction of the energy flux at the harmonic frequencies in the upper ocean. As explained above, the energy around the higher harmonics may partly originate from higher latitudes. There are many reports of the occurrence of 2f and higher harmonics in numerical simulations of near-inertial forcing by hurricanes (e.g., Price 1983; Niwa and Hibiya 1997; Meroni et al. 2016) with clear indications of wave–wave interactions. Other generation mechanisms (Danioux and Klein 2008; Danioux et al. 2008; Zedler 2009) involve the interaction of a near-inertial oscillation with the vorticity of a background eddying flow.

Fig. 9

Same as Fig. 7 but for the fluxes R(F_uτ) (yellow) and −R(F_wp), the latter taken just below MLD (red dots for downward flux and blue dots for upward flux). Upper row for winter, lower row for summer. Black dashed lines indicate the power laws (ω/f)^− 4 and (ω/f)^− 6. The MLB flux is an order of magnitude smaller than the surface flux. Peaks at f of the surface flux and the MLB flux are visible, the latter one has downward direction as expected

7 Evaluating the wave radiation flux

In this section, we evaluate the wave radiation flux Φ of the analytical model, comparing with the diagnosed frequency spectra of the numerical model discussed in the preceding section. It was shown in Section 5 that the frequency spectrum of Φ can be expressed as:

$$ \begin{array}{@{}rcl@{}} {\Phi}(\omega) = -T_{\Phi}(\omega) F^\tau(\omega) \end{array} $$

(43)

where T_Φ denotes the transfer function given by Eq. 41 while F^τ describes the spectrum of the specific windstress forcing for the radiation flux, given by Eq. 39. Remember that the transfer function in the relation (43) relies on the specific wavenumber spectrum of the windstress model and introduction of an unknown wavenumber scale k_a for the spectral peak.

We use the spectra diagnosed from the numerical model described in Section 6 to validate the frequency response functions T_Φ(ω) for the flux Φ(ω). The cross-spectrum F_wp(ω) is determined for each position and level in the numerical model as described in Section 6, and a model-related response function is computed by forming T_wp(ω) = −F_wp(ω)/F^τ(ω), using the wind product that drives the numerical model to estimate F^τ(ω) (see Fig. 7). Figure 10 compares T_wp(ω) with T_Φ(ω), evaluated for the respective MLD d and the local N diagnosed from the numerical model as described in Section 6. For T_Φ(ω) from Eq. 41, we also need to specify the peak wavenumber k_a of the windstress spectrum.

Fig. 10

Frequency response functions $T_{wp}\sqrt {\omega ^{2} - f^{2}}$ and $T_{\Phi }\sqrt {\omega ^{2} - f^{2}}$ for the radiation flux (see text for details) diagnosed from the numerical model (red dots: downward flux, blue dots: upward flux), and evaluated for the analytical model (black lines) as function of ω/f, evaluated at two latitudes: left: at 30^∘N, right: at 50^∘N. Upper row for winter, lower row for summer. The square root factor eliminates the peak at f from the analytical transfer function. The analytical model uses the winter, respectively summer, mean MLD d and N from the numerical model and peak wavelengths 2π/k_a = 50,100, and 300km (the overall level of the analytical transfer functions decreases with increasing wavelength). The vertical purple line denotes the upper bound until which internal gravity wave motions are expected to be permitted by the numerical model (see text)

There are reports on length scales for windstress fields from several studies. Levine and Zervakis (1995) report wavelengths of size 550 to 900km in observed traveling fronts. Simmons and Alford (2012) and Rimac et al. (2016) find wavelengths around 200 to 500km in their forcing fields. A completely different range of wavelengths is considered in the work of Rubenstein (1994). He uses meso-scale wind fields and has the aim to reproduce the Garrett-Munk spectrum (Garrett and Munk 1975; Munk 1981) of the internal wave field. The windstress spectra are synthesized from data and parameterized using the model spectrum by Overland and Wilson (1984) which we use here as well. They have a peak at wavelength ranging from 10 to 100km. As pointed out by D’Asaro (1989), the larger wavelengths around 500km are associated with large-scale atmospheric weather systems, while the smaller ones occur for small-scale hurricanes and moving fronts.

The response function T_Φ(ω) is shown in Fig. 10 for three values of k_a corresponding to wave lengths of 50, 100, and 300 km. Note that T_wp is shown for the depth level immediately below the diagnosed mean MLD d. For the Brunt-Väisälä frequency N, we also use the diagnosed value at the respective depth. Concerning the overall spectral level and slope the analytical model fits reasonably well to the numerical solution. Note that the analytical response function and its numerical counterpart decays by some decades in the internal wave range which is captured by the theoretical model. The correspondence is best at higher latitudes and less so at low latitudes, and generally better for the summer conditions with the lower forcing and energy levels compared to the winter conditions. The likely reason is that the analytical model takes into account only the locally forced waves while in the numerical model also remotely forced waves show up and this increasingly at lower latitudes due to predominant southward propagation (see also Section 6). Moreover, the values of the numerical model in the range of upward fluxes (blue dots) cannot agree with the analytical model because the nonlinear processes responsible for the peaks at 2f and higher harmonics are not captured by the linear model. Note that the values predicted by the analytical model mostly stay below these additional peaks, marking there the level of the wind-driven linear response in the background of the upward propagating super-harmonics. The results appear relatively insensitive to the peak wavenumber k_a.

In Appendix B, we derive the integrated (over frequency and wavenumber) radiation flux predicted by the analytical model. Here, we discuss the parameter dependencies and compare again with the numerical model. We have to specify the windstress frequency spectrum which we prescribe simply as a power law, assuming $F^{\tau }(\omega )= F_{0}^{\tau } (f/\omega )^{2}$, as detailed in Appendix B. The value $F_{0}^{\tau }$ of the spectral level of the windstress frequency spectrum at near-inertial frequencies will enter the wave flux as multiplicative factor. A typical value for mid-latitude oceans is $F_{0}^{\tau } = 10^{-6} \text {m}^{4} s^{-4}/ \text {s}^{-1}$ (used in the following figures) but higher values are still acceptable for higher latitudes (see Fig. 7 of this study and e.g. Figure 5 of Willebrand (1978)). For the forcing used by the numerical model of Section 6 we find $F_{0}^{\tau }= 1\times 10^{-7} \text {m}^{4} s^{-4}/ \text {s}^{-1}$ at 20^∘ N in summer and $F_{0}^{\tau }= 3\times 10^{-6} \text {m}^{4} s^{-4}/ \text {s}^{-1}$ at 50^∘ N in winter.

The total pumping flux $R=|{{\int \limits }_{f}^{N}} {\Phi } d \omega |$ is a complicated function of k_adN/f, given by Eq. 55 in Appendix B. Figure 11 shows the dependency of the total pumping flux R on the MLD d, the Brunt-Väisälä frequency N, and wavelength 2π/k_a. The net flux R decays with increasing d. A rough approximation is $R \simeq 4\pi F^{\tau }_{0}/d$ (see Appendix B) and a typical value for the flux is a few times 10^− 8m³s^− 3 ≡ 0.01mW m^− 2 at moderate mixed layer depths.

Fig. 11

Radiation flux R (logarithm) in m³s^− 3 as function of the MLD d and Brunt-Väisälä frequency N for two values of the windstress peak wavelength 2π/k_a = [10,100]km (left and right panel, respectively) at 30^∘ latitude. For further parameters, see text

8 Evaluating the surface flux

In Section 5, we have derived the theoretical prediction for the frequency spectrum of the surface energy flux in the form:

$$ \begin{array}{@{}rcl@{}} \mathcal{F}(\omega) = T_\mathcal{F}(\omega) F^\tau(\omega) \end{array} $$

(44)

where $T_{{\mathcal{F}}}=T_{\Phi } + T_{diss}$. The transfer function T_Φ and T_diss are given by Eqs. 41 and 42, respectively. In Fig. 12, we compare this theoretical model with data from the numerical model. The cross-spectrum $F_{u\tau }(\omega ) = (\langle \textbf { u}_{0} \cdot \boldsymbol \tau _{0}^{*}\rangle +c.c.)/2$ is computed as described in Section 6.2. The figure shows F_uτ(ω)/F^τ(ω) together with $T_{{\mathcal{F}}}(\omega )$. The only parameter in this expression is α. The comparison is convincing with α independent of frequency. The overall dependence of the transfer function on frequency is well captured though the proper value of α differs substantially with respect to season and latitude. For the winter profiles, the α parameter is 0.02 to 0.05, and for summer 0.5 to > 2.5, larger values corresponding to lower latitudes. Note that the depth scale αd does not have such a large range: the MLD at the analyzed positions is 212 m at 30^∘ and 334 m at 50^∘, hence $\alpha d \sim 20 {\cdots } 35 \text {m}$. The summer values are 13 m at 30^∘ and 51 m at 50^∘, hence $\alpha d \sim 25 {\cdots } 30 \text {m}$.

Fig. 12

Frequency response functions of the surface flux. Comparison of F_uτ(ω)/F^τ(ω) (black dots) with $T_{{\mathcal{F}}}(\omega )$ (full) where the dissipation scale α is set to 0.02 (blue), 0.1 (purple, 0.5 (red), and 2.5 (yellow). Upper row winter, lower row summer at latitudes 30^∘ N and 50^∘ N

Integration of ${\mathcal{F}}(\omega )$ over frequencies meets the problem that the part from dissipated flux is singular without addition of friction. The problem is mediated when Rayleigh friction is added to the horizontal momentum balance in Eq. 12. We find (57) in Appendix C to replace T_diss from Eq. 42 and the dissipated flux integrates approximately to $F_{diss} \simeq 2 F_{0}^{\tau }(\alpha d) \ln {(f/2r)}$ (the exact expression is in Appendix C). It depends logarithmically on the friction coefficient r and is typically less than 10^− 6m³s^− 3 ≡ 1mW m^− 2 at moderate mixed layer depths and with r = 0.01. Dissipation of energy in the ML thus exceeds radiation at the MLB by an order of magnitude (see further discussion in the summary Section 9).

Figure 13 finally displays the ratio of radiative to surface flux, R/(R + F_diss) for the theoretical model. Because of the linear nature of theoretical model, the ratio is independent of the overall level of the windstress spectrum; however, there is a dependence on the dominant wavenumber k_a as we use here the full expression for R as given in Appendix B. Note that the ratio is typically 0.1, decreasing with d and, to a less degree, also with N.

Fig. 13

Ratio R/(R + F_diss) as function of the MLD d and Brunt-Väisälä frequency N for two values of the windstress peak wavelength 2π/k_a = [10,100]km (left and right panel, respectively) at 30^∘ latitude. Further parameters: α = 0.1,r/f = 0.01

9 Summary and discussion

Internal gravity waves are subject to plentiful forcing processes (see e.g. Munk 1981; Olbers 1983; Polzin and Lvov 2011). A prominent mechanism in the upper ocean is the generation by inertial pumping, induced by windstress acting in the mixed layer (Gill 1984). It is the focus of the present study. We have developed a new analytical model of wind-induced wave radiation from the surface mixed layer. Central in the study is the energy balance of the mixed layer, including the wind-induced surface energy flux and the resulting radiation flux carried by internal waves below the mixed layer base. We discuss the parametrical dependence of the energetics on the mixed layer depth and the local Brunt-Väisälä frequency, as well as the spatial and temporal scales of the windstress and derive analytical expressions for the fluxes of energy induced by the windstress at the surface and by wave radiation at the mixed layer base.

A short critical assessment of friction in the Pollard-Millard slab model (Pollard and Millard 1970) starts our analysis of the wind-driven surface layer. We regard our model as extension of the slab model, abandoning, however, the Achilles heel of the nonphysical Rayleigh friction. The new model resolves the pressure gradient in the momentum balance of the mixed layer. The main assumptions of the new model are:

• The model is linear, local on a f -plane, and the mixed layer depth is considered temporally constant.

• Statistical homogeneity in the horizontal is assumed for the energy balance of the mixed layer; wave propagation in a β-effect environment is absent.

• The stress profile in the mixed layer is regarded as given (not parameterized).

• Vertical pumping is induced by spatial variations of the applied windstress and the forcing in the frequency range of the local f to the local N goes into the internal wave field.

• A long-wave approximation is made, essentially kd ≪ 1, where k is the wavenumber and d the mixed layer depth.

• The windstress covariance spectrum is modeled by an isotropic vector process and the associated wavenumber-frequency spectrum is taken factorized.

By this construction, the model captures only the locally forced wave field. On this base, the radiation flux spectrum is derived as function of wavenumber and frequency and the system parameters (mixed layer depth, Brunt-Väisälä frequency below the mixed layer, scales of the windstress spectrum) and integrated over the Fourier space resulting in the flux in physical space. While the radiation flux is expressed by the spectrum of the windstress fluctuations as a consequence of the long-wave approximation, the surface flux of energy, directly induced by the wind, is found to depend on interior stress properties, namely the stress divergence at the surface. It could therefore only be studied with further approximation to come up with a closed view of the energetics of the mixed layer. A comparison of the radiation flux and the surface flux in the spectral frequency domain, predicted from the analytical model, with a numerical circulation model of the North Atlantic with state-of-the-art mixed-layer physics is found by and large successful, given the simplicity and restriction of the analytical model—linear physics and horizontal homogeneity, both of which are not realized in the numerical model.

A central property of the wind-forced motions in the upper ocean is the amount of energy that escapes dissipation in the mixed layer and is able to radiate from the mixed layer base as internal gravity waves into the stratified ocean. Figure 14 displays maps of the ratio of the wind-induced radiative flux and the surface flux of energy for winter and summer, diagnosed from the numerical model of the North Atlantic, described in Section 6. Note that the color scale is linear and includes positive and negative values, the latter coming about from a negative net radiation flux at the respective position. Such features could not occur in the theoretical model with the energy balance in Section 4 appropriate for horizontal homogeneity. In the numerical model, however, we have horizontal variations and in particular, the horizontal wave-induced energy flux must be present in the energy balance of the mixed layer with a non-vanishing divergence. Hence, we face in some regions an access in the up- or downward radiation of energy compared to the local supply by the surface flux. Figure 14 shows that this mainly occurs in the tropics and in the Gulf Stream and its extension but also in smaller spots in the open ocean. The maps shows positive ratios around 0.1 for summer and winter in the open ocean regions where the spectral analysis is made, in good agreement with the theoretical model, but large values exceeding 1 occur in summer in more southerly regions. For winter, we are facing a mixed behavior with large positive and negative ratios south of 40^∘ N where near-inertial waves, generated more northward, may accumulate and induce divergences of horizontal energy flux.

Fig. 14

Maps of the ratio of radiation flux to surface flux from the numerical model described in Section 6. Only power at super-inertial frequencies has been included to compute the net fluxes. Left: mean over the 6 winter months of 2003 to 2005. Right: mean over 6 summer months of 2004 and 2005. The scale is linear

The flux of energy from the wind to the surface mixed layer has been computed by many models of varying complexity (e.g., Pollard and Millard 1970; Alford 2001; Plueddemann and Farrar 2006; Scott and Xu 2009). The flux of energy at the mixed layer base and the ratio of radiation flux to surface flux, though eminently more important for the ocean circulation, is much less known. Previous estimates come from numerical experiments and mostly refer to global-scale means (Furuichi et al. 2008; Zhai et al. 2009). Rimac et al. (2016), using a global model with 1/10^∘ resolution, reports 11.4% for the ratio on average with a decreasing tendency with respect to increasing mixed layer depth: for shallow values the ratio is about 0.2 and for deep layers the ratio decreases to about 0.04. In our North Atlantic model, we find that in total 83% of the wind-generated flux (in super-inertial frequencies and winter) is dissipated in the mixed layer and 17% is radiated as waves into the interior ocean.

Appendices

Appendix A: Green’s function of the $\hat w$-problem

The solution of Eq. 17 with the rigid lid and the inhomogeneous boundary condition (20) is written in term of the Green’s function for the problem with homogeneous boundary conditions, defining the modified vertical velocity

$$ \begin{array}{@{}rcl@{}} \tilde w = \hat w + \hat W z/d \quad\text{with}\quad \tilde w = 0 \quad\quad\text{at}\quad\quad z=0, -d \end{array} $$

(45)

The solution is $\tilde w = \tilde w_{1} +\tilde w_{2}$ where:

$$ \begin{array}{@{}rcl@{}} \tilde w_1 &=& {\int}_{-d}^0 d\xi G(z,\xi) \frac{\partial \mathcal{W}}{\partial \xi} \quad\quad\text{with}\quad\quad \mathcal{W} = \frac{\boldsymbol{\Omega}\cdot \partial_z{\hat { \boldsymbol\tau}} }{\omega^2-f^2} \\ \tilde w_2 &=& - \frac{q^2 \hat W}{d}{\int}_{-d}^0 d\xi G(z,\xi) \xi = \hat W \left( \frac{z}{d} -\frac{\sinh qz}{\sinh qd}\right) \end{array} $$

(46)

with the Green’s function:

$$ \begin{array}{@{}rcl@{}} G(z,\xi) &=& \frac{1}{q\sinh qd} \left\{ \begin{array}{cc} \sinh q (z+d) \sinh q\xi \quad\quad\text{for}\quad\quad -d\le z\le \xi \le 0 \\ \sinh q z \sinh q (\xi+d) \quad\quad\text{for}\quad\quad -d \le \xi \le z\le 0 \end{array} \right. \end{array} $$

(47)

Hence

$$ \begin{array}{@{}rcl@{}} \hat w= - \hat W \frac{\sinh qz}{\sinh qd} + \tilde w_1 \end{array} $$

(48)

The condition (21) then results in

$$ \begin{array}{@{}rcl@{}} \hat W &=& -\frac{1}{\cosh q d + i m/ q \sinh qd} {\int}_{-d}^0 d\xi (\mathcal{W} - \mathcal{W}_d) \cosh q \xi \end{array} $$

(49)

Here ${\mathcal{W}}_{d}={\mathcal{W}}(z=-d)$. The derivative of $\hat w = \tilde w - \hat W z/d$ at the surface z = 0, needed to evaluate the horizontal surface velocity from Eq. 18, becomes

$$ \begin{array}{@{}rcl@{}} \partial_z {\hat w}\Big|_{z=0} &=& - \frac{q}{\sinh qd} {\int}_{-d}^0 d\xi (\mathcal{W} - \mathcal{W}_d) \cosh q (\xi + d)\\ &&+ \mathcal{W}_0-\mathcal{W}_d - \frac{ q \hat W}{\sinh qd} \end{array} $$

(50)

with ${\mathcal{W}}_{0} = {\mathcal{W}}(z=0)$. In both integrals in Eqs. 49 and 50 we have performed a partial integration. We arrive at

$$ \begin{array}{@{}rcl@{}} \hat W &=& - K\mathcal{C} \end{array} $$

(51)

$$ \begin{array}{@{}rcl@{}} \partial_z{\hat w}|_{z=0} &=& \frac{q}{\sinh qd}\left( K \mathcal{C} - \tilde{\mathcal{C}}\right) + \mathcal{W}_0 - \mathcal{W}_d \end{array} $$

(52)

with

$$ \begin{array}{@{}rcl@{}} \mathcal{C} &=& {\int}_{-d}^0 d\xi \left( \mathcal{W} -\mathcal{W}_d \right)\cosh q\xi \quad\quad\text{and}\\\tilde{\mathcal{C}} &=& {\int}_{-d}^0 d\xi \left( \mathcal{W} -\mathcal{W}_d \right)\cosh q(\xi+d)\\ K &=& \frac{1}{\cosh qd + i(m/q) \sinh qd} \end{array} $$

(53)

The solutions presented in Section 4 follow as lowest order expansion in terms of qd.

Appendix B: Total radiation flux

The total flux is obtained by integration over frequency in the internal wave range f ≤ ω ≤ N. We specify the windstress frequency spectrum simply as a power law, assuming

$$ \begin{array}{@{}rcl@{}} F^\tau(\omega) = F_0^\tau (f/\omega)^{2} \end{array} $$

(54)

oriented at the findings of the forcing of the numerical model (see Fig. 7). With a minor approximation the total radiation flux can be expressed in analytical form. Since the power of the windstress fluctuations is at low frequencies we may replace N² − ω² in Eq. 41 by N². Then

$$ \begin{array}{@{}rcl@{}} R &=& \Big|{\int}_f^N \mathrm{d}\omega {\Phi}(\omega) \Big| \\ &=& \frac{2 \pi k_a F^\tau_0 N/f }{(1 + \eta^2)^2} \left( \frac{2}{\eta} - \frac{3\pi}{4}\left( 1-\eta^2\right) - 2\eta \ln \eta - \frac{(\pi/2) (1 - 3 \eta^2) + 4\eta\ln \eta}{1 + \eta^2} \right) \end{array} $$

(55)

with η = k_adN/f. We have assumed here for simplicity that N/f ≫ 1,k_ad ≪ 1. It is found that Eq. 55 is valid for a large range of parameter values. It is well approximated by the first term in the brackets, hence:

$$ \begin{array}{@{}rcl@{}} R \simeq \frac{4 \pi k_a F^\tau_0 N/f }{\eta (1 + \eta^2)^2} = \frac{4 \pi F^\tau_0 }{d (1 + \eta^2)^2} \simeq \frac{4 \pi F^\tau_0 }{d} \end{array} $$

(56)

Appendix C: Total dissipated flux

The transfer function for the dissipated flux must be expanded to include friction. Then:

$$ \begin{array}{@{}rcl@{}} T_{diss} = \frac{4}{\alpha d} \frac{\omega(\omega^2-f^2)}{(\omega^2 - f^2)^2 + 4r^2f^2} \end{array} $$

(57)

where r is a Rayleigh friction coefficient. One finds for the total dissipated flux

$$ \begin{array}{@{}rcl@{}} F_{diss} &=& {\int}_f^N T_{diss}(\omega) F^\tau(\omega) d \omega \\ &=& \frac{F_0^\tau}{\alpha d} \frac{f^2}{4 r^2 + f^2} \left( \ln \frac{(N^2 -f^2 )^2 +4 r^2 f^2}{N^4} + 4 \frac{r}{f} \arctan \frac{N^2 -f^2}{2 r f} + \ln \frac{f^2}{4 r^2} \right) \end{array} $$

(58)

The last term in the brackets is by far the largest contribution if N² ≫ f² ≫ r², hence

$$ \begin{array}{@{}rcl@{}} F_{diss} \simeq \frac{2 F_0^\tau}{\alpha d} \ln \frac{f}{2 r} \end{array} $$

(59)