## Abstract

A large-eddy simulation (LES) model is configured to investigate the effect of the horizontal (northward) component of Earth’s rotation on upper-ocean turbulence. The focus is on the variability of the effect with latitude/hemisphere in the presence of surface gravity waves and when capped by a stable stratification beneath the surface layer. When is included, the mean flow, turbulence, and vertical mixing depend on the wind direction. The value and effect of are the largest in the tropics and decrease with increasing latitudes. The variability in turbulent flows to wind direction is different at different latitudes and in opposite hemispheres. When limited by stable stratification, the variability in turbulence intensity to wind direction reduces, but the entrainment rate changes with wind direction. In wave-driven Langmuir turbulence, the variability in mean current to wind direction is reduced, but the variability of turbulence to wind direction is evident. When there is wind-following swell, the variability in the mean current to wind direction is further reduced. When there is strong wind-opposing swell so that the total wave forcing is opposite to the wind, the variability in the mean current to wind direction is reduced, but the variability of turbulence to wind direction is enhanced, compared to in Ekman turbulence. The profiles of eddy viscosity, including its shape and its value, show a strong wind direction dependence for both stratified wind-driven and wave-driven Langmuir turbulence. Our study demonstrates that wind direction is an important parameter to upper-ocean mixing, though it is overlooked in existing ocean models.

## 1. Introduction

Turbulent flows in the oceanic surface boundary layer (OSBL) are important to Earth’s climate and marine environment. They control the air–sea exchange of heat, momentum, and trace materials. They also determine the transport, dispersion, and transformation of materials such as nutrients, larvae, and pollutants in the upper ocean. OSBL turbulence is driven by surface atmospheric (or meteorological) conditions including winds, waves, and heat fluxes and is modulated by stratification and the rotation of Earth. While the effect of different surface meteorological conditions on OSBL turbulence has been relatively well studied through a series of observation programs (e.g., Price et al. 1978; Weller and Price 1988; Smith 1992; Plueddemann et al. 1996; Gargett et al. 2004; Fer 2006; Sanford et al. 2011; D’Asaro 2014; D’Asaro et al. 2014) and numerical studies (e.g., Skyllingstad and Denbo 1995; McWilliams et al. 1997; Li et al. 2005; Belcher et al. 2012), the effect of Earth’s rotation, particularly its horizontal component, is relatively less studied and is insufficiently understood.

Studying the effect of planetary rotation on OSBL dates back to the seminal work by Ekman (1905), who derived analytic solutions for the mean horizontal current in OSBL, known as the Ekman spiral, under the influence of a steady uniform surface wind and the vertical component of Earth’s rotation ( where is the angular speed of Earth’s rotation, and *φ* is the latitude). The Ekman theory successfully explains the wind-driven transport at the surface; however, recent studies (e.g., Price and Sundermeyer 1999) have demonstrated discrepancies between observed velocity profiles with Ekman’s solutions.

The horizontal component () of planetary rotation is neglected in the momentum equation in circulation models because it is small compared to other terms, such as the Coriolis force and pressure gradient term in the momentum budgets (White and Bromley 1995). However, recent studies (Etling and Wippermann 1975; Leibovich and Lele 1985; Garwood et al. 1985a,b; Galperin et al. 1989; Kantha et al. 1989; Coleman et al. 1990; Hassid and Galperin 1994; Wang et al. 1996; Zikanov et al. 2003; McWilliams and Huckle 2006; Wang 2006) reveal that alters the mean horizontal velocity as well as vertical mixing in the OSBL.

Through a series of linear stability analyses, Etling and Wippermann (1975) and Leibovich and Lele (1985) showed that the horizontal component of planetary rotation can act as a destabilizing agent in the planetary Ekman layer by enlarging the bands of unstable wavenumbers and thereby reducing the critical Reynolds number.

Using a bulk mixed layer model, Garwood et al. (1985a,b) suggested that the interaction between Reynolds stress and horizontal component of planetary rotation influences the conversion of turbulent kinetic energy (TKE) between horizontal and vertical components. They hypothesized that the horizontal component of Earth’s rotation, together with wind and buoyancy, induces a deeper mixed layer in the western tropical Pacific and a shallower mixed layer in the eastern tropical Pacific.

Using the level 2.5 second-order turbulence closure model (Mellor and Yamada 1982), Galperin et al. (1989) and Kantha et al. (1989) examined the impact of the horizontal component of Earth’s rotation on turbulent mixing in stably stratified, shear-driven turbulent flows near the equator and concluded that the effect of on mixed layer depth (MLD) is small due to the stable stratification. However, Hassid and Galperin (1994) found that the effects of Earth’s rotation (both horizontal and vertical components) can be significant under neutral and unstable stratification.

More recently, the effect on OSBL turbulence was studied using turbulence-resolving models. In a direct numerical simulation (DNS) study of the turbulent atmospheric Ekman layer with neutral stratification, Coleman et al. (1990) showed that the wind direction led to 6% variations in the surface friction velocity and 9° difference in the maximum and minimum of the angle between the wall shear stress and the free-stream velocity at 45°N. Using a large-eddy simulation (LES) model, Wang et al. (1996) extended the bulk mixed layer model to include entrainment and found that the effect of on OSBL turbulent mixing at the equator is small when entraining heat flux is significant. In a subsequent study, Wang (2006) found that the effect of on turbulent convection at 60°N is also small.

Zikanov et al. (2003) revisited the classic Ekman theory and investigated the flow dependence on latitude and wind direction by considering a turbulent shear flow driven by a steady surface wind using an LES model. They showed that both the mean current and the turbulence change with wind direction due to , and the effect increases with decreasing latitudes. They demonstrated that wind direction exerts a nonnegligible impact on wind-driven Ekman turbulence due to the horizontal (northward) component of Earth’s rotation: turbulence is enhanced when the wind has a westward and southward component and is suppressed when the wind has an eastward and northward component in the absence of other processes in the Northern Hemisphere.

