Parameterization of Wave-Induced Mixing Using the Large Eddy Simulation (LES) (I)

: Turbulent motions in the thin ocean surface boundary layer control exchanges of momentum, heat and trace gases between the atmosphere and ocean. However, present parametric equations of turbulent motions that are applied to global climate models result in systematic or substantial errors in the ocean surface boundary layer. Signiﬁcant mixing caused by surface wave processes is missed in most parametric equations. A Large Eddy Simulation model is applied to investigate the wave-induced mixed layer structure. In the wave-averaged equations, wave e ﬀ ects are calculated as Stokes forces and breaking waves. To examine the e ﬀ ects of wave parameters on mixing, a series of wave conditions with varying wavelengths and heights are used to drive the model, resulting in a variety of Langmuir turbulence and wave breaking outcomes. These experiments suggest that wave-induced mixing is more sensitive to wave heights than to the wavelength. A series of numerical experiments with di ﬀ erent wind intensities-induced Stokes drifts are also conducted to investigate wave-induced mixing. As the wind speed increases, the inﬂuence depth of Langmuir circulation deepens. Additionally, it is observed that breaking waves could destroy Langmuir cells mainly at the sea surface, rather than at deeper layers.


Introduction
The ocean mixed layer is central to heat and momentum exchange processes between the upper ocean and the atmosphere. Waves play an integral part as they modulate those exchanges across the air-sea interface and consequently, their influence on mixed layer dynamics and structure cannot be ignored [1]. Specifically, as wind blows over the ocean surface, momentum is transported into the ocean, resulting in enhancement of both mixing and turbulent kinetic energy. Stokes drift also plays a role in modulating upper-ocean turbulence through large-scale Coriolis-Stokes and small-scale Stokes-vortex forces. It is the interaction of Stokes drift with wind-driven surface shear flow that induces Langmuir circulation (LC) formation [2]. As one of the most significant features of the mixed layer, LC appears as an array of alternating horizontal roll vortices, whose axes are roughly parallel to the wind direction [3][4][5]. LC could enhance the vertical mixing and alter vertical profiles of temperature and velocity [6][7][8][9], making it crucial in upper ocean studies. Upper ocean turbulence

