## Abstract

The scaling of turbulent kinetic energy (TKE) and its vertical component (VKE) in the upper ocean boundary layer, forced by realistic wind stress and surface waves including the effects of Langmuir circulations, is investigated using large-eddy simulations (LESs). The interaction of waves and turbulence is modeled by the Craik–Leibovich vortex force. Horizontally uniform surface stress *τ*_{0} and Stokes drift profiles **u**^{S}(*z*) are specified from the 10-m wind speed *U*_{10} and the surface wave age *C _{P}*/

*U*

_{10}, where

*C*is the spectral peak phase speed, using an empirical surface wave spectra and an associated wave age–dependent neutral drag coefficient

_{P}*C*. Wave-breaking effects are not otherwise included. Mixed layer depths

_{D}*H*

_{ML}vary from 30 to 120 m, with 0.6 ≤

*C*/

_{P}*U*

_{10}≤ 1.2 and 8 m s

^{−1}<

*U*

_{10}< 70 m s

^{−1}, thereby addressing most possible oceanic conditions where TKE production is dominated by wind and wave forcing.

The mixed layer–averaged “bulk” VKE 〈*w*^{2}〉/*u**^{2} is equally sensitive to the nondimensional Stokes *e*-folding depth *D**_{S}/*H*_{ML} and to the turbulent Langmuir number La_{t} = *u**/*U _{S}*, where

*u** = |

*τ*_{0}|/

*ρ*in water density

_{w}*ρ*and

_{w}*U*= |

_{S}**u**

^{S}|

_{z=0}. Use of a

*D**

_{S}scale-equivalent monochromatic wave does not accurately reproduce the results using a full-surface wave spectrum with the same

*e*-folding depth. The bulk VKE for both monochromatic and broadband spectra is accurately predicted using a surface layer (SL) Langmuir number La

_{SL}=

*u**/〈

*u*

^{S}〉

_{SL}, where 〈

*u*

^{S}〉

_{SL}is the average Stokes drift in a surface layer 0 >

*z*> − 0.2

*H*

_{ML}relative to that near the bottom of the mixed layer. In the wave-dominated limit La

_{SL}→ 0, turbulent vertical velocity scales as

*w*

_{rms}∼

*u**La

^{−2/3}

_{SL}. The mean profile (

*z*) of VKE is characterized by a subsurface peak, the depth of which increases with

*D**

_{S}/

*H*

_{ML}to a maximum near 0.22

*H*

_{ML}as its relative magnitude /〈

*w*

^{2}〉 decreases. Modestly accurate scalings for these variations are presented. The magnitude of the crosswind velocity convergence scales differently from VKE. These results predict that for pure wind seas and

*H*

_{ML}≅ 50 m, 〈

*w*

^{2}〉/

*u**

^{2}varies from less than 1 for young waves at

*U*

_{10}= 10 m s

^{−1}to about 2 for mature seas at winds greater than

*U*

_{10}= 30 m s

^{−1}. Preliminary comparisons with Lagrangian float data account for invariance in 〈

*w*

^{2}〉/

*u**

^{2}measurements as resulting from an inverse relationship between

*U*

_{10}and

*C*/

_{P}*U*

_{10}in observed regimes.

## 1. Introduction

Exchanges of heat, water, momentum, and chemical species between the atmosphere and the ocean interior are mediated by mixing within the upper ocean boundary layer. This study seeks to quantify the role of surface waves in setting the level of turbulent kinetic energy (TKE) in this layer. This TKE level figures prominently in many ocean boundary layer models, including turbulence closure schemes of Mellor and Yamada (1982) and the *K*-profle parameterization (KPP; Large et al. 1994). These models ignore any surface wave effects and set the level of TKE to that found in solid-wall boundary layers. In contrast, observations (D’Asaro 2001; Tseng and D’Asaro 2004) show that the vertical turbulent kinetic energy (VKE) magnitude is higher in the ocean boundary layer than in solid-wall boundary layers with the same applied stress. We hypothesize that Langmuir circulations, generated through the Craik–Leibovich (CL) mechanism (Craik and Leibovich 1976), can quantitatively explain the enhanced level of VKE. Here, large-eddy simulations (LESs) are used to develop accurate scalings for the VKE enhancement under realistic wind and wave forcing.

### a. The Craik–Leibovich mechanism and Langmuir turbulence

In the CL mechanism, Langmuir circulations arise from the interaction of the Stokes drift **u**^{S} of surface waves and wave-averaged currents driven by a surface stress *τ*_{0} = |*τ*_{0}| = *ρ _{w}*

*u**

^{2}, where

*u** is the friction velocity in water of density

*ρ*(see Thorpe 2004 for a review). The CL interaction is represented in the momentum equations by a “vortex force,”

_{w}perpendicular to both the Stokes drift **u**^{S} and the vorticity *ζ** ^{E}* =

**∇**×

**u**

*of the Eulerian mean flow*

^{E}**u**

*= [*

^{E}*u*,

^{E}*υ*,

^{E}*w*]. While the Eulerian mean is obtained by averaging fluid velocity over the wave phase in a fixed reference frame, the Lagrangian mean velocity

^{E}**u**

*≅*

^{L}**u**

*+*

^{E}**u**

*is the average obtained in a wave-following Lagrangian reference frame.*

^{S}Early work (Craik and Leibovich 1976; Leibovich 1980, 1983) envisioned Langmuir circulations as coherent structures composed of counterrotating pairs of surface roll vortices aligned downwind and separated by a crosswind convergence zone along a downwind surface jet. Subsequent studies (Skyllingstad and Denbo 1995; McWilliams et al. 1997) using high-resolution nonhydrostatic LES methods have introduced the concept of Langmuir turbulence, in which the CL vortex force drives a turbulent boundary layer flow, where coherent Langmuir circulation structures form and dissipate locally and episodically. Vortex-force TKE production was shown to increase equilibrium mixed layer TKE levels with a decreasing turbulent Langmuir number La_{t} = *u**/*U _{S}*, where

*u** = |

*τ*_{0}|/

*ρ*and

_{w}*U*= |

_{S}**u**

^{S}|

_{z=0}. For wind- and wave-dominated regimes and typical open ocean values of La

_{t}, the CL mechanism acts to increase levels of mean VKE significantly over simulations without Stokes drift. (Overbars indicate averages over a homogeneous ensemble.) The General Lagrangian Mean (GLM) theory of Andrews and McIntyre (1978) provides an exact treatment of the wave–current interaction in terms of the pseudomomentum

**p**of the waves. The CL mechanism is an approximation of the inviscid GLM forces that are valid to

*O*(

*ξ*

^{3}) for small values of

*ξ*=

*A*/

*λ*, where

*A*is the wave amplitude and

*λ*its wavelength, in which Stokes drift components

*u*

^{S}

_{i}=

*u*

^{L}

_{i}−

*u*

^{E}

_{i}approximate

*p*when the waves are irrotational. In this study we further approximate the Stokes drift to

_{i}*O*(

*ξ*

^{2}).

Prior numerical work has focused primarily on simulating the response of mixed layer turbulence to a surface stress aligned with a monochromatic deep-water wave of frequency *ω* and wavenumber **k**. The Stokes drift at level *z* is

