## Abstract

The present study performs a wave-resolving simulation of wind-driven currents under monochromatic surface gravity waves using the latest nonhydrostatic free-surface numerical model. Here, phase speed of the waves is set much greater than the current speed. Roll structures very similar to observed Langmuir circulations (LCs) appear in the simulation only when both waves and down-wave surface currents are present, demonstrating that the rolls are driven by the wave–current interaction. A vorticity analysis of simulated mean flow reveals that the rolls are driven by the torque associated with wave motion, which arises from a correlation between wave-induced vorticity fluctuation and the wave motion itself. Furthermore, it is confirmed that the wave-induced torque is very well represented by the curl of the vortex force (VF), that is, the vector product of mean vorticity and Stokes drift velocity. Therefore, it is concluded that the simulated rolls are LCs and that the wave effects are well represented by the VF expression in the present simulation. The present study further revisits the scaling assumptions made by previous studies that derived VF formulation and shows that there is disagreement among the previous studies regarding the applicability of VF formulation when the wave orbital velocity (proportional to the amplitude times the frequency) is much smaller than the mean flow velocity. The result from the present simulation shows that the VF expression is still valid even with such small wave amplitudes, as long as phase speed of the waves is much greater than the current speed.

## 1. Introduction

Among various processes that occur in the oceanic surface boundary layer, Langmuir circulations (LCs) receive attention because they are believed to modulate the air–sea heat/material exchange through an enhancement of turbulent mixing (Langmuir 1938; Smith 1992; Belcher et al. 2012; D’Asaro 2014; Li et al. 2016). In the formation of LCs, the interaction between surface waves and mean flow is believed to play a central role. Craik and Leibovich (1976, hereinafter CL76) derived an expression for the wave–current interaction called the vortex force (VF) and showed that the VF induces roll structures similar to observed LCs. Today, a number of studies investigate the effects of LCs on mixed layer turbulence using large-eddy simulations (LES) in which the wave–current interaction is represented by the VF in the momentum equation (e.g., Skyllingstad and Denbo 1995; McWilliams et al. 1997). The VF representation is popular because it enables us to simulate wave effects under the rigid-lid upper boundary condition, with which the difficulty in simulating nonhydrostatic surface waves is avoided.

A minimal form of the wave-averaged momentum equation based on CL76 is

where *t* is time, **u** = (*u*, *υ*, *w*) is velocity, *p* is pressure, **u**^{St} is the Stokes drift velocity associated with surface waves, is vorticity, *ρ* is the density of water, and *ν* is the eddy viscosity coefficient. An overbar indicates Eulerian averaging over a wave period. Wave effects on mean flow are represented as external forces, namely, the gradient of the Bernoulli head Π and the vortex force term . Since Bernoulli head gradient cannot affect vorticity, the VF is the key term for LCs. In many applications of this formulation, **u**^{St} is prescribed using statistical parameters of waves (e.g., amplitude and wavenumber) and the irrotational linear solution of surface waves. Hereinafter, we refer to the wave-averaged momentum equation that includes the wave effect as the VF as the Craik–Leibovich (CL) equation.

When wind-driven surface current is aligned with the wave propagation direction (i.e., down-wave surface current is present), the VF causes dynamic instability, and the instability forms rolls aligned with the wind-wave direction (Craik 1977; Leibovich 1977, 1983). This driving mechanism is called the CL2 mechanism and is widely accepted as a major driving mechanism of LCs.

The CL equation has been used in various contexts to describe the effects of surface waves on mean flow. As written above, LES based on the CL equation has become a common approach to study the mixing effects of LCs and to develop or verify mixed layer turbulence models that account for the impact of LCs (e.g., Skyllingstad and Denbo 1995; McWilliams et al. 1997, 2012, 2014; Noh et al. 2004; Li et al. 2005; Harcourt and D’Asaro 2008; Van Roekel et al. 2012; Tejada-Martinez et al. 2012; Harcourt 2013, 2015; Reichl et al. 2016), and it is reported that some of these parameterizations represent better the observed spatial distribution of mixed layer depth (e.g., Li et al. 2016). Uchiyama et al. (2009, 2010); Weir et al. (2011) used a generalized form of the CL equation derived by McWilliams et al. (2004, hereinafter MRL04) to simulate littoral currents and rip currents. Hamlington et al. (2014) and Suzuki et al. (2016) investigated the effect of surface waves on submesoscale phenomena using a variant of the CL equation.

However, there remains an important issue: neither observations, laboratory experiments, nor numerical simulations have yet provided direct verification of the CL equation, a formulation arrived at variously through complex mathematical approaches (e.g., CL76; Andrews and McIntyre 1978; Holm 1996), and the CL2 mechanism, upon which the prevailing understanding of LCs is based. The lack of direct laboratory or numerical validation of the CL equation has resulted in a series of discussions regarding its validity, and a consensus view on this issue has not been reached yet (Mellor 2003, 2016; Ardhuin et al. 2008a, 2017; Aiki and Greatbatch 2012, 2013, 2014). Since this argument is mostly over the theoretical aspect of the wave–current interaction, experimental approaches, such as observations, laboratory experiments, or numerical simulations, that do not invoke the approximations and assumptions behind the CL equation can provide good evidence to support or refute the CL theory and eventually deepen our understanding of the wave–current interaction, the driving mechanism of LCs, and, therefore, ocean surface mixing.

