## Abstract

The Reynolds stress equation is modified to include the Craik–Leibovich vortex force, arising from the interaction of the phase-averaged surface wave Stokes drift with upper-ocean turbulence. An algebraic second-moment closure of the Reynolds stress equation yields an algebraic Reynolds stress model (ARSM) that requires a component of the vertical momentum flux to be directed down the gradient of the Stokes drift, in addition to the conventional component down the gradient of the ensemble-averaged Eulerian velocity. For vertical and horizontal component fluctuations, the momentum flux must be closed using the form , where the coefficient is generally distinct from the eddy viscosity or eddy diffusivity . Rational expressions for the stability functions , , and are derived for use in second-moment closure models where the turbulent velocity and length scales are dynamically modeled by prognostic equations for and . The resulting second-moment closure (SMC) includes the significant effects of the vortex force in the stability functions, in addition to source terms contributing to the and equations. Additional changes are made to the way in which is limited by proximity to boundaries or by stratification. The new SMC model is tuned to, and compared with, a suite of steady-state large-eddy simulation (LES) solutions representing a wide range of oceanic wind and wave forcing conditions. Comparisons with LES show the modified SMC captures important processes of Langmuir turbulence, but not without notable defects that may limit model generality.

## 1. Introduction

Upper-ocean mixing models without explicit representations of surface waves may implicitly represent their impact when tuned to oceanic observations because of the natural correlation between wind and wave forcing. However, such models may be inaccurate if dimensional scales of surface waves do not scale simply with the wind, as is the case for variations in sea state or wave age at a given wind speed, or for variations in the relative strength of wave versus wind effects with the upper ocean mixed layer depth.

Mixed layer models have primarily sought to explicitly articulate surface wave effects in two generally distinct ways. One approach, after Craig and Banner (1994), accounts for the loss of energy from waves into the phase-averaged turbulent velocity fluctuations relative to an ensemble mean as a source term at the surface for turbulent kinetic energy (TKE) . This is added as a surface flux or boundary condition in boundary layer turbulence closures such as “*k*–*ɛ*” models, where a second equation may predict dissipation *ɛ*, or in analogous “*k–kl*” models (a.k.a., q2–q2l) such as Mellor and Yamada (1982) and Kantha and Clayson (1994, hereafter KC94). Enhanced near-surface TKE and associated impacts on the turbulence length scale serve to increase the vertical eddy viscosity and diffusivity , where stability functions , are determined from a second-moment algebraic closure of Reynolds stress and flux equations for an equilibrium state.

Another way waves can alter mixed layer models is by accounting for the dynamical effects of surface wave Stokes drift , which appears in the equation for momentum through the Craik–Leibovich (CL) vortex force of Craik and Leibovich (1976) and through an associated Bernoulli pressure term. These dynamical effects can be included through modifications of second-moment turbulence equations or through other closure assumptions. The modified second-moment closures of D’Alessio et al. (1998) and Kantha and Clayson (2004, hereafter KC04), incorporate the corresponding CL vortex production TKE source term of TKE production into the prognostic equation for *k* or *q*^{2} and infer an analogous modification in the prognostic equation for *ɛ* or *q*^{2}*l*.

Including the dynamics associated with surface waves in upper-ocean turbulence models appears warranted. Turbulence observations in a wide variety of ocean regimes find that vertical TKE below wave-bounded mixed layers is significantly elevated above levels in wall-bounded layers (D’Asaro 2001). The CL vortex force appears most likely responsible for this broadly observed elevation in TKE: large-eddy simulations (LESs) of mixed layer turbulence that include the CL vortex force after McWilliams et al. (1997), and carried out for model forcing that covers a natural range of wind seas, predicts profiles that are consistent with these elevated levels observed by Lagrangian floats (Harcourt and D’Asaro 2008, hereafter HD08). This circumstance motivates the further development here of a turbulent mixing model that is focused on the effects of the CL vortex force on the largest mixed layer eddies of Langmuir turbulence. The omission of wave breaking from this picture assumes this larger surface energy flux is injected to small scales and dissipates rapidly through a forward turbulent cascade near the surface, without strongly impacting the dynamics of larger Langmuir structures at the scale of the mixed layer depth. Stokes–breaker interactions modeled in Sullivan et al. (2004a,b) stand as a notable hypothesis counter to this assumption.

The model formulation presented here was initially motivated by the observation that the stability functions and used in KC04 were derived in KC94 from an algebraic Reynolds stress model (ARSM) that includes only the local forcing effects of stratification and shear, but not the CL vortex force. If the additional vortex force TKE production is included in the prognostic model equations predicting and , it is inconsistent to omit them from the ARSM predicting and . The remedy requires a fundamentally different closure assumption, namely that a component of momentum flux is proportional to the Stokes shear .

## 2. Second-moment closure with Craik–Leibovich vortex forcing

### a. Reynolds equations