to *O*(*A*^{2}), or 2*ω***k***η*^{2}_{rms}*e*^{2kz} for rms surface displacement *η*_{rms}, where *k* = *ω*^{2}/*g* and *g* is gravity. This study includes some such simulations of monochromatic waves, which may be interpreted as an approximation to CL forcing due to waves with broadly distributed energy spectra, or as more accurate representations of forcing due to narrowband swell. The use of monochromatic waves with wavenumber *k** to approximate broadband CL forcing is predicated on the dynamic equivalence of Stokes drift profiles with the same *e*-folding depth scale *D**_{S} = 1/2*k** (Li and Garrett 1993).

This study focuses primarily on accurately simulating the response to pure wind seas accross a realistic range of oceanic conditions, without invoking the monochromatic approximation. Pure wind seas have broadband “equilibrium” displacement spectra *F*(*ω*, *θ*) that represent the cumulative effect over some fetch or duration of surface stress due to local wind conditions. To *O*(*A*^{2}), such broadband wave spectra have Stokes drift

composed of a Stokes drift spectrum with elements *δ***U**^{S}/*δ**ω*, which each decay exponentially below the surface (Kenyon 1969).

### b. Previous work: Scalings from LES models

Previous LES-based studies have suggested various scalings for the TKE components as a function of the Langmuir number. A modification of the vertical velocity scale “W” in the KPP mixed layer parameterization was proposed in McWilliams and Sullivan (2000), based on LES results at La_{t} = 0.3 and observational evidence that suggested *w*_{rms} ∼ *U*_{S} for La_{t} ≪ 1. This proposed dependence of vertical velocity on La_{t} implies here that

although McWilliams and Sullivan note as well that the dependence of *w* on *U*_{S} entailed by this in the absence of surface stress is not plausible. Smyth et al. (2002) examine modifications to the KPP mixed layer model based on Eq. (4) and extensions to it.

Min and Noh (2004) compare LES results for the near-surface crosswind velocity *υ*_{rms}|_{z≅0} in rotating Langmuir turbulence at La_{t} = [0.23, 0.25, 0.35, 0.45] with a set of scaling predictions:

for *m* = [1, 2, 3], and characterize these as more consistent with *υ*_{rms}|_{z≅0} ∝ *u**^{2/3}*U*^{1/3}_{S} (*m* = 3) than they are with *υ*_{rms}|_{z≅0} ∝ *u**^{1/2}*U _{S}*

^{1/2}(

*m*= 2) or

*υ*

_{rms}|

_{z≅0}∝

*U*(

_{S}*m*= 1). However, these results are somewhat ambiguous because they are based on a set of simulation cases with two different

*D**

_{S}/

*H*

_{ML}values. Moreover, the proportional scaling form Eq. (5) is relevant to La

_{t}≪ 1, but may not be the best comparison for their larger (La

_{t}= 0.45) Langmuir number case.

Li et al. (2005) simulate 13 model cases at *D**_{S}/*H*_{ML} = 0.12 with *H*_{0} = 0, and graphically present (their Fig. 4) LES results for the equilibrium level of bulk VKE 〈*w*^{2}〉, obtained by averaging *w*^{2} first over the horizontal domain and the mixed layer depth *z* > −*H*_{ML}(*t*) and then in time over a statistically steady-state period. We find that

in Fig. 1 roughly fits their LES data (courtesy of M. Li 2006, personal communication). Also shown in Fig. 1 is a curve inferred from the modification of the KPP vertical velocity scale suggested by McWilliams and Sullivan (2000) [Eq. (4)], adopting the same zero-wave La_{t} → ∞ limit as Eq. (6). They do not agree at low La_{t} (i.e., when the CL force is important).

Several LES studies have examined the sensitivity of TKE to the nondimensional ratio *D**_{S}/*H*_{ML}, where *H*_{ML} is the mixed layer depth. Skyllingstad et al. (2000) presents simulations with different Stokes depth scales that result in different rates of deepening. McWilliams and Sullivan (2000) mention sensitivity tests that suggest a TKE dependence on the relative Stokes depth for *D**_{S}/*H*_{ML }*greater than* the value of 0.13 in their simulations, at odds with subsequent reports. Li et al. (2005) compare three cases at *D**_{S}/*H*_{ML} = [0.24, 0.12, 0.06] at a fixed Langmuir number La_{t} = 0.34. They report little difference in the profile of VKE components between the two largest values of *D**_{S}/*H*_{ML} but a significant reduction in for *D**_{S}/*H*_{ML} = 0.06. While Li et al. suggest a sensitivity at small *D**_{S}/*H*_{ML}, the concurrent impact on the horizontal components of TKE from variations in *D**_{S}/*H*_{ML} appears uniformly minor. Min and Noh (2004) vary *D**_{S}/*H*_{ML} = 0.127 and 0.064 at La_{t} = 0.35 (their cases “A” and “B”), *D**_{S}/*H*_{ML} = 0.127 at La_{t} = 0.25 (“E”), and *D**_{S}/*H*_{ML} = 0.064 at La_{t} = 0.23 (“C”). Their Fig. 3 shows little change in *υ*_{rms}|_{z≅0} between A and B, but a significant difference between C and E, thereby supporting a strong *D**_{S}/*H*_{ML} dependence. As a whole, these test cases are sparse and inconsistent.

### c. Previous work: Scalings from data

Upper ocean Doppler measurements by Smith (2001) suggest *υ*_{rms}|_{z≅0} ∼ *U _{S}* for near-surface crosswind convergence velocity, but this behavior appears only episodically, and with a varying offset in the magnitude of

*υ*

_{rms}|

_{z≅0}/

*u**. Pluddemann et al. (1996) suggests

*υ*

_{rms}|

_{z≅0}∼

*u**

^{1/2}

*U*

^{1/2}

_{S}. There appears to be no consensus here.

D’Asaro (2001) and Tseng and D’Asaro (2004) report 〈*w*^{2}_{rms}〉 ∝ *u**^{2} from mixed layer Lagrangian float measurements, consistently across a wide variety of ocean-forcing regimes, including recent observations from hurricanes (D’Asaro 2003; D’Asaro and McNeil 2006). The strength of this result and the surprising lack of dependence on surface wave properties provide the major motivation for the modeling studies presented here.

### d. Wave breaking

Observations of the dissipation rate ɛ of TKE near the ocean surface show an increase above the “law-of-the-wall” expectation ɛ = *u**^{3}/*κ**z* by up to two orders of magnitude (Agrawal et al. 1992; Anis and Moum 1995; Drennan et al. 1996; Melville 1996; Terray et al. 1996). This increase is generally attributed to surface wave breaking, either at the macroscale associated with white caps (Rapp and Melville 1990; Duncan 1981) or at the microscale (Jessup et al. 1997). An increase is also expected due to additional TKE production from the vortex force. This additional contribution from the CL mechanism modifies the mixed layer TKE budget,

adding a contribution

because of vortex-force production, to the standard contributions from shear production *P*_{Shear} = −∂/∂*x _{j}*, buoyant production

*P*

_{Buoy}, and the divergence

*D*

_{Tran.}of TKE transport. Reported LES results indicate that vortex-force TKE production increases near-surface dissipation severalfold over the law-of-the-wall scaling. This, however, falls far short of the observed near-surface increase in TKE dissipation. Other TKE production processes, aside from CL vortex-force production, must be responsible for balancing the observed increase in dissipation. However, it is a separate question to ask whether these other processes, such as wave breaking, should also be expected to significantly increase bulk mixed layer TKE

*levels*. Answering this question by model–data comparison entails including a model of wave breaking either explicitly in the resolved equations or as a source of unresolved turbulent energy, and then comparing the resulting model turbulence levels to observations with a clear correspondence to model-forcing conditions.

Our primary goal is to predict VKE levels over a wide range of oceanic conditions so that these may be compared to observations to test underlying theories of the forcing processes. Several wave-breaking parameterizations have been proposed and included in mixed layer turbulence simulations with varying results (Noh et al. 2004; Sullivan et al. 2004a, b, 2005). However, the correct relative distribution of near-surface TKE production by breaking waves over resolved and unresolved model scales is not yet well established, nor is the correspondence between such model forcing and varying sea-state conditions. We therefore chose to forgo the additional uncertainties that would be introduced by including parameterizations of wave breaking in model forcing, thereby adopting as a working hypothesis that TKE production by wave breaking does not significantly impact bulk VKE levels or entrainment rates. Ultimately, comparisons with ocean data will be required to establish the significance of omitting contributions to turbulence production from wave breaking or other surface processes.

Section 2 reviews the monochromatic-forcing paradigm and describes a new spectral-forcing approach. Section 3 describes our numerical techniques. Section 4 describes a set of simulation cases spanning a wide range of wave age and wind speed, using both monochromatic- and spectral-forcing paradigms. Analysis of these results in section 5 yields a unified scaling approach that accounts for variability in the layer-averaged bulk VKE due to changes in the Langmuir number and the vertical distribution of Stokes drift. Section 5 also examines how the shape of the mean-VKE profile changes as a function of these parameters. Results are reformulated in terms of wind speed, wave age, and mixed layer depth, then followed by a summary and conclusions in section 6.

## 2. Estimation of Stokes drift profile and wind stress

Quantitative model–data comparisons designed to test wave–current interaction theories require a clear correspondence between ocean and model forcing due to **u**^{S}(*z*) and *τ*_{0}. Without a clear correspondence, the meaning of any agreement or discrepancy between predicted and observed turbulence statistics becomes less clear. Displacement spectra used to compute Stokes drift can be decomposed into

into a scalar component Φ(*ω*) and a normalized spreading function *G _{N}*(

*ω*,

*θ*). Stokes drift calculations are sensitive to certain details of these functions. In particular, the Phillips (1958) scaling Φ(

*ω*) ∼

*ω*

^{−5}for the scalar spectrum at high frequencies contributes much less to the surface Stokes drift than the scaling Φ(

*ω*) ∼

*ω*

^{−4}of Toba (1973), which, for a frequency-independent

*G*, requires truncation at the capillary limit for the Stokes drift integral [Eq. (3)] to converge at

_{N}*z*= 0. Attention to this behavior of the spectral tail is therefore critical to quantitative model–data comparisons that depend on

*U*.

_{S}### a. The monochromatic-forcing paradigm

Previous studies of Langmuir turbulence use a monochromatic downwind Stokes drift profile:

Li and Garrett (1993) assume the dynamic equivalence of Stokes drift profiles **u**^{S}(*z*) with the same *e*-folding depth scale *D**_{S} = 1/2*k**, based on Pierson–Moskowitz (PM) spectrum for fully developed seas, where the phase speed *C _{p}* =

*ω*/