However, to the best of our knowledge, no study has systematically investigated the influence of three possibly important factors: hemisphere, stable stratification below the mixed layer, and nonbreaking surface gravity waves. Hemisphere is a nonnegligible factor when studying the effect of rotation on OSBL turbulence. It affects the direction of the Coriolis force, which limits the evolution of turbulence, since the vertical rotation is reversed in opposite hemispheres. At low to midlatitudes where is relatively large, the mixed layer is relatively shallow due to strong surface heating. Stable stratification below the mixed layer suppresses the development of turbulence and traps the Ekman transport within a relatively shallow mixed layer, which also limits the exchange of materials inside and outside the boundary layer. Surface gravity waves are ubiquitous over the global ocean. They drive coherent Langmuir circulations (LCs; or Langmuir turbulence) that modify the mean currents, mixing, and dissipation profiles (e.g., Sullivan et al. 2007; Belcher et al. 2012). LCs were first observed by Langmuir (1938), later theoretically modeled by Craik and Leibovich (1976), and more recently studied numerically via the LES of the Craik–Leibovich equations (e.g., McWilliams et al. 1997; Li et al. 2005; Tejada-Martínez and Grosch 2007; Harcourt and D’Asaro 2008; Pearson et al. 2015).

The objectives of this study are 1) to expand the analysis of turbulent Ekman flow under the influence of , 2) to study the variability of the effect for varying latitude and hemisphere, and 3) to investigate how the addition of nonbreaking surface gravity waves and stratification modify the effect of . The objectives are achieved by analyzing a series of solutions from a large-eddy simulation model. The rest of the paper is organized as follows. Section 2 presents the theoretical background necessary for understanding of the effect of on upper-ocean turbulence. Section 3 describes the LES model and model configuration. Section 4 presents the numerical solutions and explores the effect of on the depth of OSBL, eddy viscosity profile, TKE production, TKE dissipation, and residence time of TKE in the absence of stratification and surface gravity waves, and it also describes the latitude/hemisphere dependence of the effect and the modification of the effect due to stratification, Langmuir turbulence, swell, and the combination of stratification and waves. Section 5 discusses the effect of on eddy viscosity profiles. Section 6 is a discussion and summary.

## 2. Background theory

The effect of on turbulence can be understood by the TKE and the momentum flux budget equations. Following Zikanov et al. (2003), the equations are

where *u*, *υ*, and *w* are the components of velocities in *x* (downwind), *y* (crosswind), and *z* (positive upward) directions, respectively (i.e., the equations are in the wind coordinate; see Fig. 1). A horizontal mean is denoted by angle brackets , and a deviation from the horizontal mean is denoted by prime (for a variable *X*, ). Other variables that appear in the above equations are the horizontal component of Earth’s rotation vector projected in *x* () and *y* () directions, where *θ* is the angle measured from the northward direction (the direction of ) to the downwind direction in a clockwise sense (e.g., northward wind for ). The dots on the right-hand side (RHS) of Eqs. (1)–(5) represent other terms, including shear production, turbulent transport, the Coriolis effects associated with the vertical component of planetary rotation, buoyancy, production/destruction due to Stokes vortex force, and pressure transport (McWilliams et al. 2012) that are identical to when .

The horizontal component of Earth’s rotation affects OSBL turbulence in two ways: 1) through the redistribution of turbulent kinetic energy between horizontal and vertical directions and 2) through the modification of the vertical momentum fluxes ( and ) (e.g., Zikanov et al. 2003). While the sign of terms in the budget equation of a variable does not imply the increase or decrease of the variable, Zikanov et al. (2003) showed that for an unstratified Ekman layer, is more negative, is more positive, and volume-averaged TKE increases when both and are negative (). In contrast, is less negative, is less positive, and volume-averaged TKE decreases when both and are positive () in the Northern Hemisphere (see Figs. 8a,b, 9a in Zikanov et al. 2003).

## 3. Model description and configuration

### a. Model description

The model used in the study is the National Center for Atmospheric Research (NCAR) LES model for a wave-driven oceanic boundary layer (McWilliams et al. 1997; Sullivan and McWilliams 2010), with the addition of the horizontal component of planetary rotation. The filtered (i.e., LES) Craik–Leibovich equations for momentum in a Cartesian coordinate system with Earth’s rotation are (e.g., McWilliams et al. 1997; Moeng and Sullivan 2002; Suzuki and Fox-Kemper 2016)

where filtered (resolved) variables are denoted by lower cases. Parameter *π* is the generalized pressure field (e.g., Sullivan et al. 2007): with resolved pressure *p*, reference density , resolved velocity , and Stokes drift . Term *e* is the subgrid-scale (SGS) TKE calculated following Sullivan et al. (2007); the SGS viscosity is calculated without the near-surface correction by Sullivan et al. (1994); are the resolved vorticities; are the SGS fluxes as evaluated in the LES parameterization model; *f* is the Coriolis parameter; and *ρ* is the water density. An upwind scheme is used for scalars in the vertical advection.

The model has been used to study the oceanic surface boundary layer turbulence driven by a variety of surface and lateral boundary conditions (e.g., Sullivan and McWilliams 2010; Kukulka et al. 2011; McWilliams et al. 2012; Sullivan et al. 2012; Van Roekel et al. 2012; Hamlington et al. 2014; McWilliams et al. 2014). It has also been successfully applied to the study of tracers including marine debris (Brunner et al. 2015; Kukulka and Brunner 2015; Kukulka et al. 2016), gas bubbles (Liang et al. 2011, 2012, 2017), and other tracers (Smith et al. 2016).

