## Abstract

A new mechanism of instability leading to development of Langmuir circulations is proposed and studied based on the hypothesis that the turbulence and, correspondingly, the eddy viscosity are reduced in regions of higher than average contaminant concentration. Here, bubbles are considered as the contaminant, although it is known that surfactants and some particles also are capable of turbulence reduction. The analysis shows that only a very small local decrease in eddy viscosity is needed to initiate the instability. Simplifications to the momentum and bubble turbulence models, as well as neglect of vertical advection, make it possible to analytically solve the perturbation equations and determine the characteristic scale with the maximum growth rate. The scale of the fastest-growing Langmuir circulations is found to be a function of the concentration of bubbles in the near-surface layer, the surface current shear (wind shear), the Stokes drift created by the surface waves, and the eddy viscosity. In contrast to the results of earlier models, the analysis predicts that the maximum change in the current velocity along the direction of the wind (the so-called jets and wakes) is at the surface, not below it. The ratio of perturbation of along-wind surface current and vertical velocity generated by circulations (the pitch) and the aspect ratio of the Langmuir rolls are in reasonable agreement with the experimental data.

## 1. Introduction

The phenomenon of Langmuir circulations (LC) is one of the more important processes in the upper ocean because, to a large degree, it is the process responsible for atmospheric and oceanic interaction, including the exchange of heat, momentum, and gas in the surface layer (Langmuir 1938; Leibovich 1983; Smith 1992, 1998; Thorpe 2004; Sullivan and McWilliams 2010). LC appear at the sea surface as a system of streaks aligned roughly with the direction of the wind. The surface streaks usually contain foam, seaweeds, and organic films and may visually modify the spectra of the short surface waves.

The streaks are produced by paired vortices, which are generated by an instability of the wind-driven shear surface current interacting with the surface waves (Leibovich 1983). The streaks correspond to the convergence regions of the flow produced by adjacent vortices, as the component of surface velocity transverse to the streak is directed toward the streak. The component of the surface velocity along the streak, which is parallel to the wind, is higher in the streak than in the surrounding areas. The areas of high and low surface velocity, which coincide with areas of convergence and divergence of the circulations, are often called surface jets and surface wakes, respectively. The spanwise characteristic scale of LC is usually most pronounced at the onset of the circulations, often appearing as an organized grid of streaks along the wind with a well-defined spacing between them (Langmuir 1938; Leibovich 1983; Smith 1992; Weller and Price 1988; Zedel and D. Farmer 1991). At later times, the streaks at the surface appear more random and are usually described as Langmuir turbulence (Sullivan and McWilliams 2010; McWilliams et al. 1997).

There is an enormous body of experimental and theoretical research work devoted to the study of LC. Included are well-known papers that present a picture of the evolution of the understanding of the physics of LC over the last several decades (Leibovich 1983; Thorpe 2004; Sullivan and McWilliams 2010). Most ocean experiments and observations are related to LC with characteristic scales of several meters or larger. In recent years, small-scale LC with characteristic scales of centimeters have been studied in laboratory experiments (Melville et al. 1998) and later were observed at sea during gusts of wind under calm weather conditions (Veron and Melville 2001). The focus of this paper, however, is on the large-scale LC observed in the ocean.

The current theoretical description of LC is based on the Craik–Leibovich model (Craik and Leibovich 1976; Leibovich 1977, 1980; Leibovich and Paolucci 1981; Leibovich 1983; Phillips 2002, 2005). The Craik–Leibovich model introduced the vortex force, which is the average force imposed by surface waves on a nonuniform current (Craik and Leibovich 1976; Garrett 1976; Leibovich 1977; Andrews and McIntyre 1978). The model is used to study the inviscid instability of a wind-driven surface current in the presence of surface waves. This instability leads to formation of circulations, also identified as rolls, along the direction of the current similar to LC. Unfortunately, the instability does not produce the scales characteristic of the most unstable perturbations, as observed at sea and in laboratory experiments. A nonstationary model was developed to take into account the change in the current with time (Leibovich and Paolucci 1981), and a preferred scale of neutral stability was identified. In that case, the instability was studied in an asymptotic sense, and the growth rate of the unstable mode was not clearly defined.

Later, the Craik–Leibovich model was modified to include the effect of a nonuniform current and strong shear currents on surface waves. This led to the development of the generalized Craik–Leibovich model (Phillips 2002, 2005). A new mode of instability was identified, which had a better-defined characteristic scale that compared favorably to the results of both laboratory conditions and ocean experiments (Veron and Melville 2001; Smith 1992). Nevertheless, another major problem with each of these models remains: the structure of unstable perturbations described by the models does not correspond to the physically occurring LC. Specifically, the theories predict that for small viscosity the maximum change in the current velocity along the direction of the wind (the so-called jets and wakes) is not at the surface, but below the surface (Leibovich 1977; Phillips 2005; Li and Garrett 1993). The velocity of the current at the surface is hardly changed by these perturbations according to the current models. This problem is a reflection of the fundamental nature of Craik–Leibovich-type models for which the velocity of the initial current is restructured by the advection induced by the developing vortices. The prediction by these models that the maximum change in the downwind velocity occurs, not at the surface, but below it, is also a limitation of large-eddy simulation (LES)-based efforts to simulate LC (Skyllingstad and Denbo 1995; Sullivan et al. 2007; Harcourt and D’Asaro 2008; Grant and Belcher 2009; Kukulka et al. 2009).