The CL vortex force due to surface wave Stokes drift is incorporated after McWilliams et al. (1997) into the Navier–Stokes equations for wave-phase-averaged Eulerian velocity as

where and are gravitational and Coriolis components, and is viscosity, with incompressible flow . Nonhydrostatic pressure is scaled to on reference density and modified by Stokes drift Bernoulli terms to . Eulerian advection of a thermodynamically active scalar with expansion coefficient and diffusivity is modified by the Stokes advection

Ensemble fluctuations , , , and of buoyancy are defined relative to ensemble averages , , and . Fluctuations in Stokes drift are assumed here to have no significant coherence with turbulent fluctuations or . Mean and are governed by

and

In Eqs. (3) and (4), the Lagrangian derivative is redefined as following [from following in Eqs. (1) and (2)], and incompressibility is now applied separately to both ensemble mean and fluctuating velocities.

Following the development of the equilibrium model in KC94, we begin by modifying the second-moment equations for the impact of the CL vortex force. Using , the equations for Reynolds stress and flux may be obtained from Eqs. (3) and (4) as

and

Note in Eq. (5) that the form of the Reynolds stress shear production (second term on right, first in brackets) and that of the new contribution from CL vortex production (last right term in brackets) differ subtly in the arrangement of indices, and only become trivially equivalent in the trace of these two terms contributing to TKE production.

### b. Closure

Standard closure assumptions invoked in KC94 are used here with two minor changes, generalizing the deformation of turbulence by shear to formally include Stokes drift effects. The two generalizations are made ad hoc, and introduce two new model constants, and , in addition to the others that reappear here following the notation in KC94: *A*_{1}, *A*_{2}, *B*_{1}, *B*_{2}, *C*_{1}, *C*_{2}, *C*_{3}, and *S _{q}*.

Assumptions unchanged from KC94 include those for dissipation and the dissipation (or master) turbulence length scale :

Pressure–strain rate correlations are modeled in KC94 to include a return-to-isotropy term after Rotta (1951) and a term for distortion by shear after Crow (1968). Here, the latter is generalized to include a component due to Stokes shear subject to a separate constant :

The pressure–scalar correlations modeled in KC94 after Moeng and Wyngaard (1986) are extended by the same similarity assumption, subject to the second additional closure constant :

While a wide variety of closure assumptions are made in about as many second-moment closure (SMC) models, the treatment considered here is constrained to the form of the closures used in KC94 and KC04. Closure assumptions for the third moment and for the pressure–velocity transport terms are not altered as follows:

In the hierarchy of Mellor and Yamada (1974), closure is obtained at different “levels” by retaining terms to a given order in an expansion of Eq. (5) in terms of a small nondimensional parameter controlling the departure of turbulence from isotropy. At a level qualitatively referred to as “2½,” Eq. (5) is replaced by a prognostic energy equation

for the trace of the Reynolds stress tensor, and the individual components are determined by algebraic expressions,

The closure of KC94 follows the expansion procedure of Galperin et al. (1988) in using the equilibrium solution of Eq. (12) to substitute for terms in Eq. (5) before truncating the expansion at a given order in nondimensional anisotropy. The resulting “quasi-equilibrium” model is sometimes referred to as the “level 2¼” closure to distinguish it from the level 2½ version. When the CL vortex force is included in this model, Eq. (12) remains as in level 2½ and the algebraic Reynolds stress model (ARSM) becomes

For covariance fluxes at this level, the algebraic relation is

and, following KC94, scalar variances are assumed to follow

Turbulence closure still requires assumptions for the form of fluxes and for the dissipation length scale.

### c. The ARSM, flux forms, and stability functions

The form of turbulent fluxes and the associated stability functions arise as solutions of the ARSM under boundary layer assumptions that disregard horizontal variability in the mean fields. Coriolis terms may be included in the ARSM but are discarded following KC94 for the wind-driven boundary layers considered here, with large natural Rossby numbers for wind-driven turbulence with friction velocity and mixed layer depths . However, the omitted Coriolis terms may be significant in different parameter regimes, such as deep convection. For these assumptions, CL vortex force terms are now included in the ARSM of the level 2¼ quasi-equilibrium model with the following components:

The terms and , appearing after distribution on the right side of Eqs. (17d,e) for vertical momentum fluxes and , preclude using a down-gradient momentum flux assumption to obtain algebraic closure for this equilibrium state because there are no nontrivial substitutions from within the set of Eq. (17) that will render these terms tractably proportional to the Eulerian shear . Incorporating CL vortex production into the equilibrium model therefore forces the assumption that momentum flux has the form

incorporating a new stability function for the component down the Stokes drift gradient, while scalar fluxes retain the form

