Influence of Coriolis Parameter Variation on Langmuir Turbulence in the Ocean Upper Mixed Layer with Large Eddy Simulation

The ocean upper mixed layer, caused by turbulent mixing, plays an important role in the weather and climate forecasts because it is the link between the oceans and atmosphere and directly controls the air–sea exchange of momentum, heat, and gases. Langmuir turbulence (LT), produced by wave–current interactions (Craik and Leibovich, 1976; Leibovich, 1983), is one of the key processes in creating and maintaining the upper mixed layer over the global ocean (Langmuir, 1938; Smith, 1998; McWilliams and Sullivan, 2000; de Boyer et al., 2004; Thorpe, 2004; Belcher et al., 2012; Liang et al., 2012; Kukulka et al., 2013; Wijesekera, et al., 2013; Basovich, 2014; D’Asaro, 2014; McWilliams et al., 2014; Yang et al., 2014; Noh et al., 2016; Smyth et al., 2017).

A significant number of studies on LT characteristics have been conducted using the large eddy simulation model, pioneered by Skyllingstad and Denbo (1995), McWilliams et al. (1997), and Noh et al. (2004). The main goals of LT studies are to reveal the complex characteristics of LT and develop an accurate LT parameterization (McWilliams and Sullivan, 2000; Smyth et al., 2002; Grant and Belcher, 2009; Noh et al., 2011, 2016; Harcourt, 2013; Sutherland et al., 2014; Furuichi and Hibiya, 2015; Harcourt, 2015; Li et al., 2017) that can be used to interpret the effects of LT on the upper mixed layer and improve the prediction ability of the oceanic general circulation models (OGCMs) for the upper mixed layer variation (Noh et al., 2009, 2016; Belcher et al., 2012; Fan and Griffies, 2014). Previous studies on LT have focused mainly on the variations in the wind and wave fields, including when the wind and wave fields are aligned and balanced (Skyllingstad and Denbo, 1995; McWilliams et al., 1997; Sullivan et al., 2007; Grant and Belcher, 2009), when the wind and wave fields are misaligned (van Roekel et al., 2012; McWilliams et al., 2014; Xuan et al., 2019, 2020) and unbalanced (Noh et al., 2004, 2010; Li et al., 2005; Polton et al., 2005; Polton and Belcher, 2007; Grant and Belcher, 2009; Sutherland et al., 2014; McWilliams et al., 2014; Pearson et al., 2015), how the LT characteristics vary in special wind event conditions (Skyllingstad et al., 1999; Li et al., 2009; Kukulka et al., 2010; Sullivan et al., 2012; Hoecker-Martínez et al., 2016), and how the wave breaking influences LT characteristics (Noh et al., 2004; Sullivan et al., 2007). These previous studies have provided a substantial amount of important information for improving LT parameterization.

Furthermore, influence of the Coriolis parameter ( f ) on LT characteristics has also been recognized [e.g., McWilliams et al. (1997) and Polton et al. (2005)]. The Coriolis–Stokes force [${ - f}{{u}_{{\rm{s}}} }$ , where ${{u}_{\rm{s}}}$ is the Stokes drift velocity], through the anti-Stokes Eulerian flow (${u} \approx { - }{{u}_{{\rm{s}}} }$) opposing the Stokes drift velocity, can modify the shape of the vertical profile of the horizontal velocity and the relative magnitude between the horizontal velocity components over the upper mixed layer depth, when the wind forcing remains unchanged and the enlarged Stokes drift velocity is set in the different ideal simulation cases (Polton et al., 2005; McWilliams et al., 2014). The entrainment depth is proportional to the Ekman depth scale (${{{{{u}}_{{*}}}}/ {{f}}}$, where ${{{u}}_{{*}}}$ is the sea surface friction velocity induced by wind) when the thermocline appears below the upper mixed layer (Noh and Choi, 2018), and the penetration depth of LT can reach the Ekman depth scale in an unstratified mixed layer (Polton and Belcher, 2007). The magnitude and curvature of the vertical momentum flux are associated with the ratio of the upper mixed layer depth (h) to the Ekman depth scale (e.g., ${{{{fh}}}/ {{{{u}}_{{*}}}}}$) (Grant and Belcher, 2009). However, previous studies on the impact of Coriolis parameter on LT characteristics were often based on a fixed value of Coriolis parameter that focuses mainly on the middle latitudes (Polton et al., 2008; Grant and Belcher, 2009; Noh et al., 2011; McWilliams et al., 2012, 2014; Li and Fox-Kemper, 2017). Insufficient variations of the Coriolis parameter were used to display the influence of the Coriolis parameter on LT characteristics (Min and Noh, 2004; Polton and Belcher, 2007; Noh et al., 2010; van Roekel et al., 2012; Pearson et al., 2015). These previous studies illustrate that the impact of Coriolis parameter on LT characteristics is significantly important, but it has not been investigated systematically how LT characteristics change as latitude varies from low to high. At present, the variation of LT characteristics with change of latitude has not been studied. This may be an important reason that the upper mixed layer variations simulated by the OGCMs including LT parameterization are significantly different with the observations in the high latitudes (Fan and Griffies, 2014; Noh et al., 2016; Li and Fox-Kemper, 2017). Hence, in this study, we focus mainly on how LT characteristics vary with change of latitude.

The remainder of this paper is organized as follows. In section 2, the large eddy simulation model and ideal simulation cases are described. The results of the mean and turbulent statistics of LT and comparison with previous studies are presented in section 3. The dynamic analysis and discussion about influence mechanism of the Coriolis parameter on the mean and turbulent statistics of LT are given in section 4. Section 5 presents a summary of the findings.

The large eddy simulation model used in this research solves incompressible, rotating Boussinesq, and filtered Craik–Leibovich equations developed by McWilliams et al. (1997), Sullivan et al. (2007), Sullivan et al. (2012), and McWilliams et al. (2014). Turbulent sub-grid-scale fluxes are calculated by a modified prognostic turbulent kinetic energy equation (McWilliams et al., 1997). The Craik–Leibovich momentum equations compute LT dynamics by a vortex force (${{u}_{\rm{s}}} {\omega }$) (Craik, 1977; Leibovich, 1977, 1983) that contains the Stokes drift velocity (${{u}_{\rm{s}}}$) resultant from surface wave fields and the vertical vorticity induced by horizontal current (${\omega }{\text{ = }}\nabla {{u}_{{{\rm{hc}}}}}$, ${{u}_{{{\rm{hc}}}}}$ is the horizontal current) generated by the wind force. The details of the governing equations are provided in the appendix.

The sea surface friction velocity (${{u}_{{*}}}$) is calculated based on the steady wind velocity (Liu et al., 1979) given by

where ${{{{\tau}}}_{{\rm{a}}}}$ is the wind stress at the sea surface, ${\rho _{{\rm{o}}}}$ is the sea water density, ${\rho _{{\rm{a}}}}$ is the air density, ${{U}_{{\rm{a}}}}$ is the wind velocity (z = 10 m), and ${{C}_{{\rm{d}}}}$ = 10⁻³ (0.79+0.0509${{U}_{{\rm{a}}}}$) is the drag coefficient for ${{U}_{{\rm{a}}}}$ ≥ 10 m s⁻¹.

The Stokes drift velocity (${{u}_{\rm{s}}}$) based on a full wave spectrum for a steady wind (Kenyon, 1969; McWilliams and Restrepo, 1999; Sullivan et al., 2007; Li et al., 2017) is computed by