The perturbation to the along-wind vertical velocity distribution is rarely discussed when the experimental and theoretical results are compared. It was discussed in detail by Li and Garrett (1993), who take into account that turbulent viscosity changes the velocity at the surface, relative to an inviscid case, and creates a nonuniform current structure. They describe this change as follows: the water is constantly accelerated as it converges on the jets and the downwelling there reduces the vertical gradient of the downwind flow, reducing the viscous (turbulent) counter to the wind drag. However, at the initial stage of the process even with the viscosity-induced change to the downwind current, the greatest change of the horizontal velocity still occurs below the surface, which can be seen from comparison of Figs. 1b(2) and 1a(2) in the paper by Li and Garrett (1993). Moreover, it has also been noted by Li and Garrett (1993) that the theoretically described ratio of the downwind jet velocity to the vertical velocity of the LC-produced perturbations is much smaller than the observed value. It is known from observations that the greatest change in the shear current takes place near the surface. Most measurements of the downwind velocity related to LC have been conducted in the near-surface region (Smith 1992, 1998; Melville et al. 1998), and these observations confirm that the significant change in the current takes place near the surface. However, the Craik–Leibovich model predicts only small changes in the horizontal velocity at the surface (Leibovich 1977; Phillips 2005; Li and Garrett 1993).

In addition, there are other important questions about the physics of the LC that are not answered by current theoretical models. One of them concerns the conditions for the development of the instability that generates LC. The controlling parameter of the Craik–Leibovich model is a so-called Langmuir number (Leibovich 1983), which expresses a balance between viscid and inertial forces; it also has the usual interpretation of an inverse Reynolds number. The predicted Langmuir number can differ from experimental values by orders of magnitude (Melville et al. 1998). In fact, in experiments by Smith (1992), a dramatic generation of the LC occurs while the Langmuir number remains essentially unchanged. The discrepancy between theoretical predictions and experimental data indicates that some relevant physical mechanism is missing from the model. It is the goal of the current study to identify the mechanism and to include it in a model of the instability that leads to the generation of LC.

In short, although the Craik–Leibovich model correctly describes some circulations caused by the interaction of surface waves with wind-driven current, these circulations lack some of the basic qualitative features of the observed phenomenon known as Langmuir circulations. Therefore, in the present paper we do not use the Craik–Leibovich model as a starting point in our analysis. In the following, we take as our working hypothesis that bubbles, particles, surfactants, and other suspensions present in the upper layer of the ocean can suppress turbulence, thereby reducing the local “eddy viscosity” and changing the evolution of the wind-driven surface current. We suppose that the contaminants reduce the eddy viscosity in some local region, so that the wind stress at the surface will drive water faster in that region and generate surface jets. The interaction of those jets with the surface waves through the vortex force will introduce circulations with regions of convergence along the jets. The circulations will, in turn, concentrate contaminants into the convergence regions, which will further reduce the eddy viscosity and provide a positive feedback. We note that it has been observed experimentally that the bubble concentration is much higher in the regions where the LC converge and that these regions coincide with higher velocity of the wind-driven surface currents (Smith 1992; Zedel and Farmer 1991; Thorpe et al. 2003a,b).

While we clearly recognize that the state of the art in turbulence modeling has long moved beyond eddy viscosity and understand that a rigorous solution of this problem will require a sophisticated approach to the coupled bubble and flow field turbulence, here we take a significantly simpler approach so as not to obscure the essence of the proposed mechanism. We assume a single eddy viscosity for the unperturbed surface layer and allow for spatial–temporal changes to its value as a result of a local accumulation of bubbles. An increase in bubble concentration might be, for example, the result of wave breaking as a consequence of an increase in wind speed. In addition, we assume that the thickness of the surface layer containing the bubbles is small compared to the characteristic scales of the flow, which allows description of the bubble concentration using a two-dimensional model. This assumption is supported by the experimental data (Smith 1992). Finally, we will show, through a comparison of our theoretical results with the experimental data, that the reduction in the eddy viscosity, due to an introduction of bubbles into the turbulent upper layer, required for development of the LC is very small.