### b. Model configuration

The model is configured on a rectangular domain of 300 m × 300 m × 300 m with 200 × 200 × 128 grids. The horizontal grids are evenly spaced ( m) and the vertical grids are stretched with the finest resolution ( m) at the surface, and sensitivity tests confirm that further decreases in grid spacing and grid anisotropy factor do not qualitatively change the result (not shown). Four groups of simulations, differing in surface meteorological conditions, are summarized in Table 1. Case *E* is experiments driven solely by wind; Case *L* is experiments driven by both wind and wave in equilibrium with wind; Case *S* is experiments driven by wind and wind-following swell; and Case *M* is experiments driven by wind and wind-opposing swell. These selected swell conditions are extreme; the effect of under a realistic weaker swell can be inferred from the combination of all four cases. Superscript *s* denotes experiments with stratification, and subscript *h* denotes experiments with both horizontal and vertical components of planetary rotation. Eight wind directions with constant wind speed (from northward wind clockwise to northwestward wind, with 45° interval) are included. The time step is dynamically determined by the Courant–Friedrichs–Lewy (CFL) condition (Sullivan et al. 1996). Each simulation was run for approximately three inertial periods following a spinup simulation for about one inertial period. In each experiment, all the statistics are averaged over two inertial periods after turbulence becomes fully developed, and a sensitivity test shows that a longer averaging period will not affect the statistics. In this paper, we choose 5 m s^{−1} as the wind speed at 10 m above the surface. The corresponding wind stress *τ* is 0.040 N m^{−2}, and the waterside friction velocity is 0.0063 m s^{−1}, using the drag law by Liu et al. (1979). For the case with stratification, the initial mixed layer is 40 m deep. Below the mixed layer, a stable stratification of , corresponding to a Brunt–Väisälä frequency of 0.0099 s^{−1}, is prescribed. The radiation boundary condition is applied at the lower boundary of the computational domain (Klemp and Durran 1983). We also choose 15°N as the latitude for controlled experiments (). Additional experiments at 15°S and 45°N were also carried out to investigate how the effect changes with latitude and hemisphere. The turbulent Langmuir number [, where is the Stokes drift at the center of the first vertical grid cell] is 0.35 for the wave-driven Langmuir case, and the Stokes drift follows the formula by McWilliams and Restrepo (1999) and Sullivan et al. (2007): , where *g* is the gravitational acceleration, is the radial frequency, and *F* is the wave spectrum proposed by Donelan et al. (1985) and Alves et al. (2003).

## 4. Model results

### a. Unstratified turbulent Ekman layer

In this section, we first revisit the turbulent Ekman layer by Zikanov et al. (2003) and expand the analysis of the role of in the turbulent Ekman layer. Our discussion is based on a series of numerical experiments at 15°N in the absence of stratification and nonbreaking surface gravity waves (experiments and *E*). For the convenience of subsequent discussion, we use the following notation: for a given variable *X*, denotes the horizontal average over the boundary (mixed) layer (layer-averaged variable); that is, , where is the boundary (mixed) layer depth, denotes vertically integrated *X* over the boundary (mixed) layer (i.e., ), and denotes the ensemble average.

An instantaneous snapshot of vertical velocity fluctuation can be used to illustrate the effect on the turbulent flow structure. Eddies are larger and more coherent horizontally when the wind blows southward to westward (; i.e., from southward wind to westward wind in a clockwise sense) than when , while eddies are smaller and more intermittent in horizontal dimensions when the wind blows northward to eastward (; Fig. 2). To quantify this, we use the velocity correlation function (Davidson 2015)

with , then calculate the shortest length where . Eddies that are larger and more coherent horizontally will have larger , and vice versa (see Fig. 3g). The instantaneous *w* patterns of different wind directions are not noticeably different close to the surface (not shown) and are evidently different in the middle of the boundary layer. This is due to the role of in rotating the horizontal velocity into vertical velocity. Although the equilibrium value of vertical velocity variance is not directly determined by the sign of the terms with (see section 2), the enhancement or reduction of vertical velocity variance is associated with the signs of , , , and [see Eq. (3)]. Vertical velocity variance increases when both and are positive. This occurs, for instance, when the wind is southwestward () in the Northern Hemisphere, when both and are negative, is negative, and is positive. Comparisons of instantaneous snapshots for *w* and show that the scale of eddies and vertical TKE are larger than when . In this regime, vertical turbulent mixing is strengthened. In contrast, vertical velocity variance decreases when both and are negative. This occurs, for example, when the wind is northeastward () in the Northern Hemisphere, when both and are positive, is negative, and is positive. Comparisons of instantaneous snapshots for *w* and show that the scale of eddies and vertical TKE are smaller than when . Vertical turbulent mixing is also weakened in this regime. The enhanced (weakened) vertical turbulent mixing also affects the vertical extent of the boundary layer. For instance, 60 m is inside the boundary layer for a southward wind, while it is outside the boundary layer for a northward wind.