where ${F}\left( {\sigma } \right)\, =\, \left( {{{{{a}_{{\rm{n}}}}}/ {{{\sigma }^5}}}} \right)\,\exp \left[ { \,- {{b}_{{\rm{n}}}}{{\left( {{{g} /{{{U}_{{\rm{a}}}}{\sigma }}}} \right)}^{n}}} \right]$, ${{a}_{{\rm{n}}}} = \left( {{{{{f}_{{\rm{o}}}}} / {2\pi }}} \right) \times {\left( {2{\pi }{{\nu }_{{\rm{o}}}}} \right)^{n}}$, ${{b}_{{\rm{n}}}} = \left( {{5/ {n}}} \right){\left( {2{\pi }{{\nu }_{{\rm{o}}}}} \right)^{n}}$, ${{f}_{{\rm{o}}}}$= 0.275 and ${{\nu }_{{\rm{o}}}}$=1.40 (Kenyon, 1969; Sullivan et al., 2007), and n = 2 (McWilliams and Restrepo, 1999).

In order to illustrate variations in LT characteristics with a change of latitude (LAT), we simulate the ideal cases of LAT = 20°N, 30°N, 40°N, 50°N, 60°N, 70°N, and 80°N, corresponding to Coriolis parameter f = 4.988 × 10⁻⁵ s⁻¹, 7.292 × 10⁻⁵ s⁻¹, 9.3744 × 10⁻⁵ s⁻¹, 1.1172 × 10⁻⁴ s⁻¹, 1.2630 × 10⁻⁴ s⁻¹, 1.3704 × 10⁻⁴ s⁻¹, and 1.4362 × 10⁻⁴ s⁻¹, and inertial period t_ip= 2π/f = 35.0 h, 23.9 h, 18.6 h, 15.6 h, 13.8 h, 12.7 h, and 12.2 h, respectively. The initial upper mixed layer depth is h = 45 m, and neutral stratification is assumed from the sea surface to the upper mixed layer base (h). The thermocline is stably stratified at a rate of ${{{{\rm{d}}}{T}} /{{{\rm{d}}}{z}}}$= 0.2 K m⁻¹; the strong stratification in the thermocline (McWilliams et al, 2014) clearly shows the variation of the upper mixed layer depth and sea surface temperature with latitude. The sea surface friction velocity is $ {u_*} $= 0.0189 m s⁻¹, which corresponds to the wind velocity ${{U}_{{\rm{a}}}}$ = 15.1 m s⁻¹ (z = 10 m). The sea surface Stokes drift velocity ${{u}_{{{\rm{os}}}}}$= 0.21 m s⁻¹. Hence, the turbulence Langmuir number ${{\rm{L}}}{{{\rm{a}}}_{{\rm{t}}}} = \sqrt {{{{{u}_{*}}} \mathord{\left/ {\vphantom {{{{u}_{*}}} {{{u}_{{os}}}}}} \right. } {{{u}_{{{\rm{os}}}}}}}}$= 0.30 (McWilliams et al., 1997) in the LT regime (Li et al., 2005; Sutherland et al., 2014). The directions of the wind and wave fields are aligned in the all-idealized simulation cases. The sea surface friction velocity (${{{u}}_{{*}}}$) and the upper mixed layer depth (h) are used to normalize the LT simulation results. To normalize the entrainment flux, a constant, small heat flux of ${{{Q}}_{{*}}}$ = 1.2 × 10⁻⁶ K m s⁻¹ is imposed on the sea surface (the positive heat flux is the heat flux out of the ocean), which implies a Monin–Obukhov length of L = –7226 m. Hence, h/L = –0.0062 is very small, and the upper mixed layer is in a very weak convection turbulence region, but in a significantly strong Langmuir turbulence region.

The computational domain is 256 m × 256 m in the horizontal directions. The depth of the domain is 128 m for all simulation cases. The number of grid points is 512 × 512 in the horizontal directions, and the corresponding grid spacing is 0.5 m. The number of grid points is 512 in the vertical direction, and the corresponding grid spacing is 0.25 m in the vertical direction. The flux- and energy-carrying large-scale eddies in the upper mixed layer can be adequately resolved by this spatial resolution (McWilliams et al., 1997, 2014; Noh et al., 2004; Sullivan et al., 2007; Grant and Belcher, 2009; Li et al., 2009; Li and Fox-Kemper, 2017). The lateral boundary condition is periodic, the bottom boundary condition is radiation and free slip due to the simulated domain being located in the open and deep ocean, and the top boundary condition is free slip with the imposed wind shear.

The initial fields are static for all simulation cases (Grant and Belcher, 2009; Li and Fox-Kemper, 2017). The mean and turbulent statistics presented in this paper are obtained by the spatial average in the x–y plane for every integral time step except for the vertical velocity fields (McWilliams et al., 1997; Sullivan et al., 2007). The solution reaches a statistical equilibrium state after 25 h, as found by Sullivan et al. (2007). The run time of all our simulation cases is t_tot = 74.4 h. The mean and turbulent statistics presented in this paper are averaged over 35 h ≤ t ≤ 70 h to smooth the relative profiles, as suggested by McWilliams et al. (2014). The horizontally and temporally averaged values of the mean and turbulent statistics presented in this paper are denoted by the angle bracket $\left\langle \varphi \right\rangle $, where $\varphi $ represents an arbitrary physical quantity.

The vertical velocity fields can directly impact the vertical transport of mass and energy in the upper mixed layer and the variation of the upper mixed layer depth (Basovich, 2014; McWilliams et al., 2014). Figure 1 shows a snapshot of the three-dimensional vertical velocity field of LT. It is characterized by the coherent structure of Langmuir cells and many Y-joints, which are consistent with the observations in the real ocean (Thorpe, 2004; Yang et al., 2014) and the computations using large eddy simulation in previous studies (Noh et al., 2004; Sullivan et al., 2012; Sullivan and McWilliams, 2019; Xuan et al., 2019). Moreover, the downwelling jets are narrow and strong compared to the upwelling velocities, consistent with the previous studies of Sullivan et al. (2007) and McWilliams et al. (2014). The penetration depth of the strong downwelling jets for LAT = 20°N (Fig. 1a) is much deeper than that for LAT = 80°N (Fig. 1b). This is because the penetration depth of the strong downwelling jets is proportional to the Ekman depth scale (${{{{{u}}_{{*}}}}/ {{f}}}$) (Polton and Belcher, 2007; Noh and Choi, 2018).

Figure 1.  Snapshot of three-dimensional vertical velocity (w) fields for (a) latitude = 20°N and for (b) latitude = 80°N.