The underlying assumption of the model of LC described here is very difficult to validate by direct measurements. However, the suppression of turbulence by surfactants (Procaccia and L’vov 2008), bubbles (Madavan et al.1985; Mohanarangam et al. 2009), and some particles (Zhao et al. 2010), leading to a reduction in boundary layer drag, has been studied for decades. Following an argument first advanced by Lumley (1977), any additive that increases the molecular viscosity away from a solid surface, while leaving the viscosity at the surface unchanged, will reduce drag. For clarity in the presentation, we choose to focus on bubbles. For a low bubble concentration, we would expect the dynamic viscosity to increase through an Einstein-type relation, while the density will decrease, increasing the kinematic viscosity and reducing the effective turbulent Reynolds number. Although the described mechanism of reductions of eddy viscosity was mainly studied in relation to drag reduction in the boundary layer, it does not depend on the presence of a rigid surface. In the ocean, the effect of such a change in the concentration of bubbles on the eddy viscosity seems plausible.

Since the model presented here contains some unmeasured parameters, comparisons of the theoretical results with experimental data are often qualitative. Nevertheless, quantitative estimates can be made for relationships that do not include these unmeasured parameters. One example is the relationship between the characteristic scale of fastest-growing circulations and rate of their growth. Here, we find that the growth rate is proportional to the eddy viscosity and inversely proportional to a square of the characteristic scale with a coefficient of proportionality near one. As there are no unknown parameters in this relationship, it can be checked against known experimental data in cases where the spatial and temporal scales of growing LC are observed. Estimate of the eddy viscosity from the equation should be in a range of generally accepted values and, fortunately, the relationship between the rate of growth and characteristic scale of the fastest-growing circulations provides reasonable estimates of the eddy viscosity for all known cases. Where comparison with experimental data is possible, we focus on the comprehensive experimental study of LC in the ocean by Smith (1992), as that study provides the most detailed measurement of the relevant physical parameters.

Finally, we note that since the system is linear to both the Craik–Leibovich instability and the contaminant-related instability proposed here, they may develop simultaneously. An important difference between them is that the instability due to the effect of the bubbles on the eddy viscosity leads to exponential growth of the circulations, while perturbations due to the Craik–Leibovich instability grow algebraically with time. In any event, numerical simulations are required to study the simultaneous effect of both mechanisms.

## 2. Statement of the problem and formulation of the model

We consider the instability of wind-driven surface currents resulting from their interaction with surface waves, which is governed by (Landau and Lifshitz 1959)

and

Here, is the water velocity, is density, is pressure, is eddy viscosity, and is the vortex force. The indices in (1)–(2) are 1, 2, or 3 and correspond to the , , and coordinates, whose respective velocity components are , , and . The directions of the wind, surface current, and surface waves are considered to be the same and coincide with the *x* axis. The *y* axis is horizontal and transverse to the direction of the current (spanwise direction), and the *z* axis is vertical and directed upward. The unperturbed eddy viscosity is assumed to be known. We note that by using the incompressible forms of (1) and (2), we must neglect the change in density caused by the presence of the bubbles, which will play a role in the turbulence reduction. For the purpose of this analysis, where we have in fact assumed a very simple “turbulence model,” the details of the mechanism of drag reduction are not needed to show the strong influence of the local change in eddy viscosity on the development and evolution of the LC.

The vortex force is given as (Leibovich 1977, 1980)

Here, is the Stokes drift velocity. In the simple case of a monochromatic surface wave in deep water, the Stokes drift velocity is given as (Phillips 1977)

Here, is the amplitude of the surface wave with its wavenumber, and is the acceleration of gravity. The boundary conditions for the velocity at the surface with the wind stress applied in the direction are (Phillips 1977)

where is the wind stress at the surface.

We further assume that the eddy viscosity is linearly proportional to the concentration of bubbles, which decreases exponentially with depth:

Here, is the dynamic eddy viscosity due to the turbulence but without bubbles; is the concentration of bubbles near the surface, which depends on the horizontal coordinates; is inverse characteristic scale of the vertical change in the concentration; and is some unknown coefficient. Reduction in the eddy viscosity due to the presence of the bubbles corresponds to a negative value of . We assume that the second term in expression (5) is small:

We also assume that the thickness of the layer containing the bubbles is much smaller than the characteristic vertical scale of the initial current and of the characteristic horizontal and vertical scales of the growing circulations:

It will be shown that the horizontal and vertical scales of the developing circulations are of the same order and thus are described by one value . The assumption expressed by condition (7) is justified, since the characteristic scale of the layer with bubbles is usually on the order of a meter, while the characteristic scale of the initial current is on the order of meters and the characteristic scale of LC is on the order of tens of meters (Smith 1992).

When condition (7) is satisfied, then the concentration of the bubbles in the near-surface layer can be described by the following simplified equation (Landau and Lifshitz 1959):