These experimental approaches have not been commonly taken because each approach has different issues. In field observations, the measurement of waves and wave-averaged flow is not sufficiently accurate for comparing the actual wave effect on the averaged flow and the wave-related terms in the CL equation. As for laboratory experiments, the resulting flow in small tanks cannot be free from the effect of the side and bottom walls of wave tanks, which leads to complications for comparing results with LCs in the open ocean. Melville et al. (1998) and Veron and Melville (2001) conducted a series of careful laboratory experiments with a large wind-wave channel to examine mechanisms of LCs with negligibly small side/bottom effects, but the velocity scale of the resulting flow was comparable to the phase speed of waves. Such a situation usually only occurs for waves with *O*(0.1) m wavelengths, so the direct comparison with LCs in the typical open ocean condition is again difficult. Numerical simulations that explicitly simulate wave motions, namely, wave-resolving simulations (WRS), can be a possible approach, but numerical models that can explicitly simulate three-dimensional turbulence and deep-water waves have not been developed until recently. Takagaki et al. (2015) reported that they simulated LC-like structures in a direct numerical simulation of air–water multiphase flow. Tsai et al. (2017) investigated the LC-like flow obtained from a direct numerical simulation of free-propagating surface waves and concluded that the flow was driven by the CL2 mechanism. Tsai et al.’s (2017) conclusion is based on a comparison of the simulated results with the linear stability analysis results from the streamwise-homogeneous CL equation. However, the consistency between the simulated wave force and the CL formulation (i.e., VF and Bernoulli head), which is a controversial issue in today’s community (e.g., Mellor 2016; Ardhuin et al. 2017), still remains uninvestigated.

There is another issue: Given that the CL equation is valid, is the equation always an appropriate way to describe the wave–current interaction in the upper ocean? Although it is used in a variety of situations (e.g., mixed layer turbulence, coastal circulations, submesoscale phenomena), less care is given to the assumptions under which the equation is derived. In fact, its derivations include, explicitly or implicitly, some assumptions about spatiotemporal scales of waves and mean flow. For example, CL76 assumed that the characteristic spatial scale of mean flow is comparable to that of waves (wavelength) and that the characteristic velocity scale of mean flow is comparable to Stokes drift velocity. Evidently, this scaling is not satisfied by all of the various phenomena mentioned above.

The present study intends to address these questions. First, using a nonhydrostatic free-surface numerical model developed recently, a wave-resolving simulation of surface gravity waves and mean flow is conducted to clarify if the wave effect is consistent with the CL equation. Since the mathematical argument starting from Mellor (2003) is still going on, we consider it beneficial to support or refute (whichever the result is) the validity of the CL formulation from a different perspective: an experimental approach. Detailed dynamical analyses are conducted to study the driving mechanism of simulated flow, and it is explained from the Eulerian point of view (in contrast to the CL theory, which is based on the Lagrangian point of view). Consistency between the simulation results and the CL theory is found. Then, in order to clarify the scales of waves and mean flow for which the CL equation has been derived, a systematic description of spatiotemporal scaling assumptions in previous studies is provided. It is shown that there is disagreement among the previous studies regarding the applicability of the VF formulation when the amplitude of waves is small. Using the WRS results, it is demonstrated that the VF formulation is still valid even when the wave orbital velocity is much smaller than the characteristic velocity scale of mean flow as long as the latter is much smaller than the wave phase speed.

Throughout this paper, we focus on the situations where phase speed of the waves is much greater than the current velocity scale, and hence the wave motions are approximately linear and irrotational, as most applications of the CL equation to LCs in the open ocean assume. These assumptions are common because in many cases ocean surface waves are well described with the linear irrotational theory. Furthermore, in such cases, the Stokes drift profile can be prescribed using the well-known linear analytic solution. Note that interaction of currents with rotational waves, which is relevant to waves with smaller phase speed (wavelength) and plays important roles in local air–sea interaction processes (e.g., Craik 1982; Melville et al. 1998; Veron and Melville 2001; Phillips 2005), is not the scope of the present study.

First, simulation settings are described in section 2, and the results and analyses are shown in section 3. Then, in section 4, scaling assumptions in previous studies are systematically examined to clarify the regimes (in terms of relative spatial and temporal scales of waves and mean flow) in which the VF has been derived, and the validity of one of these assumptions is examined using the WRS results. Finally, the findings of the present study are summarized in section 5.

## 2. Problem setup

To directly investigate the driving mechanism of LCs, we designed a wave-resolving simulation of surface waves and currents in a configuration similar to a wave tank with a wave generator in one end. For this purpose, we used a nonhydrostatic free-surface version of the so-called kinaco model (Matsumura and Hasumi 2008) in which nonhydrostatic pressure and free-surface elevation are solved simultaneously as in Casulli and Zanolli (2002). This scheme and the efficient Poisson–Helmholtz equation solver (Matsumura and Hasumi 2008) make it possible to correctly simulate the interaction of deep-water waves and currents with reasonably limited computational resources ( appendix A).

To contrast the role of waves, experiments with and without waves are conducted. The control run is the basic experiment, in which surface gravity waves, as well as wind-driven currents, are explicitly simulated. The no-wave run is the same as the control run except that surface waves are not forced. In both experiments, system rotation (Earth rotation) is not considered, and an incompressible fluid of uniform density is considered in order to focus on LC dynamics. For the same reason, uniform and isotropic eddy viscosity is assumed.

Consider a fluid in a rectangular water tank of streamwise, spanwise, and vertical dimensions of 16*λ*_{10}, 4*λ*_{10}, and ≥2*λ*_{10}, respectively. Here, *λ*_{10} indicates the wavelength of the deep-water wave with a period of 10 s, which is approximately 156 m. Cartesian coordinates (*x*, *y*, *z*) are assigned, where the *x*, *y*, and *z* axes are aligned with the streamwise, spanwise, and upward directions, respectively (Fig. 1). The reference level (*z* = 0) is set so that the bottom of the tank is at *z* = −2*λ*_{10}.