Experimental Setup
The three-dimensional primitive equations of momentum, temperature, humidity, and other scalars are solved by LES models. In these models, two turbulence scales are employed to represent turbulence processes. Scales that are larger than a certain filter width are directly resolved, and those that are smaller than the model grid sizes are parameterized by a subgrid-scale turbulence model [49]. It is upon the separation of scales that LES models are built. In this study, simulations are carried out using the parallelized Large-Eddy Simulation Model (PALM) [49] which can simultaneously account for LC and BW effects. Craik and Leibovich [2] indicate that the vortex force drives the formation of LC. The vortex force (u s × ω) is the interaction between Stokes drift (u s ) and vorticity (ω), and the interaction between Stokes drift and planetary vorticity is termed the Coriolis-Stokes force (u s × f , f is the Coriolis parameter.), the vortex force represents the generation of LC in PALM, the Coriolis-Stokes force influences the structure of the ocean surface boundary layer. F i is the random forcing representing BW-induced small scale turbulence (assuming that the random forcing exists only at the sea surface (z = 0)); the random forcing only produces turbulence with the integral length l 0 and the time scale t 0 = 0.1l 0 /(αu * ), where α is a constant of proportionality and u * is the frictional velocity in the ocean. The modified filtered momentum equation is introduced to the PALM as follows [7,49]: where u i is the velocity component, t is the time, u sj is the Stokes velocity component, which represents the value u s in Equation (3), p is pressure, ρ 0 = 1 kg/m 3 is the density of dry air at the sea surface, ε is the Levi-Civita symbol, ω k is the vorticity, τ ij is the subgrid-scale Reynolds stress. Subscripts i, j, k ∈ {1, 2, 3} and x is the grid location. G(0; 1) is the Gaussian distribution [50], δ is the Kronecker delta. For simplicity, both the wind stress and wave fields are assumed to be in the x direction. The model depends on the non-hydrostatic, filtered, incompressible Navier-Stokes equations in Boussinesq approximation form, it is based on the Cartesian grid. In PALM, there are two versions, atmospheric version and ocean version, respectively. There is an option to turn on the ocean version. A variety of boundary conditions are provided by PALM. Dirichlet or Neumann boundary conditions could be used for the velocity, the temperature and the perturbation pressure. Neumann conditions are used in our model setting. The Neumann condition, which yields a free slip condition, is chosen as boundary conditions in the experiments, and Neumann conditions are also applied in the subgrid-scale turbulent kinetic energy. The vertical coordinate of the model points upper ward, and the original point is set at the sea surface, so that the depth z in the ocean is negative value. The model domain is 300 and 60 m in the horizontal and vertical directions, respectively, with 1 m grid resolution. Starting from rest, surface cooling (0.47 × 10 -4 K m/s) is applied during the beginning period of 900 s to initiate turbulent motions. An equilibrium is reached after 8 h [7]; the simulation time is 13 h in our experiments. The wave field is assumed to be steady and monochromatic, allowing for a multitude of outcomes for different wave conditions with varying wavelengths and wave heights to be studied (Section 3.1). The associated Stokes velocity is given by the following Equation (3), where g is the gravitational acceleration, H is the wave height and λ is the wave length. z is the depth and U s is the Stokes velocity at the ocean surface. The wavelength and wave height are set at 40 m and 1 m, respectively. The friction velocity in the ocean is set at 0.01 m/s derived from the wind stress (u * = τ/ρ w ) [43] which represents a wind speed of nearly 9.5 m/s and a Langmuir number less than 1, ρ w = 1003 kg/m 3 is the density of water, τ is the wind stress. Here, this friction velocity is used to calculate the random forcing and the Langmuir number in PALM. The random forcing F i represents BW, which is only at the surface. The wind stress is applied as the boundary condition, which is an initialization parameter in PALM. The turbulence length scale at the surface is set at 1 m and the random forcing time scale (t 0 = 0.1l 0 /(αu * )) is set at 3.34 [7]. The turbulent Langmuir number, quantifies the relative influence of wind-driven shear instability and the Stokes drift Coriolis force in the turbulent mixing layer. Observed values of La t are less than 1 [51], to guarantee the formation of the Langmuir Circulation we consider a specific range of Langmuir number, which is less than 1. In this range, turbulent eddies begin to resemble Langmuir cells. To examine the effects of wave parameters on mixing, wave conditions with varying wavelengths and wave heights are used to drive the model into a variety of Langmuir turbulence outcomes. Table 1 displays the values of the turbulent Langmuir number for different combinations of wavelengths and wave heights. It can be noted that the Langmuir number decreases with an increase in the wave height, and increases with an increase in wavelengths.  (1) and (2), this study introduces the Stokes velocity Equation (5) proposed by McWilliams and Restrepo [52] and examines the response of the upper mixed layer to a variety of LC and BW conditions under different wind speeds (Section 3.2). The equation τ = ρ 0 C D U 2 = ρ w u 2 * is for wind stress, and C D is the drag coefficient, calculated as shown in Equation (6) [53]. U is the 10-m wind speed (vector quantity) and W = |U| is the absolute wind speed (scalar).
To examine the effects of LC and BW in the upper mixed layer, a series of numerical experiments were conducted. Here, the B and the L, respectively, denote the effects of BW and LC. O represents the experiments without the effects of LC and BW. EXP and EXP WIND were used to distinguish Stokes velocities calculated with Equations (3) and (5), respectively.