Here, and are horizontal components of the velocity of the flow at the surface, and is the eddy diffusivity of the bubbles. We assume here that the characteristic time to restore the vertical structure of the bubble distribution is small compared to the characteristic time of the change in the bubble concentration caused by horizontal diffusion and by near-surface current, as described by (8). The exponential distribution of the bubbles is formed due to wave breaking and dissolution on the one hand and turbulent diffusion and the rise of the bubbles under buoyancy on the other. While the downwelling current created by LC could affect the distribution of the bubbles and there is experimental evidence that bubbles are carried down by strong circulations, we focus here on the initial (linear) stage of the instability when all perturbations are small and the buoyancy can restore the vertical distribution of bubble concentration near the surface. More realistic models of bubble distribution are still being developed [see discussion by Sullivan and McWilliams (2010)].

Now consider small perturbations in the initial current, directed along the direction. The velocity of the horizontal current is a function of the vertical coordinate, of the pressure, and of the initial concentration of bubbles. The latter is independent of the horizontal coordinates. Under the assumption that the characteristic development time of the perturbations resulting from the instability is much smaller than the time to change the wind-driven current, the initial values can be considered time independent. A similar assumption was used in the early work on the Craik–Leibovich instability (Leibovich 1977). Thus, we can introduce perturbations to (1), (2), and (8) with boundary conditions (4) in the following form:

Here, is the velocity of the main current, is initial pressure, and is the initial concentration of bubbles. In expressions (9) and (10), the small perturbations in velocity and pressure are functions of the spanwise coordinate , the depth , and the time as , , , and . The concentration of bubbles near the surface is a function of and : . As a result of the perturbation to the bubble concentration, the turbulent eddy viscosity will also be perturbed:

Here, is the unperturbed dynamic eddy viscosity in the presence of the bubbles.

Now substitute the expressions (9)–(11) into (1), (2), and (8) and boundary condition (5), subtract the mean flow equations, and keep only the linear terms with respect to small perturbations as follows:

Corresponding boundary conditions are

Of course, all perturbations should be zero at large values of depth. In (12) and (17), and are first and second derivatives of the initial current velocity with respect to the vertical coordinate, and is the derivative of the initial dynamic eddy viscosity. The last terms in (13) and (14) are the and components of the vortex force.

To simplify (12)–(14) and boundary conditions (17), based on condition (6), we will neglect all the terms in the equations containing the products , , , and their derivatives, as well as all the terms containing . We note that formal justification for neglect of these terms requires the solution of the perturbation equations. With the solution in hand, it is straightforward to show that these terms are small as a consequence of inequality (6) (this is demonstrated in the appendix). The simplified system of (12)–(14) is

and

The corresponding boundary conditions [(17)] are

Now, further simplify this system of equations by introducing the streamfunction of the flow in the plane :

That allows us, in a standard way, to eliminate (15) and to combine (19) and (20) into one equation that does not contain terms with the pressure perturbation .

When , this system returns to the Craik–Leibovich model, as expected. In that model, changes in the initial horizontal current are mainly caused by the advective term , produced by the growing circulations. In principle then, both mechanisms for instability might operate simultaneously. In the case considered here, the mechanism of instability caused by the effect of the bubbles on near-surface turbulence is prevalent and an analytical solution of the instability problem can be obtained and examined. To obtain such a solution, the advection term in (18) must be neglected. While it is shown in the appendix, based on data from Smith (1992), that for some relative changes of the eddy viscosity due to the presence of the bubbles (for the higher end of the obtained estimates), this term is small compared to the last term in (18); in general, the terms are of the same order. Nevertheless, the changes in horizontal velocity caused by the vertical advection term occur well below the surface, at the depth of about (as estimated in the appendix), where is the characteristic vertical scale of the initial current. Therefore, the vortex force in (19) and (20) that drives the developing circulations is significantly smaller for changes in due to the vertical advection term than for changes caused by effect of the bubbles that occur near the surface, since the Stokes drift (and the vortex force) decreases exponentially with depth. For example, at the depth on the order of several meters and wavelength of the surface waves of 30 m [an estimate from the paper by Smith (1992)], the Stokes drift is lower by more than an order of magnitude compared to its value at the surface. For shorter surface waves, the decrease in the vortex force is even more dramatic. Therefore, in the current analysis we neglect the advection term in (18). In addition, the term (which can be estimated as ) is small compared to in the brackets in (18) under condition (7) and is neglected as well. Finally, a simplified system of equations, with boundary conditions at the surface that describes the evolution of perturbations leading to the generation of LC is obtained:

The value of in (23) is taken at the surface , since it changes little within the thin layer containing the bubbles. The last term in (24) is a combination of derivatives of the components of the vortex force in (19) and (20). The value is the derivative of the Stokes drift velocity with respect to and is derived from expression (3):

## 3. Approximate solution for growing perturbations

To obtain the relationship between the growth rate and the characteristic scale of the perturbations, to find the scale that corresponds to the fastest growth rate, and to examine the structure of the evolving circulations and their dependence on depth, we introduce