Under the influence of *f*, the mean surface currents (, ) are deflected to the right of the wind and rotate clockwise with depth in the Northern Hemisphere (Ekman 1905), and the velocity spiral is independent of wind direction. In the presence of , both the deflection and the magnitude of the mean horizontal current vary with wind direction (Fig. 3a): when the wind blows northward to southeastward (), the mean surface current is deflected more to the right of the wind, as indicated by a larger angle between wind stress and surface current. When the wind blows southward to northwestward (), the deflection is less (cf. Fig. 7 in Zikanov et al. 2003).^{1} For instance, the mean surface current is 31.1° to the right of the wind when the wind blows northward (), while it is only 18.7° to the right when the wind blows southward (). The change in mean horizontal velocity by is due to the modulation of the vertical momentum fluxes (,) by (Zikanov et al. 2003). When the wind blows northward to southeastward (), the near-surface gradient for increases due to , resulting in a more negative . When the wind blows in the opposite directions (), the near-surface gradient for decreases due to , resulting in a less negative . The change in is less dramatic than the change in because of the smaller change in the near-surface gradient than in the near-surface gradient. The mean surface horizontal current ( and ) is particularly important in determining the transport and dispersion of many pollutants, such as spilled oil and marine debris, that are positively buoyant and are trapped near the surface. It is also important to determine the wind work to the ocean. The maximum time rate of wind work (, where is water density) when the wind blows eastward () is 30% larger than the minimum when the wind blows westward () at 15°N (not shown). The magnitude of the mean current decays faster when the wind blows northward to eastward than the opposite wind directions (Fig. 3a). The effect of on the mean horizontal Ekman current in the absence of stratification and surface gravity waves is consistent with the previous studies (Coleman et al. 1990; Zikanov et al. 2003; McWilliams and Huckle 2006).

Boundary layer depth is also significantly modified by (Fig. 3b). Here, the boundary layer depth is defined as the depth at which the horizontally and time-averaged momentum flux is 10% of the surface momentum flux . Sensitivity tests indicate that a smaller criterion (less than 10%) of surface momentum flux used in the definition of influences the magnitude of , but does not change the dependence of on wind direction (not shown). Compared with the case (dashed line), the boundary layer is deeper for southeastward to westward winds () and is shallower for other directions. The maximum is twice as large as the minimum .

The depth of the equilibrium OSBL is a function of the characteristic length scale of the turbulent flow. The commonly used length scales to characterize boundary layers are the Ekman depth (Rossby and Montgomery 1935) under the sole influence of wind stress and the Obukhov length (Monin and Obukhov 1954) with the combined effect of wind stress and a destabilizing buoyancy flux. However, neither length scale includes the wind direction dependence of boundary layer depth under the influence of .

We find that the wind direction dependence of is the same as the wind direction dependence of the layer-averaged vertical velocity variance (Fig. 3c), indicating the variation of is controlled by the vertical TKE. The is the average taken over the boundary layer, ranging from 0.51 to 0.72. The turbulent kinetic energy is redistributed from the vertical to horizontal component when both and are positive () and from the horizontal to vertical direction when both and are negative (). The redistribution of TKE can be inferred from the ratio between vertical and horizontal variances (; Fig. 3d). The ratio is larger than the case when and smaller than the case when . This is associated with larger for southward to westward winds as well. In this case, turbulence with a larger can penetrate deeper, which induces a deeper .

The vertically integrated mean kinetic energy (MKE) over the boundary layer [] is shown as a function of wind direction *θ* in Fig. 3e. It is clear that increases for cases , compared with the case, and decreases for cases . The variation of is affected by both the variation of momentum flux that can affect the shear production of turbulence and the variation of wind work at the water surface.

The vertically integrated turbulent kinetic energy over the boundary layer {}, normalized by , is shown in Fig. 3f. In particular, is smaller than the case for , which is consistent with Zikanov et al. (2003, their Fig. 9a). Note that the *f* plane used in Zikanov et al. (2003) is chosen at 90°N, which has a stronger suppression on turbulence due to a larger *f*.

The variability in to wind direction at the vertical depth (), where the instantaneous snapshot of is taken, is shown in Fig. 3g. The distance is smaller when the wind blows northward to eastward than the case. In contrast, is larger when the wind blows southward to westward than the case. This is consistent with the instantaneous snapshot (Fig. 2).

Contrary to the conjecture by Zikanov et al. (2003) that the difference in TKE production () among wind directions governs the difference in the magnitudes of TKE, the wind direction dependence of total TKE production (Fig. 4a) is different from that of (Fig. 3f). For instance, the total production of TKE is the largest (smallest) for eastward (westward) wind (Fig. 4a). The production of TKE is maximum at the surface where the mean current shear (, ) is the largest. It therefore has the same wind direction dependence as the surface mean current (Fig. 3a). Although can influence the profiles of momentum flux (,), its effect is larger in the middle of the boundary layer and is negligible close to the surface where shear is the largest (cf. Fig. 8 in Zikanov et al. 2003).

The wind direction dependence of layer-averaged eddy turnover time (the ratio of total TKE to dissipation; Fig. 4b) dominates the wind direction dependence of . It is also larger than the case for cases and smaller than the case for cases . The largest eddy turnover time is approximately twice as large as the smallest value in unstratified Ekman turbulence.

The wind direction dependence of eddy turnover time can also be interpreted by examining the dissipation length scale , defined as

where is the turbulent kinetic energy, and is the dissipation rate (Tennekes and Lumley 1972; Grant and Belcher 2009). The averaged dissipation length scale over the boundary layer is shown in Fig. 4c. Larger corresponds to larger , meaning that larger eddies with more energy spend more time dissipating energy into smaller eddies.

### b. Modification of the effect due to latitude and hemisphere

Compared with 15°N, the magnitude of the mean surface current (, ) decreases due to increasing *f* at 45°N (Fig. 5a). The variation of the angle between wind stress and surface current also decreases as latitude increases: the difference between the largest and smallest angle decreases from 12.6° at 15°N to 5.7° at 45°N. The boundary layer depth increases as *θ* increases when , which qualitatively agrees with the result that growth rate for the most unstable mode increases as the Ekman-layer surface-current direction increases based on a series of linear instability analyses (cf. Fig. 19 in Leibovich and Lele 1985). Moreover, in the same hemisphere, the variabilities in , , , , and , defined as the ratio of the difference between the maximum and minimum values to the minimum value, decrease from 120%, 43%, 79%, 98%, and 76% at 15°N (Figs. 3b–f) to 38%, 16%, 30%, 44%, and 19% at 45°N (Figs. 5b–f), respectively.