_{p}*k*of spectral peak waves at

_{p}*k*=

_{p}*ω*

^{2}

_{p}/

*g*reaches

*C*= 1.2

_{p}*U*

_{10}. For a range of empirical parameters, Li and Garrett estimate the Stokes wind ratio as

and the *e*-folding Stokes depth scale as

for wind speed *U _{W}*. Adjusting for an apparently different wind reference height

^{1},

*D**

_{S}= 0.14

*U*

^{2}

_{10}/

*g*or 0.19 (2

*k*)

_{p}^{−1}, and

*U*/

_{S}*U*

_{10}= (0.015 to 0.016). Using

*C*from Large and Pond (1981) gives

_{D}For *U _{S}*/

*U*

_{10}= 1.55%, La

_{t}ranges from 0.28 at

*U*

_{10}= 10 m s

^{−1}to 0.33 at

*U*

_{10}= 30 m s

^{−1}, while

*D**

_{S}increases from 1.4 to 13 m, as illustrated in Fig. 2.

Figure 2 also plots the choices of *D**_{S} and La_{t} from the subsequent LES studies of Skyllingstad and Denbo (1995), McWilliams et al. (1997), Min and Noh (2004), Noh et al. (2004), and Li et al. (2005), restricted to model cases with negligible surface heat flux. The range of simulated Langmuir numbers corresponds predominately to values 0.2 < La_{t} < 0.5 reported for field measurements (Smith 1992). The focus on depth scales 3.2 m < *D**_{S} < 6.4 m would imply fully developed seas at 15 m s^{−1} < *U*_{10} < 21 m s^{−1}. Many of these combinations are unlikely to be realized, while many possible combinations are not sampled. We will attempt a more systematic approach.

### b. A spectral-forcing paradigm

A modified version of the spectrum *F*(*ω*, *θ*) reported in Eqs. (9.1)–(9.4) of Donelan et al. (1985) is used to compute *u ^{S}*(

*z*) as a function of varying wind speed

*U*

_{10}and wave age

*C*/

_{p}*U*

_{10}. This empirical spectrum, based on measurements at

*ω*< 4

*ω*, combines

_{p}with a nearly normalized spreading function

where the shape parameters *α*, *γ*, and *σ* depend on ϖ = *ω*/*ω _{p}* and on wave age

*C*/

_{p}*U*

_{10}. Banner (1990) suggests that the Φ(

*ω*) ∼

*ω*

^{−4}behavior of the wind wave spectrum in Eq. (14) should transition to Φ(

*ω*) ∼

*ω*

^{−5}at higher frequencies, and it achieves this by replacing

in the Donelan et al. (1985) argument of sech^{2}, with

The transition from *ω*^{−4} to *ω*^{−5} scaling is then effectively due to changes in the ratio

between the radial norms of the two spreading functions [Χ(ϖ) = 1 for ϖ ≤ 1.6]. In a revision of this formulation, we first resolve the ambiguity between Φ(*ω*) and *F*(*ω*, *θ*) by normalizing the modified spreading function

exactly to one, with *β* = *β*_{Ban} substituted at ϖ ≤ 1.6. The transition to Φ(*ω*) ∼ *ω*^{−5} is then accomplished by modifying the frequency spectrum by Χ(ϖ):

Wave energy *η*^{2}_{rms} = ∫^{∞}_{0} Φ′(*ω*) *d**ω* for this spectrum is not significantly changed from that in Donelan et al. (1985), *η*^{2}_{rms}*g*^{2}/*U*^{4}_{10} = (6.365 × 10^{−3})[*U*_{10}/(2*π**C _{p}*)]

^{−3.3}: for wave ages

*C*/

_{p}*U*

_{10}= [0.6, 0.8, 1.0, 1.2], finescale discrete integration gives nondimensional energy levels

*η*

^{2}

_{rms}

*g*

^{2}/

*U*

^{4}

_{10}= [0.525, 1.36, 2.63, 5.21] × 10

^{−3}, respectively. Figure 3a compares the Donelan et al. (1985) scalar energy spectrum with the adjusted Donelan–Banner (DB) spectrum and the wave age–dependent Joint North Sea Wave Project (JONSWAP) spectrum (Hasselmann et al. 1973), as reformulated by Lewis and Allos (1990) at

*C*/

_{p}*U*

_{10}= 1.17 and

*C*/

_{p}*U*

_{10}= 0.6. The above adjustment to the high frequencies is minor in the energy spectrum, while large differences are apparent between the DB spectrum and JONSWAP, particularly for

*C*/

_{p}*U*

_{10}= 1.17, where JONSWAP and PM spectra are nearly the same.

Stokes drift profiles are found by substituting this Donelan–Banner spectrum *F*′(ϖ, *θ*) = Φ′(ϖ)*G _{N}*(ϖ,

*θ*) into Eq. (3). Scaling the DB Stokes drift on

*U*

_{10}gives

The integral is evaluated by discrete approximation over 0 < ϖ_{i} < ϖ_{cut} as **u**^{S} ≅ Σ(*δ***U*** ^{S}*/

*δ*ϖ)

_{ϖi}

*e*

^{2kiz}Δϖ

_{i}with

*i*= 1, . . . ,

*N*

_{ϖ}frequencies. This is combined with an analytical tail contribution from above ϖ

_{cut}> 20, approximated by

Figure 3b shows how the DB Stokes drift spectra differ significantly from the JONSWAP form, both for young and mature seas. Note that Stokes drift spectra for young seas are substantially narrower than those of mature seas. The effect of the adjustment at high frequencies in the DB spectrum is also apparent here, where it serves to keep the Stokes drift integral Eq. (21) finite.

Figure 3c plots the surface Stokes wind ratio *U*^{S}/*U*_{10} evaluated from Eq. (21) at *z* = 0. It depends only on *C _{p}*/

*U*

_{10}. A finescale approximation was evaluated with

*N*

_{ϖ}= 6000 and Δϖ/

*ω*= 0.001. Adding a small tail contribution above ϖ

_{cut}= 100 gives the accurate values plotted in Fig. 3c (solid curve). For the LES model forcing, a slightly less accurate estimation used

*N*

_{ϖ}= 50 discrete frequencies separated by a larger interval Δϖ/ϖ = 0.1 and cut off ϖ

_{cut}= 0.25 (1 + Δϖ/ϖ)

^{Nϖ−1}. Plotted in Fig. 3c (circles) and listed below in Tables 1 –3, these values are not significantly different. The

*e*-folding depth scale

*D**

_{S}was determined iteratively from the Stokes drift computed at

*z*. Figure 3d shows 2

*D**

_{S}

*k*versus

_{p}*C*/

_{p}*U*

_{10}, scaling the Stokes depth on the (significantly larger) Stokes depth (2

*k*)

_{p}^{−1}of the dominant wavenumber

*k*=

_{p}*g*/

*C*

^{2}

_{p}. For typical ocean wave ages

*C*/

_{p}*U*

_{10}> 0.6,

*e*-folding depths for the DB spectrum are 26%–37% larger than the 2

*D**

_{S}

*k*= 0.19 value obtained by Li and Garrett (1993) for the fully developed PM spectrum.

_{p}The above DB model has a slightly larger surface Stokes wind ratio *U*^{S}/*U*_{10} = 1.75% at *C _{p}*/

*U*

_{10}= 1.2 (Fig. 3c) than the value

*U*

^{S}/

*U*

_{10}= (1.5 to 1.6)% entailed by the analysis of Li and Garrett (1993) using the PM scalar spectrum. Field measurements in Smith (2006) obtain Stokes drift from the difference between near-surface Lagrangian advection of bubble cloud features and their mean Eulerian Doppler velocity. Fitting the observed correlation between the wave drift of bubbles and winds around