into (23)–(25) under boundary conditions (26) with the understanding that the perturbations vanish at large values of depth. Here, is growth rate of the instability, is the spanwise wavenumber, and are complex functions of , and is a complex number. Equation (23) for with the boundary condition at the surface for given , takes the form:

Here, . We will show that for the fastest-growing scale , and, correspondingly, is roughly . This notion is used to obtain an approximate solution of (23)–(26).

Under the condition specified in (7), where the layer with bubbles is thin compared to the horizontal and vertical scales of the perturbations, we have . The first term in the brackets of expression (31) can be neglected and the denominator simplified to give an approximate expression for , as a function of , as

For a reduction in eddy viscosity, is negative, and (32) shows that the velocity perturbation along the wind-driven current increases linearly with bubble concentration.

From (25) for a given value of , we find the simple expression for concentration:

where

Here, is the derivative of at the surface, which, according to (22), is the horizontal spanwise surface velocity, and is a normalized coefficient of diffusion relative to one. Expression (33) indicates that perturbations of transverse velocity and of concentrations of the bubbles are shifted relative to each other by a quarter of the spanwise period ; their phases are in quadrature ().

with boundary conditions

Here,

and

Equations (34)–(37) describe the structure of the developing transverse circulations. The value , in expression (37), has the dimensions of a wavenumber, and we will show that this value defines the characteristic scale and, correspondingly, the rate of growth of LC.

The right-hand side of expression (34) for the streamfunction is stated now as a function of its own derivative at the surface, which is the transverse velocity [see (22) and (36)]. Assuming that the wavelength of the surface waves is much greater than the thickness of the layer containing the bubbles (), we neglect the second term in the exponent in the right-hand side of (34). The solution of simplified (34), which satisfies boundary condition (35) at the surface, is

For and , an approximate expression for the streamfunction, which can be used in the region away from the surface, is

Here, the last term in expression (38) is neglected as it is small in comparison to the other two terms everywhere except in the surface layer to the depth . The derivative of the streamfunction represents the spanwise velocity of the perturbations and can be found from (39):

It should be noted that this approximate expression [(40)] is valid even at the surface, as the inaccuracy introduced by neglecting the last term in (38) is not significant. That term is only required to satisfy the boundary condition (35) exactly, which requires the second derivative of the streamfunction to be equal to zero at the surface. In fact, the approximation given in expression (39) satisfies the boundary condition when the amplitude of the streamfunction is equal to zero at the surface (). This approximate expression [(40)] will be used to study the structure of the circulation’s streamlines later in the text.

Here, we note that (33), (39), and (40) confirm that the perturbations are indeed the circulations. Specifically, the horizontal component of the velocity, described in expression (40) as a function of , changes sign at some depth. It is also clear, from expression (33), that the amplitude of the spanwise velocity of the perturbation at the surface is related to the amplitude of the change in the concentration of bubbles. Specifically, at the surface, the transverse velocity as a function of [see (28)] changes its sign at the point of maximum bubble concentration, since the phases of these values are in quadrature. Thus, the circulations described by expression (39) converge at the region of maximum bubble concentration and, correspondingly, around the jet of the perturbed current in the direction.

By using expressions (36) and (40), a relationship between the growth rate of the developing circulations and their characteristic spanwise scale can be obtained. Substituting (36) and (37) into (40), making equal to zero, and canceling on both sides, gives

The value is a function of both the LC wavenumber and its growth rate according to the definition presented in the text below expression (30):

Thus, the solution of equations (41) and (42) provides the relationship between and . Our interest is in determining whether or not the algebraic equations (41) and (42) provide for a maximum circulation growth rate , with a characteristic wavenumber , and, if so, what the values are. As a closed-form analytical solution is not available, we introduce a parameter , such that , , and

From (41) we have

while from (42), we obtain

Clearly, the range of interest is , since that corresponds to a positive rate of circulation growth . Substituting (44) into (45), we can write

Expressions (44) and (46) represent both and as functions of and the three parameters , , and . Extrema of the function correspond to the points where its derivative equals zero, which may be written with respect to the parameter as

For the derivative (47) to be zero, either its numerator has to be zero or its denominator must go to infinity. As the magnitude of the coefficient , introduced in expression (33), can be only in the range and with , the denominator of (47) [the derivative of (43) with respect to ] does not have singularities and condition (47), for the extrema of , becomes

The solution of (49) provides the extrema for , as a function of the parameter . It can be verified numerically that (49) has only one real root in the interval for all values . Correspondingly, the rate of growth of the circulations as a function of has only one extremum, a maximum, since according to (46) the value is positive except at an of one or infinity, where it is zero. Expression (44) demonstrates that the wavenumber of the growing circulations is a monotonic function of . That means that the function , represented parametrically by expressions (44) and (46), also has only one maximum as shown in Fig. 1. In the figure, the normalized values of the growth rate and wavenumber are used, where