Vertical Velocity, Wavelength and Wave Height
In order to analyze the influence of BW and LC on the mixed layer, vertical velocity is examined in this subsection with four experiments being performed. Experiment EXPO considered vertical velocity profiles without the influence of BW and LC, EXPB included only BW, EXPL included only LC, and EXPLB considered both LC and BW. The three-dimensional instantaneous vertical velocity fields (from the surface to a depth of 51 m) from the four experiments are plotted in Figure 1a (EXPO), c (EXPB), e (EXPL) and g (EXPLB), respectively. Corresponding vertical velocity fields at 51 m are plotted in Figure 1b  From Figure 1c,g (EXPB and EXPLB, respectively), it can be observed that small-scale turbulence as generated by BW appeared only near the surface, enhancing the velocity field intensity. From Figure 1a-d, there were minimal differences between EXPB and EXPO. When BW was considered, strong small-scale turbulence was enhanced near the surface (Figure 1a,c). In Figure 2d, stripes, indicative of LC, could be observed below the near-surface, confirming that the influence of BW was limited to the near-surface. Downward velocity in EXPLB (Figure 2c) was slightly weaker than in EXPL (Figure 2d), demonstrating that BW may disturb the formation of Langmuir cells at the sea surface.
When LC was considered (EXPL and EXPLB), vertical velocity fields were significantly affected (Figure 1e,g). In the presence of LC (Figures 1e and 2c,d), stripes were parallel to one another, indicating the formation of Langmuir cells. Moreover, the downwelling velocity field, which increased up to 0.03 m/s, was narrower and more intense than the upwelling velocity field. In Figures 1e and 2c, the stripes widened, and their number decreased as depth increased. Furthermore, it can also be observed that the distance between stripes increased with depth. From Figure 1f,h, the downwelling zone was fragmented and became dispersed, showing that the influence of LC was restricted by depth. The comparison of Figure 1b,d,f,h illustrates that the vertical velocity with LC was larger than that without LC. The value of vertical velocity was nearly unchanged when only considering BW (EXPB), which means the influence of BW was limited.
In summary, from the analysis of the above four experiments, it can be concluded that the inclusion of the LC generated stripes of strong vertical velocities down to a certain depth, the inclusion of the BW generated small-scale turbulence and weakened the near-surface LC, and the LC could influence ocean mixing to a greater depth than BW. sea surface.
When LC was considered (EXPL and EXPLB), vertical velocity fields were significantly affected (Figure 1e,g). In the presence of LC (Figures 1e and 2c,d), stripes were parallel to one another, indicating the formation of Langmuir cells. Moreover, the downwelling velocity field, which increased up to 0.03 m/s, was narrower and more intense than the upwelling velocity field. In Figures  1e and 2c, the stripes widened, and their number decreased as depth increased. Furthermore, it can also be observed that the distance between stripes increased with depth.
From Figure 1f, and h, the downwelling zone was fragmented and became dispersed, showing that the influence of LC was restricted by depth. The comparison of Figure 1b,d,f,h illustrates that the vertical velocity with LC was larger than that without LC. The value of vertical velocity was nearly unchanged when only considering BW (EXPB), which means the influence of BW was limited.
In summary, from the analysis of the above four experiments, it can be concluded that the inclusion of the LC generated stripes of strong vertical velocities down to a certain depth, the inclusion of the BW generated small-scale turbulence and weakened the near-surface LC, and the LC could influence ocean mixing to a greater depth than BW.  To test the sensitivity of wavelengths and wave heights on LC in the mixed layer, a series of wave conditions featured by varying wavelengths and wave heights were used to drive the LES model. In Figure 3, a horizontal cross-section of vertical velocity (z = 8 m) is displayed under different wavelengths and wave heights. The effects of LC and BW were included in all the experiments. Results show that when the wavelength was constant, the magnitude of the vertical velocity increased as wave height increased. Furthermore, it is observed that when wave heights increased, the Langmuir number decreased (Table 1)  To test the sensitivity of wavelengths and wave heights on LC in the mixed layer, a series of wave conditions featured by varying wavelengths and wave heights were used to drive the LES model. In Figure 3, a horizontal cross-section of vertical velocity (z = 8 m) is displayed under different wavelengths and wave heights. The effects of LC and BW were included in all the experiments. Results show that when the wavelength was constant, the magnitude of the vertical velocity increased as wave height increased. Furthermore, it is observed that when wave heights increased, the Langmuir number decreased (Table 1) and the vertical velocity increased. The stripes became more evident and more coherent. Alternatively, when the wave height was constant, the Langmuir number increased as wavelength increased (Table 1), but the value of the vertical velocity remained almost unchanged. The above results suggest that LC was more sensitive to changes in wave heights rather than changes in wavelengths.

