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 uS(z) are specified from the 10-m wind speed U10 and the surface wave age CP/U10, where CP is the spectral peak phase speed, using an empirical surface wave spectra and an associated wave age–dependent neutral drag coefficient CD. Wave-breaking effects are not otherwise included. Mixed layer depths HML vary from 30 to 120 m, with 0.6 ≤ CP/U10 ≤ 1.2 and 8 m s−1 < U10 < 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 〈w2〉/u*2 is equally sensitive to the nondimensional Stokes e-folding depth D*S/HML and to the turbulent Langmuir number Lat = u*/US, where u* = |τ0|/ρw in water density ρw and US = |uS|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 LaSL = u*/〈uS〉SL, where 〈uS〉SL is the average Stokes drift in a surface layer 0 > z > − 0.2HML relative to that near the bottom of the mixed layer. In the wave-dominated limit LaSL → 0, turbulent vertical velocity scales as wrms ∼ u*La−2/3SL. The mean profile (z) of VKE is characterized by a subsurface peak, the depth of which increases with D*S/HML to a maximum near 0.22HML as its relative magnitude /〈w2〉 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 HML ≅ 50 m, 〈w2〉/u*2 varies from less than 1 for young waves at U10 = 10 m s−1 to about 2 for mature seas at winds greater than U10 = 30 m s−1. Preliminary comparisons with Lagrangian float data account for invariance in 〈w2〉/u*2 measurements as resulting from an inverse relationship between U10 and CP/U10 in observed regimes.
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 uS of surface waves and wave-averaged currents driven by a surface stress τ0 = |τ0| = ρwu*2, where u* is the friction velocity in water of density ρw (see Thorpe 2004 for a review). The CL interaction is represented in the momentum equations by a “vortex force,”
perpendicular to both the Stokes drift uS and the vorticity ζE = ∇ × uE of the Eulerian mean flow uE = [uE, υE, wE]. While the Eulerian mean is obtained by averaging fluid velocity over the wave phase in a fixed reference frame, the Lagrangian mean velocity uL ≅ uE + uS is the average obtained in a wave-following Lagrangian reference frame.
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 Lat = u*/US, where u* = |τ0|/ρw and US = |uS|z=0. For wind- and wave-dominated regimes and typical open ocean values of Lat, 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 uSi = uLi − uEi approximate pi when the waves are irrotational. In this study we further approximate the Stokes drift to 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(A2), or 2ωkη2rmse2kz 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/2k* (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(A2), such broadband wave spectra have Stokes drift
composed of a Stokes drift spectrum with elements δUS/δω, 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 Lat = 0.3 and observational evidence that suggested wrms ∼ US for Lat ≪ 1. This proposed dependence of vertical velocity on Lat implies here that
although McWilliams and Sullivan note as well that the dependence of w on US 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 Lat = [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/3U1/3S (m = 3) than they are with υrms|z≅0 ∝ u*1/2US1/2 (m = 2) or υrms|z≅0 ∝ US (m = 1). However, these results are somewhat ambiguous because they are based on a set of simulation cases with two different D*S/HML values. Moreover, the proportional scaling form Eq. (5) is relevant to Lat ≪ 1, but may not be the best comparison for their larger (Lat = 0.45) Langmuir number case.
Li et al. (2005) simulate 13 model cases at D*S/HML = 0.12 with H0 = 0, and graphically present (their Fig. 4) LES results for the equilibrium level of bulk VKE 〈w2〉, obtained by averaging w2 first over the horizontal domain and the mixed layer depth z > −HML(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 Lat → ∞ limit as Eq. (6). They do not agree at low Lat (i.e., when the CL force is important).
Several LES studies have examined the sensitivity of TKE to the nondimensional ratio D*S/HML, where HML 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/HML 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/HML = [0.24, 0.12, 0.06] at a fixed Langmuir number Lat = 0.34. They report little difference in the profile of VKE components between the two largest values of D*S/HML but a significant reduction in for D*S/HML = 0.06. While Li et al. suggest a sensitivity at small D*S/HML, the concurrent impact on the horizontal components of TKE from variations in D*S/HML appears uniformly minor. Min and Noh (2004) vary D*S/HML = 0.127 and 0.064 at Lat = 0.35 (their cases “A” and “B”), D*S/HML = 0.127 at Lat = 0.25 (“E”), and D*S/HML = 0.064 at Lat = 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/HML 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 ∼ US 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/2U1/2S. There appears to be no consensus here.
D’Asaro (2001) and Tseng and D’Asaro (2004) report 〈w2rms〉 ∝ 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 PShear = −∂/∂xj, buoyant production PBuoy, and the divergence DTran. 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 uS(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 GN(ω, θ). 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 GN, requires truncation at the capillary limit for the Stokes drift integral [Eq. (3)] to converge at z = 0. Attention to this behavior of the spectral tail is therefore critical to quantitative model–data comparisons that depend on US.
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 uS(z) with the same e-folding depth scale D*S = 1/2k*, based on Pierson–Moskowitz (PM) spectrum for fully developed seas, where the phase speed Cp = ωp/kp of spectral peak waves at kp = ω2p/g reaches Cp = 1.2U10. For a range of empirical parameters, Li and Garrett estimate the Stokes wind ratio as
and the e-folding Stokes depth scale as
For US/U10 = 1.55%, Lat ranges from 0.28 at U10 = 10 m s−1 to 0.33 at U10 = 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 Lat 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 < Lat < 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 < U10 < 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 uS(z) as a function of varying wind speed U10 and wave age Cp/U10. This empirical spectrum, based on measurements at ω < 4ωp, combines
with a nearly normalized spreading function
where the shape parameters α, γ, and σ depend on ϖ = ω/ωp and on wave age Cp/U10. 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 sech2, 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 η2rms = ∫∞0 Φ′(ω) dω for this spectrum is not significantly changed from that in Donelan et al. (1985), η2rmsg2/U410 = (6.365 × 10−3)[U10/(2πCp)]−3.3: for wave ages Cp/U10 = [0.6, 0.8, 1.0, 1.2], finescale discrete integration gives nondimensional energy levels η2rmsg2/U410 = [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 Cp/U10 = 1.17 and Cp/U10 = 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 Cp/U10 = 1.17, where JONSWAP and PM spectra are nearly the same.
Stokes drift profiles are found by substituting this Donelan–Banner spectrum F′(ϖ, θ) = Φ′(ϖ)GN(ϖ, θ) into Eq. (3). Scaling the DB Stokes drift on U10 gives
The integral is evaluated by discrete approximation over 0 < ϖi < ϖcut as uS ≅ Σ(δUS/δϖ)ϖie2kizΔϖ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 US/U10 evaluated from Eq. (21) at z = 0. It depends only on Cp/U10. 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 2D*Skp versus Cp/U10, scaling the Stokes depth on the (significantly larger) Stokes depth (2kp)−1 of the dominant wavenumber kp = g/C2p. For typical ocean wave ages Cp/U10 > 0.6, e-folding depths for the DB spectrum are 26%–37% larger than the 2D*Skp = 0.19 value obtained by Li and Garrett (1993) for the fully developed PM spectrum.
The above DB model has a slightly larger surface Stokes wind ratio US/U10 = 1.75% at Cp/U10 = 1.2 (Fig. 3c) than the value US/U10 = (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 U10 = 10 m s−1, Smith obtains US/U10 = 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 US/U10 = 1.53%. This corresponds to less than fully developed waves, Cp/U10 = 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 k2F(ω, θ) 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 ξI(100ωp) are still small. As a function of wave age, ξI(100ωp) ranges from 0.055 for Cp/U10 = 1.2 to 0.066 for Cp/U10 = 0.6 in the DB-based simulations below, while accounting for over 99% of US. 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.
A wind stress consistent with the wave-age dependence of the DB model Stokes drift is the neutral drag coefficient CD of Donelan et al. (1993), which is a function of Cp/U10. Surface stress τ0 = ρaCDU210 in air of density ρa is evaluated from surface roughness length z0 using
Figure 3e plots this age-dependent drag coefficient versus U10, 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 CD on Cp/U10. 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 U10 > 30 m s−1 where evidence suggests that the increase of CD with wind speed “saturates” at a maximum value.
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 uS is replaced by a grid-filtered quantity 〈uS〉Δ, 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 uS(−Δz/2) and 〈uS〉Δ at z = −Δz/2, but for the DB spectral model, at the vertical resolution (Δz = 1.42) of most model cases below, using unfiltered uS(−Δz/2) would underestimate 〈uS(−Δz/2)〉Δ by between 1% (for U10 = 32.6 m s−1 and Cp/U10 = 1.2) and 22% (for U10 = 8.3 m s−1 and Cp/U10 = 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 U10 and Cp/U10 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 CD reduce the forcing parameter space by determining Lat, D*S, US/U10, and uS/US as functions of U10 and Cp/U10. The nondimensional ratio D*S/HML depends on an average HML 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 U10 values versus four values of wave age. Surface stress and Stokes drift profiles were determined as functions of U10 and Cp/U10 using the spectral-forcing paradigm described in section 2b. Wave ages ranged over typical oceanic values of 0.6 < Cp/U10 < 1.2 for storm forcing, with a wind speed range of 8 m s−1 < U10 < 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 U10 values from Σ1 as wind speed UW 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 US/UW. Surface stress was determined using the neutral drag coefficient of Large and Pond (1981), modified to extend the linear drag relation to the lowest (U10 < 10 m s−1) values. The resulting dependence of D*S(U10) versus Lat(U10) 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 Lat 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, U10 > 20 m s−1. Set Σ3a continues the Cp/U10 = 0.6 subset of Σ1 over three hurricane-strength wind speeds, 33 m s−1 < U10 < 70 m s−1, while Σ3b uses a constant CD = 2.3 × 10−3, slightly below high-wind saturation values of (2.5 − 2.8) × 10−3 suggested by Donelan et al. (2004). Solutions at Cp/U10 = 0.6 from Σ1 are contiguous with Σ3a; Σ3b connects with Σ1 near U10 = 20 m s−1, where CD = 2.3 × 10−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 Cp/U10 = 1.2 subset of Σ1 for mature seas, but with the initial mixed layer depth HML 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−1p = C2pg−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(Cp/U10 = 1.2) distinguishes D*S/HML dependence from Δz/HML 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 U10 = 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 U10 = 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 U10 = 32.6 m s−1, Cp/U10 = 0.6 case in Σ1. Initial conditions for Σ4 were generated by vertically compressing the final states from the Cp/U10 = 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, HML, indicated for each example profile. Here, HML 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, HML 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 〈w2〉/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 w2 horizontally over the LES domain. A time series of the mixed layer average 〈w2〉 is obtained by averaging these profiles over the mixed layer 0 > z > −HML, every 50 time steps. The final steady-state average of 〈w2〉 is obtained over a terminal time interval Tavg ≫ HVKE/〈w2〉, 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/HML ≅ 0.12.
Within each case set, 〈w2〉/u*2 increases with Lat for varying wind U10 and fixed wave age CP/U10 (i.e., along the dotted lines). This is opposite of the anticipated trend (i.e., Fig. 1) for these variations in Lat alone. Conversely, 〈w2〉/u*2 decreases with decreasing wave age at fixed U10, but with a stronger dependence on Lat than found for the LES cases reported by Li et al. (2005), particularly at the lower winds speeds 8 m s−1 < U10 < 15 m s−1, where 〈w2〉/u*2 ∼ La−3t. These differences result from the strong sensitivity of 〈w2〉/u*2 to D*S/HML through its dependence on both U10 and CP/U10. These results demonstrate that within the natural parameter space of open ocean storm conditions dominated by wind and wave forcing with 0.6 < Cp/U10 < 1.2, typical variations in D*S are as significant to the level of mixed layer VKE as variations in Lat.
b. Bulk VKE scaling analysis
Figure 7 shows the excess (〈w2〉/u*2 − 0.64) of VKE over the Lat → ∞ limit of 〈w2〉/u*2 ≅ 0.64, scaled by the La−2t dependence obtained (Fig. 1) from the data of Li et al. (2005) and plotted as a function of D*S/HML. At small D*S/HML < 0.05, the dependence appears tight and linear, with
followed by a weakening dependence over 0.05 < D*S/HML < 0.1. The augmentation (〈w2〉/u*2 − 0.64) of scaled VKE by Langmuir turbulence appears proportional over D*S/HML < 0.1 to the Lat 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/HML = 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 HML. The increase of 〈w2〉/u*2 with D*S/HML 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/HML.
More generally, Fig. 7 demonstrates that the Lat → ∞ “zero-wave” limit and the D*S/HML → 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, 〈w2〉/u*2 decreases with D*S/HML toward the zero-wave limit 〈w2〉/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/HML → 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/HML → 0 limit as it does in the wave-free limit Lat → ∞. That Eq. (26) implies no extra bulk VKE as D*S/HML → 0 therefore demonstrates that not all wave-related TKE production is equally important in setting the level of 〈w2〉/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/HML among the different case sets—motivates a different turbulence scaling that is obtained by replacing Lat with an SL turbulent Langmuir number
based on a near-surface average 〈uS〉SL over the surface layer −DSL < z < 0. A reference level uSref = uS(zref) from depth zref 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 LaSL can collapse the simulation results into a single curve, depending largely on the selection of the near-surface averaging interval DSL.
We now find definitions for the averages in LaSL and a formula to fit this scaling. By varying the specification of DSL and uSref, 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 〈w2〉/u*2 ≅ 0.4 in the zero-wave limit. The exponential function [Eq. (29)] gives be ≅ 0.64, an appropriate wave-free limit, but implies that for a fixed D*S/HML and u*, 〈w2〉 approaches a finite value near 4u*2 in the Lat → 0 limit, irrespective of the divergence of vortex-force TKE production [Eq. (8)] as US → ∞. The relative depths of zref and DSL were allowed to vary in these nonlinear fits, using initial guesses of −HML and HML/10, respectively. Minimum errors were consistently found for both functional forms when LaSL was defined by a surface layer average 〈uS〉SL of the Stokes drift down to DSL ≅ 0.2HML, relative to a lower mixed layer reference drift of uSref = uS(zref) at −0.6HML > zref > −0.9HML.
A p = 4/3 power law for small LaSL is consistent with a dominant balance in the TKE budget [Eq. (7)] between vortex-force TKE production of O(HML)-scale eddies and their dissipation:
Combining a power law for small LaSL and an exponential decay for large LaSL 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 LaSL = 1, defining the surface layer averaging depth of 〈uS〉SL in Eq. (27) by DSL = 0.2HML, choosing p = 4/3, and setting the LaSL → ∞ limit to 〈w2〉/u*2 = 0.64. Varying ap and bp in Eq. (28), and the relative reference depth zref/HML for uSref in LaSL Eq. (27), gives
One point from the Σ3a set (Lat = 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 U10, highest CD, largest HML 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 uSref, 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 〈w2〉 against exp(−Δz/D*S), and Fig. 8c plots it against Δz/HML. 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 HML 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 Dmax(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 2e/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 Lat. (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/HML increases from O(0.1), the scaling of converges on /u*2 ∼ La−4/3t for Σ1 cases, consistent with the SL–Langmuir scaling [Eq. (32)]. Low-wind Σ1 cases with D*S/HML < O(0.05) exhibit a scaling /u*2 ∼ La−7/2t for constant U10 that is only marginally stronger than that found in these cases for 〈w2〉. Comparing Σ1(1.2) and Σ4(1.2) also shows that is sensitive to HML at larger values of D*S/HML. 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 /〈w2〉 of the VKE profile, complemented by its relative depth Dmax(w2)/HML.
Figures 10b,c show /〈w2〉 and Dmax(w2)/HML, 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.2HML, the depth of the surface layer. The magnitude of the subsurface peak in correspondingly decreases from 2 − 3 〈w2〉 to roughly 1.6 〈w2〉. For large D*S/HML, the shape /〈w2〉 of the profile converges on a shape that depends primarily on z/HML.
Note that DB spectral-forcing cases (Σ1, Σ3, and Σ4) with different HML and Lat are clustered together by their similar values of D*S/HML. 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 Dmax(w2) in the broadband-forcing cases, which threads the values from the large eddy and net profiles at low D*S/HML. A corresponding estimate for /〈w2〉 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/HML < 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 〈w2〉. A uniformly rescaled version of Eq. (34),
is also included in Fig. 10b as the more likely prediction for /〈w2〉 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−4tu*2. These trends at moderate wind speeds within the Σ1(1.2) and Σ2(1.2) subsets suggest that manifestations of an apparent υrms ∼ US scaling in observations of near-surface crosswind convergence may be due to the natural covariance of US and 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 〈w2〉/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.2HML. 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 〈w2〉/u*2 = 1.35 ± 0.07 for VKE scaling in 8 to 21.5 m s−1 winds. A slightly lower level, 〈w2〉/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 < Cp/U10 < 0.6) may be found in rapidly moving hurricanes. Such an inverse relationship of U10 and Cp/U10 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 ∼2u*2 occurring at a depth less than HML/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 = ρWu*2Û10 and a profile of surface wave Stokes drift uS(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 U10 and wave age CP/U10 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 〈w2〉/u*2 reaches levels as high as 1.8, several times its Lat → ∞ 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 〈w2〉/u*2 in LES results can be predicted accurately from a surface layer Langmuir number LaSL [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 〈w2〉/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/HML is small and by HML when it is large, as shown in Fig. 10c, and partially parameterized by Eq. (33).
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 〈w2〉/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.
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.
LES Model Description
The wave-averaged Lagrangian velocity is uLi = uEi + uSi, where i ∈ [1, 2, 3] indicate [x̂i] = [x̂, ŷ, ẑ], with uL = [uLi] and Eulerian uE = [uEi]. 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 ui = [u, υ, w] ≡ uLi is omitted. The evolution of mean Lagrangian momentum uses
to solve prognostically for u in terms of ζ = ∇ × u and the Stokes drift curl ζS = ∇ × uS,
equivalent to standard ∂tuE formulations with CL vortex force uS × ζE. Generalized pressure
is solved diagnostically using ∇ · u = 0. Subgrid stress contributions Γi = ∂jτij ≡ ∂τij/∂xj use an advected subgrid TKE equation. Potential temperature is computed prognostically from
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 KM, as well as thermohaline scalar fluxes in terms of resolved gradients and eddy diffusivity KH:
Shear production in the subgrid TKE equation,
with ɛ = cɛe3/2/ld, combines (Eulerian) stress with Lagrangian shear: −τij∂jui = −τij(∂juEi + ∂juSi).
Length scale ld is model resolution scale δ, and cɛ and cK are constants, except when stratification or proximity to a boundary reduces the turbulence length scale below δ. 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 ∂zuEx,y at z = −Δz is adjusted from ΔuEx,y/Δz, assuming ΔuEx,y is a finite difference between grid volume averages of a logarithmic function of z.
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〉Δ − 〈uS〉Δ)], where “∇Δ × ” is the discrete curl. Equivalently, uS × ζE would become 〈uS〉Δ × ζE. Grid-filtered Stokes drift,
is computed from a discrete Stokes spectrum δUS/δω|ωi of Nω elements, plus a “tail” contribution 〈uS|ω>ωcut〉Δ. Assuming Φ(ω) ∼ ω−5 above, ωcut = ωNω + ΔωNω/2,
where a± = −2(z ± Δz/2)kcut and the integral
are evaluated using standard numerical implementations of the complementary error function erfc.
Corresponding author address: R. R. Harcourt, Applied Physics Laboratory, University of Washington, 1013 NE 40th St., Seattle, WA 98105. Email: email@example.com
The original PM reference height 19.5 m may be inferred for UW from Eq. (15) in Li and Garrett (1993).