It is now clear that the model developed in this paper, unlike previous models, provides the particular characteristic scale of the LC that correspond to the highest rate of growth of the circulations . The values of and are functions of parameter and can be easily found numerically from (44) and (46) using the function , which in turn is found from (49). The functions , , and are plotted in Fig. 2. Figure 2 shows that both the most unstable wavenumber and highest growth rate of circulations are greater for smaller values of and, correspondingly, for smaller coefficients of bubble diffusion . The value of is not known. As the horizontal “diffusion” of bubbles is caused by the turbulence in the same way as is the diffusion of momentum (eddy viscosity), we assume that and are equal and thus the normalized coefficient of diffusion relative to one is

All further considerations use condition (51), although they could be easily modified for any value of . For the case described by (51), the value of is found from (49) as

The values of and in this case are described according to (44) and (45) by the following expressions:

Here, and . Expressions (53) will be used for comparison of the predictions of the model with the experimental results.

The streamlines of the circulations in the plane are shown in Fig. 3, where and are normalized horizontal spanwise and vertical coordinates. The streamlines in Fig. 3 are plotted based on expression (39) multiplied by the function that describes the change of perturbations in the spanwise direction. The plot shows only the structure of the streamlines (contour lines of streamfunction), since, for simplicity, the streamfunction is normalized such that

In addition, only one-half of the paired circulations is shown. The depth of the center of the circulations , where the velocity is zero, can be found by using the definition of the streamfunction (22) and by setting the derivative of expression (54) with respect to the vertical coordinate to zero as

The maximum vertical velocity of the circulations, that is, the derivative of expression (54) with respect to the horizontal coordinate, also occurs at the depth .

An important parameter observed experimentally, but not previously described theoretically, is the aspect ratio of the Langmuir rolls. This parameter is defined as a ratio of the half-spacing of the LC to the depth of the rolls. For the case of LC in nonstratified water, this ratio is reported to be close to one (Leibovich 1983; Smith 1992). The circulations presented in Fig. 3 correspond to the case of no stratification and clearly have the aspect ratio of about one. It should be noted that the aspect ratio of developing circulations were observed in numerical simulations based on the Craik–Leibovich model (e.g., Leibovich and Paolucci 1981; Kukulka et al. 2009; Tejada-Martinez et al. 2009) and this ratio can be on the order of one. Nevertheless, the aspect ratio in these simulations depends on the parameters used and, more importantly, changes significantly as the circulations develop. According to the model presented in this paper, the fastest-growing LC in nonstratified fluid always have the same aspect ratio of one.

Another parameter of LC found experimentally (Weller and Price 1988; Smith et al. 1987) is the pitch , which is a ratio of the surface downwind jet velocity to the maximum downwelling velocity. This parameter is discussed in detail by Li and Garrett (1993), who estimated its value from available experimental data to be about 1.5 or higher. At the same time, they pointed out that the maximum value of the pitch predicted by the Craik–Leibovich model is less than 0.5. Thus, according to the currently accepted theoretical model, the downwind current strength is significantly smaller than the downwelling velocity, which is in contradiction to experimental results.

The value can be easily found from the results obtained in the current paper. The surface jet velocity is defined as the difference between the surface velocity in the wind direction at the regions of convergence (the jet) and the regions of divergence (the wake) (Li and Garrett 1993). Therefore, this value is twice the perturbation amplitude of the surface velocity in the current model, which is given by expression (32). The pitch can be expressed as follows:

where the subscripts represent maximal values of the perturbation velocities and . This ratio can be found by using expression (32) for and expressions (28) and (39) to derive as the derivative of streamfunction with respect to spanwise coordinate . Taking into account expressions (36) and (37), the pitch can be expressed as follows:

From (57), it is clear that the ratio of the downwind jet to downwelling velocity depends on the eddy viscosity, the thickness of the upper layer containing the bubbles, and the derivative of the Stokes drift at the surface, which is .

## 4. Discussion of theoretical predictions and experimental data

The theoretical results derived from the model of LC developed in this paper qualitatively explain many of the observed features of the developing circulations. As a first test of the validity of the model, we consider its predictions based on the experimentally observed relationship between the spacing of the circulations and the rate of their growth. To do this, we turn to the experimental study by Smith (1992) in which details of the birth and evolution of the LC were observed. As an initial state, Smith (1992) reports surface waves but no LC under a relatively strong wind of about 8 m s^{−1} that had been blowing steadily for several hours. When the wind speed abruptly changes from 8 to 13 m s^{−1} (over 1 min), the LC develop in approximately 10 min from the initial increase in wind speed. Coincidentally and independently, a significant increase in the intensity of the surface waves occurs at the same time as the increase in the wind speed: strong long surface waves come to the area of the observations. The spacing of the LC that is first observed is about 15 m, though within 2 h, the LC spacing increases gradually to roughly 100 m. During this period, the depth of the mixed layer created by LC increases at roughly half of the rate of increase of the spacing.