Wind Speed
However, although Equation (3) was used to calculate Stokes drift, and in the above experiments (see Figure 3), with different wavelengths and wave heights, was used to validate the effects that different monochromatic wave-generated LC has on upper-ocean mixed layer vertical velocity, the real ocean is polychromatic, rather than monochromatic. Consequently, Equation (5) was better at representing the effects LC had on wave-induced mixing. Using Equation (5), Figure 4 presents horizontal cross-sections of instantaneous vertical velocity fields from two new experiments (where EXP WIND L considered just LC and EXP WIND LB considered both LC and BW) with increasing depths (z = 1, 8, 18, and 42 m). It can be observed that stripes appeared in both experiments but the phenomenon was more apparent in EXP WIND L than in EXP WIND LB, especially near the sea surface. It can be highlighted that as compared to Figure 1c (EXPL), the instantaneous vertical velocity profiles of Figure 4 were stronger. When z = 42 m, the downwelling zone's striping pattern could still be observed, suggesting that the depth of LC influence was deeper than in EXPL (Figure 1f). In the presence of LC, near-surface velocity fields exhibited striping patterns. In contrast, in the presence of BW (EXP WIND LB), near-surface velocity fields were dominated by small-scale turbulence and striping patterns appeared several meters below the surface. Comparing EXP WIND LB and EXP WIND L, it can be inferred that BW weakened LC. Leibovich [54] points out that when wind speeds were on the order of 8-18 m/s, phenomena such as turbulence and submesoscale processes may have intermittently occurred. Therefore, wind speeds from 8-16 m/s were chosen to investigate the response of the upper-ocean mixed layer to wind speeds under different Langmuir number settings ( Table 2).