Figure 2a shows variations in the temperature of the upper mixed layer and the structure of the thermocline inversion. Figure 2b indicates a change in the entrainment flux ($\left\langle {{w'T'}} \right\rangle $, ${w'}$ represents the vertical fluction velociy). The increase in the entrainment flux at the smaller ${{{{fh}}} /{{{{u}}_{{*}}}}}$ (lower latitudes) (Fig. 2b) causes a stronger cooling of the upper mixed layer (Fig. 2a). This can also be observed by the flat temperature gradient in the entrainment layer (Fig. 2a). Hence, these results suggest that a variation in Coriolis parameter can modulate the entrainment flux and the upper mixed layer temperature. The decrease in the entrainment flux with ${{{{fh}}} /{{{{u}}_{{*}}}}}$ (latitude) changing from 0.119 (20°N) to 0.342 (80°N) (Fig. 2b) is partly due to the fact that the Coriolis parameter may also modify the gradient of horizontal current (shown later in Fig. 4) near the upper mixed layer base, that is, the velocity shear can also affect the entrainment flux (Sullivan et al., 2007). Furthermore, the positive entrainment flux near the surface layer represents the cold-water sink.

Figure 2.  Vertical profiles of the normalized (a) temperature ($\left\langle {T} \right\rangle $) (For reference, the initial mixed-layer sounding is shown as a black fine line. The reference temperature T_r = 302 K) and (b) entrainment flux ($\left\langle {{w'T'}} \right\rangle $) with a variation of ${{{fh}} \mathord{\left/ {\vphantom {{{fh}} {{{u}_{*}}}}} \right. } {{{u}_{*}}}}$.

Figure 4.  Vertical profiles of the normalized (a) downwind Lagrange velocity ($\left\langle {u} \right\rangle + {{u}_{{\rm{s}}}}$) (the vertical solid line indicates zero value) and (b) crosswind velocity ($\left\langle {v} \right\rangle $) with a variation of ${{{fh}} \mathord{\left/ {\vphantom {{{fh}} {{{u}_{*}}}}} \right. } {{{u}_{*}}}}$.

The changes of the upper mixed layer and the entrainment flux with a variation of $ {{fh} \mathord{\left/ {\vphantom {{fh} {{u_*}}}} \right. } {{u_*}}} $ are illustrated in this paragraph. Figure 3 shows the changes in the sea surface temperature ($\left\langle {{{\rm{SST}}}} \right\rangle$), the upper mixed layer depth ($\left\langle {{{h}_{{{\rm{LAT}}}}}} \right\rangle$), and the entrainment flux extremum [min ($\left\langle {{w'T'}} \right\rangle $)] in the thermocline, as suggested by previous studies (Grant and Belcher, 2009; McWilliams et al., 2014; Li and Fox-Kemper, 2017). A comparison of the variations between the $\left\langle {{{\rm{SST}}}} \right\rangle $ (Fig. 3a), the $\left\langle {{{h}_{{{\rm{LAT}}}}}} \right\rangle$ (Fig. 3b), and the min ($\left\langle {{w'T'}} \right\rangle $) (Fig. 3c) indicates that their variations are evident for ${{{{fh}}} \mathord{\left/ {\vphantom {{{{fh}}} {{{{u}}_{{*}}}}}} \right. } {{{{u}}_{{*}}}}}$≤ 0.266 (LAT ≤ 50°N) but relatively weak for ${{{{fh}}} \mathord{\left/ {\vphantom {{{{fh}}} {{{{u}}_{{*}}}}}} \right. } {{{{u}}_{{*}}}}}$≥ 0.266 (LAT ≥ 50°N). This result indicates that ${{{{fh}}} \mathord{\left/ {\vphantom {{{{fh}}} {{{{u}}_{{*}}}}}} \right. } {{{{u}}_{{*}}}}}$= 0.266 (LAT = 50°N) is a key value (latitude), that is, the rate of change in the upper mixed layer quantities caused by Coriolis parameter (i.e., Ekman effect) is clearly different between${{{{fh}}} \mathord{\left/ {\vphantom {{{{fh}}} {{{{u}}_{{*}}}}}} \right. } {{{{u}}_{{*}}}}}$≤ 0.266 (LAT ≤ 50°N) and ${{{{fh}}} \mathord{\left/ {\vphantom {{{{fh}}} {{{{u}}_{{*}}}}}} \right. } {{{{u}}_{{*}}}}}$≥ 0.266 (LAT ≥ 50°N).

Figure 3.  Variation of the normalized (a) sea surface temperature ($\left\langle {{{\rm{SST}}}} \right\rangle$) (the reference temperature T_r = 302 K), (b) upper mixed layer depth ($\left\langle {{{h}_{{{\rm{LAT}}}}}} \right\rangle$) (the reference upper mixed layer depth h = 45 m), and (c) entrainment flux extremum in the entrainment layer (min ($\left\langle {{w'T'}} \right\rangle $)) as a function of ${{{fh}} \mathord{\left/ {\vphantom {{{fh}} {{{u}_{*}}}}} \right. } {{{u}_{*}}}}$. The horizontal black lines indicate ${{{fh}} \mathord{\left/ {\vphantom {{{fh}} {{{u}_{*}}}}} \right. } {{{u}_{*}}}}$= 0.266 (latitude = 50°N).

Figure 4 shows the influence of Coriolis parameter on the vertical profile of the horizontal velocity. As plotted in Fig. 4a, when ${{{{fh}}} \mathord{\left/ {\vphantom {{{{fh}}} {{{{u}}_{{*}}}}}} \right. } {{{{u}}_{{*}}}}}$ (latitude) changes from 0.119 (20°N) to 0.342 (80°N), the region of the positive downwind velocity ($\left\langle {u} \right\rangle + {{u}_{s}}$) shrinks gradually and the vertical distribution of the downwind velocity becomes more uniform. Diminution in the crosswind velocity ($\left\langle {v} \right\rangle $) is distinct (Fig. 4b), whereas the variation in the vertical distribution of the crosswind velocity is not remarkable above z / h = –1 and it is notable below z / h = –1 with ${{{{fh}}} \mathord{\left/ {\vphantom {{{{fh}}} {{{{u}}_{{*}}}}}} \right. } {{{{u}}_{{*}}}}}$ (latitude) changing from 0.119 (20°N) to 0.342 (80°N). The reason is that the Coriolis–Stokes force (${ - f}{{u}_{s} }$) modulates the relative magnitude and vertical distribution of the horizontal velocity by the anti-Stokes Eulerian flow (${u} \approx { - }{{u}_{s} }$) opposing the Stokes drift velocity (Polton et al., 2005; McWilliams et al., 2014); that is, both the enlargement in the negative downwind velocity and the reduction in the crosswind velocity with the Coriolis parameter changing from weak to strong presented here is similar to that with Stokes drift velocity varying from small to large under an unchanged wind forcing condition, as discussed in the previous studies of Polton et al. (2005) and McWilliams et al. (2014). Furthermore, a decrease in the penetration depth of the horizontal velocity is evident when ${{{{fh}}} \mathord{\left/ {\vphantom {{{{fh}}} {{{{u}}_{{*}}}}}} \right. } {{{{u}}_{{*}}}}}$ (latitude) changes from 0.119 (20°N) to 0.266 (50°N), while the penetration depth is almost unchanged as ${{{{fh}}} \mathord{\left/ {\vphantom {{{{fh}}} {{{{u}}_{{*}}}}}} \right. } {{{{u}}_{{*}}}}}$ (latitude) varies from 0.266 (50°N) to 0.342 (80°N). This result is consistent with the variation of the upper mixed layer depth influenced by Coriolis parameter as shown in Fig. 3b.

Variations in the values of the depth-averaged downwind velocity ($\overline {\left\langle {u} \right\rangle + {{u}_{s}}} $, hereafter the overbar ($\overline {\left\langle \varphi \right\rangle } $, where the $\left\langle \varphi \right\rangle $ represents the arbitrary horizontally and temporally averaged physical quantity) indicates averaging over the upper mixed layer (–h ≤ z ≤ 0) as calculated by Polton et al. (2008) and McWilliams et al. (2014), the depth-averaged crosswind velocity ($\overline {\left\langle {v} \right\rangle } $), and the ratio of $\overline {\left\langle {u} \right\rangle + {{u}_{s}}} $ to $\overline {\left\langle {v} \right\rangle } $ are shown in Figs. 5a, 5b and 5c, respectively. As observed in Figs. 5a and 5b, both the enhancement in $\overline {\left\langle {u} \right\rangle + {{u}_{s}}} $ and the reduction in $\overline {\left\langle {v} \right\rangle } $ are evident when ${{{{fh}}} \mathord{\left/ {\vphantom {{{{fh}}} {{{{u}}_{{*}}}}}} \right. } {{{{u}}_{{*}}}}}$ (latitude) changes from 0.119 (20°N) to 0.266 (50°N), while their changes are very weak when ${{{{fh}}} \mathord{\left/ {\vphantom {{{{fh}}} {{{{u}}_{{*}}}}}} \right. } {{{{u}}_{{*}}}}}$ (latitude) varies from 0.266 (50°N) to 0.342 (80°N). These results suggest that the influence of Coriolis parameter by Coriolis–Stokes force on horizontal velocity changes from intense to moderate across ${{{{fh}}} \mathord{\left/ {\vphantom {{{{fh}}} {{{{u}}_{{*}}}}}} \right. } {{{{u}}_{{*}}}}}$= 0.266 (LAT = 50°N). In addition, the increase in ${{\left( {\overline {\left\langle {u} \right\rangle + {{u}_{s}}} } \right)} \mathord{\left/ {\vphantom {{\left( {\overline {\left\langle {u} \right\rangle + {{u}_{s}}} } \right)} {\overline {\left\langle {v} \right\rangle } }}} \right. } {\overline {\left\langle {v} \right\rangle } }}$is remarkable for ${{{{fh}}} \mathord{\left/ {\vphantom {{{{fh}}} {{{{u}}_{{*}}}}}} \right. } {{{{u}}_{{*}}}}}$≤ 0.266 (LAT ≤ 50°N) as compared to that for $ {{fh} \mathord{\left/ {\vphantom {{fh} {{u_*}}}} \right. } {{u_*}}} $≥ 0.266 (LAT ≥ 50°N) (Fig. 5b) as ${{{{fh}}} \mathord{\left/ {\vphantom {{{{fh}}} {{{{u}}_{{*}}}}}} \right. } {{{{u}}_{{*}}}}}$ (latitude) changes from 0.119 (20°N) to 0.342 (80°N), revealing that the reduction in the transition from $ \overline {\left\langle u \right\rangle + {u_s}} $ to $\overline {\left\langle {v} \right\rangle } $ by Lagrange Ekman effect (${f}{{u}_{L}}$, ${{u}_{L}}{\text{ = }}{u} + {{u}_{s}}$ is the Lagrange downwind velocity) (McWilliams et al., 2014) is noticeable in the former as compared to that in the latter. In addition, $\overline {\left\langle {v} \right\rangle } $ is always larger than $\overline {\left\langle {u} \right\rangle + {{u}_{s}}} $ in our simulation cases, which indicates that the basic characteristics of LT (Li et al., 2005; Sutherland et al., 2014; Li and Fox-Kemper, 2017), that is, the velocity in the crosswind direction is larger than that in the downwind direction, is not altered by Coriolis parameter.

Figure 5.  Variation of normalized depth-averaged (a) downwind Lagrange velocity ($\overline {\left\langle {u} \right\rangle + {{u}_{{\rm{s}}}}}$) and (b) crosswind velocity ($\overline {\left\langle {v} \right\rangle } $) as a function of latitude. (c) Ratio of $\overline {\left\langle {u} \right\rangle + {{u}_{{\rm{s}}}}}$ to $\overline {\left\langle {v} \right\rangle } $ as a function of ${{{fh}} \mathord{\left/ {\vphantom {{{fh}} {{{u}_{*}}}}} \right. } {{{u}_{*}}}}$. The horizontal black lines indicate ${{{fh}} \mathord{\left/ {\vphantom {{{fh}} {{{u}_{*}}}}} \right. } {{{u}_{*}}}}$= 0.266 (latitude = 50°N).

In the present simulation cases where the wind and wave fields are aligned, the horizontal flow is homogeneous and the flow is static as z→−∞. The momentum budget for the steady flow (McWilliams et al., 1997; Sullivan et al., 2007) is calculated by

Vertical profiles of the vertical momentum flux are shown in Fig. 6. When $ {{fh} \mathord{\left/ {\vphantom {{fh} {{u_*}}}} \right. } {{u_*}}} $ (latitude) changes from 0.119 (20°N) to 0.342 (80°N), the downwind vertical momentum transport ($\left\langle {{u'w'}} \right\rangle $) reduces (Fig. 6a), while the crosswind vertical momentum transport ($\left\langle {{v'w'}} \right\rangle $) increases (Fig. 6b), displaying that some of the vertical exchange of the momentum shifts from the downwind direction to the crosswind direction. In addition, the curvature of the $\left\langle {{u'w'}} \right\rangle $ profile grows with an increase of the upper mixed layer depth when the Ekman depth scale is constant in the ideal simulation cases, i.e., the value of ${{{{fh}}} \mathord{\left/ {\vphantom {{{{fh}}} {{{{u}}_{{*}}}}}} \right. } {{{{u}}_{{*}}}}}$ increases, as analyzed by Pearson et al. (2015). The simulation cases designed here show that the enhancement in the curvature of the $\left\langle {{u'w'}} \right\rangle $ profile directly associated with an enhancement in Coriolis parameter is quick for ${{{{fh}}} \mathord{\left/ {\vphantom {{{{fh}}} {{{{u}}_{{*}}}}}} \right. } {{{{u}}_{{*}}}}}$≤ 0.266 (LAT ≤ 50°N) but very slow for ${{{{fh}}} \mathord{\left/ {\vphantom {{{{fh}}} {{{{u}}_{{*}}}}}} \right. } {{{{u}}_{{*}}}}}$≥ 0.266 (LAT ≥ 50°N). This result demonstrates that, at ${{{{fh}}} \mathord{\left/ {\vphantom {{{{fh}}} {{{{u}}_{{*}}}}}} \right. } {{{{u}}_{{*}}}}}$= 0.266 (LAT = 50°N), an enlargement in the curvature of the $\left\langle {{u'w'}} \right\rangle $ profile changes from sensitive to insensitive as ${{{{fh}}} \mathord{\left/ {\vphantom {{{{fh}}} {{{{u}}_{{*}}}}}} \right. } {{{{u}}_{{*}}}}}$ (latitude) varies from 0.119 (20°N) to 0.342 (80°N). We also note that the enlargement in $\left\langle {{v'w'}} \right\rangle $ by an enhancement of Coriolis parameter with ${{{{fh}}} \mathord{\left/ {\vphantom {{{{fh}}} {{{{u}}_{{*}}}}}} \right. } {{{{u}}_{{*}}}}}$ (latitude) varying from 0.119 (20°N) to 0.342 (80°N) presented here is similar to the enhancement in $\left\langle {{v'w'}} \right\rangle $ with the upper mixed layer depth varying from shallow to deep and the sea surface buoyancy changing from strong to weak, as analyzed by Grant and Belcher (2009) and Pearson et al. (2015), respectively.

Figure 6.  Vertical profiles of the normalized (a) downwind vertical momentum flux ($\left\langle {{u'w'}} \right\rangle $) and (b) crosswind vertical momentum flux ($\left\langle {{v'w'}} \right\rangle $) with a variation of ${{{fh}} \mathord{\left/ {\vphantom {{{fh}} {{{u}_{*}}}}} \right. } {{{u}_{*}}}}$.

The changes in the depth-averaged values of the downwind vertical momentum flux ($\overline {\left\langle {{u'w'}} \right\rangle } $), the crosswind vertical momentum flux ($\overline {\left\langle {{v'w'}} \right\rangle } $), and the ratio of $\overline {\left\langle {{u'w'}} \right\rangle } $ to $ \overline {\left\langle {v'w'} \right\rangle } $ with a variation of ${{{fh}} \mathord{\left/ {\vphantom {{{fh}} {{{u}_{{*}}}}}} \right. } {{{u}_{{*}}}}}$ are shown in Fig. 7. Decay in $\overline {\left\langle {{u'w'}} \right\rangle } $ is more evident as ${{{fh}} \mathord{\left/ {\vphantom {{{fh}} {{{u}_{{*}}}}}} \right. } {{{u}_{{*}}}}}$ (latitude) varies from 0.119 (20°N) to 0.266 (50°N) compared to from 0.266 (50°N) to 0.342 (80°N) (Fig. 7a), indicating that the modulation effect of Coriolis parameter on $\overline {\left\langle {{u'w'}} \right\rangle } $ for $ {{fh} \mathord{\left/ {\vphantom {{fh} {{u_*}}}} \right. } {{u_*}}} $≤ 0.266 (LAT ≤ 50°N) is stronger than for ${{{fh}} \mathord{\left/ {\vphantom {{{fh}} {{{u}_{{*}}}}}} \right. } {{{u}_{{*}}}}}$≥ 0.266 (LAT ≥ 50°N). The rate of change of $\overline {\left\langle {{v'w'}} \right\rangle } $ does not have a clear variation as $ {{fh} \mathord{\left/ {\vphantom {{fh} {{u_*}}}} \right. } {{u_*}}} $ (latitude) varies from 0.119 (20°N) to 0.342 (80°N) (Fig. 7b), i.e., the variation of ${{{fh}} \mathord{\left/ {\vphantom {{{fh}} {{{u}_{{*}}}}}} \right. } {{{u}_{{*}}}}}$ (latitude) in the influence of Coriolis parameter on the rate of change of $\overline {\left\langle {{v'w'}} \right\rangle } $ is not evident.

Figure 7.  Normalized depth-averaged (a) downwind vertical momentum flux ($\overline {\left\langle {{u'w'}} \right\rangle } $) and (b) crosswind vertical momentum flux ($\overline {\left\langle {{v'w'}} \right\rangle } $) with a change of latitude. (c) Ratio of $\overline {\left\langle {{u'w'}} \right\rangle } $ to $\overline {\left\langle {{v'w'}} \right\rangle } $ as a function of ${{{fh}} \mathord{\left/ {\vphantom {{{fh}} {{{u}_{*}}}}} \right. } {{{u}_{*}}}}$. The horizontal black lines indicate ${{{fh}} \mathord{\left/ {\vphantom {{{fh}} {{{u}_{*}}}}} \right. } {{{u}_{*}}}}$= 0.266 (latitude = 50°N).

A decrease in ${{\overline {\left\langle {{u'w'}} \right\rangle } } \mathord{\left/ {\vphantom {{\overline {\left\langle {{u'w'}} \right\rangle } } {\overline {\left\langle {{v'w'}} \right\rangle } }}} \right. } {\overline {\left\langle {{v'w'}} \right\rangle } }}$ is notable when${{{{fh}}}/{{{{u}}_{{*}}}}}$ (latitude) changes from 0.119 (20°N) to 0.342 (80°N) as seen in Fig. 7c, which is induced by a reduction in $\overline {\left\langle {{u'w'}} \right\rangle } $ (Fig. 7a) and an enhancement in $\overline {\left\langle {{v'w'}} \right\rangle } $ (Fig. 7b). It is also found that diminution in ${{\overline {\left\langle {{u'w'}} \right\rangle } } \mathord{\left/ {\vphantom {{\overline {\left\langle {{u'w'}} \right\rangle } } {\overline {\left\langle {{v'w'}} \right\rangle } }}} \right. } {\overline {\left\langle {{v'w'}} \right\rangle } }}$ is quicker for ${{{{fh}}} /{{{{u}}_{{*}}}}}$ (latitude) varying from 0.119 (20°N) to 0.266 (50°N) relative to from 0.266 (50°N) to 0.342 (80°N) (Fig. 7c), which is attributed to the fact that a decay in $\overline {\left\langle {{u'w'}} \right\rangle } $ is more vigorous for ${{{{fh}}}/{{{{u}}_{{*}}}}}$≤ 0.266 (LAT ≤ 50°N) than for ${{{{fh}}} /{{{{u}}_{{*}}}}}$≥ 0.266 (LAT ≥ 50°N) (Fig. 7a). The change in the rate of enlargement in $\overline {\left\langle {{v'w'}} \right\rangle } $ is not obvious (Fig. 7b).

The turbulent kinetic energy budget is fundamental to the analysis of LT (McWilliams et al., 1997; Polton and Belcher, 2007; Grant and Belcher, 2009; Noh et al., 2011; Wang et al., 2018). It can be written as

where ${t}$ is time, $\left\langle E \right\rangle = \dfrac{1}{2}\left( {\left\langle {{{u'}^2}} \right\rangle + \left\langle {{{v'}^2}} \right\rangle + \left\langle {{{w'}^2}} \right\rangle } \right)$ is the total turbulent kinetic energy, $\left\langle {{{P}_{{\rm{s}}}}} \right\rangle = - \left\langle {{u'w'}} \right\rangle {{\partial \left\langle {u} \right\rangle }/ {\partial {z} - \left\langle {{v'w'}} \right\rangle }}{\partial \left\langle {v} \right\rangle }/ {\partial {z}}$ is the shear production, $\left\langle {{{P}_{{\rm{L}}}}} \right\rangle = - \left\langle {{u'w'}} \right\rangle \partial {{{{u}_{{\rm{s}}}}} \mathord{\left/ {\vphantom {{{{u}_{s}}} {\partial {z}}}} \right. } {\partial {z}}}$ is the Langmuir production, $\left\langle {{{P}_{{\rm{b}}}}} \right\rangle = {\alpha g}\left\langle {{w'T'}} \right\rangle$ (${T'}$ is the temperature perturbation) is the buoyancy production, $\left\langle {F} \right\rangle = - \partial \left( \left\langle {{w'P}} \right\rangle + \left\langle {{w'E'}} \right\rangle \right) / {\partial {z}}$ (where ${P}$ is the modified pressure and ${E'}$ is the total turbulent kinetic energy perturbation) is the transport production, and $\left\langle {\varepsilon } \right\rangle $ is the dissipation rate of turbulent kinetic energy.

Vertical profiles of the turbulent kinetic energy budget terms are plotted in Fig. 8. As shown in Fig. 8, the diminution in every term of the turbulent kinetic energy budget becomes weak as ${{{fh}} \mathord{\left/ {\vphantom {{{fh}} {{{u}_{{*}}}}}} \right. } {{{u}_{{*}}}}}$ (latitude) changes from 0.119 (20°N) to 0.342 (50°N) in the upper mixed layer. This result indicates that the effect of the Coriolis parameter on the turbulent kinetic energy budget terms weakens as the ${{{fh}} \mathord{\left/ {\vphantom {{{fh}} {{{u}_{{*}}}}}} \right. } {{{u}_{{*}}}}}$ increases. In addition, as the${{{fh}} \mathord{\left/ {\vphantom {{{fh}} {{{u}_{{*}}}}}} \right. } {{{u}_{{*}}}}}$ (latitude) changes from 0.119 (20°N) to 0.342 (80°N), the difference of the turbulent kinetic energy budget terms between the different cases in the lower part of the upper mixed layer is larger than that in the upper part of the upper mixed layer. The negative value of $\left\langle {{{P}_{{\rm{s}}}}} \right\rangle$ (Fig. 8b) is induced by the change in vertical gradient of the downwind Euler current (not shown), which is the sink for the turbulent kinetic energy budget. The positive value of $\left\langle {{{P}_{{\rm{b}}}}} \right\rangle$ near the surface layer (Fig. 8d) indicates the transformation of the available potential energy to turbulent kinetic energy, which is the source for the turbulent kinetic energy budget. Reduction in the transport production term $\left\langle {F} \right\rangle $ over the upper mixed layer depth is evident, as observed in Fig. 8e, reflecting that, as compared to the other terms, the effect of Coriolis parameter on $\left\langle {F} \right\rangle $ (Grant and Belcher, 2009) is considerable for modifying its vertical distribution and relative magnitude. These results suggest that the Coriolis parameter can influence the growth of the upper mixed layer depth by modulating the relative magnitude of turbulent kinetic energy budget terms in the entrainment layer, as shown by comparing Figs. 2 and 8.

Figure 8.  Vertical profiles of the normalized (a) total kinetic energy ($\left\langle {E} \right\rangle $), (b) shear production ($\left\langle {{{P}_{{\rm{s}}}}} \right\rangle$), (c) Langmuir production ($\left\langle {{{P}_{{\rm{L}}}}} \right\rangle$), (d) buoyancy production ($\left\langle {{{P}_{{\rm{b}}}}} \right\rangle$), (e) transport production ($\left\langle {F} \right\rangle $), and (d) dissipation rate ($\left\langle \varepsilon \right\rangle $) with a variation of ${{{fh}} \mathord{\left/ {\vphantom {{{fh}} {{{u}_{*}}}}} \right. } {{{u}_{*}}}}$.

To further investigate the differences of the influence of Coriolis parameter on the turbulent kinetic energy budget, we plot in Fig. 9 the depth-averaged values of the turbulent kinetic energy budget terms as a function of ${{{fh}} \mathord{\left/ {\vphantom {{{fh}} {{{u}_{{*}}}}}} \right. } {{{u}_{{*}}}}}$. Decreases in the depth-averaged turbulent kinetic energy ($\overline {\left\langle {E} \right\rangle } $) (Fig. 9a), shear production ($ \overline {\left\langle {{{P}_{s}}} \right\rangle } $) (Fig. 9b), Langmuir production ($ \overline {\left\langle {{{P}_{L}}} \right\rangle } $) (Fig. 9c), and dissipation rate ($\overline {\left\langle {\varepsilon } \right\rangle } $) (Fig. 9f) are evident when $ {{fh} \mathord{\left/ {\vphantom {{fh} {{u_*}}}} \right. } {{u_*}}} $ (latitude) changes from 0.119 (20°N) to 0.266 (50°N) as compared to when ${{{fh}} \mathord{\left/ {\vphantom {{{fh}} {{{u}_{{*}}}}}} \right. } {{{u}_{{*}}}}}$ (latitude) varies from 0.266 (50°N) to 0.342 (80°N), which reveals that the impact of the Coriolis parameter on these terms transforms from strong to weak at ${{{fh}} \mathord{\left/ {\vphantom {{{fh}} {{{u}_{{*}}}}}} \right. } {{{u}_{{*}}}}}$= 0.266 (LAT = 50°N). On the other hand, the rate of reduction of the depth-averaged buoyancy production ($ \overline {\left\langle {{{P}_{b}}} \right\rangle } $) (Fig. 9d) and transport production ($ \overline {\left\langle {F} \right\rangle } $) (Fig. 9e) across ${{{fh}} \mathord{\left/ {\vphantom {{{fh}} {{{u}_{{*}}}}}} \right. } {{{u}_{{*}}}}}$= 0.266 (LAT = 50°N) is not very pronounced.

Figure 9.  Normalized depth-averaged value of (a) the total turbulent kinetic energy ($\overline {\left\langle {E} \right\rangle } $), (b) the shear production ($\overline {\left\langle {{{P}_{{\rm{s}}}}} \right\rangle }$), (c) the Langmuir production ($\overline {\left\langle {{{P}_{{\rm{L}}}}} \right\rangle }$), (d) the buoyancy production ($\overline {\left\langle {{{P}_{{\rm{b}}}}} \right\rangle }$), (e) the transport production ($\overline {\left\langle {F} \right\rangle } $), and (f) the dissipation rate ($\overline {\left\langle \varepsilon \right\rangle } $) as a function of ${{{fh}} \mathord{\left/ {\vphantom {{{fh}} {{{u}_{*}}}}} \right. } {{{u}_{*}}}}$. The horizontal black lines indicate ${{{fh}} \mathord{\left/ {\vphantom {{{fh}} {{{u}_{*}}}}} \right. } {{{u}_{*}}}}$= 0.266 (latitude = 50°N).

The results in the preceding sections indicate that the influence of Coriolis parameter on the mean and turbulent statistics of Langmuir turbulence is evident for ${{{fh}} \mathord{\left/ {\vphantom {{{fh}} {{{u}_{{*}}}}}} \right. } {{{u}_{{*}}}}}$≤ 0.266 (LAT ≤ 50°N) but weak for ${{{fh}} \mathord{\left/ {\vphantom {{{fh}} {{{u}_{{*}}}}}} \right. } {{{u}_{{*}}}}}$≥ 0.266 (LAT ≥ 50°N). The f, h, and ${{u}_{*}}$ variables represent the Coriolis parameter, the initial upper mixed layer depth, and the friction velocity of sea surface induced by wind, respectively. This can be attributed to the different rates of change of the Coriolis parameter ( f and f ^1/2 ), Ekman depth scale (H_EK = ${{{{u}_{{*}}}}/ {{f}}}$) (Pearson et al., 2015; Noh and Choi, 2018), and inertial period (t_ip = ${{2\pi } / {{f}}}$) at different latitudes. Their changes for a variation of latitude from 20°N to 80°N are shown in Fig. 10.

Figure 10.  Variation of (a) Coriolis parameter ( f ), (b) square root of Coriolis parameter ( f ^1/2 ), (c) Ekman depth scale (H_EK), and (d) inertial period (t_ip) as a function of latitude.

A comparison of Fig. 3b with Fig. 10c reveals that the evolution between the upper mixed layer depth (h_LAT) (Fig. 3b) and the Ekman depth scale (H_EK) (Fig. 10c) is similar when LAT ≤ 50°N (${{{fh}}/ {{{u}_{{*}}}}}$≤ 0.266), whereas a variation in the upper mixed layer depth is evidently slow relative to Ekman depth scale when LAT ≥ 50°N (${{{fh}}/{{{u}_{{*}}}}}$≥ 0.266). In addition, a comparison between the variations with latitude of h_LAT (Fig. 3b) and f ^1/2 (Fig. 10b) shows similar differences as variations between h_LAT and H_EK. Previous studies based on ideal simulation cases have shown that although the penetration depth of Langmuir turbulence can reach the Ekman depth scale in an unstratified mixed layer (Polton and Belcher, 2007), the penetration depth is also suppressed by the pycnocline [Grant and Belcher (2009); Pearson et al. (2015); Noh and Choi (2018)]. Hence, under conditions in which the Ekman depth scale for LAT ≥ 50°N (${{{fh}} \mathord{\left/ {\vphantom {{{fh}} {{{u}_{{*}}}}}} \right. } {{{u}_{{*}}}}}$≥ 0.266) has a remarkably weak decrease and is shallow relative to that for LAT ≤ 50°N (${{{fh}} \mathord{\left/ {\vphantom {{{fh}} {{{u}_{{*}}}}}} \right. } {{{u}_{{*}}}}}$≤ 0.266), the suppression effect of the pycnocline on the penetration depth of Langmuir turbulence becomes more significant when LAT ≥ 50°N (${{{fh}} \mathord{\left/ {\vphantom {{{fh}} {{{u}_{{*}}}}}} \right. } {{{u}_{{*}}}}}$≥ 0.266), which causes the rates of decay of the upper mixed layer depth and the Ekman depth scale to be dissimilar when LAT ≥ 50°N (${{{fh}} \mathord{\left/ {\vphantom {{{fh}} {{{u}_{{*}}}}}} \right. } {{{u}_{{*}}}}}$≥ 0.266). This is also the reason that variations in the penetration depth of the horizontal velocity (Fig. 4), the vertical momentum flux (Fig. 6), and the turbulent kinetic energy budget terms (Fig. 8) significantly decrease when ${{{fh}} \mathord{\left/ {\vphantom {{{fh}} {{{u}_{{*}}}}}} \right. } {{{u}_{{*}}}}}$≥ 0.266 (LAT ≥ 50°N), as compared to when ${{{fh}} \mathord{\left/ {\vphantom {{{fh}} {{{u}_{{*}}}}}} \right. } {{{u}_{{*}}}}}$≤ 0.266 (LAT ≤ 50°N).

An enlargement in the curvature of the downwind vertical momentum flux ($ \left\langle {u'w'} \right\rangle $) profile changes from sensitive to insensitive at ${{{fh}} \mathord{\left/ {\vphantom {{{fh}} {{{u}_{{*}}}}}} \right. } {{{u}_{{*}}}}}$= 0.266 (LAT = 50°N), as shown in Fig. 6a. The curvature of the $\left\langle {{u'w'}} \right\rangle $ profile is associated with the parameter${{{fh}} \mathord{\left/ {\vphantom {{{fh}} {{{u}_{{*}}}}}} \right. } {{{u}_{{*}}}}}$, that is, the ratio of the upper mixed layer depth to the Ekman depth scale, as analyzed by Pearson et al. (2015). Here, the combined effect of the increased Coriolis parameter (Fig. 10a) and the decreased upper mixed layer depth (Fig. 3b) leads to the result that the enlargement in the curvature of the $\left\langle {{u'w'}} \right\rangle $ profile varies from quick to slow at ${{{fh}} \mathord{\left/ {\vphantom {{{fh}} {{{u}_{{*}}}}}} \right. } {{{u}_{{*}}}}}$= 0.266 (LAT = 50°N) (Fig. 3b).

In addition, the rate of change of the depth-averaged crosswind vertical momentum flux ($\overline {\left\langle {{v'w'}} \right\rangle } $) does not have a clear variation with ${{{fh}} \mathord{\left/ {\vphantom {{{fh}} {{{u}_{{*}}}}}} \right. } {{{u}_{{*}}}}}$, as shown in Fig. 7b. The crosswind vertical momentum flux ($\left\langle {{v'w'}} \right\rangle $) weakens as the upper mixed layer depth becomes shallow, as suggested by the previous study of Grant and Belcher (2009), and enhances with an increase of Coriolis parameter, as presented here (Fig. 6b). Hence, the combined effect of a decreased upper mixed layer depth (Fig. 3b) and an increased Coriolis parameter (Fig. 10a) with a variation of latitude (${{{fh}} \mathord{\left/ {\vphantom {{{fh}} {{{u}_{{*}}}}}} \right. } {{{u}_{{*}}}}}$) from 20°N (0.119) to 80°N (0.342) explains the nearly unchanged rate of change of the $\overline {\left\langle {{v'w'}} \right\rangle } $ with a variation of ${{{fh}} \mathord{\left/ {\vphantom {{{fh}} {{{u}_{{*}}}}}} \right. } {{{u}_{{*}}}}}$ (Fig. 7b).

Variation in the depth-averaged value of the horizontal velocity ($\overline {\left\langle {u} \right\rangle + {{u}_{s}}} $ and $\overline {\left\langle {v} \right\rangle } $) (Fig. 5) is more evident with a change of ${{{fh}} \mathord{\left/ {\vphantom {{{fh}} {{{u}_{{*}}}}}} \right. } {{{u}_{{*}}}}}$ (latitude) when ${{{fh}} \mathord{\left/ {\vphantom {{{fh}} {{{u}_{{*}}}}}} \right. } {{{u}_{{*}}}}}$≤ 0.266 (LAT ≤ 50°N), relative to when ${{{fh}} \mathord{\left/ {\vphantom {{{fh}} {{{u}_{{*}}}}}} \right. } {{{u}_{{*}}}}}$≥ 0.266 (LAT ≥ 50°N). The mechanism is that the modulation effect of the Coriolis–Stokes force ($- {{f}}{{u}_{{\rm{s}}} }$, ${{u}_{{\rm{s}}}}$ is the Stokes drift velocity) by anti-Stokes Eulerian flow (${u} \approx { - }{{u}_{{\rm{s}}}}$) opposing the Stokes drift velocity (Polton et al., 2005) and Lagrange Ekman effect (${f}{{u}_{{{L}}}}$, ${{u}_{{\rm{s}}}} = {u} + {{u}_{{\rm{s}}}}$ is the Lagrange downwind velocity) (McWilliams et al., 2014) on the horizontal velocity changes from strong to weak with Coriolis parameter varying from quick to slow at LAT = 50°N (Fig. 10a).

The inhibiting effect of the Coriolis parameter on the turbulent kinetic energy budget terms varies from intense to weak with a change of ${{{fh}} \mathord{\left/ {\vphantom {{{fh}} {{{u}_{\text{*}}}}}} \right. } {{{u}_{\text{*}}}}}$ from 0.119 to 0.342 at ${{{fh}} \mathord{\left/ {\vphantom {{{fh}} {{{u}_{\text{*}}}}}} \right. } {{{u}_{\text{*}}}}}$= 0.266 (LAT = 50°N) (Fig. 9). The reason is that the wind acts to remove energy from the flow after a quarter to a half inertial period (Skyllingstad et al., 2000), that is, the time interval of the wind removing the energy from the flow shortens rapidly for LAT ≥ 50°N and slightly for LAT ≤ 50°N (Fig. 10d).

The combination of the Coriolis–Stokes force, the Lagrange Ekman effect, and the ratio of the upper mixed layer depth to the Ekman depth scale can modify the mean and turbulent statistics of Langmuir turbulence by the nonlinear form. Variation in the depth-averaged value of the vertical momentum flux ($\overline {\left\langle {{u'w'}} \right\rangle } $ and $\overline {\left\langle {{v'w'}} \right\rangle } $) (Fig. 7), total turbulent kinetic energy ($\overline {\left\langle {E} \right\rangle } $) (Fig. 9a), and dissipation rate ($\overline {\left\langle \varepsilon \right\rangle } $) (Fig. 9f) is closed to linear with a change of ${{{fh}} \mathord{\left/ {\vphantom {{{fh}} {{{u}_{{*}}}}}} \right. } {{{u}_{{*}}}}}$, though the rate of change has distinct difference between ${{{fh}} \mathord{\left/ {\vphantom {{{fh}} {{{u}_{{*}}}}}} \right. } {{{u}_{{*}}}}}$≤ 0.266 (LAT ≤ 50°N) and ${{{fh}} \mathord{\left/ {\vphantom {{{fh}} {{{u}_{{*}}}}}} \right. } {{{u}_{{*}}}}}$≥ 0.266 (LAT ≥ 50°N). This result indicates that the nonlinear variation in Coriolis parameter with latitude (Fig. 10a) induces an approximately linear change in these terms. In addition, changes in the sea surface temperature (Fig. 3a), upper mixed layer depth (Fig. 3b), and strongest heat flux (Fig. 3c), and the depth-averaged value of the downwind velocity ($\overline {\left\langle {u} \right\rangle + {{u}_{{\rm{s}}}}}$) (Fig. 5a), crosswind velocity ($\overline {\left\langle {v} \right\rangle } $) (Fig. 5b), shear production ($\overline {\left\langle {{{P}_{{\rm{s}}}}} \right\rangle }$) (Fig. 9b), buoyancy production ($\overline {\left\langle {{{P}_{{\rm{b}}}}} \right\rangle }$) (Fig. 9d), and transport production ($ \overline {\left\langle {F} \right\rangle } $) (Fig. 9e) are nonlinear. This suggests that the nonlinear evolution in Coriolis parameter with latitude (Fig. 10a) produces the nonlinear change in these quantities. These results demonstrate that the influence of the Coriolis parameter on the characteristics of Langmuir turbulence can also be divided into two categories, linear and nonlinear, respectively, with a nonlinear change of Coriolis parameter.

The effects of the added swell waves on the mean and turbulent statistics of Langmuir turbulence are similar to the enhanced effects of the Coriolis parameter, which is due to the Coriolis–Stokes force (${ - f}{{u}_{{\rm{s}}} }$) containing both the Coriolis parameter ( f ) and the Stokes drift velocity (${{u}_{s}}$) (Polton et al., 2005), that is, the enlarged Stokes drift velocity also enhances anti-Stokes Eulerian flow (${u} \approx { - }{{u}_{s}}$) opposing the Stokes drift velocity. The added swell waves aligned with the wind direction may amplify the mean Lagrange current (McWilliams et al., 2014), while the added swell waves misaligned with the local wind and wave fields can cause the vertical profiles of the mean Lagrange current and the vertical momentum flux to have a fatter Ekman spiral and rotate toward the swell-wave direction, respectively (Sullivan et al., 2012; McWilliams et al., 2014).

The Coriolis parameter is negative (f < 0) in the southern hemisphere, representing that the horizontal current rotates to its left, induced by the Ekman effect. Hence, the crosswind velocity of Langmuir turbulence caused by the Ekman effect is positive in the southern hemisphere, when the direction of the wind and wave fields aligns with the positive x-direction. The influence of the Ekman effect and Coriolis–Stokes force on variation of the mean horizontal velocity with a change of latitude in the southern hemisphere is consistent with that in the northern hemisphere. The vertical momentum flux associated with the Coriolis parameter and the horizontal velocity [Eqs. (4a) and (4b)] in the southern hemisphere has consistent variation with change of latitude to that in the northern hemisphere. The variations of the entrainment depth and the turbulent kinetic energy budget terms with the latitude changing from low to high are associated with a change of the Ekman depth scale (${{{{u}_{*}}} \mathord{\left/ {\vphantom {{{{u}_{*}}} {\left| {f} \right|}}} \right. } {\left| {f} \right|}}$) (Polton and Belcher, 2007; Noh and Choi, 2018) and the inertial period (${{2\pi } \mathord{\left/ {\vphantom {{2\pi } {\left| {f} \right|}}} \right. } {\left| {f} \right|}}$) (Skyllingstad et al., 2000), respectively, which cannot be changed by positive or negative values of Coriolis parameter.

We have used a large eddy simulation model to explore the variation of Langmuir turbulence characteristics with the change of Coriolis parameter as latitude (LAT) varies from 20°N to 80°N. The rate of change of the mean and turbulent statistics of Langmuir turbulence and the influence of Langmuir turbulence on the upper mixed layer are significant for ${{{fh}} \mathord{\left/ {\vphantom {{{fh}} {{{u}_{{*}}}}}} \right. } {{{u}_{{*}}}}}$≤ 0.266 (LAT ≤ 50°N) conditions, while it is obviously weak for ${{{fh}} \mathord{\left/ {\vphantom {{{fh}} {{{u}_{{*}}}}}} \right. } {{{u}_{{*}}}}}$$ {{fh} \mathord{\left/ {\vphantom {{fh} {{u_*}}}} \right. } {{u_*}}} $≥ 0.266 (LAT ≥ 50°N). The f, h, and ${{u}_{*}}$ variables represent the Coriolis parameter, the initial upper mixed layer depth, and the friction velocity of sea surface induced by wind, respectively. In particular, the rate of change of the crosswind vertical momentum flux does not have a significant variation between ${{{fh}} \mathord{\left/ {\vphantom {{{fh}} {{{u}_{{*}}}}}} \right. } {{{u}_{{*}}}}}$≤ 0.266 (LAT ≤ 50°N) and ${{{fh}} \mathord{\left/ {\vphantom {{{fh}} {{{u}_{{*}}}}}} \right. } {{{u}_{{*}}}}}$≥ 0.266 (LAT ≥ 50°N).

The results in the present paper demonstrate that the influence of Coriolis parameter by the Coriolis–Stokes force, the Ekman effect, and the ratio of the upper mixed layer depth to the Ekman depth scale on the mean and turbulent statistics of Langmuir turbulence is complex with a change of latitude and need to be taken into account for the improvement of Langmuir turbulence parameterization.

Acknowledgements. The research of Guojing LI, Dongxiao WANG, and Yeqiang SHU was supported by the National Key Research and Development Program of China (Grant No. 2018YFC1405701), the National Natural Science Foundation of China (Grant Nos. 92158204, 41506001, 42076019, 42076026 and 41876017), and the Project supported by Southern Marine Science and Engineering Guangdong Laboratory (Guangzhou) (Grant No. GML2019ZD0304). Lian SHEN acknowledges the support by University of Minnesota. The large eddy simulation model is provided by National Center for Atmospheric Research. All numerical calculations were carried out at the High Performance Computing Center (HPCC) of the South China Sea Institute of Oceanology, Chinese Academy of Sciences.

The large eddy simulation model of Langmuir turbulence satisfies the following equations:

where ${{x}_{i}}$ (i = 1, 2, 3) are the Cartesian coordinates, ${{u}_{i}}$ (i = 1, 2, 3) are the resolved velocity components in the ${{x}_{i}}$ directions, ${t}$ is the time, ${{f}_{k}}$ is the Coriolis parameter, ${\omega _{i}}$ is the vorticity component, ${{u}_{{s}{j}}}$ is the Stokes drift velocity, ${P} = {p}/ {{\rho }_{{\rm{o}}}} + {{2{e}} / 3} + {1/ {2\left[ {{{\left( {{{u}_{i}} + {{u}_{{s}{i}}}} \right)}^2} - {{u}_{i}}{{u}_{i}}} \right]}}$ is the modified pressure, ${p}$ is the pressure, ${e}$ is the turbulent kinetic energy of sub-grid-scale, ${\tau _{{ij}}} = - {{\nu }_{t}}{{S}_{{ij}}}$ is the momentum flux of sub-grid-scale, $ {v_t} $ is the turbulent eddy viscosity, ${{S}_{{ij}}} = {1 / 2}\left( {\partial {{u}_{i}}}/ \partial {{x}_{j}} + {\partial {{u}_{j}}} / {\partial {{x}_{i}}} \right)$ is the strain rate tensor of resolved velocities, ${{\xi }_{{ijk}}}$ is the standard antisymmetric tensor, ${\delta _{{i3}}}$ is the Kronecker delta, $\rho = {\rho _{{\rm{o}}}}\left( {1 - \alpha T} \right)$ is the density, ${\rho _{{\rm{o}}}}$ is the reference density, $\alpha $ is the thermal expansion coefficient, g is the gravity acceleration, ${{S}_{{{\rm{SGS}}}}} = - {\tau _{{ij}}}{{S}_{{ij}}}$ is the shear production of sub-grid-scale, ${{B}_{{{\rm{SGS}}}}} = {{{g}{\tau _{{{\rm{T}}}{k}}}} \mathord{\left/ {\vphantom {{{g}{\tau _{{T}{k}}}} {{{T}_0}}}} \right. } {{{T}_0}}}$ is the buoyancy production of sub-grid-scale, T₀ is the reference temperature, $\varepsilon = {{0.93{e^{{3 \mathord{\left/ {\vphantom {3 2}} \right. } 2}}}} \mathord{\left/ {\vphantom {{0.93{e^{{3 \mathord{\left/ {\vphantom {3 2}} \right. } 2}}}} \Delta }} \right. } \Delta }$ ($\Delta = \sqrt[3]{{\Delta {x}\Delta {y}\Delta {z}}}$, where$\Delta {x}$, $\Delta {y}$, $\Delta {z}$ are the grid spacings) is the dissipation rate, ${{D}_{{{\rm{SGS}}}}} = {{\partial \left( {2{\nu _{t}}{{\partial {e}} \mathord{\left/ {\vphantom {{\partial {e}} {\partial {{x}_{i}}}}} \right. } {\partial {{x}_{i}}}}} \right)} \mathord{\left/ {\vphantom {{\partial \left( {2{\nu _{t}}{{\partial {e}} \mathord{\left/ {\vphantom {{\partial {e}} {\partial {{x}_{i}}}}} \right. } {\partial {{x}_{i}}}}} \right)} {\partial {{x}_{i}}}}} \right. } {\partial {{x}_{i}}}}$ is the diffusion production of sub-grid-scale, ${\nu _{t}}$ is the turbulent eddy viscosity, ${\tau _{{j}{{\rm{T}}}}} = - {\nu _{T}}{{\partial {T}} \mathord{\left/ {\vphantom {{\partial {T}} {\partial {{x}_{i}}}}} \right. } {\partial {{x}_{i}}}}$ is the heat flux of sub-grid-scale, T is the resolved temperature and ${\nu _{T}}$is the turbulent eddy diffusivity. The ${\nu _{{\rm{t}}}} = 0.1{l}{{e}^{{1 \mathord{\left/ {\vphantom {1 2}} \right. } 2}}}$ and ${\nu _{{\rm{T}}}} = \left( {1 + \left( {{{2{l}} \mathord{\left/ {\vphantom {{2{l}} \Delta }} \right. } \Delta }} \right)} \right){\nu _{t}}$ terms are suggested by Moeng (1984), where $ l = \Delta $ within the upper mixed layer and ${l} = 0.76{{e}^{{1 \mathord{\left/ {\vphantom {1 2}} \right. } 2}}}\left( {{{{g}\partial {T}} \mathord{\left/ {\vphantom {{{g}\partial {T}} {{{T}_0}\partial {z}}}} \right. } {{{T}_0}\partial {z}}}} \right)$ in the thermocline.