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

Dongxiao WANG, Guojing LI*,3, Lian SHEN, and Yeqiang SHU,3

1State Key Laboratory of Tropical Oceanography, South China Sea Institute of Oceanology,Chinese Academy of Sciences, Guangzhou 510301, China

2Department of Mechanical Engineering and St. Anthony Falls Laboratory,University of Minnesota, Minneapolis, Minnesota 55455, USA

3Southern Marine Science and Engineering Guangdong Laboratory (Guangzhou), Guangzhou 510301, China

ABSTRACT Langmuir turbulence is a complex turbulent process in the ocean upper mixed layer. The Coriolis parameter has an important effect on Langmuir turbulence through the Coriolis–Stokes force and Ekman effect, however, this effect on Langmuir turbulence has not been systematically investigated. Here, the impact of the Coriolis parameter on Langmuir turbulence with a change of latitude (LAT) from 20°N to 80°N is studied using a non-hydrostatic large eddy simulation model under an ideal condition. The results show that the ratio of the upper mixed layer depth to Ekman depth scale (RME)RME = 0.266 (LAT = 50°N) is a key value (latitude) for the modulation effect of the Coriolis parameter on the mean and turbulent statistics of Langmuir turbulence. It is found that the rate of change of the sea surface temperature, upper mixed layer depth, entrainment flux, crosswind velocity, downwind vertical momentum flux, and turbulent kinetic energy budget terms associated with Langmuir turbulence are more evident at RME ≤ 0.266 (LAT ≤ 50°N) than at RME ≥ 0.266 (LAT ≥50°N). However, the rate of change of the depth-averaged crosswind vertical momentum flux does not have a clear variation between RME ≤ 0.266 and RME ≥ 0.266. The complex changes of both Langmuir turbulence characteristics and influence of Langmuir turbulence on the upper mixed layer with latitude presented here may provide more information for further improving Langmuir turbulence parameterization.

Key words:Langmuir turbulence,Coriolis parameter,the upper mixed layer,large eddy simulation

1.Introduction

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 [? fus, where usis the Stokes drift velocity], through the anti-Stokes Eulerian flow (u ≈?us) 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.

2.Large eddy simulation model and ideal simulation cases