Wind Speed
However, although Equation (3) was used to calculate Stokes drift, and in the above experiments (see Figure 3), with different wavelengths and wave heights, was used to validate the effects that different monochromatic wave-generated LC has on upper-ocean mixed layer vertical velocity, the real ocean is polychromatic, rather than monochromatic. Consequently, Equation (5) was better at representing the effects LC had on wave-induced mixing. Using Equation (5), Figure 4 presents horizontal cross-sections of instantaneous vertical velocity fields from two new experiments (where EXP WIND L considered just LC and EXP WIND LB considered both LC and BW) with increasing depths (z = 1, 8, 18, and 42 m). It can be observed that stripes appeared in both experiments but the phenomenon was more apparent in EXP WIND L than in EXP WIND LB, especially near the sea surface. It can be highlighted that as compared to Figure 1c (EXPL), the instantaneous vertical velocity profiles of Figure 4 were stronger. When z = 42 m, the downwelling zone's striping pattern could still be observed, suggesting that the depth of LC influence was deeper than in EXPL (Figure 1f). In the presence of LC, near-surface velocity fields exhibited striping patterns. In contrast, in the presence of BW (EXP WIND LB), near-surface velocity fields were dominated by small-scale turbulence and striping patterns appeared several meters below the surface. Comparing EXP WIND LB and EXP WIND L, it can be inferred that BW weakened LC. Leibovich [54] points out that when wind speeds were on the order of 8-18 m/s, phenomena such as turbulence and submesoscale processes may have intermittently occurred. Therefore, wind speeds from 8-16 m/s were chosen to investigate the response of the upper-ocean mixed layer to wind speeds under different Langmuir number settings (Table 2).  Three representative wind speeds were chosen: 8, 10 and 12 m/s ( Figure 5). When the wind speed was 8 m/s, the vertical velocity decreased as depth increased. The stripes became less obvious, the distance between them got larger, and at 43 m depth, the stripes disappeared. When the wind speed was 10 m/s, the same phenomenon was observed as depth increased but stripes were still noticeable at 43 m. When the wind speed was 12 m/s, stripes also existed. They were more obvious than when the wind speed was less than 12 m/s. The vertical velocity increased as the wind speed increased. Stripe width increased with wind speeds. When the wind speed was 8 m/s, at z = 43 m, the velocity field was less isotropic. When the wind speed was 12 m/s, at z = 43 m, stripes could still be observed on the velocity field. Consequently, it can be surmised that different intensities of wind speed generated different LC. As wind speed increased, the influence depth of LC deepened.  Three representative wind speeds were chosen: 8, 10 and 12 m/s ( Figure 5). When the wind speed was 8 m/s, the vertical velocity decreased as depth increased. The stripes became less obvious, the distance between them got larger, and at 43 m depth, the stripes disappeared. When the wind speed was 10 m/s, the same phenomenon was observed as depth increased but stripes were still noticeable at 43 m. When the wind speed was 12 m/s, stripes also existed. They were more obvious than when the wind speed was less than 12 m/s. The vertical velocity increased as the wind speed increased. Stripe width increased with wind speeds. When the wind speed was 8 m/s, at z = 43 m, the velocity field was less isotropic. When the wind speed was 12 m/s, at z = 43 m, stripes could still be observed on the velocity field. Consequently, it can be surmised that different intensities of wind speed generated different LC. As wind speed increased, the influence depth of LC deepened.

Eddy Viscosity
The profile of mean eddy viscosity for momentum K m ( Figure 6) is deduced from Equation (7), where , ̅ is the mean horizontal velocity. Eddy viscosity is usually used to quantify the upper ocean mixing strength and mixed layer depth as a mixing coefficient [55]. In both Figure 6a,b, we can see that Langmuir cells increased eddy viscosity, congruent with results from Noh et al. [56] which suggests that LC made the eddy viscosity K m much larger within the mixed layer. When considering LC (LB and L), the eddy viscosity of polychromatic waves (Figure 6b) was much larger than that of monochromatic waves (Figure 6a). The result of the picture suggests that BW could not, alone, change the eddy viscosity to a noticeable degree. However, only in conjunction with LC, could BW change eddy viscosity profiles and even decrease eddy viscosity. Moreover, LC together with BW could influence ocean mixing strength and mixed layer depth.

Eddy Viscosity
The profile of mean eddy viscosity for momentum K m ( Figure 6) is deduced from Equation (7), where (u, v) is the mean horizontal velocity. Eddy viscosity is usually used to quantify the upper ocean mixing strength and mixed layer depth as a mixing coefficient [55]. In both Figure 6a,b, we can see that Langmuir cells increased eddy viscosity, congruent with results from Noh et al. [56] which suggests that LC made the eddy viscosity K m much larger within the mixed layer. When considering LC (LB and L), the eddy viscosity of polychromatic waves (Figure 6b) was much larger than that of monochromatic waves (Figure 6a). The result of the picture suggests that BW could not, alone, change the eddy viscosity to a noticeable degree. However, only in conjunction with LC, could BW change eddy viscosity profiles and even decrease eddy viscosity. Moreover, LC together with BW could influence ocean mixing strength and mixed layer depth.