*U*

_{10}= 10 m s

^{−1}, Smith obtains

*U*/

^{S}*U*

_{10}= 1.25%. Adjusting for a 1.5-m

*e*-folding decay depth for near-surface bubble concentration, and for other acoustic sheltering effects, Smith infers a surface value 0.92/0.75 larger than

*U*

^{S}/

*U*

_{10}= 1.53%. This corresponds to less than fully developed waves,

*C*/

_{p}*U*

_{10}= 1.0 in the DB model. The difference may indicate that DB Stokes drift is slightly high, that the estimated effects of bubble cloud thickness or sheltering are low, or that the waves measured by Smith are not fully mature. It may also simply be that the Smith observations were made in marginally less pure wind seas than the measurements of Donelan et al. (1985), and a greater directional spread in the surface wave spectrum may reduce the Stokes drift. More generally, any single surface wave spectrum will be observed to differ substantially from the empirical spectrum with associated constants above. Choosing this particular form permits an examination of how, on average, sea-state variations such as wave age may impact mixed layer turbulence. Further comparisons, such as those made with monochromatic wave–forcing cases, can lead to more general predictions of how turbulent mixing varies with arbitrary surface waves.

Is the CL mechanism an accurate approximation to the exact GLM wave–current interaction theory when applied to such empirical spectra? While it is straightforward to test that the GLM expansion parameter *ξ* = *A*/*λ* = *k**η*_{rms}/*π* is small for monochromatic waves, the most appropriate measure based on integrating the slope spectrum *k*^{2}*F*(*ω*, *θ*) cos(*θ*) diverges logarithmically for wind seas with Φ(*ω*) ∼ *ω*^{−5} at high frequencies. However, if the integrals determining both the Stokes drift [Eq. (3)] and a generalized parameter

are truncated at *ω*_{cut} = 100*ω _{p}*, almost all the Stokes drift is included while the values of

*ξ*(100

_{I}*ω*) are still small. As a function of wave age,

_{p}*ξ*(100

_{I}*ω*) ranges from 0.055 for

_{p}*C*/

_{p}*U*

_{10}= 1.2 to 0.066 for

*C*/

_{p}*U*

_{10}= 0.6 in the DB-based simulations below, while accounting for over 99% of

*U*. For practical purposes, this integral formulation of the GLM expansion parameter may therefore still be considered small. In any case, applying CL theory thus to the Stokes drift of an empirical broadband wave spectrum commits no theoretical errors not already implied by comparing ocean data to LES results with equivalent monochromatic forcing, and it corrects several gross simplifications in the monochromatic approach.

_{S}A wind stress consistent with the wave-age dependence of the DB model Stokes drift is the neutral drag coefficient *C _{D}* of Donelan et al. (1993), which is a function of

*C*/

_{p}*U*

_{10}. Surface stress

*τ*

_{0}=

*ρ*

_{a}*C*

_{D}*U*

^{2}

_{10}in air of density

*ρ*is evaluated from surface roughness length

_{a}*z*

_{0}using

Figure 3e plots this age-dependent drag coefficient versus *U*_{10}, along with the somewhat smaller wave-independent neutral drag coefficients of Smith (1988) and Large and Pond (1981). A more recent analysis of a larger surface flux dataset in Drennan et al. (2003) obtains similar dependence of *C _{D}* on

*C*/

_{p}*U*

_{10}. The nondimensionalization of LES results shown below will permit the substitution of other surface drag formulations in model–data comparisons. Eqs. (24) and (25) are not expected to hold for the very young wave fields in laboratory tank conditions, nor for very high winds

*U*

_{10}> 30 m s

^{−1}where evidence suggests that the increase of

*C*with wind speed “saturates” at a maximum value.

_{D}## 3. LES model and vortex-force formulation

The LES model used here was originally developed by Deardorff (1980) and Moeng (1984) to simulate atmospheric boundary layers. It combines pseudospectral horizontal resolution with vertical finite differencing, and it employs an advected subgrid TKE to parameterize unresolved fluxes. Similar oceanic adaptations of this LES model are used by Skyllingstad and Denbo (1995) and McWilliams et al. (1997). The version employed here was adapted for the ocean by Garwood et al. (1994) and Harcourt (1999). Lherminier et al. (2001) and Harcourt (1999) use this model without Stokes drift forcing to study deep convection in polar oceans and the response of isobaric and Lagrangian drifters. Simulations of Labrador Sea deep convection using this LES model and observed meteorological forcing were found to be in close agreement with VKE measured under the Lagrangian float data in Harcourt et al. (2002). Most recently, LES modeling in Harcourt (2005) examines thermobaric cabbeling instabilities below the wintertime Antarctic ice cover. These previous LES studies do not include surface wave effects.

The abbreviated model description provided in appendix A serves to describe the unique numerical implementation of the Craik–Leibovich force used here. A source term for the vortex-force production is included in the prognostic equation for the advected subgrid TKE *e* [Eq. (A6)], which is used to effect turbulence closure. This makes the subgrid TKE budget consistent with the resolved scales [Eq. (8)]. The Stokes drift is filtered to the grid scale, making it more formally consistent with LES variables as the representation of quantities filtered over a volume corresponding to the model resolution (see Moeng and Wyngaard 1988). Appendix B describes how the unfiltered Stokes drift **u**^{S} is replaced by a grid-filtered quantity 〈**u*** ^{S}*〉

_{Δ}, and how this step is applied to model forcing described by a discretely approximated Stokes drift spectrum.

For monochromatic Stokes drift forcing with *D**_{S} > Δ*z*, there is little difference between *u _{S}*(−Δ

*z*/2) and 〈

**u**

^{S}〉

_{Δ}at

*z*= −Δ

*z*/2, but for the DB spectral model, at the vertical resolution (Δ

*z*= 1.42) of most model cases below, using unfiltered

*u*(−Δ

_{S}*z*/2) would underestimate 〈

*u*(−Δ

_{S}*z*/2)〉

_{Δ}by between 1% (for

*U*

_{10}= 32.6 m s

^{−1}and

*C*/

_{p}*U*

_{10}= 1.2) and 22% (for

*U*

_{10}= 8.3 m s

^{−1}and

*C*/

_{p}*U*

_{10}= 0.6). By representing the Stokes drift in the equations for resolved momentum by its filtered form [Eq. (B1)], we avoid a corresponding underestimate of vortex- force TKE production in Eq. (8).

## 4. The LES model cases

All simulations were run with constant forcing to determine the steady-state dynamic scaling of low-order turbulence statistics for a wind- and wave-forced ocean mixed layer. Four major case sets were performed, identified as “Σ_{1}” through “Σ_{4}.” Each is described in turn below and in Tables 1 –4. Model forcing is varied with dynamical scales *U*_{10} and *C _{p}*/

*U*

_{10}over ranges typical of mixed layers, where TKE production is dominated by wind and wave forcing. The choices of a model for the wave spectrum and surface drag

*C*reduce the forcing parameter space by determining La

_{D}_{t},

*D**

_{S},

*U*/

_{S}*U*

_{10}, and

**u**

^{S}/

*U*as functions of

_{S}*U*

_{10}and

*C*/

_{p}*U*

_{10}. The nondimensional ratio

*D**

_{S}/

*H*

_{ML}depends on an average

*H*

_{ML}value that is only partially determined by initial conditions.

All simulation cases use a 288 m × 288 m horizontally periodic domain with radiative bottom boundary conditions (Klemp and Duran 1983). The vertical extent of the domain was varied to accommodate different mixing depths. Unless otherwise indicated for several cases repeated at a high resolution, numerical domains are resolved using grid intervals of Δ*x* = Δ*y* = 1.5 m and Δ*z* = 1.42 m. The time step Δ*t* varies among simulation cases as indicated in Tables 1 –4. Surface wind and Stokes drift vectors all point north, and planetary rotation **Ω** is fixed at latitude 48°N. Velocity components *u* and *υ* are thereby crosswind and downwind, respectively, but surface crosswind velocity statistics refer here to “*υ*” rather than “*u*,” for notational consistency with previous Langmuir turbulence studies.

### a. Primary simulation case set Σ_{1}

In the primary set Σ_{1} (Table 1), 28 simulation cases are composed on a 8 × 4 grid of seven *U*_{10} values versus four values of wave age. Surface stress and Stokes drift profiles were determined as functions of *U*_{10} and *C _{p}*/

*U*