Based on LES solutions of Eq. (1), several studies (McWilliams and Sullivan 2000; Smyth et al. 2002) have inferred the need for such a fundamentally different assumption for the momentum flux term in the *K*-profile parameterization (KPP; Large et al. 1994), making in part proportional to the Stokes drift gradient . Recently, McWilliams et al. (2012) explored the implications for KPP of assuming that momentum is directed down the gradient of Lagrangian momentum, that is, . Here, a component of momentum directed down arises as a direct consequence of properly accounting in the algebraic closure for the contributions to second-moment production from the CL vortex force of Eq. (1). These contributions, leading through the last right term in brackets of Eq. (5) to terms proportional to in Eqs. (17d,e) of the Reynolds stress, entail instead an eddy coefficient that is generally distinct from the eddy viscosity .

Substitution for terms on the right in Eqs. (17d,e) and (17i) yields three equations relating , , and :

The fluxes in Eqs. (18) and (19) are then substituted into Eq. (20) and the cross terms are rearranged to give

This yields three equations

that produce rational expressions for the stability functions flux forms Eqs. (18) and (19) in terms of new nondimensional local forcing functions

in addition to those already appearing in second-moment closures without ARSM CL vortex forcing

The resulting stability functions are

where

As in the Galperin et al. (1988) equilibrium model used in KC94 and extended here, the stability functions do not depend explicitly on nondimensional shear forcing .

Such expressions for stability functions [Eqs. (23)–(28)] are typically subject to “realizabilty constraints” when invoked in the context of a SMC model where dynamic conditions may at any time be far from the equilibrium state they represent. These limit the permitted ranges of some nondimensional forcing functions to avoid producing impossible states such as *G _{M}* < 0, or with improbably small levels of . Otherwise, a new wave-modified SMC model of upper-ocean turbulence is obtained by combining these stability functions Eqs. (18)–(28), the prediction of

*q*

^{2}by Eq. (12), and the prediction of the dissipation or “master” length scale

*l*to determine contributions to vertical fluxes down the local gradients of temperature, salinity, momentum and Stokes drift. The stability functions are restated in the appendix as ratios of polynomials, along with a summary of SMC model constants introduced in sections 3 and 4.

## 3. LES solutions for SMC comparison

### a. LES forcing case sets

To tune the new model, SMC predictions are compared here with steady-state solutions from LES in HD08, where the Craik–Leibovich vortex force models the interaction of waves and turbulence. Steady-state forcing cases in HD08 are specified from the 10-m wind speed and the surface wave age , where is the spectral peak phase speed, using an empirical surface wave spectra and an associated wave-age-dependent neutral drag coefficient *C _{D}*. The set of LES simulations is composed of several subsets of forcing cases labeled as “

**Σ**

_{i},” where

**Σ**

_{1}is a matrix of forcing cases varying wave age over four values in for each of eight values in , with mean mixed layer depths of 61–82 m. Set Σ

_{2}repeats the eight fully developed cases of

**Σ**

_{1}using a monochromatic approximation from Li and Garrett (1993) to specify . Set

**Σ**

_{4}repeats the same eight case components of

**Σ**

_{1}with reduced by one-half. In HD08, case sets

**Σ**

_{3a}and

**Σ**

_{3b}provide a continuation of the subset of

**Σ**

_{1}to hurricane-strength winds using a drag coefficient that continues to increase with in

**Σ**

_{3a}, and one that remains constant for in

**Σ**

_{3b}. Mixing dynamics in these high-wind, young-sea case sets are governed more strongly by entrainment zone shear. These case sets

**Σ**

_{1},

**Σ**

_{2},

**Σ**

_{3ab}, and

**Σ**

_{4}containing 55 LES steady-state solutions for wind and wave forcing that are described in greater detail in HD08 and are listed in correspondingly numbered tables therein. They are supplemented here by three additional wave-free forcing cases corresponding to

**Σ**

_{1}(

*U*

_{10}= 8.3 m s

^{−1}),

**Σ**

_{2}(

*U*

_{10}= 8.3 m s

^{−1}),

**Σ**

_{3}(

*U*

_{10}= 32.6 m s

^{−1}), and

**Σ**

_{4}(

*U*

_{10}= 8.3 m s

^{−1}), for a total of 59 cases. Figures 1–3 present dynamically scaled overviews of relevant mean LES profiles. Profiles from a subset of cases

**Σ**

_{1},

**Σ**

_{2}, and

**Σ**

_{4}are shown in Figs. 1 and 2, and profiles from high wind cases of

**Σ**

_{3ab}are in Fig. 3.

### b. Scaling for bulk and near-surface TKE components

We reported (HD08) that the scaling of bulk, mixed-layer-averaged vertical kinetic energy (VKE) on friction velocity is predicted for aligned wind and wave forcing by a surface layer (SL) Langmuir number that is a variant of the McWilliams et al. (1997) turbulent Langmuir number , where the surface Stokes drift is replaced by a near-surface average over the upper 1/5th of the mixed layer, relative to a reference value in the lower layer:

This scaling effectively absorbs variations in the Stokes drift *e*-folding depth scale relative to the mixed layer depth, as well as higher order effects owing to variations in the shape of between monochromatic and broadband surface wave spectra. Figure 1 compares the performance of the two different Langmuir numbers in scaling the profiles, rather than the bulk averages, of TKE components with , that is, on either . Crosswind , downwind , and vertical TKE component profiles are shown for LES forcing sets **Σ**_{2} (Figs. 1a–c) and **Σ**_{4} (Figs. 1j–l), and for the mature *C _{p}*/

*U*

_{10}= 1.2 (Figs. 1d–f) and young

*C*/

_{p}*U*

_{10}= 0.6 (Figs. 1g–i) sea components of

**Σ**

_{1}. These TKE component profiles are shown for high wind cases

**Σ**

_{3ab}in Figs. 3a and 3b. As noted in HD08, is constant for the empirical spectra at fixed wave age, so La

_{t}varies only with changes in

*C*with

_{D}*U*

_{10}at constant

*C*/

_{p}*U*

_{10}.

This causes the scaling of TKE components on La_{t} to appear more effective within some Fig. 1 plots of LES case subsets at fixed *C _{p}*/

*U*

_{10}than it does between differing

*C*/

_{p}*U*

_{10}plots. The LES model used in HD08 advects , a dynamically simulated subgrid turbulent kinetic energy (TKE), used in combination with a subgrid length scale , set by either the model resolution scale or by the length scale of smaller unresolved turbulence, to set the local nonlinear subgrid viscosity, that is, . The LES profiles (Figs. 1a–c and 2b) combine explicitly resolved TKE with a contribution from the advected subgrid TKE used for LES closure. It is apportioned equally into the TKE components except near the surface, where an adjustment is made for presumed subgrid anisotropy (Harcourt et al. 2002; HD08) to produce more consistent profiles and bulk averages across variations in grid resolution. The overall scaling comparison shows that, although scales the layer-averaged bulk VKE studied in HD08 accurately and provides a better interior mixed layer scaling for net and component TKE profiles, provides a better scaling for near-surface TKE levels and for corresponding boundary values in the SMC.

### c. Near-surface scaling of dissipation length scale l with depth

To evaluate the SMC prediction of and the relevance of the underlying equilibrium model, the dissipation length scale is diagnosed from the LES solutions and averaged over a period of steady-state evolution. The dissipation of explicitly resolved TKE is the shear-production source term in the conservative prognostic equation for advected *e*_{SG}, and that includes parameterized subgrid buoyancy fluxes and dissipation , where both the coefficient and the subgrid length scale may depend on and stratification in the pycnocline. To compare with in a SMC or an equilibrium model, the dissipation scale in the LES may be diagnosed as

using a traditional KC94 value for the coefficient in Eqs. (7a) and (12). The LES closure is based on a conservative prognostic equation for advected subgrid TKE , where the major production source is dissipation of resolved TKE and where dissipation . Because mean subgrid dissipation profiles were not routinely saved during the LES runs, the diagnosis of here relies on a combination of offline computations of profiles for select cases and an estimate using mean stability and subgrid TKE profiles and .

Figure 2 shows the diagnosis of the dissipation length scale based on and the estimate of for the forcing sets **Σ**_{2} and **Σ**_{4}, and for the young *C _{p}*/

*U*

_{10}= 0.6 and mature

*C*/

_{p}*U*

_{10}= 1.2 sea components of

**Σ**

_{1}, arrayed as in Fig. 1; Figs. 3d–f provide the corresponding profiles for high-wind cases

**Σ**

_{3ab}. Profiles of estimated are shown scaled by both and , and profiles of , scaled by both and , are shown in Figs. 2 and 3.

In Figs. 2and 3 profiles are also shown for estimated versus , for the same LES case subsets along with profiles indicating the buoyancy length scale , using buoyancy frequency . These profiles are given for the LES cases in Fig. 1, and an additional profile for each set shows a wave-free, low-wind case. In a near-surface region, the estimated LES dissipation length-scale estimate is constrained approximately to , a scaling significantly higher than law of the wall predictions of turbulence dissipation and length scale for nonslip shear boundary layers without wave effects. Even the profiles for wave-free low-wind cases well exceed the scaling expected for nonslip wall boundaries, though they are ultimately much smaller in the interior and lower layers than the other profiles with Langmuir turbulence.

Figure 4 compares profiles of computed offline from saved model fields for three LES forcing cases, to length-scale profiles estimated using (as in Fig. 2). These LES cases are identified in HD08 as **Σ**_{1} (*U*_{10} = 14.8 m s^{−1}, *C _{p}*/

*U*

_{10}= 0.8) for a young sea, Σ

_{1}(

*U*

_{10}= 18.1 m s

^{−1},

*C*/

_{p}*U*

_{10}= 1.2) for a mature sea, and for a high-wind case with very young waves and strong pycnocline shear,

**Σ**

_{3b}(

*U*

_{10}= 44.5 m s

^{−1},

*C*/

_{p}*U*

_{10}= 0.6).