Using these values, the magnitude of the eddy viscosity can be estimated. If the model that produces the relationship between the growth rate and spacing of LC [(45)] is valid, then the magnitude of turbulent viscosity should be in the range that corresponds to the experimental conditions. In the case described by Smith (1992), the observed spacing, , of dominant LC was 15 m with characteristic time of the development of the LC (inverse of the growth rate ) of about 10 min from the onset of the strong wind to about 100 m for a time of LC development less than 2 h. We may find the eddy viscosity from expression (45):

Using the value of from (52) and the values of and estimated from the experimental data, we find that at the beginning of the generation of LC, the eddy viscosity is about , and at the end of the considered time interval it is about . We note that these values of the eddy viscosity are reasonable for the relatively strong storm conditions of the observations. If the coefficient of diffusion of the bubbles is higher than the eddy viscosity (we assumed that they are equal), then ; from Fig. 2, the magnitude of is larger than that given by expression (52), and the estimates of viscosity obtained above can be significantly lower. Some indications that is larger than 1 can be inferred from the experimental data (Smith 1998).

The increase of the eddy viscosity with time is a result of the intensive wave breaking induced by the significant increase of the wind speed. The increase of turbulence in the upper ocean due to the wave breaking has been extensively studied [see, e.g., the discussion in the paper of Sullivan and McWilliams (2010)], but a full description of this phenomenon is not yet available. A second feature of the wave breaking is the creation of bubbles. In our model, it is the concentration of these bubbles in the convergence (jet) regions of the flow and the subsequent local reduction in the eddy viscosity that triggers the LC instability.

The increase in eddy viscosity due to the wave breaking explains the observed growth of LC spacing during the interval of observation. The evolution of LC from smaller to larger scales after the initial stage of generation is commonly observed (Leibovich 1983). An initial stage follows a strong increase of the wind that causes a sharp increase in the bubble concentration. Thus, relatively small-scale LC with a high growth rate are generated initially. With time, both the concentration of the bubbles and the eddy viscosity increase gradually. If we assume that they increase at comparable rates, then according to expression (37) the spacing of LC, which is inversely proportional to the characteristic wavenumber , will also increase. Indeed, the value of is proportional to , , and . The last value is inversely proportional to , as shown below. Thus, K will decrease with time and the spacing of LC will increase. We may also estimate the parameter to find the change of eddy viscosity due to the presence of bubbles that is required in the model from expression (37), which we modify as

Here, the value in (37) is replaced by its more general expression: the derivative of the Stokes drift velocity with respect to the vertical coordinate . We do this because an estimate of this parameter is given by Smith (1992) from the experimental data: . The derivative of the wind-driven current can be estimated from the boundary condition [(4)] , where is the wind stress applied to the surface. For the wind speed of 15 m s^{−1} the value is about 0.3 N m^{−2} (see, e.g., Phillips 1977). The parameter can be expressed from (53) as . At the beginning of the observation, we estimate the wavenumber of the fastest-growing LC for a spacing of 15 m and an eddy viscosity of and a characteristic thickness of the layer containing the bubbles of 1 m. For these values, expression (59) gives . At the end of the observation, we take a spacing of 100 m and the eddy viscosity of and get . Thus, the required change of eddy viscosity due to the effect of bubbles on the turbulence is very small, less than a percent.

We note that the value of pitch, estimated for the same parameters in (57), is about . This value is in good agreement with the generally observed values (Li and Garrett 1993). In addition, the aspect ratio for the LC computed in this paper for a nonstratified case is close to 1, which also corresponds to the value experimentally observed (Smith 1992). A crucial role of the bubbles in the development of LC is supported by another experimental study conducted by Smith (1998). A detailed analysis of this experiment will be the subject for an additional paper. Here, we mention one observed phenomenon that can be best explained with the model presented in this paper. It is known that, when the wind is slowly changing its direction, the Langmuir circulations do not exactly align with the wind. The effect observed by Smith (1998) is that during the change in the wind direction there is a difference between the directions of the rolls and the bands of bubbles. This difference cannot be explained if the bubbles are considered as simple tracers, as in the previous models, and not considered as an active part of a mechanism of the generation of LC.

## 5. Conclusions

We have examined the onset and evolution of Langmuir circulations using the hypothesis that the turbulent eddy viscosity may be locally reduced by an accumulation of contaminants such as bubbles, surfactants, or particles. For simplicity, we have focused here on bubbles. The reduction in eddy viscosity due to the presence of bubbles has been studied for decades as a means of reducing drag in wall-bounded shear flows. Accumulation of bubbles near the surface may be a result of wave breaking, due to an increase in wind speed. Regions with a higher concentration of bubbles will be driven faster by the wind, becoming surface jets. Interaction of these jets with surface waves produces a converging circulation around the jet, further concentrating bubbles and initiating the instability.