_{10}using the spectral-forcing paradigm described in section 2b. Wave ages ranged over typical oceanic values of 0.6 <

*C*/

_{p}*U*

_{10}< 1.2 for storm forcing, with a wind speed range of 8 m s

^{−1}<

*U*

_{10}< 33 m s

^{−1}. Figure 4 plots

*D**

_{S}as a function of the Langmuir number for simulation set Σ

_{1}. All domains of Σ

_{1}extended vertically to

*z*= −113.92 m, and mean mixed layer depths for steady-state turbulence statistics varied from 61 to 82 m.

### b. Monochromatic simulation case set Σ_{2}

For comparison with Σ_{1}, steady-state solutions were obtained for an auxiliary case set Σ_{2} (Table 2), with Stokes drift wave forcing given by a monochromatic wave with depth scale *D**_{S} as described in section 2a. The domains of Σ_{2} extended vertically to *z* = −113.92 m as in Σ_{1}, accommodating mean mixed layer depths between 60 and 77 m. Stokes drift profiles for Σ_{2} were calculated using *U*_{10} values from Σ_{1} as wind speed *U _{W}* in the PM-based expressions of Li and Garrett [(1993); Eqs. (18) and (20) therein]. The midpoint of 1.45% is chosen from the range 1.4%–1.5% suggested by Li and Garrett for

*U*/

_{S}*U*. Surface stress was determined using the neutral drag coefficient of Large and Pond (1981), modified to extend the linear drag relation to the lowest (

_{W}*U*

_{10}< 10 m s

^{−1}) values. The resulting dependence of

*D**

_{S}(

*U*

_{10}) versus La

_{t}(

*U*

_{10}) for set Σ

_{2}is also plotted in Fig. 4. A belated realization that the formulas of Li and Garrett are apparently based on 19.5-m winds means that the reference height corresponding to the calculated wave forcing is different from that in Σ

_{1}. However, this also means that values of

*u** and La

_{t}are additionally low for neutral conditions in this case set Σ

_{2}by 5%–8% and 3%–4%, respectively, because 19.5-m winds were thereby used in a drag formula for winds referenced to 10 m. This set is nevertheless useful for distinguishing the effects of forcing from wave drift profiles due to monochromatic versus broadband wave spectra.

### c. High-wind case set Σ_{3}

A third case set is composed of two subsets, Σ_{3a} and Σ_{3b}, distinguished by different drag coefficients at high winds, *U*_{10} > 20 m s^{−1}. Set Σ_{3a} continues the *C _{p}*/

*U*

_{10}= 0.6 subset of Σ

_{1}over three hurricane-strength wind speeds, 33 m s

^{−1}<

*U*

_{10}< 70 m s

^{−1}, while Σ

_{3b}uses a constant

*C*= 2.3 × 10

_{D}^{−3}, slightly below high-wind saturation values of (2.5 − 2.8) × 10

^{−3}suggested by Donelan et al. (2004). Solutions at

*C*/

_{p}*U*

_{10}= 0.6 from Σ

_{1}are contiguous with Σ

_{3a}; Σ

_{3b}connects with Σ

_{1}near

*U*

_{10}= 20 m s

^{−1}, where

*C*= 2.3 × 10

_{D}^{−3}for both. The forcing scales for Σ

_{3}are detailed in Tables 3a, b. Turbulence statistics were obtained for mean mixed layer depths between 81 and 118 m in domains between 170.87 and 227.85 m deep.

### d. Reduced mixed layer depth set Σ_{4}

A fourth case set repeats the *C _{p}*/

*U*

_{10}= 1.2 subset of Σ

_{1}for mature seas, but with the initial mixed layer depth

*H*

_{ML}cut in half by compressing the initial temperature and salinity profiles. The forcing scales for Σ

_{4}are detailed in Table 4. Average final mixed layer depths of 30–56 m in Σ

_{4}and 61–82 m in Σ

_{1}can be determined from the tabulated ratios using

*k*

^{−1}

_{p}=

*C*

^{2}

_{p}g

^{−1}. The vertical domains extended to 56.96 and 85.44 m for mixed layer depths averaging between 30 and 56 m. A comparison between this set and Σ

_{1}(

*C*/

_{p}*U*

_{10}= 1.2) distinguishes

*D**

_{S}/

*H*

_{ML}dependence from Δ

*z*/

*H*

_{ML}effects.

### e. Initial and final conditions

All the numerical simulations of the primary set Σ_{1} used as their initial condition a quasi-steady-state solution drawn from the *U*_{10} = 21.3 m s^{−1} case of the PM-based set Σ_{2}. Initial conditions for set Σ_{2} were similarly spun up using the surface stress of *U*_{10} = 11 m s^{−1} and a Stokes drift profile derived from observations of incidental origin. Initial conditions for the high-wind cases in Σ_{3} were drawn from the final state of the *U*_{10} = 32.6 m s^{−1}, *C _{p}*/

*U*

_{10}= 0.6 case in Σ

_{1}. Initial conditions for Σ

_{4}were generated by vertically compressing the final states from the

*C*/

_{p}*U*

_{10}= 1.2 set of Σ

_{1}in half, while adjusting velocity fields for incompressibility. These choices were made primarily for the sake of expediency, and the initial conditions of turbulence are not deemed to influence the final state beyond altering its mean hydrography.

Figure 5 shows final profiles of mean surface-relative density *ρ*_{0}(*z*) − *ρ*_{0}(−Δ*z*/2) from several example cases. Variations in both forcing strengths and initial states for each case lead to different final mixed layer depths, *H*_{ML}, indicated for each example profile. Here, *H*_{ML} is defined as the greatest depth above the actively entraining pycnocline, where the stability *N* ^{2} of horizontally averaged profiles is less than 10% of its maximum value within the pycnocline. This depth tends to scale with the penetration depth of TKE from the surface layer (SL), while largely excluding the (smaller) TKE levels produced by shear within the pycnocline. Here, *H*_{ML} is approximately 0.8 of that depth, where both mean stability and the magnitude of entrainment fluxes are largest.

## 5. Simulation results

### a. Bulk VKE averages

Figure 6 plots the *u**-scaled, time and bulk mixed layer depth–averaged vertical TKE 〈*w*^{2}〉/*u**^{2} (VKE) from all simulation sets as a function of the Langmuir number, grouped by wave age. Mean profiles of VKE for quasi-steady-state entraining solutions were obtained by averaging *w*^{2} horizontally over the LES domain. A time series of the mixed layer average 〈*w*^{2}〉 is obtained by averaging these profiles over the mixed layer 0 > *z* > −*H*_{ML}, every 50 time steps. The final steady-state average of 〈*w*^{2}〉 is obtained over a terminal time interval *T*_{avg} ≫ *H*_{VKE}/〈*w*^{2}〉, the efficacy of which was evaluated by inspection. Our fit [Eq. (6); Fig. 1] to the LES data of Li et al. (2005) is included for comparison. This curve passes near-current results at similar depth scales *D**_{S}/*H*_{ML} ≅ 0.12.

Within each case set, 〈*w*^{2}〉/*u**^{2} increases with La_{t} for varying wind *U*_{10} and fixed wave age *C _{P}*/

*U*

_{10}(i.e., along the dotted lines). This is

*opposite*of the anticipated trend (i.e., Fig. 1) for these variations in La

_{t}alone. Conversely, 〈

*w*

^{2}〉/

*u**

^{2}decreases with decreasing wave age at fixed

*U*

_{10}, but with a stronger dependence on La

_{t}than found for the LES cases reported by Li et al. (2005), particularly at the lower winds speeds 8 m s

^{−1}<

*U*

_{10}< 15 m s

^{−1}, where 〈

*w*

^{2}〉/

*u**

^{2}∼ La

^{−3}

_{t}. These differences result from the strong sensitivity of 〈

*w*

^{2}〉/

*u**

^{2}to

*D**

_{S}/

*H*

_{ML}through its dependence on both

*U*

_{10}and

*C*/

_{P}*U*

_{10}. These results demonstrate that within the natural parameter space of open ocean storm conditions dominated by wind and wave forcing with 0.6 <

*C*/

_{p}*U*

_{10}< 1.2, typical variations in

*D**

_{S}are as significant to the level of mixed layer VKE as variations in La

_{t}.

### b. Bulk VKE scaling analysis