The estimates of are shown to be about 1.2 times the more accurately calculated profiles above the pycnocline. This difference is due to the replacement of by in , and therefore corresponds to a scaling for the leading omitted Taylor expansion term of . However, the dissipation length scale still well exceeds for the cases with strong Langmuir turbulence.

For comparison, profiles of a related Langmuir turbulence dissipation length scale are shown in Grant and Belcher (2009, hereafter GB09, their Fig. 6), based on the dissipation of their resolved TKE . Adopting again the traditional value for the dissipation constant , their reported results for Langmuir turbulence have approximately between and until the level is reached at middepths. These results of GB09 for Langmuir turbulence are less in violation of the expectation, transplanted from solid-wall boundary layers, that . However, the profile provided in Fig. 6 of GB09 from a comparison case with no Stokes drift translates to a dissipation scale that exceeds .

Further development and tuning of a SMC model assumes that in the near-surface region of strong vortex force TKE production, a corresponding growth in the dissipation length scale occurs, with the result that both it and the vertical TKE component may be about double their values in the nonwave case. This expectation is in line with the qualitative understanding that the presence of Langmuir circulation structures, embedded in the turbulent boundary layer flow, entails an increase in the energy injection rate and energy levels at the larger *O*(*H*_{ML}) scales characterizing the separation between jets. As a result there should be a corresponding near-surface increase in the turbulence decay time scale , and enforcing conformity to the expectation of would defeat this additional effect of the Craik–Liebovich vortex force over the very depth range where its effect should be strongest. However, given what disagreements there are between the LES results for dissipation length in this and other studies, subsequent model adjustments may be required.

## 4. The second-moment closure model

Equation (12) for *q*^{2} and the stability functions [Eqs. (25)–(28)] are extended into a SMC model with the prognostic determination of *l* after KC04 through

Standard closure constants are retained as possible following KC94 and KC04, with unaltered values for {*A*_{1} = 0.92, *A*_{2} = 0.74, *B*_{1} = 16.6, *B*_{2} = 10.1, *C*_{1} = 0.08, *C*_{2} = 0.7, *C*_{3} = 0.2, *S*_{q2} = 0.41*S _{H}*,

*S*

_{l}= 0.41

*S*,

_{H}*E*

_{2}= 1.0}, and new constants for Stokes effects in third moment closures are taken here to be and . The surface boundary condition for energy is informed here by HD08 near-surface LES results for

*q*

^{2}(Fig. 2, middle column, dashed lines), and using the more appropriate Langmuir number scaling assumption it is set at . This gives under typical oceanic conditions, and reverts to the KC94 level in wave-free cases. Note that neither this SMC surface boundary condition nor the forcing of the HD08 LES model cases contain any additional TKE from wave breaking; additional changes in this boundary condition would be necessitated by directionally misaligned wind and waves. The large increase required in KC04 for the length-scale diffusion coefficient

*S*

_{l}over its KC94 value has not been found necessary. Indeed

*S*

_{l}=

*S*

_{q2}appears necessary to produce broadly well-behaved profiles in the lower mixed layer. Several more discontinuous departures from standard SMC implementations are described below.

Near-surface LES comparisons motivated a modification of the wall damping function with a dependence on La_{t} in , a change from that is equivalent within the wall function dissipation term to an increase in the Von Kármán constant from 0.4 to 0.7 for a typical oceanic case with . Simulations without Stokes drift but with strong entrainment shears (not shown) can also have larger near-surface scalings for . If these LES results are correct, it suggests further modifications of this wall damping term may be necessary, but these are not implemented in the SMC model version presented here, where the focus is typical open ocean conditions with waves.

The pattern of *l* values diagnosed from LES solutions in the lower boundary layer are better replicated by increasing the buoyancy forcing coefficient to *E*_{3} = 5.0, in line with the generally larger values suggested by Burchard (2001) as an alternative to restricting *l* values used to compute the stability function so that it not fall below the Ozmidov scale . However, because tuning for entrainment fluxes still benefits from , that is, after KC94, this limitation was retained. Trial and error with SMC stability suggests it is also necessary to restrict *l* values on the upper end in the stability functions to keep , in analogy to the standard realizability condition of KC94, which is also retained. Here, these limitations, and an additional requirement for stability that , were applied to *l* values throughout its use in determining the stability functions, but not directly to the prediction of .

Furthermore, the blunt shape of some interior profiles (e.g., over for the higher wind cases in Fig. 2f) are reproduced better by adding a dependence on the projections

and

of the Eulerian and Stokes shear direction into the Lagrangian , through modifying the coefficient of shear production from 1.8 to , and vortex force production by . The underlying rationale is that neither shear nor vortex force processes can effect an increase in the integral length scale unless the relevant shears are in line with the Lagrangian one, a condition suggested by results of Van Roekel et al. (2012) for the scaling of energy when wind and wave directions are not aligned.