The computational domain is divided into three subdomains, namely, the forcing domain at 0 ≤ *x* ≤ 2*λ*_{10}, the free propagation domain at 2*λ*_{10} ≤ *x* ≤ 12*λ*_{10}, and the damping domain at 12*λ*_{10} ≤ *x* ≤ 16*λ*_{10}. Waves are generated in the forcing domain using external forcing terms, propagate along the *x* axis through the free-propagation domain, and are finally damped in the damping domain using external forcing terms.

The equation of motion and the continuity equation are

Here, **g** = (0, 0, *g*) (*g* = 9.8 m s^{−2}) is the gravitational acceleration, and **F** represents external force that is used to artificially generate and damp wave motion. The reference density *ρ* is 1020 kg m^{−3}, and the eddy viscosity coefficient *ν* is 0.01 m^{2} s^{−1}. This choice of eddy viscosity coefficient suppresses small-scale turbulence and allows the separation of mean and wave motion, and thereby we can focus on the wave–current interaction processes. As will be seen later (e.g., Fig. 4), the Reynolds number of the resulting mean flow is roughly *O*(10^{−3}).

The top of the computational domain is a free surface *η*, which follows

Here, *F*_{η} is a forcing term for *η*, which corresponds to water flux. Vector quantities with a subscript of *H* indicate the projection onto the horizontal plane. The atmospheric pressure is constant. Wind stress is imposed at the surface. It is (*τ*_{x}, *τ*_{y}) = (0.02, 0) N m^{−2} throughout most of the domain, except in the forcing domain, where small noise is added to *τ*_{y} to seed instability. The bottom of the domain is flat and free-slip. The side boundaries (*y* = 0 and 4*λ*_{10} are periodic. Boundary conditions at the upstream (*x* = 0) and downstream (*x* = 16*λ*_{10}) ends are related to the wave forcing described below.

In the forcing domain, the external forces **F** and *F*_{η} are imposed to force wave motion. These forcings serve as a wave generator in the wave tank. (Note that wave motion is forced not by wind but by external forcing in the present study.) They are given by the following:

These forces restore the velocity and the free surface to the orbital motion of linear surface gravity waves, , and . The coefficients, and , control the strength of the forcing. They vary in the *x* direction as shown in Fig. 1.

The forced wave motion is of the deep-water wave with an amplitude of and a period of 10 s. Therefore,

Here, *σ*_{10} = 0.63 s^{−1} and *k*_{10} = 2*π*/*λ*_{10} = 0.040 m^{−1} are the angular frequency and the wavenumber of the forced wave, respectively, which are calculated from the dispersion relation of linear deep-water waves. The amplitude of the forced wave *a*_{force} is 0.5 m in the control run and 0 m in the no-wave run.

At the upstream end of the domain, *x* = 0, the velocity component normal to the boundary, is set to match the orbital motion forced in the forcing domain:

For *η*, a free-end condition is imposed there:

Although this boundary condition does not match the orbital motion of the deep-water waves, the resulting noise is damped as waves pass through the forcing domain. Inside the free-propagation domain, the leading-order motion of surface waves is the freely propagating deep-water wave with a period of 10 s. In the no-wave run, *x* = 0 is effectively a solid wall.

In the damping domain, nonzero external force is imposed so that the surface waves will be artificially damped to avoid reflection at the downwind end. Here, the external force restores the velocity and the free surface to the state of rest:

The magnitudes of *χ*_{u} and *χ*_{η} in the damping domain (Fig. 1) are tuned so that either reflection at the downstream end (which is a free-slip solid wall) or reflection by the damping force will hardly occur. Note that *χ*_{u} and *χ*_{η} are zero throughout the free propagation domain.

### Numerical configuration

The numbers of grid points in the (*x*, *y*, *z*) direction are (1024, 256, 128), respectively. The grid spacing is uniform (*λ*_{10}/64) in the *x* and *y* directions. The vertical grid spacing is uniform (*λ*_{10}/64) except for the top layer, where it changes according to the surface elevation (*λ*_{10}/64 + *η*).

The simulations are started from rest (**u** = 0) and a flat surface (*η* = 0) and integrated until the rms of (defined later) at the *y*–*z* section that we use in the analysis becomes quasi steady (*t* = 36 h). In the control run, the time step of integration is 0.5 s, which resolves one wave period with 20 steps. In the no-wave experiment, the time step of 2 s is used because we need not resolve wave motion.

## 3. Results

### a. Characteristic temporal scales of simulated flows

To clarify what are “mean flow” and “waves” in the WRS result, we first look into the characteristic temporal scales of the flow. To visualize the temporal scales, the normalized energy density spectrum of *w*, computed from a time series of *w* at a fixed point, is shown in the lower-right corner of Fig. 2. It can be seen that there is a clear gap between the temporal scales of low-frequency flow (continuous spectrum in Hz) and high-frequency flow (spikes at 0.1, 0.2, 0.3 … Hz). The former corresponds to the roll-like flow described later, and the latter corresponds to 10-s-period wave orbital motion and its harmonics.

The separation in time scales motivates us to filter low-frequency and high-frequency motions by introducing a wave-period Eulerian average :

Here, *T*_{10} = 10 s is the period of the forced waves, and *t*_{a} represents the time at which the analysis is conducted (in this case, *t*_{a} is 36 h). The integral is obtained from the sum of the integrand at every time step during the average period. Hereinafter we refer to this 10-s-average as “wave average” or “wave mean.”

### b. Overview of simulated flows

A three-dimensional sketch of the flow simulated in the control run is shown in Fig. 2. The colored lines are the streamlines of the wave-averaged flow at time *t* = 36 h (computed by simulating trajectories through a frozen field). Multiple streamwise bands of downwelling jets are observed. On both sides of these bands, the streamlines form helixes, showing the presence of roll cells. Superposed above these streamlines, the transparent surface shows the instantaneous surface elevation *η*.

Figure 3 shows the horizontal section of velocity () near the surface (*z* ≈ −0.15*λ*_{10}) and the distribution of 1024 surface trapped passive particles uniformly released at *t* = 33 h in the unforced region (2*λ*_{10} ≤ *z* ≤ 12*λ*_{10}) in the control run. Note that wave orbital motion is filtered by wave averaging.

In Fig. 3, some features of the observed LCs (see, e.g., CL76) can be identified. There is an array of downwelling jets aligned with the wave-propagating direction. Downstream surface current velocity is strong over the downwelling jets. Moreover, the position of surface convergence corresponds to that of the downstream current maximum and of the downwelling jets. In the upstream part of the free-propagation domain (), flow patterns are nearly two dimensional (i.e., uniform in *x*). Downwelling jets form a pattern aligned with the wind and the waves, with spanwise intervals of ~0.5*λ*_{10}. As the rolls develop while drifting downstream, they merge, start to bifurcate, and become less uniform in *x*.

The horizontal section of the no-wave run result is also shown in Fig. 3. In the absence of waves, no roll structure is seen (flow is nearly uniform in the *y* direction). It should also be noted that the roll structure is not found in the experiment with wind forcing with the same magnitude as and in the opposite direction to that of the control run (not shown). These results demonstrate that the presence of both waves and down-wave surface currents is crucial for the LC-like rolls in the control experiment.

Figure 4 shows the *y*–*z* sections of the flow of the control run at *x* = 6*λ*_{10}. Including this figure, all the *y*–*z* slices shown below are smoothed in *x* direction over one wavelength to remove small phase-dependent noise, unless noted explicitly. The flow structure in Fig. 4 again shows the features of observed LCs: the downstream velocity maximum is closely associated with downwelling jets, and the maximum downwelling velocity (~0.02 m s^{−1}) is comparable to the mean downstream drift speed at the surface (~0.04 m s^{−1}). Furthermore, clear roll structures can be seen in the cross section. The aspect ratio of the roll cells is about unity. The magnitude of spanwise velocity is comparable to that of vertical velocity, in marked contrast to purely shear-driven turbulence.

Hereinafter, we shall call the rolls observed in the control run “LCs”, since 1) the downwind waves and wind-driven currents are crucial for generating the rolls and 2) the rolls have multiple features that are commonly considered to be those of LCs.

### c. Driving mechanism of simulated LCs

To investigate the driving mechanism of the LCs in the control run, the mean vorticity equation is examined.

In Eq. (2) in the free-propagation domain where the forcing terms are absent, the wave effects on mean flow appear in the nonlinear term. In order to extract the wave effects on mean flow, we consider a separation of the nonlinear term using the wave average:

Temporal averaging is defined in Eq. (5), and indicates deviation from the temporal average. This deviation is dominated by wave orbital motion and contains nearly no turbulent fluctuation as shown in Fig. 2 and thus can be called “wave fluctuation.” Therefore, contains wave effects only, so it can be regarded as the rotational component of the Reynolds stress divergence associated with waves. [Note that the flow is laminar in the no-wave case (Fig. 3), so all high-frequency fluctuations and associated stress in the control run can be attributed to wave effects.] Likewise, is the nonlinear term arising from the mean flow.

Taking the curl of the averaged Eq. (2) after omitting the forcing term leads to the vorticity equation for mean flow:

where the decomposition in Eq. (6) is applied.

Figure 5 shows and terms in the *x* component of Eq. (7), along the vertical slice at *x* = 6*λ*_{10}. The lhs term of Eq. (7) is evaluated as the sum of the rhs terms. Since the time derivative term (the lhs term; Fig. 5b) is very small, the mean flow nonlinear term (the first term of the rhs; Fig. 5c), the wave nonlinear term (the second term; Fig. 5d), and the viscosity term (the third term; Fig. 5e) are nearly balanced. Figure 6 shows and some important terms in the *z* component of Eq. (7), at the same slice as Fig. 5. We found that the rhs terms of the *z* component of Eq. (7) were again nearly balanced, and only major two terms, the mean flow nonlinear term (Fig. 6b) and the viscosity term (Fig. 6c), are shown. It is noteworthy that has a positive value where is also positive (Figs. 5a and 6a), reflecting the common structure of LCs (Fig. 4). We first look into the *x*-component terms of Eq. (7) in detail and then the *z*-component terms.

In the downwelling jet region, the mean flow nonlinear term , where indicates the *x* component of the vector quantity, shows a pattern that has the opposite sign to the pattern of streamwise vorticity near the surface and the same sign at the subsurface levels (Figs. 5a,c). This can be interpreted as the vertical transport of vorticity due to the downwelling jets of the LCs. The wave nonlinear term shows a pattern that has the same sign as vorticity (Figs. 5a,d), indicating that it reinforces the circulation. Its magnitude is largest near the surface, where wave motion is strongest. The viscosity term shows a pattern that has the opposite sign to vorticity (Figs. 5a,e), which corresponds to the viscous diffusion of vorticity. Therefore, the simulated LCs are driven by torque imposed by the wave-associated Reynolds stress. Hereinafter we call the term “wave torque.”

The wave torque term can be rewritten as

Here we have used the incompressibility of wave component velocity . Among the six terms on the rhs of Eq. (8), it is found that the three terms , , and explain most of the variation of the wave torque (Figs. 5d,f). The sum of the three terms and the total wave torque coincide with a correlation coefficient of 0.98 and a regression line slope of 0.94, over the slice shown in Fig. 5.

To understand how these three terms force mean streamwise vorticity, vorticity fluctuation ** ω**′ is examined. The temporal evolution of the vorticity fluctuation component is described by the following equation:

We found that the rhs of Eq. (9) is dominated by the first term in the parentheses, that is, advection/stretching/tilting of the mean vorticity by wave motion ( appendix B).

This and the phase relation of deep-water waves together lead to the following interpretation of the driving terms, , , and (Fig. 7). They can all be understood as the correlation between the leading-order wave orbital motion and wave-induced vorticity fluctuation ** ω**′, and they all supply positive when and vice versa, which is consistent with the vorticity patterns seen in Figs. 5 and 6.

(Fig. 7a): The fluctuation of streamwise vorticity that is in phase with arises through the term , that is, tilting of the mean vertical vorticity by wave motion.

^{1}When , the correlation between and leads to a downward flux of the positive streamwise vorticity. Since the magnitude of the wave motion decays with depth, the flux converges, resulting in the mean supply of positive .(Fig. 7b): In the same way as in (i), the variation in phase with is induced by the tilting of the mean vertical vorticity. When , the forward-tilted vorticity () is stretched by wave motion, while the backward-tilted vorticity () is squeezed. This leads to a net stretching of the streamwise vorticity, resulting in the mean supply of positive .

(Fig. 7c): The variation of in phase with arises from the term , that is, the stretching of mean vertical vorticity. When , the stretched vertical vorticity () is tilted forward by wave motion, while squeezed vertical vorticity () is tilted backward. This leads to a net forward-tilting of vertical vorticity, resulting in the mean supply of positive .

As shown above, positive (negative) streamwise wave torque arises where mean vertical vorticity is also positive (negative), and the WRS result shows a pattern of vertical vorticity favorable for reinforcing preexisting rolls. Then how does the vertical vorticity maintain such a spatial pattern? The rhs terms of the *z*-component vorticity in Eq. (7) (Figs. 6b,c) show that the mean nonlinear term supplies in such a pattern and is balanced with viscous diffusion. The mean nonlinear term is expanded as follows:

The *y*-to-*z* tilting term , the stretching term , and the flux term are shown in Fig. 6. Tilting and stretching contribute to intensify near the surface (Figs. 6d,e), while the flux term represents downward transport of by downwelling jets (Fig. 6f). The *x*-to-*z* tilting term is insignificant and not shown. Since the only source of vorticity is the wind-driven shear , is induced by tilting of , intensified by vertical stretching, and advected downward by downwelling jets.

### d. Comparison with the CL theory

Are the WRS results and their interpretation consistent with the prevailing understanding of LCs based on the CL theory?

The CL equation states that the evolution of mean flow vorticity is described by

In other words, it states that the wave torque is expressed as the curl of the VF, .

In Fig. 5, at the same slice (*x* = 6*λ*_{10}) is shown. Here, is estimated using the linear analytic expression with the wave parameters of the forced wave:

The analytic expression of Stokes drift velocity approximates well the directly evaluated Stokes drift velocity (Lagrangian minus Eulerian mean velocity) as shown in appendix A. Evidently, the curl of the wave nonlinear term is represented very well by the curl of the VF. The wave torque and the VF curl agreed with a correlation coefficient greater than 0.99 and a regression line slope of 1.14. They also agree well at other positions in *x* (not shown). These results demonstrate that the CL equation (VF formulation) is valid in the present configuration.

In fact, the flow simulated using the CL equation (with VF forcing) with the wave parameters in the control run instead of using the full Eqs. (2)–(4) produced qualitatively and quantitatively similar flow to the LCs in the control run (not shown).

The expansion of with the prescribed form of yields

Figure 5 shows that the latter component is dominant. Therefore, the forcing mechanism of streamwise vorticity can also be interpreted as the down-wave tilting of the mean vertical vorticity by Stokes drift velocity, which is consistent with the CL2 mechanism (Craik 1977; Leibovich 1977, 1983). This can be recognized as an interpretation based on the Lagrangian point of view, in contrast to the Eulerian interpretation (i.e., the torque induced by the wave-associated Reynolds stress) presented in the previous subsection.

As for vertical vorticity , the WRS result shows that spanwise to vertical tilting of vorticity plays an important role in inducing aligned with and hence in developing LCs (Fig. 6). We have confirmed that the tilting is dominant over the other terms in the initial developing phase of LCs, as the CL2 instability theory suggests. In addition, the present result shows that vertical stretching of vorticity is also important in maintaining when the LCs have fully developed.

In his criticism of VF formulation, Mellor (2016) argued that the VF is much smaller than other wave stress terms and thus negligible, because VF terms are scale,^{2} the wave stress terms are , and is much smaller than *c*. Here, *a*, *k*, *c*, and *L* are amplitude, wavenumber, phase speed, and length scale of wave-averaged flow, respectively.

However, this does not mean that the VF terms are negligible. The terms that are retained by Mellor are nonrotational components of the vector field, and it is the terms that drive rotational motion. In fact, the WRS result shows that the wave torque scales as the terms neglected by Mellor. In the WRS result, *a* = 0.5 m, *k* = 0.040 m^{−1}, *c* = 15.6 m s^{−1}, *L* ≈ *k*^{−1} = 25.0 m, and m s^{−1}. The measured wave torque is *O*(10^{−7}–10^{−6}) s^{−2} (Fig. 5d). The VF curl scales as s^{−2}; whereas the terms retained by Mellor scale as s^{−2}.^{3}

## 4. Scaling assumptions of the CL equation

### a. Scaling assumptions by preceding studies

In this section, we first revisit the scaling assumptions made by previous studies that derived the CL equation in order to clarify situations in which the CL equation is known to be valid. As noted before, we shall restrict our attention to the cases where waves are linear and irrotational. To describe the scales of the wave and mean flow quantities, we first introduce a set of notations as in CL76, MRL04, and Aiki and Greatbatch (2014). To begin with, suppose that velocity field is decomposed into the “mean flow” component and the “wave” component .^{4}

Let the characteristic wavenumber, angular frequency, and wave slope (the product of wave amplitude and wavenumber) of the surface waves be denoted by *k*, *σ*, and *ε*, respectively. Assuming deep-water waves, the scales of the phase speed, amplitude, and orbital velocity follow *σk*^{−1}, *εk*^{−1}, and *εσk*^{−1}, respectively. The temporal scale of waves is *σ*^{−1}, and the horizontal and vertical spatial scales of wave quantities are both *k*^{−1} (because of the characteristics of deep-water waves).

Next, we introduce *γ*, *δ*_{H}, and *δ*_{V} as nondimensional numbers that represent the ratio of the temporal/horizontal/vertical scales of the mean flow and waves:

Then, the temporal, horizontal, and vertical scales of the mean flow quantities are expressed as , , and , respectively.

Let *α*_{H} and *α*_{V} be the ratios of the horizontal and vertical velocity scales of mean flow to the phase speed of waves, respectively. They can be considered as the Froude numbers based on the horizontal/vertical velocity of mean flow. The horizontal and vertical velocity scales of mean flow are expressed as *α*_{H}*σk*^{−1} and *α*_{V}*σk*^{−1}, respectively. The continuity equation leads to

Scaling the temporal change rate of mean flow with advection by mean flow itself leads to

Now the characteristic scales of the wave and mean flow quantities can all be described with four nondimensional numbers (*ε*, *γ*, *δ*_{H}, *δ*_{V}), the angular frequency *σ*, and the wavenumber of surface waves *k*. Among the four nondimensional numbers, *ε* characterizes the nonlinearity of waves, and the other three represent the ratio of the spatiotemporal scale between mean flow and waves. The notation introduced here is listed in Table 1.

Since the temporal scale of mean flow is longer than wave period, the following relationship is required for all the derivations of the VF:

#### 1) Applicability of irrotational linear wave solution

First, let us consider under which conditions the wave motion is linear and irrotational. For waves to be linear, . Scaling each term using our notation leads to *ε* ≪ 1. The momentum equation linearized about the mean flow field is

Here is a pressure field that corresponds to , and viscosity is assumed to be sufficiently small.

The leading-order motion of waves is irrotational when the dominant balance of Eq. (16) is between the time derivative term (the first term on the lhs) and the pressure gradient term (the rhs). Therefore, for the *x* component of Eq. (16),

Rewriting , and leads to the irrotationality conditions (together with the linearity condition):

#### 2) CL76

Assuming LCs as , CL76 suggested the following scale assumptions:

In CL76, the condition in Eq. (18a) is implicitly assumed by evaluating the spatial derivative of the wave and mean flow quantities using the same scale. Substitution of these relationships into Eq. (14) yields *γ* = *ε*^{2}. Note that the same scaling is assumed by Aiki and Greatbatch (2014), who showed that the CL equation can be derived in the vertically Lagrangian, horizontally Eulerian coordinate system (Mellor 2003, 2008; Aiki and Greatbatch 2012).

#### 3) L80

Leibovich (1980, hereinafter L80) demonstrated that mean flow equation derived from the generalized Lagrangian mean (GLM) theory (Andrews and McIntyre 1978) is equivalent to the CL equation of CL76 under certain approximations. L80 states that, when the following requirements are satisfied, the CL equation is an appropriate approximation to GLM equation:

When the second condition is satisfied, the pseudomomentum that appears in the GLM equation is identical to the Stokes drift velocity, reducing the equation to the CL equation. L80 requires the third and fourth conditions to ensure the irrotational wave field to hold. Using our notation, the fourth requirement is rewritten as

#### 4) MRL04

MRL04 derived a set of equations that comprehensively describes the interaction among coastal circulations, “long” gravity waves, and wind waves. In this system, wind-wave effects on coastal circulations are expressed in terms of the Bernoulli head gradient and VF. They assumed the following scale relationships between coastal circulations (corresponds to ) and wind waves :

#### 5) ARB08

Ardhuin et al. (2008b, hereinafter ARB08) made use of GLM theory to derive a set of equations for mean flow over varying bottom topography. In their derivation, they made various assumptions to justify their use of the irrotational linear solution as the leading-order wave motion. Although their assumptions are more specific, we will represent them with the irrotational linear solution applicability [Eqs. (17a)–(17d)] to simplify the comparison with other studies. For the same reason, their assumptions regarding topography are not considered here. In addition to these conditions, they assumed mean flow to be hydrostatic. Therefore, their assumptions can be expressed as

#### 6) Scaling diagram

The spatiotemporal scaling assumptions of the previous studies are illustrated in Fig. 8. Here, *γ*^{−1} and are chosen for the horizontal and vertical axes, among the four nondimensional parameters that characterize the relationship between wave and mean flow; and scaling diagrams for different choices of the characteristic wave slope *ε* are shown. Assuming LCs as the characteristic mean flow, *δ*_{V} = 1 is chosen here (i.e., the vertical scale of mean flow is comparable to wavelength). The orange markers indicate the scale assumed in CL76 and the scale assumed in MRL04. The area shaded with red satisfies the requirements of L80, while the blue area satisfies the requirements of ARB08. The irrotational linear solution of surface waves is applicable in the lower right area bounded by the gray broken line. In these diagrams, the relationship *γ*^{−1} ≫ 1 is plotted as *γ*^{−1} > 10. The same rule is applied for *γ*/*δ*_{H} ≪ 1 and so on. Figure 8 indicates that L80 and ARB08 cover a wide range of regimes. This is an advantage of the GLM theory; that is, no explicit choice of spatiotemporal scale is needed.

### b. Vortex force in small-amplitude cases

As shown in Fig. 8, the assumed spatiotemporal scales, defined in terms of wave slope *ε*, change with wave amplitude. Notably, the area that meets the requirements of L80 shrinks as the amplitude becomes small, because L80 requires much greater orbital velocity than mean flow velocity for irrotational wave field to hold. In contrast, the conditions in Eq. (17) suggest that the irrotationality of wave field is not restricted by small wave amplitude, as long as the Froude number (*γ*/*δ*_{H,V}) is small. This implies that the requirements by L80 [Eqs. (19a) and (19b)] are too strict especially when wave amplitude is small. In addition, ARB08 does not require mean flow velocity to be smaller than orbital velocity. Therefore, according to ARB08, the CL equation remains applicable even when wave amplitude is small.

Here arises a disagreement regarding the applicability of the CL equation when wave orbital velocity scale is smaller than or comparable to characteristic mean flow velocity scale. It is beneficial to investigate this disagreement and clarify whether the CL equation (especially the VF formulation) is valid in such cases because such conditions can be seen in waves interacting with frontal or coastal currents, whose speed can be as large as wave orbital velocity, while current speed is much smaller than wave phase speed.

We will take an experimental approach to investigate the above disagreement since our argument leading to Eq. (17) is based on simple scale comparison of terms and thus rather crude. To achieve this, a simulation (hereinafter referred to as the small-amplitude experiment) is carried out with an orbital velocity scale smaller than the characteristic velocity scale of mean flow. The experiment is initialized with the flow field at *t* = 36 h of the control run in the previous section, where LCs are active. The same forcing as the control run is imposed, except that *a*_{force} = 0.005 m, which is one-hundredth of the control run. The characteristic velocity scale of mean flow is m s^{−1} (section 3), the orbital velocity scale of the forced wave is 0.003 m s^{−1}, and the phase speed of waves is 15.6 m s^{−1}. Therefore, here the conditions of Eqs. (17b) and (17c) are satisfied, while the requirements of L80 [Eqs. (19a) and (19b)] are not. Time integration is performed for 15 min, and the last 10 s of the integration period is used for time averaging. The LCs do not rapidly decay during the 15 min, so the mean flow velocity scale is still much greater than the orbital velocity scale. The large amplitude waves of the control run are damped during the 15 min, leaving no effect on the analysis result.

To investigate the applicability of the VF expression, the same dynamical analysis as in the previous section is performed. In the small-amplitude case, a “true” wave nonlinear term cannot be obtained by temporal filtering only, because most of the temporal variation is attributed to “trend” or “drift” over the average window, and the small wave motion is overwhelmed. To extract the effects of very weak wave motion , a Fourier analysis is attempted. After subtracting the trend component (assumed to be linear in time), the complex Fourier coefficient for the 10-s-period motion is computed. The wave torque associated with 10-s-period waves is obtained as follows:

where denotes the real part, and superscript asterisk indicates the complex conjugate.

The wave torque obtained from the above analysis is shown in Fig. 9. It agrees very well with the curl of the VF . Therefore, it is demonstrated that the VF expression is applicable even when the orbital velocity scale is not much greater than the characteristic velocity scale of the mean flow, as long as the Froude number is sufficiently small as in the present situation. As noted previously, such situations may be found in frontal or coastal regions. For example, a typical wave orbital velocity scale is *O*(1) m s^{−1}, which is comparable to current speed in strong frontal regions like ACC region, where LCs are suggested to deepen the mixed layer considerably (Belcher et al. 2012).

## 5. Conclusions

In this study, we first performed a wave-resolving simulation of surface waves, wind-induced current, and their interactions using a nonhydrostatic free-surface numerical model. When waves and wind-induced down-wave currents were present, roll structures of mean flow very similar to observed LCs appeared, and they did not appear if either surface waves or down-wave surface currents were absent. From these results, we concluded that the roll structures observed in the present simulation were LCs. Dynamical analysis of the mean streamwise vorticity equation showed that the streamwise vorticity of simulated LCs was forced by the torque induced by the Reynolds stress associated with wave motion, and it was shown how wave motion induced mean rotational forces as waves interacted with the mean vorticity field from the Eulerian point of view. For wave torque to enhance preexisting streamwise vorticity, the spatial patterns of the streamwise and vertical mean vorticity must be aligned in the same location in the spanwise direction. Analysis of the mean vertical vorticity equation clarified that not only spanwise to vertical tilting but also vertical stretching played a central role in maintaining such vorticity patterns. It was also confirmed that the rotational component of the wave-associated Reynolds stress is well represented by the vortex force, demonstrating the validity of the CL equation (VF formulation).

Then, we reviewed the spatiotemporal scaling assumptions made by the studies that had derived the CL equation. Among them, the derivations by Leibovich (1980) and Ardhuin et al. (2008b) that utilized GLM strategy covered a wide range of scales. However, the limitation indicated by Leibovich (1980) could be relaxed (at least for the simulated condition), as demonstrated in the present WRS result. With this relaxation, it was confirmed that VF formulation was possible in conditions where the mean flow was stronger than the wave orbital motion. Such situations may include examples of wave effects on strong frontal flow.

The original derivation of the VF by Craik and Leibovich (1976) showed that, if some scaling assumptions (section 4) were satisfied, the vorticity fluctuation ** ω**′ arises through advection/stretching/tilting of the mean vorticity, and the wave torque could be substituted by the VF curl . These arguments were experimentally confirmed in the dynamical analysis in section 3. Although CL76’s scaling assumptions were roughly satisfied in the control run, it was also confirmed that these still hold in some other scaling regimes, such as in the small-amplitude case and a case with a large

*δ*value that was realized by setting small

*ν*(not shown).

The present study investigated LCs driven by a particular set of wave parameters (period and amplitude) and wind stress. Various sets of wave parameters and wind stress need to be investigated for an extensive understanding of LCs in the real ocean. Notably, further investigation is needed for situations where the mean flow is strong (Froude number is large). Wave parameters in the present study (period of 10 s and amplitude of 0.5 m) led to the phase speed of 15.6 m s^{−1} and the wave slope of 0.02, which resulted in rather fast swells passing over rather weak LCs [~0.02 m s^{−1}, i.e., *α* ~ *O*(10^{−3})]. In reality, the Froude number can be much greater depending on the sea states, which may lead to the inapplicability of irrotational linear wave solutions. In such cases, the mean flow effects on waves must be explicitly taken into account, and the VF formulation may have to be modified, as in Craik (1982) and Phillips (2005). Meanwhile, active LCs tend to reduce shear in the mean field, which may affect the irrotationality of waves. In addition, the present study only simulated monochromatic waves. Presence of multiple waves might affect the dynamics of LCs, for example, through the CL1 mechanism (Leibovich 1983). As in the present study, wave-resolving simulations will be a useful approach to investigate wave–current interaction in such high–Froude number or multiple-wave situations. The interaction in these situations will be investigated in our future study.

## Acknowledgments

The authors appreciate the insightful comments and suggestions from the editor and the two reviewers. This research is partially supported by Grant-in-Aid for Scientific Research on Innovative Areas (MEXT KAKENHI Grants JP15H05824 and JP15H05825), Grant-in-Aid for JSPS Research Fellow (JSPS KAKENHI Grant JP17J07923), and Initiative on Promotion of Supercomputing for Young or Women Researchers, Information Technology Center, The University of Tokyo.

### APPENDIX A

#### Deep-Water Waves and Stokes Drift Reproduced in Model

Here, we present the capability of our numerical model to properly reproduce deep-water waves and Stokes drift.

First, the dispersion relation of the simulated waves is evaluated. The two-dimensional domain in the *x*–*z* plane is assumed, with a fixed depth of *H* and periodic horizontal boundaries at *x* = 0 and *λ*. As an initial condition, the surface elevation of *η* = *a* cos(2*π*/*λ*)*x* is added to the state of rest (*a* is wave amplitude and much smaller than *λ* and *H*). Periods of resulting standing waves *T* are evaluated at the antinodes, and experiments are conducted with various *λ* values.

Figure A1 shows the measured wave angular frequency (circles), *σ* = 2*π*/*T*, and the dispersion curve of linear surface waves (solid line), , where *k* = 2*π*/*λ* is the wavenumber. It shows that the model is capable of simulating nonhydrostatic surface waves from the shallow to deep-water regimes.

Then, in order to directly measure Stokes drift, particle tracking is performed. Here, the flow field from the control run is used in order to demonstrate the model capability to reproduce Stokes drift and the validity of using an analytical expression of Stokes drift in Fig. 5. When LCs are fully developed in the control run, passive particles are deployed throughout the water column. Then the particles are advected by current velocity for one wave period. The Lagrangian mean velocity of a particle is

where **x**_{p}(*t*) is particle position at time *t* and is its mean position. The Lagrangian mean velocity field is interpolated to the positions where the Eulerian mean velocity is defined. Finally, Stokes drift velocity is obtained as the difference between the Lagrangian and Eulerian mean velocities.

Figure A1 shows the vertical profile of the simulated Stokes drift velocity (horizontally averaged on the *x* = 6*λ*_{10} plane) and the Stokes drift velocity calculated analytically from linear theory [Eq. (12)]. The Stokes drift velocity experimentally evaluated and analytically calculated agree well with each other.

A detailed implementation of our nonhydrostatic free-surface model and further analysis of simulated deep-water waves will be discussed elsewhere.

### APPENDIX B

#### Fourier Analysis of Fluctuation Vorticity

To investigate how the wave-associated vorticity fluctuation arises, a Fourier analysis of the control run result is performed. We let denote the 10-s-period complex Fourier coefficient of vorticity. Figure B1 shows the pattern of and , respectively, at the same slice as in Fig. 5. Figure B1 shows variables that are not smoothed in *x* direction.

If the temporal evolution of wave-associated vorticity fluctuation is dominated by the first term on the rhs of Eq. (9), the Fourier coefficient of vorticity fluctuation follows

Here, *i* is the imaginary unit, and is the 10-s-period Fourier coefficient of **u**. To distinguish from the directly measured , we write in Eq. (B1) as .

The distribution of and (Fig. B1) shows that the measured temporal variation of is well described by . Both real and imaginary parts of both *x* and *z* components coincide with correlation coefficients greater than 0.99 and regression line slopes between 0.95 and 0.97. Therefore, and in the driving terms of the simulated LCs in the control run are induced by advection, stretching, and tilting of the mean vorticity by wave motion.

## REFERENCES

## Footnotes

^{a}

Research Fellow, Japan Society for the Promotion of Science, Tokyo, Japan.

^{b}

Current affiliation: Atmosphere and Ocean Research Institute, The University of Tokyo, Tokyo, Japan.

© 2018 American Meteorological Society. For information regarding reuse of this content and general copyright information, consult the AMS Copyright Policy (www.ametsoc.org/PUBSReuseLicenses).

^{1}

Other terms induce fluctuation that is 90° out of phase with *w*′.

^{2}

Mellor (2016) assumed , but here we retain the original form.

^{3}

The denominators are *L*^{2} because the terms in the vorticity equation is investigated here.

^{4}

This decomposition is normally considered as an Eulerian temporal filtering.