Although is the same at the same latitude in opposite hemispheres, the dependence of turbulence on wind direction is different because of the reversed vertical rotation in opposite hemispheres. An instantaneous snapshot of vertical velocity fluctuation at 15°S is used to show the major difference in instantaneous flow structure in the middle of the water column, where the effect is obvious (Fig. 6). Turbulent flow is more coherent and is of larger horizontal scale when the wind blows westward to northward (; i.e., from a westward wind to a northward wind in a clockwise sense) than when , while turbulence is more intermittent and is of smaller spatial space when the wind blows eastward to southward ().

Under the influence of negative *f*, mean surface current (, ) is deflected to the left of the wind and rotates counterclockwise with depth in the Southern Hemisphere (Ekman 1905). The mean horizontal velocity has a different wind direction dependence at 15°S (Fig. 7a) than at 15°N: the mean surface current is deflected more to the left of the wind than without the influence of , indicated by a larger angle between wind stress and surface current when the wind blows eastward to southward (), and is deflected less when the wind blows westward to northward ().

Different from 15°N, is larger for westward to northward wind () than when at 15°S (Fig. 7b). In particular, the boundary layer is the shallowest for a southward wind () and the deepest for a northward wind (). Similar to 15°N, the maximum is twice as large as the minimum at 15°S.

Similar to the variation in , the values of (Fig. 7c), (Fig. 7d), and (Fig. 7f) are larger than the case when the wind blows westward to northward and smaller than the case when the wind blows eastward to southward. It is obvious that increases for cases , compared with the case, and decreases for cases (Fig. 7e).

The magnitude of the mean surface current decreases at 45°S (Fig. 8a) due to larger magnitude of *f*, compared to at 15°S. The variabilities in , , , , and to wind direction also reduce as latitude increases in the Southern Hemisphere from 113%, 39%, 66%, 96%, and 77% at 15°S (Figs. 7b–f) to 52%, 14%, 39%, 44%, and 26% at 45°S (Figs. 8b–f), respectively. The hemisphere dependence of the effect is associated with the signs of , , , and as well: turbulence is enhanced when the wind has a westward and southward component () in the Northern Hemisphere and strengthened when the wind has a westward and northward wind () in the Southern Hemisphere, both of which have positive values for both and (within the boundary layer, is negative in both hemispheres, and is positive in the Northern Hemisphere and negative in the Southern Hemisphere). Turbulence is reduced when the wind has an eastward and northward component () in the Northern Hemisphere and is weakened when the wind has an eastward and southward component () in the Southern Hemisphere, both of which have negative values for both and .

### c. Modification of the effect due to stratification

The mixed layer is generally shallow and is capped by stable stratification at low to midlatitudes. The modification of the effect on Ekman turbulence by stable stratification is studied with experiments and (Table 1). Similar to unstratified Ekman turbulence (Fig. 2), the horizontal scale of convergence regions for a southwestward wind is larger than for a northeastward wind inside the boundary layer (not shown), although the difference is less obvious than in the unstratified Ekman turbulence. Similar to the unstratified Ekman turbulence, the effect in the vertical velocity fluctuation is negligible at the surface in stratified Ekman turbulence. In the lower part of the boundary layer, the wind direction dependence of flow structure is insignificant because of the reduced .

Compared with unstratified Ekman cases (Cases and *E*), Fig. 9a shows that the variability in the mean flow to wind direction is reduced by stratification. The stratification reduces the difference in the angle *β* between surface mean current and wind stress significantly (Fig. 9b). The difference in *β* between maximum and minimum reduces from 12.6° in unstratified Ekman cases to 3.3° in stratified Ekman cases. The variability in mean surface current is similar with and without stratification: the largest mean surface speed for an eastward wind () is approximately 35% larger than the smallest mean surface speed for a westward wind ().

Stratification reduces the variability in vertical velocity variance (Fig. 9c): the largest is only 14% larger than the smallest. The variability in reduces from 79% in unstratified Ekman turbulence (Fig. 3d) to 25% in stratified Ekman turbulence (Fig. 9d). Stratification reduces the magnitude and the variability of to wind direction (Fig. 9e). The variability in decreases from 98% in unstratified Ekman turbulence to 37% in stratified Ekman turbulence. The variability in to wind direction in stratified Ekman turbulence (Fig. 9f) also reduces and is different from unstratified Ekman turbulence: the variability in in stratified Ekman turbulence is only 15%, which is much smaller than the variation (76%) in unstratified Ekman turbulence. The magnitude of for the case with stratification reduces by 45%, compared with the case without stratification. The dependence of wind work, as well as TKE dissipation rate (not shown), on the wind direction is enhanced in the presence of stratification. The largest wind work when the wind blows eastward in stratified Ekman turbulence is 36% larger than the smallest when the wind blows westward, compared with the variability (30%) in unstratified Ekman turbulence.

The average entrainment rate is calculated as (Sullivan et al. 1998), where is the depth where the potential temperature is 0.25 K smaller than the surface value. For stratified cases, the boundary layer turbulence is limited by a shallow mixed layer depth , which is the temporal and horizontal average of . The variability in to wind direction during the relatively short integration is small (not shown). In unstratified Ekman turbulence, is strongly dependent on wind direction. In the regions with a shallow mixed layer, is limited by the stratification and is equal to (about −45 m at 15°N). The largest entrainment rate ( m day^{−1}) is about 116% of the smallest ( m day^{−1}; Fig. 9g). The stronger entrainment for southward to westward winds () also explains the reduced for these wind directions (Fig. 9d). When entrainment strengthens, horizontal velocity variances increase within the entrainment layer (e.g., Wang et al. 1996), leading to a lower vertical to horizontal velocity variance ratio.