Figure 7 shows the excess (〈*w*^{2}〉/*u**^{2} − 0.64) of VKE over the La_{t} → ∞ limit of 〈*w*^{2}〉/*u**^{2} ≅ 0.64, scaled by the La^{−2}_{t} dependence obtained (Fig. 1) from the data of Li et al. (2005) and plotted as a function of *D**_{S}/*H*_{ML}. At small *D**_{S}/*H*_{ML} < 0.05, the dependence appears tight and linear, with

followed by a weakening dependence over 0.05 < *D**_{S}/*H*_{ML} < 0.1. The augmentation (〈*w*^{2}〉/*u**^{2} − 0.64) of scaled VKE by Langmuir turbulence appears proportional over *D**_{S}/*H*_{ML} < 0.1 to the La_{t} dependence obtained by fitting data of Li et al. (2005), but with a small offset: our fit [Eq. (6)] to the Li et al. (2005) result corresponds to a level (0.098 at *D**_{S}/*H*_{ML} = 0.12 in Fig. 7) 18% below the equivalent Σ_{2} monochromatic result. This difference is in part due to the inclusion here of a subgrid contribution and to differences in the definition of *H*_{ML}. The increase of 〈*w*^{2}〉/*u**^{2} with *D**_{S}/*H*_{ML} saturates at values above 0.1, but with significantly greater scatter that appears to separate the case subsets by their Stokes drift spectrum bandwidth. This entails a sensitivity to details of the Stokes drift profile, which is not accounted for by its nondimensional *e*-folding scale *D**_{S}/*H*_{ML}.

More generally, Fig. 7 demonstrates that the La_{t} → ∞ “zero-wave” limit and the *D**_{S}/*H*_{ML} → 0 “short-wave limit are equivalent for this bulk nondimensional measure of VKE. As *D**_{S} becomes smaller than approximately 10% of the mixed layer, 〈*w*^{2}〉/*u**^{2} decreases with *D**_{S}/*H*_{ML} toward the zero-wave limit 〈*w*^{2}〉/*u**^{2} ≅ 0.64. For fixed wind and wave forcing in progressively deeper layers, this means that Langmuir turbulence ceases to penetrate a significant fraction of the mixed layer as *D**_{S}/*H*_{ML} → 0, eliminating its impact on the bulk VKE average.

The vertical integration of vortex-force TKE production [Eq. (8)] does not, however, vanish in the *D**_{S}/*H*_{ML} → 0 limit as it does in the wave-free limit La_{t} → ∞. That Eq. (26) implies no extra bulk VKE as *D**_{S}/*H*_{ML} → 0 therefore demonstrates that not all wave-related TKE production is equally important in setting the level of 〈*w*^{2}〉/*u**^{2}. Instead, it is the longest-lived contribution to the mixed layer TKE budget that governs the bulk turbulence level. We infer that it is the CL production associated with the mixed layer–scale eddies that is most important, particularly the eddies generated by the strong surface convergences characteristic of Langmuir circulations. This—and the sensitivity to Stokes bandwidth in the departure from Eq. (26) with increasing *D**_{S}/*H*_{ML} among the different case sets—motivates a different turbulence scaling that is obtained by replacing La_{t} with an SL turbulent Langmuir number

based on a near-surface average 〈*u*^{S}〉_{SL} over the surface layer −*D*_{SL} < *z* < 0. A reference level *u*^{S}_{ref} = *u*^{S}(*z*_{ref}) from depth *z*_{ref} within the lower mixed layer is subtracted because vortex-force production [Eq. (8)], and therefore any augmentation in TKE level, must vanish within a mixed layer with uniform Stokes drift. It was found that plotting VKE intensity versus La_{SL} can collapse the simulation results into a single curve, depending largely on the selection of the near-surface averaging interval *D*_{SL}.

We now find definitions for the averages in La_{SL} and a formula to fit this scaling. By varying the specification of *D*_{SL} and *u*^{S}_{ref}, similarly small errors can be obtained for fitting the model data to either a general power law,

or an exponential function,

A nonlinear fit to Eq. (28) suggests *p* = 1.30 (nearly 4/3) in the wave-dominated small Langmuir number limit. However, this fit extrapolates to an unacceptably low value 〈*w*^{2}〉/*u**^{2} ≅ 0.4 in the zero-wave limit. The exponential function [Eq. (29)] gives *b _{e}* ≅ 0.64, an appropriate wave-free limit, but implies that for a fixed

*D**

_{S}/

*H*

_{ML}and

*u**, 〈

*w*

^{2}〉 approaches a finite value near 4

*u**

^{2}in the La

_{t}→ 0 limit, irrespective of the divergence of vortex-force TKE production [Eq. (8)] as

*U*→ ∞. The relative depths of

_{S}*z*

_{ref}and

*D*

_{SL}were allowed to vary in these nonlinear fits, using initial guesses of −

*H*

_{ML}and

*H*

_{ML}/10, respectively. Minimum errors were consistently found for both functional forms when La

_{SL}was defined by a surface layer average 〈

*u*

^{S}〉

_{SL}of the Stokes drift down to

*D*

_{SL}≅ 0.2

*H*

_{ML}, relative to a lower mixed layer reference drift of

*u*

^{S}

_{ref}=

*u*

^{S}(

*z*

_{ref}) at −0.6

*H*

_{ML}>

*z*

_{ref}> −0.9

*H*

_{ML}.

A *p* = 4/3 power law for small La_{SL} is consistent with a dominant balance in the TKE budget [Eq. (7)] between vortex-force TKE production of *O*(*H*_{ML})-scale eddies and their dissipation:

Combining a power law for small La_{SL} and an exponential decay for large La_{SL} appears to give the best result consistent with both this scaling of TKE production and the zero-wave limit for turbulence intensity. There is also good physical reason for such a shift between two regimes corresponding to Eqs. (28) and (29), because vertical mixing by strong Langmuir circulations largely removes the shear profile responsible for turbulence production in wave-free boundary layers. With some trial and error, a good compromise between simplicity and accuracy was obtained by splitting the fit among continuous functions at La_{SL} = 1, defining the surface layer averaging depth of 〈*u*^{S}〉_{SL} in Eq. (27) by *D*_{SL} = 0.2*H*_{ML}, choosing *p* = 4/3, and setting the La_{SL} → ∞ limit to 〈*w*^{2}〉/*u**^{2} = 0.64. Varying *a _{p}* and

*b*in Eq. (28), and the relative reference depth

_{p}*z*

_{ref}/

*H*

_{ML}for

*u*

^{S}

_{ref}in La

_{SL}Eq. (27), gives

One point from the Σ_{3a} set (La_{t} = 0.435 in Fig. 6) is anomalously high, due to either relatively much larger TKE production in the pycnocline with rapid entrainment or the excitation of a box mode. This highest *U*_{10}, highest *C _{D}*, largest

*H*

_{ML}case was excluded from the fit [Eq. (32)]. Both the low- and high-SL–Langmuir number parameterizations are plotted in Fig. 8a. The fits are not very sensitive to the precise reference level for

*u*

^{S}

_{ref}, but all fits were marginally improved by subtracting a lower mixed layer reference drift.

Equation (32) is the first major result of this study. Given surface wave spectra and drag coefficients, it predicts bulk VKE as a function of wind, waves, and mixed layer depth in the upper ocean boundary layer, assuming Craik–Leibovich dynamics. To demonstrate this result is not simply an artifact of finite-grid scale: Figs. 8b,c scale both the resolved (large eddy) and net subgrid-inclusive bulk VKE from all LES cases on the predictions of Eq. (32) and plot them against nondimensional functions of resolution. Figure 8b plots scaled 〈*w*^{2}〉 against exp(−Δ*z*/*D**_{S}), and Fig. 8c plots it against Δ*z*/*H*_{ML}. In both figures, the extrapolated line from resolved statistics is consistent with the mean (1.00) of the normalized net statistics. This demonstrates that resolution effects on resolved-scale statistics are largely removed by including a subgrid contribution in the net statistics. Results are also shown for three LES cases repeated at uniformly twice the resolution. The trends within the standard resolution cases (due to *D**_{S} and *H*_{ML} variations) are consistent with the convergence toward one with increased resolution in the repeated cases. Together, these tests indicate that limited grid resolution alone does not significantly impact the result in Eq. (32), based on net bulk statistics, beyond an approximately ±5% uncertainty level that is due as well to other error sources, such as domain size and geometry.

### c. The shape of VKE profiles

The profile of (*z*)/*u**^{2} averaged in time and horizontally at each model depth level has a subsurface maximum at a depth *D*_{max(w2)}, which varies with the Stokes depth scale, as noted in Min and Noh (2004) and Li et al. (2005). The location and level of this peak can be sensitive to inclusion of subgrid VKE in (*z*)/*u**^{2}, as illustrated for several Σ_{1} profiles in Fig. 9. This contribution is 2*e*/3 when subgrid TKE is isotropic. Adjacent to the surface, this fraction must be adjusted to reflect changes in the subgrid length scale and anisotropy. The contribution is modified near the surface, where a subgrid anisotropy of / = min(1, −*z*/2*δ*)^{2/3} is assumed. This form of near-surface isotropy is derived from laboratory observations of wall-bounded flow scaling, but the coefficient of *z*/*δ* is chosen to minimize near-surface sensitivity to the grid scale in net VKE profiles. While bulk TKE results are not usually sensitive to this detail, greater circumspection is warranted in the interpretation of near-surface profile shapes from LES models. Maximum net levels of horizontal and VKE components (e.g., Fig. 9, left) may also vary significantly with grid resolution, in part because the downgradient flux approximation in subgrid TKE transport [second from right in Eq. (A6)] is an inaccurate representation of unresolved Langmuir turbulence.

Figure 10a plots /*u**^{2}, the scaled magnitude of the subsurface peak of resolved large-eddy VKE, against La_{t}. (Overbars denote horizontal or temporal mean at a fixed depth, as opposed to the vertically averaged bulk mixed layer statistic.) Comparing Fig. 10a to Fig. 6, /*u**^{2} scales similarly to the bulk VKE. As *D**_{S}/*H*_{ML} increases from *O*(0.1), the scaling of converges on /*u**^{2} ∼ La^{−4/3}_{t} for Σ_{1} cases, consistent with the SL–Langmuir scaling [Eq. (32)]. Low-wind Σ_{1} cases with *D**_{S}/*H*_{ML} < *O*(0.05) exhibit a scaling /*u**^{2} ∼ La^{−7/2}_{t} for constant *U*_{10} that is only marginally stronger than that found in these cases for 〈*w*^{2}〉. Comparing Σ_{1}(1.2) and Σ_{4}(1.2) also shows that is sensitive to *H*_{ML} at larger values of *D**_{S}/*H*_{ML}. This means that remains largely coupled to the bulk TKE budget over our dimensional parameter range. For this reason, we analyze the scaling of the subsurface peak in terms of a shape index /〈*w*^{2}〉 of the VKE profile, complemented by its relative depth *D*_{max(w2)}/*H*_{ML}.

Figures 10b,c show /〈*w*^{2}〉 and *D*_{max(w2)}/*H*_{ML}, evaluated both from the resolved fields and from the net profiles that include a subgrid contribution. Although there is more scatter here than in bulk VKE, the trends are clear: as the Stokes depth scale increases from zero, the depth of the VKE maximum increases until it reaches 0.2*H*_{ML}, the depth of the surface layer. The magnitude of the subsurface peak in correspondingly decreases from 2 − 3 〈*w*^{2}〉 to roughly 1.6 〈*w*^{2}〉. For large *D**_{S}/*H*_{ML}, the shape /〈*w*^{2}〉 of the profile converges on a shape that depends primarily on *z*/*H*_{ML}.