Toward the bottom of the mixed layer a limitation is imposed on the relative vertical decay with depth of each eddy coefficient to not exceed while , where is the maximum value of within the layer where and , and where is the lower layer stratification at the depth where . Applied at grid level as , this only impacts the tail of relatively small eddy coefficients in the lower mixed layer and is restricted to on the presumption that it represents the effect of internal gravity wave dynamics in the pycnocline (cf. Polton et al. 2008). The resulting eddy coefficient “tails” have a significant impact on the steady-state entrainment rates when pycnocline shear is not large because without them buoyancy flux is extinguished immediately below the entrainment minimum, resulting in the growth of excessively concentrated stratification at that depth, by comparison to the LES. As an additional measure to prevent the growth of instabilities in transient solutions, the value of in the pycnocline is also not permitted to exceed .

Given these model features governing the prediction of *q ^{2}l* and eddy diffusivities, the coefficient

*E*

_{6}of the vortex TKE production is left to be determined by tuning SMC predictions to match LES results. Here, this is done for the ensemble of LES cases on the basis of and the mixed layer average of . However, it is notably more difficult to tune SMC performance to a large ensemble of LES results than it would be to a small handful; Figs. 5 and 6 illustrate that choosing a value for

*E*

_{6}presents a dilemma between tuning to minimize error in predicting (Figs. 5b,d) with

*E*

_{6}= 4.0 versus predicting energy levels (Figs. 5a,c) using

*E*

_{6}= 7.0. To predict mean vertical TKE (Figs. 5e,f) by inference from the SMC equilibrium model [Eq. (17c)] for comparison with Lagrangian float measurements (D’Asaro 2001; Harcourt and D’Asaro 2010), an intermediate value of

*E*

_{6}= 5 would be better. Figure 6 shows SMC to LES comparisons on the (power) rate of energy conversion into work done against gravity by the buoyancy flux due to entrainment, computed both by integrating from the depth of the minimum to the surface (Figs. 6a,b) and as an integral over the full model domain (Figs. 6c,d). These comparisons for two different -scaled metrics of entrainment (Figs. 6a–d) favor the larger value

*E*

_{6}= 7.0 when, as is usually the case, predicting vertical buoyancy flux is the primary consideration. Some outlying values are due to excessive entrainment in LES case set

**Σ**

_{4}owing to inadequate vertical domain size and excessive interaction with the radiative bottom boundary conditions. Still, even after discounting these, the tuning to entrainment is not strongly compelling; there is a large scatter (Figs. 6a,c) and the poor dissipation length-scale predictions for

*E*

_{6}= 7.0 (Fig. 5a) belie the model’s purportedly more accurate articulation of the underlying physics. The equivocation here on the value of

*E*

_{6}echoes different values reported by KC04 (

*E*

_{6}= 4.0) and Carniel et al. (2005, p. 36), which recharacterizes the KC04 value to

*E*

_{6}= 7.2. Kantha et al. (2010) reconfirm that the value in KC04 was reported incorrectly and that the larger value should be used. However, it is unclear if the ambiguity here stemming from tuning priorities is directly relevant to the similarly different choices for

*E*

_{6}reported for these prior studies.

### SMC model profiles

Profiles of SMC eddy coefficients and the underlying stability functions are shown in Fig. 7 for the three example cases of Fig. 4, with the more accurate calculation of . The profiles for these cases demonstrate how the shape and order of these mixing coefficients changes with *E*_{6} and with the relative rate of vortex force TKE production, switching from in the lower mixed layer to under stronger Stokes shear and toward the surface. Large changes in eddy coefficients correspond as much, if not more, to variations in stability functions than they do to changes in . This reflects both the increased ARSM contributions to from vortex force production [Eq. (17c)], and the transfer of that increase into the leading turbulent flux terms, that is, the first right side terms in Eqs. (17e,d,i). Neither of these processes is included in the KC04 model.