The variability in velocity variance to wind direction is reduced by stratification. For stratified Ekman turbulence with , reduces by 4% at the surface, compared with the unstratified Ekman case. The maximum variance of vertical velocity fluctuation decreases by 15% due to the limitation of stratification on turbulent mixing. When is considered, the difference reduces by 13% when the wind is northeastward and increases by 22% when the wind is southwestward at the surface. The ordering of velocity variances is above for all the wind directions. Both and are positive in the whole boundary layer, while decreases rapidly below (Fig. 10). The large and near the base of the mixed layer () when the wind blows southward to westward indicate that the horizontal turbulence is strengthened when the entrainment rate is strong (e.g., Kukulka et al. 2010).

Stratification also decreases both the magnitude of eddy turnover time and its variability to wind direction (not shown). The eddy turnover time in the case is decreased by 41% due to stratification. When is considered, the variability in eddy turnover time reduces from 86% in the unstratified Ekman turbulence to 52% in the stratified Ekman turbulence. For unstratified cases, larger represents more time for turbulence dissipating energy and, hence, larger TKE contained in the water column, and vice versa. For stratified cases, stratification rather than eddy turnover time may have more impacts on TKE, which means TKE is more likely to be determined by the depth that turbulence can reach than eddy turnover time in the presence of stratification as production changes less.

### d. Modification of the effect by Langmuir turbulence

Nonbreaking surface gravity waves have two primary effects because of phase-averaged Stokes drift: the mean horizontal current profile is modified by Coriolis–Stokes forcing (e.g., Polton et al. 2005), and turbulent kinetic energy is enhanced through Stokes production (e.g., McWilliams et al. 1997). Also, wave–current interactions generate coherent streamwise vortices, which are signatures of Langmuir turbulence (e.g., Belcher et al. 2012). We first investigate the impact of on OSBL turbulence for wind-following waves. In this regime, Stokes drift tilts vertical vorticity perturbation into horizontal vortices that align with the wind stress. Cases and *L* represent simulations driven by wind-following, nonbreaking surface gravity waves.

For unstratified Langmuir turbulence and , the downwelling branches of the LCs are oriented roughly downwind with relatively small horizontal scales near the surface; LCs induce streamwise streaky patterns in the vertical velocity (not shown) and narrow forward-looking Y junctions (Farmer and Li 1995; Sullivan et al. 2007). With increasing depth, the streamwise streaks are more coherent and rotate clockwise, consistent with previous studies (e.g., McWilliams et al. 1997). In the lower part of the boundary layer (e.g., ), the streak-like patterns are more fragmented. When is considered, there are only slight differences at the surface for varying wind directions. With increasing depth, the streamwise streaks become more coherent and orient obliquely to the wind direction when the wind blows southward to westward. In contrast, the streak-like pattern becomes less organized, with mostly small scales when the wind blows northward to eastward. At the same time, the asymmetry between stronger downwelling and weaker upwelling branches of the LCs is enhanced for a southwestward wind and is reduced for a northeastward wind. The profile of conditional averaged vertical velocity (McWilliams et al. 1997; Van Roekel et al. 2012) indicates that has a minimal effect on the coherent structure of Langmuir turbulence close to the surface (not shown).

The amplitude of the Eulerian current decreases dramatically due to wind-following waves (inset of Fig. 11a), consistent with the Stokes–Ekman theory (e.g., McWilliams et al. 1997). However, the difference in the surface Lagrangian current is slight for cases , compared with the case. The hodographs with deviate from the case, compared with cases (Fig. 11a). The difference of maximum and minimum in the angle between surface Eulerian mean current and wind stress is 19.3°, while the difference is 9.5° in the angle between surface Lagrangian mean current and wind stress (not shown).

The boundary layer is deeper in Langmuir turbulence than in Ekman turbulence (Fig. 3b), consistent with previous studies (e.g., McWilliams et al. 2012). The difference in between Ekman and Langmuir simulations for the same wind direction ranges from 0 to 35 m. The largest is 87% larger than the smallest (Fig. 11b). The largest is 21% larger than the smallest (Fig. 11c). In contrast, the largest is 83% larger than the smallest (Fig. 11d). The variabilities in and reduce to 81% and 46%, respectively (Figs. 11e,f).

Shear production decreases by 73% in unstratified Langmuir cases with , compared to their Ekman counterpart, due to the reduced shear of the Eulerian mean flow. The largest shear production is about 138% that of the smallest. The reduction in mean velocity (shear) is also due to the Stokes–Coriolis forcing. Although Stokes production () can be nearly 3 times as large as shear production, the differences of Stokes production among different wind directions are negligible. The variability in TKE production with wind direction is mainly determined by Eulerian shear production in unstratified Langmuir cases (Fig. 12). Wind-following waves increase the magnitude of but lower its variability to wind direction (not shown).

### e. Modification of the effect by swell