Discussion and Conclusions
As previously highlighted, most traditional turbulence models do not consider the effects of LC and BW in their simulations. Furthermore, as model resolutions increase, existing parameterizations need to be more refined. Consequently, in our research, the LC and BW are added to the large eddy simulation to examine the wave-induced mixing and its parameterization. Six kinds of experiments are carried out so that the effects of BW and LC could be examined singularly, in addition to their joint effects. Through these experiments, it is found that: (1) The presence of LC could generate stripes of strong vertical velocities to a certain depth, BW could produce small-scale turbulences on the sea surface and could weaken LC near the surface, and LC can influence the mixing to a greater depth than BW. Therefore, it is proposed that the vortex force is a much more effective mechanism in deepening the mixed layer than BW. (2) It is also identified that as wave heights increase, stripes become more evident. When wavelengths increase, by contrast, vertical velocities are not modified significantly, thus suggesting that the wave-induced mixing is more sensitive to the wave height than to the wavelength. in accordance with previous findings [56]. Additionally, although it is identified that BW could not influence eddy viscosity alone, BW working in tandem with LC leads to reductions in eddy viscosity. The eddy viscosity of polychromatic waves is larger than that of monochromatic waves. Wave-induced mixing could influence the eddy viscosity and enhance the mixing efficiency.
Based on the above conclusions, it is estimated that a mean eddy viscosity assumption in the presence of Langmuir cells is not enough. The wave-induced mixing should be considered in the eddy viscosity assumption. A new parametric relationship in the upper layer should be derived to describe this phenomenon more accurately. As only a few studies which cover this topic exist [55], parameterized results still contain errors. This problem will be addressed in a future publication in which a new parameterization which relates to the Langmuir number, the Stokes drift and the mixed layer depth will be derived. Stratification in the mixing layer and the application of the new parameterization will also be analyzed.

Discussion and Conclusions
As previously highlighted, most traditional turbulence models do not consider the effects of LC and BW in their simulations. Furthermore, as model resolutions increase, existing parameterizations need to be more refined. Consequently, in our research, the LC and BW are added to the large eddy simulation to examine the wave-induced mixing and its parameterization. Six kinds of experiments are carried out so that the effects of BW and LC could be examined singularly, in addition to their joint effects. Through these experiments, it is found that: (1) The presence of LC could generate stripes of strong vertical velocities to a certain depth, BW could produce small-scale turbulences on the sea surface and could weaken LC near the surface, and LC can influence the mixing to a greater depth than BW. Therefore, it is proposed that the vortex force is a much more effective mechanism in deepening the mixed layer than BW. (2) It is also identified that as wave heights increase, stripes become more evident. When wavelengths increase, by contrast, vertical velocities are not modified significantly, thus suggesting that the wave-induced mixing is more sensitive to the wave height than to the wavelength. (3) The Stokes drift velocity and the vertical velocity increase with the wind speed with stripes becoming more evident as well. The influence of LC increases with wind speed. (4) With increasing depth, the number of stripes declines and stripes become more coherent. (5) Model results and calculations indicate that LC increases eddy viscosity within the mixed layer, in accordance with previous findings [56]. Additionally, although it is identified that BW could not influence eddy viscosity alone, BW working in tandem with LC leads to reductions in eddy viscosity. The eddy viscosity of polychromatic waves is larger than that of monochromatic waves. Wave-induced mixing could influence the eddy viscosity and enhance the mixing efficiency.
Based on the above conclusions, it is estimated that a mean eddy viscosity assumption in the presence of Langmuir cells is not enough. The wave-induced mixing should be considered in the eddy viscosity assumption. A new parametric relationship in the upper layer should be derived to describe this phenomenon more accurately. As only a few studies which cover this topic exist [55], parameterized results still contain errors. This problem will be addressed in a future publication in which a new parameterization which relates to the Langmuir number, the Stokes drift and the mixed layer depth will be derived. Stratification in the mixing layer and the application of the new parameterization will also be analyzed.