The shape of the diffusivity profiles differs between the coefficients and varies with the relative penetration of Langmuir turbulence into the layer, with the order of buoyancy and momentum coefficients reverting to at much lower levels when . Under strong Langmuir forcing (Figs. 7b,c), profiles of these coefficients are more similar to those prescribed for the more empirical *K*-profile parameterization of Large et al. (1994) than they are to those in KC04 (their Fig. 1), due in large part to the higher value of *E*_{3} = 5.0 used for the effect of buoyancy flux on through Eq. (31). This increased sink of in response to entrainment of heavy water brings the dissipation length-scale profiles of Fig. 8 more in line with LES results in the lower mixed layer. Energy profiles in Fig. 9 show interior SMC predictions of to be more accurate in conjunction with the overpredicted values when *E*_{6} = 7.0, but sometimes more consistent near the surface when the *E*_{6} = 4.0 SMC values are close to the LES profile. While LES and SMC profiles are significantly elevated over the SMC case throughout most of the mixed layer in the case **Σ**_{1} (*U*_{10} = 18.1 m s^{−1}, *C _{p}*/

*U*

_{10}= 1.2) of Fig. 9b with relatively deep penetration of Stokes shear, interior SMC TKE levels differ little at depths where Stokes shear is small. In these cases (Figs. 9a,c) where Stokes shear is either restricted within the surface layer or small compared to Eulerian shear, the most significant change in the strength of mixing then comes primarily from the larger new stability functions, and to a lesser extent from increases in . Both and diverge from LES profiles in the pycnocline where the dynamics of internal gravity waves are fundamentally not represented by the SMC closure assumptions, for example, unless the measure of is restricted to scales within the inertial subrange. At depths below about 0.9, 1.0, and 1.1 times (Figs. 7a,b,c, respectively), the limitation on the decay rate with depth of the eddy coefficients serves to mitigate to some extent the cessation of mixing with SMC-predicted and when the closure ceases to be as relevant. The effect on buoyancy flux profiles in Fig. 10 allows the SMC to follow the curve of the entrainment zone into the pycnocline, which in turn prevents a barrier of stratification forming too rapidly at the mixed layer bottom. This correction generally does not extend to cover the lower tail of the decay of the LES flux profile, but background eddy diffusivities representing the mixing effects of much larger-scale near-inertial internal wave climate in the ocean (not a part of LES forcing) would become comparable to SMC predictions at these depths. Even aside from these details of pycnocline mixing there are some significant discrepancies in mixed layer entrainment, such as the comparison in Fig. 10b, that are not eliminated by this SMC model tuning.

Figure 11 demonstrates one of the major outcomes of modifications introduced into the model by the component of momentum flux down the Stokes drift gradient, most markedly by comparison with SMC results where . While LES downwind Eulerian shear profiles are positive approximately over a near-surface layer, those cases with deeply penetrating present a retrograde shear below this level that is predicted with varying skill by the new SMC; it is predicted well using *E*_{6} = 7.0 for the case in Fig. 11b. For most comparison cases, however, the retrograde SCM shear continues to strengthen toward the surface and does not revert into a downwind shear in the approach to the surface boundary. Another significant difference between LES and SMC predictions is apparent in the comparison of vertical TKE profiles in Fig. 12, where the SMC profile is inferred from the equilibrium model Eq. (17c) for . While the SMC does not explicitly predict vertical kinetic energy, this underlying ARSM expression significantly impacts the stability functions because is the energy factor in the leading terms [Eqs. (17d,e) and (17i)] for these vertical fluxes. While there is a clear improvement over SMC models with , this implicit profile remains surface intensified except at , which does not exhibit the subsurface maximum of the LES profiles, and perhaps most critically it underpredicts in the lower mixed layer. The low values of in the lower layer entail underprediction of and of vertical fluxes at these depths where Stokes drift is small, making it necessary to commit compensating errors elsewhere in the closure at the ultimate expense of model generality. The effect of instead correcting the equilibrium model through improved closure assumptions, or a higher-order closure with a separate prognostic equation for , would then be to shift the concave shape of diffusivities downward, thereby increasing lower mixed layer and pycnocline diffusivity.

In Figs. 13 and 14 the equilibrium model [Eq. (17)] for the Reynolds stress tensor is evaluated to examine the self-consistency of its predictions in the context of the three example LES case steady-state solutions. The comparison in Fig. 13 of the resolved TKE components with their corresponding right side equilibrium model expressions shows that while differences for the downwind TKE component are relatively small, a significant term responsible for rotating vertical TKE into the crosswind component is missing from the closure. Further analysis, not shown here, suggests this term is missing from the pressure–strain closure [Eq. (8)] and would represent the acceleration of near-surface water into the pressure low along Langmuir jets, in response to the downwelling driven by the Craik–Leibovich force below them. Diagnosis of the missing pressure–strain closure terms is the subject of further ongoing research to improve the SMC.

Applying the same self-consistency test to the Reynolds stress tensor cross terms in Fig. 14 shows that for (Figs. 14d–f) the large retrograde SMC shear in downwind momentum corresponds to near-momentum flux levels in the equilibrium model level diagnosed from LES that well exceed the surface flux . While the ad hoc pressure–strain closure term [Eq. (8)] provides a plausible and tractable interior contribution with , the large retrograde shear in the near-surface layer may be due to the inaccuracy of this closure assumption. Indeed, it seems unlikely that the deformation of turbulence by Stokes and Eulerian shears would be equivalent as implied by Eq. (8).

On the other hand, the LES subgrid closure accounts (as in KC04) for only the additional vortex force production of subgrid TKE and does not include a subgrid momentum flux component down the Stokes gradient. The significance of this omission would increase as the ratio of subgrid to resolved TKE increases toward the surface, coincident with the discrepancy between LES and SMC Eulerian shear. As the excessive vertical momentum flux in the surface layer is similar to that component in the LES (dot-dashed in Figs. 14d–f) down the Eulerian shear, it is therefore also possible that this shear develops erroneously or excessively in the LES solutions in response to the lack of a correct subgrid momentum flux component down the Stokes drift gradient.