Waves and their associated Stokes drifts align with surface wind stress in the previous subsection. This is also the most common regime adopted in many earlier theoretical (Craik and Leibovich 1976; Leibovich 1983) and LES studies of Langmuir turbulence, assuming waves in equilibrium with the steady wind (Skyllingstad and Denbo 1995; McWilliams et al. 1997; Harcourt and D’Asaro 2008). In a realistic ocean, waves seldom align with wind because of the influence of swell and the temporal variability of wind. By analyzing LES solutions driven by systematically varying wind–wave misalignment, Van Roekel et al. (2012) showed that a projected Langmuir number can be used to quantify the strength of Langmuir turbulence with wind–wave misalignment and concluded that as the angle between wind and wave increases, the upwind cell in a Langmuir cell vortex pair becomes stronger than the downwind cell. Recent studies have focused on the relative misalignment of winds and waves, which is critically important in determining the strength of mixing by Langmuir turbulence (Webb and Fox-Kemper 2011; Belcher et al. 2012; Van Roekel et al. 2012). In the regime of wind-opposing swell, Stokes production is negative, leading to smaller turbulence variances (Sullivan et al. 2012) than in Langmuir turbulence. In this study, we consider only two idealized scenarios: the case () when swell is twice as strong as the equilibrium wave and propagates in the same direction as the wind and the case () when swell is as strong as equilibrium wave and propagates in the opposite direction as the wind, although wind and wave can be at any angle.

The effect of in Case (Fig. 13) is more similar to Case (section 4d). There is a large decrease in the amplitude of Eulerian mean current in both Case (inset of Fig. 11a) and Case (inset of Fig. 13a). The difference of maximum and minimum in the angle between surface Eulerian mean current and wind stress is 18.4°, while the difference is 4.4° in the angle between surface Lagrangian mean current and wind stress (not shown). The largest boundary layer depth is 128 m when the wind blows westward, while the smallest depth is 71 m when the wind blows northward (Fig. 13b). The relative difference between the largest and the smallest is similar to Case . The largest is 16% larger than the smallest (Fig. 13c). The variabilities in , , and are 80%, 82%, and 47%, respectively (Figs. 13d–f).

The effect of on the pattern of surface vertical velocity fluctuation in Case is negligible as well. Away from the surface (), turbulence is more coherent and is of larger spatial scale when the wind blows southward to westward than when the wind blows in other directions. The coherent structures also orient obliquely to wind (not shown), similar to Case .

Compared with the hodograph of the Eulerian mean current with in Case (inset of Fig. 11a), the Eulerian surface mean current with increases due to the negative Stokes drift in the opposite direction of wind in Case (inset of Fig. 14a). The difference of maximum and minimum in the angle between surface Eulerian mean current and wind stress is 7.9°, while the difference is 10.4° in the angle between surface Lagrangian mean current and wind stress (not shown). In comparison with the case driven by wind-following waves (Case *L*), the boundary layer in Case *M* () is 30 m shallower due to weaker vertical mixing associated with the negative Stokes drift shear (Fig. 14b). The largest is 62% larger than the smallest (Fig. 14c). The instantaneous snapshot indicates that the vertical turbulent mixing is stronger when the wind blows southward to westward, represented by coherent streaks of downward velocity (not shown) and deeper boundary layer (not shown). The largest for a southward wind is 99% larger than the smallest for a northward wind (Fig. 14d). The largest is 46% larger than the smallest (Fig. 14e). Without the effect of , is 31% weaker when a wave propagates in the opposite direction of the wind (Fig. 14f) than when a wave is aligned with wind. Under the influence of , the largest is 87% larger than the smallest. The effect of when swell is opposite to the wind with equilibrium waves (Fig. 14) is more similar to unstratified Ekman turbulence (section 4a).

### f. Modification of effect due to both stratification and waves

Under the simultaneous influence of stratification and nonbreaking surface gravity waves, the variabilities of the mean flow (Fig. 15a) and turbulent statistics to wind direction are small, compared to all other previously discussed cases (Cases , , , , and ). Wind direction has little effect on the instantaneous flow structure in the stratified Langmuir turbulence since both stratification and Langmuir circulation reduce the variability of turbulent kinetic energy to wind direction (not shown). The difference between the maximum and minimum in the angle between surface Eulerian mean current and wind stress is 2.8°, while the difference is 2.1° for the angle between surface Lagrangian mean current and wind stress.

The magnitude of maximum entrainment rate is 14% larger in stratified Langmuir turbulence (Fig. 15b) than in stratified Ekman turbulence; however, the ratio of maximum to minimum entrainment rate is similar to that of the stratified Ekman turbulence (Case ). Equilibrium wave increases by only 3% when (dashed lines in Figs. 9g, 15b), while the variability due to wind direction is similar with and without wave forcing. The boundary layer depth is limited by stratification, similar as the result in section 4c. The layer-averaged ranges from 0.68 to 0.74 (Fig. 15c), consistent with Harcourt and D’Asaro (2008) and Van Roekel et al. (2012) for simulations with similar parameter space. The wind direction dependence of (Fig. 15d), (Fig. 15e), and (Fig. 15f) is reduced and even negligible.

The effect of with stratification and wind-opposing swell (Case ) is similar to the case with stratified Ekman turbulence (Case ; see section 4c), and the effect with stratification and wind-following swell (Case ) is similar to stratified Langmuir turbulence (Case ).

## 5. Effect of on eddy viscosity profiles

The effect of turbulence is represented by eddy viscosity *K* in predictive ocean and climate models. The LES solutions are used to study the effect of on *K* in this section. Eddy viscosity is diagnosed as (McWilliams et al. 2012)

where is the Reynolds stress, and is the Lagrangian mean velocity. The subscript denotes the horizontal components. The effects of on eddy viscosity profiles in both stratified Ekman and Langmuir cases (Cases and ) at 15°N are shown in Fig. 16. Here, eddy viscosity *K* is normalized by ( 45 m). The impact of on the *K* profiles is significant, and the value of *K* changes with wind direction.