2.1.Large eddy simulation model

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 (usω) (Craik, 1977; Leibovich, 1977, 1983)that contains the Stokes drift velocity (us) resultant from surface wave fields and the vertical vorticity induced by horizontal current (ω =▽uhc, uhcis 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 τais the wind stress at the sea surface, ρois the sea water density, ρa(bǔ)is the air density, Uais the wind velocity(z = 10 m), and Cd= 10?3(0.79+0.0509Ua) is the drag coefficient for Ua≥ 10 m s?1.

The Stokes drift velocity (us) 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

2.2.Ideal simulation cases

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?5s?1,7.292 × 10?5s?1, 9.3744 × 10?5s?1, 1.1172 × 10?4s?1,1.2630 × 10?4s?1, 1.3704 × 10?4s?1, and 1.4362 × 10?4s?1,and inertial period tip= 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 dT/dz=0.2 K m?1; 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?1,which corresponds to the wind velocity Ua= 15.1 m s?1(z =10 m). The sea surface Stokes drift velocity uos= 0.21 m s?1.Hence, the turbulence Langmuir number0.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?6K m s?1is 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 largescale 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 ttot= 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 〈φ〉, where φ represents an arbitrary physical quantity.

3.Results

3.1.Vertical velocity, upper mixed layer, and entrainment flux

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 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(〈 w′T′〉, 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.

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

Fig. 2. Vertical profiles of the normalized (a) temperature (〈 T〉) (For reference, the initial mixed-layer sounding is shown as a black fine line. The reference temperature Tr = 302 K) and (b) entrainment flux (〈 w′T′〉) with a variation of fh/u?.

The changes of the upper mixed layer and the entrainment flux with a variation of fh/u?are illustrated in this paragraph. Figure 3 shows the changes in the sea surface temperature (〈 SST〉) , the upper mixed layer depth (〈 hLAT〉), and the entrainment flux extremum [min (〈 w′T′〉)] 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 〈 SST〉 (Fig. 3a),the 〈 hLAT〉 (Fig. 3b), and the min (〈 w′T′〉) (Fig. 3c) indicates that their variations are evident for fh/u?≤ 0.266 (LAT ≤50°N) but relatively weak for fh/u?≥ 0.266 (LAT ≥ 50°N).This result indicates that fh/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/u?≤ 0.266(LAT ≤ 50°N) and fh/u?≥ 0.266 (LAT ≥ 50°N).

3.2.Horizontal velocity and vertical momentum flux

Figure 4 shows the influence of Coriolis parameter on the vertical profile of the horizontal velocity. As plotted in Fig. 4a, when fh/u?(latitude) changes from 0.119 (20°N) to 0.342 (80°N), the region of the positive downwind velocity(〈 u〉+us) shrinks gradually and the vertical distribution of the downwind velocity becomes more uniform. Diminution in the crosswind velocity (〈 v〉) 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/u?(latitude) changing from 0.119(20°N) to 0.342 (80°N). The reason is that the Coriolis–Stokes force (? fus) modulates the relative magnitude and vertical distribution of the horizontal velocity by the anti-Stokes Eulerian flow (u ≈?us) 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/u?(latitude) changes from 0.119 (20°N) to 0.266(50°N), while the penetration depth is almost unchanged as fh/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.

Fig. 3. Variation of the normalized (a) sea surface temperature (〈 SST〉) (the reference temperature Tr = 302 K), (b)upper mixed layer depth (〈 hLAT〉) (the reference upper mixed layer depth h = 45 m), and (c) entrainment flux extremum in the entrainment layer (min (〈 w′T′〉 ) ) as a function of fh/u?. The horizontal black lines indicate fh/u?= 0.266 (latitude =50°N).

Fig. 4. Vertical profiles of the normalized (a) downwind Lagrange velocity (〈 u〉+us) (the vertical solid line indicates zer o value) and (b) crosswind velocity (〈 v〉) with a variation of fh/u?.

Variations in the values of the depth-averaged downwind velocity (, hereafter the overbar (, where the 〈φ〉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, and the ratio oftoare shown in Figs. 5a, 5b and 5c, respectively. As observed in Figs. 5a and 5b, both the enhancement inand the reduction inare evident when fh/u?(latitude) changes from 0.119(20°N) to 0.266 (50°N), while their changes are very weak when fh/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/u?= )0/.266 (LAT = 50°N). In addition, the increase inis remarkable for fh/u?≤ 0.266 (LAT ≤ 50°N) as compared to that for fh/u?≥ 0.266 (LAT ≥ 50°N) (Fig. 5b) as fh/u?(latitude)changes from 0.119 (20°N) to 0.342 (80°N), revealing that the reduction in the transition fromtoby Lagrange Ekman effect ( fuL, uL=u+usis the Lagrange downwind velocity) (McWilliams et al., 2014) is noticeable in the former as compared to that in the latter. In addition,is always larger thanin 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.

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/u?(latitude) changes from 0.119(20°N) to 0.342 (80°N), the downwind vertical momentum transport (〈 u′w′〉) reduces (Fig. 6a), while the crosswind vertical momentum transport (〈 v′w′〉) 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 〈u′w′〉 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/u?increases, as analyzed by Pearson et al. (2015). The simulation cases designed here show that the enhancement in the curvature of the 〈u′w′〉 profile directly associated with an enhancement in Coriolis parameter is quick for fh/u?≤ 0.266 (LAT ≤ 50°N) but very slow for fh/u?≥ 0.266 (LAT ≥ 50°N). This result demonstrates that,at fh/u?= 0.266 (LAT = 50°N), an enlargement in the curvature of the 〈 u′w′〉 profile changes from sensitive to insensitive as fh/u?(latitude) varies from 0.119 (20°N) to 0.342 (80°N). We also note that the enlargement in 〈 v′w′〉 by an enhancement of Coriolis parameter with fh/u?(latitude) varying from 0.119 (20°N) to 0.342 (80°N) presented here is similar to the enhancement in 〈v′w′〉 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.

The changes in the depth-averaged values of the downwind vertical momentum flux, the crosswind vertical momentum flux, and the ratio oftowith a variation of fh/u?are shown in Fig. 7. Decay inis more evident as fh/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 onfor fh/u?≤ 0.266 (LAT ≤50°N) is stronger than for fh/u?≥ 0.266 (LAT ≥ 50°N).The rate of change ofdoes not have a clear variation as fh/u?(latitude) varies from 0.119 (20°N) to 0.342 (80°N)(Fig. 7b), i.e., the variation of fh/u?(latitude) in the influence of Coriolis parameter on the rate of change ofis not evident.

Fig. 5. Variation of normalized depth-averaged (a) downwind Lagrange velocity ( ) and (b) crosswind velocity ( ) as a function of latitude. (c) Ratio of to as a function of fh/u?. The horizontal black lines indicate fh/u?= 0.266(latitude = 50°N).

Fig. 6. Vertical profiles of the normalized (a) downwind vertical momentum flux (〈 u′w′〉) and (b) crosswind vertical momentum flux (〈 v′w′〉) with a variation of fh/u?.

Fig. 7. Normalized depth-averaged (a) downwind vertical momentum flux ( 〈u′w′〉) and (b) crosswind vertical momentum flux ( ) with a change of latitude. (c) Ratio of to as a function of fh/u?. The horizontal black lines indicate fh/u?= 0.266 (latitude = 50°N).

3.3.Turbulent kinetic energy budget

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

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/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/u?increases.In addition, as the fh/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〈Ps〉 (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 〈Pb〉 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〈F〉 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 〈F〉 (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.

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/u?.Decreases in the depth-averaged turbulent kinetic energy(Fig. 9a), shear production(Fig. 9b), Langmuir production(Fig. 9c), and dissipation rate(Fig.9f) are evident when fh/u?(latitude) changes from 0.119(20°N) to 0.266 (50°N) as compared to when fh/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/u?= 0.266 (LAT= 50°N). On the other hand, the rate of reduction of the depth-averaged buoyancy production(Fig. 9d) and transport production(Fig. 9e) across fh/u?= 0.266(LAT = 50°N) is not very pronounced.

Fig. 8. Vertical profiles of the normalized (a) total kinetic energy (〈 E〉) , (b) shear production (〈 Ps〉), (c) Langmuir production(〈 PL〉) , (d) buoyancy production (〈 Pb〉) , (e) transport production (〈 F〉) , and (d) dissipation rate (〈 ε〉) with a variation of fh/u?.

Fig. 9. Normalized depth-averaged value of (a) the total turbulent kinetic energy ( 〈E〉), (b) the shear production ( 〈Ps〉), (c)the Langmuir production (), (d) the buoyancy production (), (e) the transport production (), and (f) the dissipation rate () as a function of fh/u?. The horizontal black lines indicate fh/u?= 0.266 (latitude = 50°N).

4.Dynamic analysis and discussion

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/u?≤0.266 (LAT ≤ 50°N) but weak for fh/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 f1/2), Ekman depth scale (HEK= u?/f)(Pearson et al., 2015; Noh and Choi, 2018), and inertial period (tip= 2π/f) at different latitudes. Their changes for a variation of latitude from 20°N to 80°N are shown in Fig. 10.

A comparison of Fig. 3b with Fig. 10c reveals that the evolution between the upper mixed layer depth (hLAT)(Fig. 3b) and the Ekman depth scale (HEK) (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 hLAT(Fig. 3b) and f1/2(Fig. 10b) shows similar differences as variations between hLATand HEK. 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/u?≥0.266) has a remarkably weak decrease and is shallow relative to that for LAT ≤ 50°N ( fh/u?≤ 0.266), the suppression effect of the pycnocline on the penetration depth of Langmuir turbulence becomes more significant when LAT ≥ 50°N( fh/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/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/u?≥ 0.266 (LAT ≥ 50°N), as compared to when fh/u?≤ 0.266 (LAT ≤ 50°N).

An enlargement in the curvature of the downwind vertical momentum flux (〈 u′w′〉) profile changes from sensitive to insensitive at fh/u?= 0.266 (LAT = 50°N), as shown in Fig. 6a. The curvature of the 〈u′w′〉 profile is associated with the parameter fh/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 〈u′w′〉 profile varies from quick to slow at fh/u?= 0.266 (LAT = 50°N)(Fig. 3b).

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

In addition, the rate of change of the depth-averaged crosswind vertical momentum fluxdoes not have a clear variation with fh/u?, as shown in Fig. 7b. The crosswind vertical momentum flux (〈 v′w′〉) 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/u?) from 20°N(0.119) to 80°N (0.342) explains the nearly unchanged rate of change of thewith a variation of fh/u?(Fig. 7b).

Variation in the depth-averaged value of the horizontal velocity (and) (Fig. 5) is more evident with a change of fh/u?(latitude) when fh/u?≤ 0.266 (LAT ≤ 50°N), relative to when fh/u?≥ 0.266 (LAT ≥ 50°N). The mechanism is that the modulation effect of the Coriolis–Stokes force (? fus, usis the Stokes drift velocity) by anti-Stokes Eulerian flow (u ≈?us) opposing the Stokes drift velocity(Polton et al., 2005) and Lagrange Ekman effect ( fuL,us=u+usis 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/u*from 0.119 to 0.342 at fh/u*=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 (and) (Fig. 7), total turbulent kinetic energy(Fig. 9a), and dissipation rate(Fig. 9f) is closed to linear with a change of fh/u?, though the rate of change has distinct difference between fh/u?≤0.266 (LAT ≤ 50°N) and fh/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(Fig. 5a), crosswind velocity(Fig. 5b), shear production(Fig. 9b),buoyancy production(Fig. 9d), and transport production(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 (? fus) containing both the Coriolis parameter ( f ) and the Stokes drift velocity (us) (Polton et al., 2005), that is, the enlarged Stokes drift velocity also enhances anti-Stokes Eulerian flow (u ≈?us) 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 xdirection. 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?/|f|) (Polton and Belcher, 2007; Noh and Choi, 2018) and the inertial period (2 π/|f|) (Skyllingstad et al., 2000), respectively,which cannot be changed by positive or negative values of Coriolis parameter.

5.Summary

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/u?≤ 0.266 (LAT ≤ 50°N) conditions, while it is obviously weak for fh/u?fh/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/u?≤ 0.266 (LAT ≤50°N) and fh/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.

APPENDIX

Large Eddy Simulation Model

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

where xi(i = 1, 2, 3) are the Cartesian coordinates, ui(i = 1,2, 3) are the resolved velocity components in the xidirections, t is the time, fkis the Coriolis parameter, ωiis the vorticity comp[onent, usjis the] Stokes drift velocity, P=p/ρo+is the modified pressure, p is the pressure, e is the turbulent kinetic energy of sub-gridscale, τij=?νtSijis the momentum flux of (sub-grid-scale,vtis the turbulent eddy viscosity, Sij=1/2?ui/?xj+?uj/?xi) is the strain rate tensor of resolved velocities, ξijkis the standard antisymmetric tensor, δi3is the Kronecker delta,ρ=ρo(1?αT) is the density, ρois the reference density, α is the thermal expansion coefficient, g is the gravity acceleration, SSGS=?τijSijis the shear production of sub-gridscale, BSGS=gτTk/T0is the buoyancy production of su/bgrid-s cale,is the reference temperature,where? x, ?y, ?z are the grid spacings) is the dissipation rate, DSGS=?(2νt?e/?xi)/?xiis the diffusion production of sub-grid-scale, νtis the turbulent eddy viscosity, τjT=?νT?T/?xiis the heat flux of sub-grid-scale, T is the resolved temperature and νTis the turbulent eddy diffusivity. The νt=0.1le1/2and νT=(1+(2l/?))νtterms are suggested by Moeng (1984), where l=? within the upper mixed layer and l =0.76e1/2(g?T/T0?z) in the thermocline.