This analysis permits, for the first time, the calculation of a dominant characteristic scale for the developing LC with the highest growth rate. The scale of the fastest-growing LC is found to be a function of the concentration of bubbles in the near-surface layer, the surface current shear (wind shear), the Stokes drift created by the surface waves, and the eddy viscosity. Exercising the model based on the observations of Smith (1992) allows us to demonstrate that the calculated eddy viscosities are reasonable for the sea state described in the experiments and that only a very small change in the eddy viscosity is needed to trigger the instability. In contrast to earlier models, we show here that the maximum perturbation of the wind-driven current occurs at the surface and not below it. The measured ratio of the perturbation of the along-wind horizontal velocity to the vertical velocity created by LC (the pitch) and the aspect ratio of the Langmuir rolls are in reasonable agreement with the experimental data.

The analytical solution of this problem required many simplifications. The vertical advection, which affects the wind-driven current, was neglected here and has to be taken into account. More complete turbulence models, for both momentum and bubble concentration, should also be used. Therefore, a study of the nonlinear stage of development of LC will likely require direct or “large eddy” numerical simulations. Although the reduction of turbulence by the presence of bubbles, surfactants, and some particles has been thoroughly documented in wall-bounded shear flows, experimental work on the influence of bubbles on the turbulence in the surface layer of the ocean is needed. In this regard, the recent work of Gemmrich (2012), which experimentally confirms the suppression of turbulence by near-surface bubbles, is significant and may provide the basis for improving the simple model used here. Also of interest is the effect of wind amplification of circulations caused by the interaction of near-surface currents, such as ship wakes, with surface waves (Basovich 2011).

## Acknowledgments

The author is extremely grateful to the President of Cortana Corporation, Mr. K. J. Moore, and to Professor Steven Deutsch for invaluable help in the improvement of the paper in both technical clarity and style of the presentation.

### APPENDIX

#### Conditions for Neglecting Small Terms in Simplified Model Equations

Some small terms were neglected in the derivation of a simplified system of equations, which describe the evolution of the perturbations that lead to the generation of Langmuir circulations. Justification for this simplification is provided herein.

- The condition of the third term on the right-hand side of (12) being much smaller than the last term is From expression (11), we have From the solution (32), we get Substituting expressions (A2) and (A3) into expression (A1) and noticing that in the right part of the obtained inequality in the narrow upper layer can be replaced with its value at the surface , we retrieve expression (6):
- There are three terms containing derivative and derivatives of the perturbation velocity components in (13) and (14) that were neglected. The conditions of neglecting these terms can be derived using the solutions found in the paper. Here, we derive the condition for the term in (14) being small compared to the component of the vortex force. From the second expression (22), the derivative of vertical velocity near the surface is The value is found from (11), and value of the can be expressed from the solution (32). Thus, the condition takes the form The condition (A7) can be further simplified. We substitute the amplitude of the perturbation of concentration of the bubbles from (33), assuming that , and express the product of as a function of the parameter from expression (37). Recalling from (53) that , we have the simple condition This condition is always satisfied based on assumptions expressed by inequalities (6) and (7). Conditions similar to (A8) can be derived for neglecting the two other terms containing in (13).
- The condition of neglecting the vertical advection term in (18) is For simplicity, we assume that the main surface current velocity decreases with the depth exponentially: Here, is the characteristic vertical scale of the surface current [see (7)]. The absolute value of the vertical velocity is A maximum value of is reached roughly at the depth . From expression (33), it follows for that Substituting (A12) into (A11) and then the result of this substitution into (A9), we get the following condition The condition (A13) means that the maximum advection term is smaller than the term related to the change of concentration of the bubbles. This condition is rather strong, considering that the right-hand side of inequality (A13) is already small according to expression (6). We estimate two sides of inequality expression (A13) using experimental data obtained by Smith (1992). For characteristic scales of developing circulations , of tens of meters and of initial current of several meters the left-hand side of the inequality (A13) is on the order of . Estimates of the value of , according to expression (59), discussed in the paper are between and , so that for the higher estimate of the condition (A13) can be satisfied. Clearly, for a higher value of and a lower value of the condition is satisfied for a smaller change in the eddy viscosity . It should be noted though that it is not a necessary condition for the described mechanism to work and is introduced only with the goal of simplifying the analysis of the model equations. Full numerical analysis of the model equations will be required to include the advective term into the solution. Nevertheless, qualitatively it is clear that taking this term into account would only enforce the instability, since the vortex force due to the flow change caused by advection also amplifies the circulations.

## REFERENCES

*The Dynamics of the Upper Ocean.*Cambridge University Press, 336 pp.