For stratified Ekman turbulence, the maximum value of eddy viscosity changes with wind direction (Fig. 16a). The value of is the largest in stratified Ekman turbulence when the wind blows westward (), while it is the smallest when the wind blows northeastward (). The shape of the *K* profile and the depth where are similar for different wind directions (). The eddy viscosity reduces to zero at the lower 10% of the mixed layer. The eddy viscosity magnitudes when are smaller than when , while those of the opposite wind directions () are larger than when .

There are two major effects of on *K* in stratified Langmuir cases. One is the amplification or reduction of the magnitude of : the maximum eddy viscosity decreases for and increases for , compared with the case. The other is deepening or shoaling the location where . In stratified Ekman cases, however, the latter effect is negligible.

For stratified Langmuir turbulence, the magnitude of *K* increases as surface gravity waves enhance turbulent mixing in the mixed layer (Fig. 16b). The value of is the largest for westward wind (), approximately 150% as large as the smallest for eastward wind (). The location of is deeper when is larger, and vice versa. The variation of eddy viscosity profiles with respect to wind direction mainly appears in the middle of the mixed layer (). The difference in *K* profiles is negligibly small (~6%) close to the surface in the stratified Langmuir case, while the largest value of *K* for a westward wind is twice as large as the smallest value for a northeastward wind in the stratified Ekman case in the upper 10% of the mixed layer. Therefore, the effects of on turbulent mixing are mainly confined to the center of the mixed layer in stratified Langmuir turbulence. The differences of viscosity with respect to wind direction become more pronounced in stratified Langmuir turbulence, compared to stratified Ekman turbulence.

In the absence of stratification, the *K* profiles deviate from the convex shape (not shown), consistent with McWilliams et al. (2012). However, there is still a wind direction dependence of eddy viscosity in both unstratified Ekman and Langmuir turbulence, similar to the stratified case (not shown).

A complete mixing parameterization includes boundary layer depth and entrainment rate , in addition to eddy viscosity *K*. Results in section 4 (Figs. 3, 9, 15) show that the two quantities ( and ) also change with wind direction when is considered. The evaluation of the existing parameterization for eddy viscosities and the inclusion of the effect in mixing parameterizations are important in improving the performance of parameterizing turbulent mixing in ocean models and will be discussed in a future work.

## 6. Summary and discussion

In this study, we use large-eddy simulations to investigate how latitude/hemisphere, surface gravity waves, and stratification influence the effect of on OSBL turbulence. In the absence of stable stratification and nonbreaking waves, both the mean and turbulent flows change substantially with wind direction and latitude under the influence of , consistent with previous studies (Coleman et al. 1990; Zikanov et al. 2003). Horizontal rotation changes the flow structures and redistributes TKE among three spatial directions. The wind direction dependence of vertically integrated TKE over boundary layer () is mainly determined by the eddy turnover time (or dissipation length scale) instead of TKE production. The variability of turbulent flows to wind direction is different at the same latitude in opposite hemispheres and is also different at different latitudes in the same hemisphere. In the Northern Hemisphere, turbulence is stronger when the wind has a westward and southward component than when the wind has an eastward and northward component. In the Southern Hemisphere, turbulence is stronger when the wind has a westward and northward component than when the wind has an eastward and southward component. When the boundary layer is capped by stable stratification, alters the entrainment rate and also redistributes TKE among three spatial directions. Both the entrainment rate and the magnitude of eddy viscosity show strong variability with wind direction in stratified Ekman turbulence at 15°N.

Both locally wind-driven wave and swells reduce the wind direction dependence of mean current. However, there are also strong variabilities in both the magnitude and shape of eddy viscosity in the presence of stratification and Langmuir circulations. Wind-opposing swell decreases the wind direction dependence of mean velocity but increases the wind direction dependence of .

While has its largest value and effect on turbulence in the tropics, Wang et al. (1996) found that the entraining boundary layer, which is the focus of this study, lacks an equilibrium solution at the equator where *f* is zero unless large-scale equatorial circulation is included. Furthermore, Wang (2006) also showed that has little effect on convective turbulence. Therefore, convective turbulence is not discussed in this study.

Accurate simulation of sea surface temperature (SST) and MLD is important to the fidelity of Earth system models and remains one of the major challenges. Recent efforts aiming to reduce existing biases in SST and MLD have focused mainly on the inclusion of parameterization for submesoscale processes (Fox-Kemper et al. 2008; Fox-Kemper and Ferrari 2008; Fox-Kemper et al. 2011) and on the improvement of parameterization for boundary layer mixing by including wave effects (e.g., Qiao et al. 2004; Harcourt 2013, 2015; Fan et al. 2014; Li et al. 2016; Li and Fox-Kemper 2017). These efforts have achieved success to different degrees, but have not fully resolved the biases. The current study demonstrates that the effect of wind direction leads to variability in mixing and associated SST and MLD. For example, in trade wind regions, the northeast wind in the Northern Hemisphere and southeast wind in the Southern Hemisphere drive stronger vertical mixing with than without . In the Bay of Bengal and the Arabian Sea, influenced by Indian monsoon, northeast wind in the winter drives stronger mixing than southwest wind in the summer does if hydrographic conditions are identical. The difference between the minimum and maximum is 0.22 in unstratified Ekman turbulence, while increases from 0.55 in unstratified Ekman turbulence to 0.76 in unstratified Langmuir turbulence at 15°N when , indicating that the variability in mixing associated with wind direction is comparable to that due to waves. Normalized maximum eddy viscosity increases from [0.041, 0.087] to [0.129, 0.190] when equilibrium wave is added, while the variabilities in are 110% and 48% for stratified Ekman and Langmuir turbulence, respectively (Figs. 16a,b). The entrainment rate is only 3% larger in wave-driven Langmuir turbulence than in Ekman turbulence. The difference in caused by equilibrium wave forcing is substantially smaller than the variability of