Note that DB spectral-forcing cases (Σ_{1}, Σ_{3}, and Σ_{4}) with different *H*_{ML} and La_{t} are clustered together by their similar values of *D**_{S}/*H*_{ML}. The trends are similar for Σ_{2}, but with significantly larger values for both shape indices. Figure 10c plots a subjective estimate,

of the relationship between *D**_{S} and *D*_{max(w2)} in the broadband-forcing cases, which threads the values from the large eddy and net profiles at low *D**_{S}/*H*_{ML}. A corresponding estimate for /〈*w*^{2}〉 dependence in Fig. 10b is less apparent, but a fit to the relative peak magnitude determined from the large eddy component suggests that

for the DB spectral forcing cases. Net profiles give significantly higher relative magnitudes for *D**_{S}/*H*_{ML} < 0.015, but this apparent divergence is most likely due to the inadequate accuracy of the subgrid parameterization at *D**_{S} < Δ*z* for these shape indices, rather than to a decoupling of the dynamics governing and 〈*w*^{2}〉. A uniformly rescaled version of Eq. (34),

is also included in Fig. 10b as the more likely prediction for /〈*w*^{2}〉 from net TKE.

Last, Fig. 11 shows scaled near-surface crosswind convergence velocity (|_{z=−Δz/2})/*u**^{2}, obtained by averaging the resolved large-eddy statistics from the uppermost model grid level. We do not adjust for a subgrid contribution, because |_{z=−Δz/2} would depend excessively on rough anisotropy assumptions for unresolved TKE. This, and the dependence of measurement level on vertical grid size in a surface-intensified quantity, makes for a capricious turbulence statistic from both modeling and observational perspectives. While the rapid decay of resolved (|_{z=−Δz/2})/*u**^{2} toward the lowest wind speeds of each case set (overstruck x) is at least partly due to these issues of resolution, Fig. 11 nevertheless demonstrates that surface crosswind covariance does not scale with either the bulk or maximum VKE. Several different apparent scaling dependencies may be obtained from limited test cases selected among those shown, including the strong Stokes dependence |_{z=−Δz/2} ∼ La^{−4}_{t}*u**^{2}. These trends at moderate wind speeds within the Σ_{1}(1.2) and Σ_{2}(1.2) subsets suggest that manifestations of an apparent *υ*_{rms} ∼ *U _{S}* scaling in observations of near-surface crosswind convergence may be due to the natural covariance of

*U*and

_{S}*D**

_{S}in mature seas.

### d. Some dimensional predictions

Figures 12a–d translate Eqs. (32), (33), and (35) into dimensional units as a function of wind speed. Figure 12a shows the boundary layer average VKE, scaled by *u**, as a function of wind speed for a 50-m mixed layer depth, for several wave ages and for two different assumptions about the drag coefficient at a high wind speed. Figure 12b shows the results for several mixed layer depths and two extreme wave ages. These results indicate that variations in 〈*w*^{2}〉/*u**^{2} due to wind speed and wave age are similar over the representative range of oceanic values, while the variations due to mixed layer depth are somewhat smaller. Figures 12c,d further show how the shape of the VKE profile, parameterized by indices of its relative magnitude and depth, is predicted to vary as a function of wind speed for several different layer depths, and for both younger and mature waves. The depth clearly increases with increasing wind, saturating at about 0.2*H*_{ML}. This is a clear signature of the increase in Stokes depth with increasing wind speed and thus a signature of surface wave forcing.

A preliminary comparison can be made here between the model predictions of bulk VKE shown in Figs. 12a,b and published measurements. D’Asaro (2001) reports Lagrangian float measurements indicating a uniform level 〈*w*^{2}〉/*u**^{2} = 1.35 ± 0.07 for VKE scaling in 8 to 21.5 m s^{−1} winds. A slightly lower level, 〈*w*^{2}〉/*u**^{2} = 1.14 ± 0.06, was obtained in Tseng and D’Asaro (2004) under a similar range of wind conditions. Bulk VKE statistics obtained at higher wind speeds in Hurricane Dennis (D’Asaro 2003) and Hurricane Frances (unpublished, but see D’Asaro and McNeil 2006) appear to be consistent with the lower wind speed results, within substantially larger experimental uncertainties. These data are consistent with the modeling results presented here if wave age decreases from near 1 at 8 m s^{−1} to about 0.5 at 40 m s^{−1}. This, however, is the expected trend: mature seas are typically encountered in light to moderate winds, while the youngest seas encountered in the open ocean (0.4 < *C _{p}*/

*U*

_{10}< 0.6) may be found in rapidly moving hurricanes. Such an inverse relationship of

*U*

_{10}and

*C*/

_{p}*U*

_{10}follows from empirical growth rates for fetch and duration-limited waves and the synoptic length and time scales of storm forcing; seas take longer to mature at high winds, thus the wave field tends to be younger at high winds. Similarly, reported subsurface peak-VKE values, typically ∼2

*u**

^{2}occurring at a depth less than

*H*

_{ML}/4, are not inconsistent with our predictions. More careful and detailed model–data intercomparisons require further analysis of the sea states concurrent with float measurements, which are the subjects of a later paper.

## 6. Summary and conclusions