Other differences in the momentum fluxes suggest that several other smaller closure contributions are missing. The covariance (Figs. 14g–i) of horizontal momentum is assumed without consequence to be zero in KC94 and indeed the small non-Stokes drift terms on the right side of Eq. (17f) may well cancel against closure expressions not included here. However, the near-surface comparisons for both (Figs. 14g–i) and (Figs. 14a–c) suggest that additional surface-intensified closure contributions are missing in these components that may be proportional to Stokes shear. While it is possible to find functions of the local nondimensional forcing that are consistent with these missing closure terms, their physical import can be unclear and the resulting equilibrium model might not yield rational expressions for stability functions. More accurate closure expressions may then necessitate a higher level of closure with more than two prognostic equations.

## 5. Summary

A new level 2¼ second-moment closure (SMC) model was developed that extends the model of KC04 to include Langmuir turbulence effects in the algebraic Reynolds stress model (ARSM) and in the resulting stability functions and turbulent flux closure. This required adding vortex force TKE production in the ARSM, as well as the introduction of a new momentum flux component that is directed down the gradient of the Stokes drift, in addition to the conventional term down the gradient of the Eulerian momentum. Relative to KC04, the new model includes changes in the momentum flux closure [Eq. (18)] and in the response of stability functions to Stokes shear [Eq. (23)–(28)] that result directly and unequivocally from the inclusion of the CL vortex force in all components of the Reynolds stress tensor equation [Eq. (5)] that is used to derive the ARSM. Additional changes in the stability functions stem from the generalization of KC94 closure assumptions for pressure–strain and pressure–scalar correlations [Eqs. (8) and (9)] and are subject to corresponding choices in two new closure constants. Several other SMC model components were modified to conform to a suite of LES simulations for mixed layers with varying degrees of Langmuir forcing. Tuning the SMC model presents a dilemma between skill at predicting the dissipation length scale versus predicting mixed layer TKE and the entrainment rate. In general, the eddy coefficients for momentum flux due to Eulerian and Stokes shear vary independently with the relative strength of Langmuir and shear-driven turbulence because of the corresponding dependence of leading closure terms on different components of the TKE in the ARSM. Analysis of the equilibrium model using LES results suggests several closure terms are still missing, notably a pressure–strain contribution responsible for transferring vertical into crosswind TKE. The new SMC model improves the prediction of momentum profiles, reproducing a retrograde Eulerian shear in mixed layer interiors, and suggesting that downwind near-surface LES shear profiles not replicated by the SMC model may be at least partly due to missing LES subgrid flux components directed down the Stokes drift gradient.

Modifications made to the wall-damping function, to the decay of eddy coefficients below the mixed layer, and the introduction of functional dependence for equation coefficients on local nondimensional forcing scales, may each have limited general validity. While these features served to pull SMC performance into line with HD08 LES results, they may not represent the most robust solutions and should be considered subject to further testing through a broader range of LES forcing cases and modeling techniques as well as empirical validation. The new SMC model serves primarily here as a vehicle to articulate the significant effects of the CL vortex force on both the momentum flux closure, and on the SMC equilibrium model. It is presented here in the hopes that practitioners of the many other approaches to turbulence closure distinct from Mellor and Yamada (1982), KC94, and KC04 may also consider corresponding adaptations to better represent the effects of Langmuir turbulence. Future studies may include considering a broader possible set of turbulence closure assumptions, examining the interaction of energy injected by wave breaking with the new momentum flux closure, and exploring the relationship between the Stokes-mediated vertical momentum flux and the evolution of surface wave spectra.

## Acknowledgments

This work was supported by the National Science Foundation (OCE0850551 and OCE0934580), the Office of Naval Research (N00014-08-1-0575), and by a grant of HPC resources from the Department of Defense High Performance Computing Modernization Program.

### APPENDIX

#### Summary of Stability Functions and SMC Model Constants

The ARSM stability functions may be written as ratios of polynomials:

where

and where and share the common denominator

ARSM closure constants unchanged from KC04 are *A*_{1} = 0.92, *A*_{2} = 0.74, *B*_{1} = 16.6, *B*_{2} = 10.1, *C*_{1} = 0.08, *C*_{2} = 0.7, and *C*_{3} = 0.2. New ARSM closure constants are assumed to have values and . The nondimensional forcing functions , , , and [Eqs. (23) and (24)] and eddy flux coefficients , , , , and are subject to limiting values. Stability functions of KC94 are recovered for . The prognostic equations for TKE [Eq. (12)] and length scale [Eq. (31)] use *S*_{q2} = 0.41*S _{H}*,

*S*

_{l}= 0.41

*S*, ,

_{H}*E*

_{2}= 1.0,

*E*

_{3}= 5.0, , (

*E*

_{5}= 0), and with , given by Eq. (32). SMC results compared above use either

*E*

_{6}= 4.0 or

*E*

_{6}= 7.0.

## REFERENCES

^{2}l equation by Mellor and Yamada, 1982

**60,**629, doi:10.1007/s10236-010-0283-5.