In his recent review of Langmuir circulations, Thorpe (2004) states, “It is now widely accepted that what is usually identified as [Langmuir circulations] at scales of 2 m–1 km on sea or lake surfaces . . . generally arises through the interaction of the Stokes drift induced by surface waves and the vertical shear in turbulent fluid, consequently producing a vortex force as described by the so-called CL2 model of Craik and Leibovich (1976).” However, little quantitative evidence supports this claim. There is little doubt that structures similar to the canonically defined Langmuir circulations exist and that the equations of motions with the addition of the CL vortex force can create similar structures. It is less clear to what extent these structures as well as the overall boundary layer turbulence statistics, fluxes, and intensities are accurately predicted by these equations. This is because both appropriate oceanic measurements and theoretical analyses with sufficient realism to be quantitatively compared to observations have been lacking. The major motivation for the work presented here is to provide more realistic theoretical predictions of turbulent kinetic energy levels in the ocean boundary layer for comparison with recent observations. The major results are 1) scaling Eq. (32) for bulk-averaged vertical kinetic energy in the ocean boundary layer forced by wind stress and the CL vortex force and 2) predictions for the shape of the vertical kinetic profile [Eqs. (33) and (35)].

A total of 58 large-eddy simulation (LES) cases of the ocean boundary layer were conducted across a range of forcing that spans realistic oceanic storm conditions. The simulations were forced by wind stress *τ*_{0} = *ρ _{W}*

*u**

^{2}

**Û**

_{10}and a profile of surface wave Stokes drift

**u**

^{S}(

*z*) interacting with the velocity field through the Craik–Leibovich vortex force. No other wave-forcing effects were included. An empirical surface wave spectral form based on Donelan et al. (1985, 1993 and Banner (1990) was used to generate consistent stress and Stokes drift forcing for a given wind speed

*U*

_{10}and wave age

*C*/

_{P}*U*

_{10}under the assumption of pure wind seas. Wind speeds (8–70 m s

^{−1}), wave ages (0.6–1.2), and mixed layer depths (30–82 m s

^{−1}) span the range of midlatitude and tropical storm-forced conditions. Results include:

For the simulated range of surface forcing, the mixed layer-averaged, scaled VKE 〈

*w*^{2}〉/*u**^{2}reaches levels as high as 1.8, several times its La_{t}→ ∞ wave-free value (≅0.64), under the action of the CL vortex force.A conventional assumption—that surface forcing may be represented by a scale-equivalent monochromatic wave—is shown to not always yield the same scaling results as empirically based, spectrally distributed surface wave forcing.

The value of 〈

*w*^{2}〉/*u**^{2}in LES results can be predicted accurately from a surface layer Langmuir number La_{SL}[Eq. (27)] as shown in Fig. 8a and parameterized by Eq. (32). This applies to waves that are either monochromatic or spectrally distributed.Variations in 〈

*w*^{2}〉/*u**^{2}due to wind speed and wave age are of roughly equal magnitude, suggesting that wave properties cannot, in general, be ignored in boundary layer models sensitive to turbulence levels.The profile of VKE (

*z*)/*u**^{2}has a subsurface maximum. The depth of this maximum is controlled by the*e*-folding depth of the Stokes drift*D**_{S}when*D**_{S}/*H*_{ML}is small and by*H*_{ML}when it is large, as shown in Fig. 10c, and partially parameterized by Eq. (33).The scaling of the subsurface TKE maximum is similar to that of 〈

*w*^{2}〉. Deviations that depend on*D**_{S}/*H*_{ML}are shown in Fig. 10b and partially parameterized by Eq. (35).The crosswind convergence velocity variance scales differently from the VKE (Fig. 11).

Model–data comparisons show that Lagrangian float observations of a nearly uniform value of 〈*w*^{2}〉/*u**^{2} can be explained by the combination of an increase of this ratio with wind speed, the decrease of this ratio with wave age, and the natural tendency for wave age to decrease with wind speed.

## Acknowledgments

This work was supported by the National Science Foundation (OCE0117411), the Office of Naval Research (N00014-00-1-0893), and by a grant of HPC resources from the Arctic Region Supercomputing Center at the University of Alaska Fairbanks as part of the Department of Defense High Performance Computing Modernization Program. Two reviewers at this journal, and two from its prior submission to *Ocean Modelling*, provided very useful critiques of the presentation.

## REFERENCES

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**.**

**,**

**.**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

### APPENDIX A

#### LES Model Description

The wave-averaged Lagrangian velocity is *u*^{L}_{i} = *u*^{E}_{i} + *u*^{S}_{i}, where *i* ∈ [1, 2, 3] indicate [**x̂**_{i}] = [**x̂**, **ŷ**, **ẑ**], with **u**^{L} = [*u*^{L}_{i}] and Eulerian **u**^{E} = [*u*^{E}_{i}]. Equivalent indexing *i* ∈ [*x*, *y*, *z*, *θ*, *S*, *b*] extends to potential temperature, salinity and buoyancy. The Lagrangian velocity is taken to be the prognostic variable, and its superscript *u _{i}* = [

*u*,

*υ*,

*w*] ≡

*u*

^{L}

_{i}is omitted. The evolution of mean Lagrangian momentum uses

to solve prognostically for **u** in terms of ** ζ** =

**∇**×

**u**and the Stokes drift curl

*ζ*^{S}=

**∇**×

**u**

*,*

^{S}equivalent to standard ∂_{t}**u**^{E} formulations with CL vortex force **u*** ^{S}* ×

*ζ**. Generalized pressure*

^{E}is solved diagnostically using **∇** · **u** = 0. Subgrid stress contributions Γ_{i} = ∂_{j}*τ _{ij}* ≡ ∂

*τ*/∂

_{ij}*x*use an advected subgrid TKE equation. Potential temperature is computed prognostically from

_{j}with corresponding equations for salinity. Unresolved subgrid stresses are modeled in terms of the resolved symmetric strain rate of the Eulerian velocity and local nonlinear eddy viscosity *K _{M}*, as well as thermohaline scalar fluxes in terms of resolved gradients and eddy diffusivity

*K*:

_{H}Shear production in the subgrid TKE equation,

with ɛ = *c*_{ɛ}*e*^{3/2}/*l _{d}*, combines (Eulerian) stress with Lagrangian shear: −

*τ*∂

_{ij}_{j}

*u*= −

_{i}*τ*(∂

_{ij}_{j}

*u*

^{E}

_{i}+ ∂

_{j}

*u*

^{S}

_{i}).

Length scale *l _{d}* is model resolution scale

*δ*, and

*c*

_{ɛ}and

*c*are constants, except when stratification or proximity to a boundary reduces the turbulence length scale below

_{K}*δ*. Model resolution

*δ*is about 1.5 times Δ = (Δ

*x*Δ

*y*Δ

*z*)

^{1/3}. The vertical grid scale is not refined toward the surface, but the derivative ∂

_{z}

*u*

^{E}

_{x,y}at

*z*= −Δ

*z*is adjusted from Δ

*u*

^{E}

_{x,y}/Δ

*z*, assuming Δ

*u*

^{E}

_{x,y}is a finite difference between grid volume averages of a logarithmic function of

*z*.

### APPENDIX B

#### Implementation of the Vortex-Force Term

As LES equations predict volume-filtered variables, Stokes drift is applied in a filtered form,

which is subtracted from filtered momentum 〈**u**〉_{Δ} to compute vorticity in Eq. (A2): 〈**u** × (** ζ** −

*ζ*^{S})〉

_{Δ}≅ 〈

**u**〉

_{Δ}× [

**∇**

_{Δ}× (〈

**u**〉

_{Δ}− 〈

**u**

^{S}〉

_{Δ})], where “

**∇**

_{Δ}× ” is the discrete curl. Equivalently,

**u**

^{S}×

*ζ**would become 〈*

^{E}**u**

*〉*

^{S}_{Δ}×

*ζ*^{E}. Grid-filtered Stokes drift,

is computed from a discrete Stokes spectrum *δ***U**_{S}/*δ**ω*|_{ωi} of *N _{ω}* elements, plus a “tail” contribution 〈

**u**

^{S}|

_{ω>ωcut}〉

_{Δ}. Assuming Φ(

*ω*) ∼

*ω*

^{−5}above,

*ω*

_{cut}=

*ω*+ Δ

_{Nω}*ω*/2,

_{Nω}where *a*_{±} = −2(*z* ± Δ*z*/2)*k*_{cut} and the integral

are evaluated using standard numerical implementations of the complementary error function erfc.

## Footnotes

*Corresponding author address:* R. R. Harcourt, Applied Physics Laboratory, University of Washington, 1013 NE 40th St., Seattle, WA 98105. Email: harcourt@apl.washington.edu

^{1}

The original PM reference height 19.5 m may be inferred for *U _{W}* from Eq. (15) in Li and Garrett (1993).