One-dimensional modelling of upper ocean mixing by turbulence due to wave orbital motion

Mixing of the upper ocean affects the sea surface temperature by bringing deeper, colder water to the surface. Because even small changes in the surface temperature can have a large impact on weather and climate, accurately determining the rate of mixing is of central importance for forecasting. Although there are several mixing mechanisms, one that has until recently been overlooked is the effect of turbulence generated by non-breaking, wind-generated surface waves. Lately there has been a lot of interest in introducing this mechanism into ocean mixing models, and real gains have been made in terms of increased fidelity to observational data. However, our knowledge of the mechanism is still incomplete. We indicate areas where we believe the existing parameterisations need refinement and propose an alternative one. We use two of the parameterisations to demonstrate the effect on the mixed layer of wave-induced turbulence by applying them to a one-dimensional mixing model and a stable temperature profile. Our modelling experiment suggests a strong effect on sea surface temperature due to non-breaking wave-induced turbulent mixing.


Introduction
For many applications ocean waves can be adequately modelled using potential wave theory, meaning the water is presumed inviscid, incompressible and the motion adiabatic.As a result, the wave orbital motion, which in potential theory is irrotational, is often presumed incapable of producing turbulence, and therefore unable to contribute to ocean mixing.(A recent paper by Beyá et al., 2012 argued, based on experiments they carried out, that the deviations from irrotational motion are entirely negligible.We discuss this paper in Appendix B.) Yet real waves by definition do produce vorticity (for example Kinsman, 1965;Landau and Lifshitz, 1987), and moreover it may be shown that even potential waves can, by stressing vortex lines, amplify pre-existing vorticity (Phillips, 1961;Kinsman, 1965;Benilov et al., 1993).Nevertheless the effect was either ignored or thought too small to produce any significant effect, and it remains excluded from operational mixing models.However, recent work suggests that the effect is significant enough that its inclusion in geophysical models could greatly improve ocean as well as sediment suspension, weather and climate forecasts.Despite the encouraging results the existing approaches are in large part guesses, undeniably useful though they are, and a parameterisation based on the underlying physical principles is still needed.
Our aim is, ultimately, to produce a useful and efficient parameterisation for the wave-induced turbulence so that it can be applied to ocean mixing models and, by extension, to ocean circulation, weather and climate models.And because the generation of turbulence is a sink for wave energy, the same parameterisation should improve wave modelling itself by better accounting for swell dissipation.This has historically been poorly predicted (Rogers, 2002), although recent advances have improved the situation considerably (Ardhuin et al., 2010).Nevertheless, the physics is not yet settled, and while it seems likely that more than one process is responsible, the precise contributions are unknown.In this paper we assume that the bulk of the swell dissipation is due to turbulence production.
We briefly run through the theoretical and experimental justification for the mechanism in Sect. 2. In Sect. 3 we review the different approaches attempted and describe their strengths and weaknesses.We also suggest another way of Published by Copernicus Publications on behalf of the European Geosciences Union & the American Geophysical Union.
introducing the effect; this is described in Sect.4, along with its application to a one-dimensional mixing model and the results.
There are of course other mechanisms of turbulence generation involving waves, such as wave breaking (see for example Chalikov and Belevich, 1993) and Langmuir circulation (Sullivan and McWilliams, 2010).We do not discuss wave breaking in detail since it is a source independent of that produced by the wave orbital motion which concerns us here.Langmuir circulation, on the other hand, is connected, as it is driven by the Stokes drift, a property of the wave orbital motion -but it is not a complete description since the Stokes drift does not encompass the entire orbital motion.We briefly discuss this in Sect.3, though we do not treat Langmuir circulation directly.

Wave-induced turbulence
As early as the first quarter of the 20th century, Jeffreys (1920) sought to explain observed eddy viscosity in the upper ocean by the wave orbital motion, which he noticed had a similar rate of depth attenuation.Later, Dobroklonskiy (1947) and Bowden (1950) produced the first useful parameterisations for eddy viscosity and turbulence production.Further theoretical work examining the order of the effect of the orbital motion on turbulence was published in the 1960s (Phillips, 1961;Kinsman, 1965).Meanwhile, observations from the 1970s (Yefimov and Khristoforov, 1971a, b) and possibly as early as the 1950s (Francis et al., 1953 observed higher than expected stresses near the surface in windless conditions and speculated that they may be due to a "hitherto unexpected complication . . .caused by the free surface") seemed to confirm the existence of the mechanism, but the phenomenon had more recently been neglected by oceanographers, though not by experimentalists.
One reason for this neglect might be because as a mechanism for wave dissipation it is much smaller than that caused by wave breaking.Indeed, Phillips (1961) (see also Kinsman, 1965 for a more expansive treatment) tried to estimate wave dissipation caused by wave-turbulence interaction, and the conclusion that it was negligible could have had some influence among oceanographers.But the treatment assumed a Reynolds number defined using wavelength as the length scale rather than amplitude, which is the variable that controls the size of the water particle orbital motion; and wave dissipation due to turbulence was demonstrated all the same.
A further reason for the neglect could be that, especially from the 1960s, many important features of nonlinear wave dynamics were found to be well described by potential theory where viscosity and vorticity were neglected (see, for example, the discussion in Babanin and Chalikov, 2012).Because turbulence by definition requires vorticity, it might have seemed that since vorticity could be safely ignored in describing wave dynamics accurately enough that turbulence generated by them could be ignored too, but such an interpretation misses three vital points: firstly, that despite being small, the vorticity is not zero, and could still produce a measurable effect on the boundary layer if the waves were large enough; secondly, that irrotational waves can interact with and amplify pre-existing turbulence; and thirdly, that while the energy involved in turbulence production might be small relative to the wave energy, or even wave breaking, it does not necessarily follow that it is small with relation to ocean mixing.
To elaborate on that second point, Benilov et al. (1993) (this argument was revisited by Benilov, 2012) demonstrated by a direct instability analysis that vorticity grew with the passage of a linear potential plane wave, and that the only component that was stable was the one directed along the wave crest.Any treatment restricted to two-dimensional motion would miss this important aspect of the flow.Furthermore, Anis and Moum (1995), by considering the Reynoldsaveraged equations of motion partitioned three ways into mean, orbital and turbulent components, differentiated between two mechanisms for turbulence production by the wave orbital motion -one that irrotational motion could account for and another that required rotational motion.They attempted to apply both to observed scalings of turbulent dissipation which the conventional mechanisms (such as wave breaking) could not adequately explain, but they did not determine the precise contribution of each term, and both scalings, though based on different mechanisms, could be used to describe the observed profiles.
Experiments such as those of Ölmez and Milgram (1992) and Milgram (1998) sought to investigate the passage of waves on a turbulent body of water, and demonstrated that the existence of turbulence, generated independently of the wave motion, led to a measurable attenuation of the surface waves.Taken together with the theoretical arguments that the waves should interact with the turbulence, and especially that they should transfer their energy to it, this finding highlights that turbulence production and swell attenuation in the absence of breaking are two sides of the same coin, a notion exploited by Bowden (1950), Babanin (2011) and by us too, in Sect.4.1.

Existing parameterisations
Several approaches to parameterising non-breaking waveinduced mixing in numerical models have already been tried, and they have all been successfully used to improve model results.They fall into three broad categories: those that attempt to parameterise Langmuir circulation; those that are scaled with wind stress; and those that, like the present one, are based on the notion that a fundamental property of the wave orbital motion is capable of generating turbulence which can lead to mixing.Those based on wind stress include the parameterisations described by Jacobs (1978) and Huang and Qiao (2010) (also Huang et al., 2011), and we discount these immediately since we are interested in capturing the effect of turbulence generated by swell in the absence of wind stress.Parameterising Langmuir circulation is usually done by scaling it with the Stokes drift.For example, Kantha and Clayson (2004) modify the Mellor-Yamada turbulence model by augmenting the turbulence production by adding the Stokes shear vectorially to the current shear in the turbulence production terms.In the kinetic energy equation this is done directly, however in the length-scale determining equation they rescale the Stokes-driven production in order to reflect the large size of the Langmuir cells in the turbulent macro length scale.Their model parameters were chosen in order to match the LES results of McWilliams et al. (1997).While studies such as these suggest that Langmuir circulation is an important mechanism in ocean mixing, we argue that it represents only a part of the dynamics.The Stokes drift is sometimes, as by Kantha and Clayson (2004) and Ardhuin and Jenkins (2006), treated as a vertical velocity shear, though it is really only one part of a motion which results in a mean transport.The rest of the motion is important: Benilov (2012) showed that vorticity was unstable for even linear waves, which generate no drift, and Anis and Moum (1995) described how irrotational wave motion could transport turbulent diffusion downward; this mechanism was also suggested by Ölmez and Milgram (1992) in order to explain their experiments.Even the analysis by Teixeira and Belcher (2002), who favoured Stokes drift as the primary mechanism for distortion of vorticity, nevertheless demonstrated that without it the waves affected the turbulence.
It is worth noting here that our proposed alternative also involves adding wave-induced production to the turbulent kinetic energy production in a two-equation model; however, we do not attempt to extract a value for the transport such as Stokes drift or mean shear in order to augment the shear velocity.Instead the production and its vertical scaling are independently determined.Of course in practice the difference in approach is notional rather than real, since any term can always be scaled by a parameter or by altering the empirical closure coefficients, but we argue that the method applied here involves quantities more directly related to the measurements and thus, since we make relatively few assumptions about the underlying physics, is likely to provide realistic scalings.
The proportion of mixing due to Langmuir circulation alone (that is, exclusive of other effects of non-breaking waves) is unknown, and while we do not know how prevalent it is in the ocean we do know that Langmuir cells are not always present.Since an efficient mixing due to wave motion as well as an efficient mechanism for wave attenuation have been observed in laboratory experiments with no reported Langmuir circulation (Dai et al., 2010;Babanin and Haus, 2009;Ölmez and Milgram, 1992), and this mechanism should in principle always exist, we make the assumption that Langmuir circulation is a minor contributor to the mixing.Thus we follow several other authors in applying a vertical distribution of turbulent kinetic energy that does not explicitly take into account Langmuir circulation, and it is to these approaches that we compare ours.On the other hand, our source of energy for mixing is derived from swell dissipation measurements, which means that even if Langmuir circulation is more prevalent than assumed here the total energy going into mixing will not be affected; in this case, our parameterisation can at least be useful as a first-order approximation.In the end, of course, neither a purely Langmuir nor a purely non-Langmuir formulation will perfectly capture the physics.Ideally, a Langmuir-based vertical scaling might be used in conjunction with the one proposed here, with the energy for turbulence production partitioned appropriately between the two mechanisms.
Previous approaches to non-breaking wave-induced turbulence have been applied to improve modelling of sea surface temperatures, sediment concentrations and climatic events.We have discovered some areas where the parameterisations could be improved upon, which we hope will lead to still more successful implementations.Below we describe two successful approaches and explain why we felt the need for devising yet another.Qiao et al. (2004) recognised that ocean mixing is often underestimated in summer and postulated an additional source of mixing due to surface waves.Most prior attempts at implementing a wave-induced effect relied on wave breaking (see for example Chalikov and Belevich, 1993), but this has not been sufficient to explain the observed mixed-layer depths and surface temperatures.Qiao et al. (2004) sought to account for the disparity between models and observations by relating the new source of mixing to the wave spectrum.The new term was then directly added to the turbulent viscosity in the mixing scheme's turbulence model.Unlike turbulence generated by wave breaking, which is injected at the scale of the wave height, the non-breaking wave-induced turbulence is vertically distributed at the scale of the wavelength, comparable to the depth of the mixed layer.Following Prandtl's mixing length hypothesis, they defined the wave-induced component of the turbulent viscosity ν T as

The B v parameterisation of Qiao et al. (2004)
(1) This was developed by analogy with Prandtl's later definition of the kinematic turbulent viscosity, where k is the turbulent kinetic energy, l is a length scale (usually interpreted as the characteristic size of the large, energy-containing eddies) and c is a proportionality coefficient, usually related to the stability of the flow though sometimes set constant (see Wilcox, 1998).
In the equation for B v the length l z is related to the size of the water particle orbits due to the wave motion, and u z is the increment in vertical wave velocity over the distance l z .Combining these definitions results in the following form for B v (a full exposition is provided by Qiao et al., 2010): where E is the wavenumber spectrum, κ is the wavenumber vector, σ is the angular frequency of the waves, z is depth and α is a tuning parameter, assumed to be of order one and set according to the observations.This formula embodies the assumption that the turbulence will be scaled by the wave motion and attenuate according to an exponential law.They applied it to different ocean mixing schemes as well as to model a more controlled laboratory experiment, reporting improved model fidelity in all cases.
The first published application was to the Princeton Ocean Model (Blumberg and Mellor, 1987), which uses the Mellor-Yamada turbulence model to drive mixing (Qiao et al., 2004, andalso Qiao et al., 2010).This is a two-equation model, which means that it describes the evolution of the turbulence by way of two coupled transport equations: one for the turbulent kinetic energy and another for a quantity that may be related somehow to the length scale (Mellor and Yamada, 1982).Crucially, the turbulent viscosity in this kind of model is a function of both of these quantities, and in turn the transport equations both depend on the turbulent viscosity, as well as the turbulent kinetic energy and the length scale.Thus introducing the wave effect by adding directly to the turbulent viscosity adds a source of viscosity without corresponding quantities relating to the creation and destruction of turbulence energy, and this effectively amounts to modifying the turbulence scheme at a fundamental level.All turbulence models are calibrated empirically and therefore are open to modification in the light of new experimental data, but no recalibration seems to have been attempted.That notwithstanding, the addition of B v did lead to a more accurate prediction of ocean mixing in the presence of waves.A more detailed description of the issue of modifying a two-equation model in this way is illustrated in Appendix A.
In a later application B v was added to the turbulent viscosity in the KPP (K Profile Parameterisation) model (Wang et al., 2010), a popular mixing scheme in operational oceanography.It does not feature the kind of interdependence of parameters as two-equation turbulence models such as Mellor-Yamada do.The model produces a vertical profile for the turbulent viscosity, and this is used to determine the turbulent flux of a quantity.Similar to Qiao et al. (2004), the wave-induced turbulent viscosity was added to the KPP turbulent viscosity profile.Because it was unclear whether a wave effect had been implicitly included in the KPP model, they tuned B v via the constant α, choosing a value which gave the most satisfactory match between modelled and observed data.Again, this pragmatic approach led to real improvements, significantly reducing errors in ocean simulations.
The "K profile" (K is the turbulent viscosity, just ν T in our notation) is, for the purposes of the KPP model, divided into three separate regimes: two outer layers -a surface and bottom layer -and the ocean interior.The viscosity profile is matched at the boundaries between these layers, which are themselves defined based on a bulk Richardson number.Forgetting for the moment the bottom layer (in deep water this will not greatly concern us as the wave influence would be expected to vanish), the surface layer and the interior are formulated very differently.The former consists of the product of a cubic shape function and a velocity scale, which is a function of the Monin-Obukhov stability function, and the latter of simplified depth-dependent turbulent viscosities for heat, salinity and momentum.While merely adding B v works, a preferable approach would be to reformulate KPP at a more fundamental level, reconciling the assumptions underpinning its development with the added influence of nonbreaking waves.This implies nothing less than replacing the Monin-Obukhov scaling with something more suitable to the ocean.Additionally, boundary conditions for the boundary layer, taking into account the wave-induced turbulence, would need to be adjusted, as would the Richardson number, which determines the boundary layer depth (although this last one would be done automatically if the Monin-Obukhov stability function were modified).
Since there does not appear to be a ready alternative to the Monin-Obukhov theory, for now we have to conclude that KPP is, strictly speaking, incompatible with the desire to model wave-induced turbulence.Indeed, this is not a new observation: D'Alessio (2002) considered Monin-Obukhov similarity theory and wave breaking, finding that when waves are breaking the Monin-Obukhov functions no longer describe the ocean boundary layer accurately.Yet he also remarked that despite this discrepancy, KPP was nevertheless successful as a mixing model.In his study, though, the deep mixing characterised by turbulence generated by the wave orbital motion was absent and, as he suggested, since wave breaking is confined to the near-surface layer, the difference could be insignificant.
A third application of the B v parameterisation was to a coupled ocean-atmosphere climate model in simulations of the South Asian Summer Monsoon (Song et al., 2012).
Here the mixing scheme used was that of Pacanowski and Philander (1981), which is much simpler than either a twoequation model or KPP.It consists of separate turbulent viscosities for momentum and scalars, which are empirically determined and dependent on a constant background value, in addition to a term based around the Richardson number.B v was added to these viscosities, and once again the model results were much closer to the observations.Being a far less sophisticated mixing model, adding the wave-induced viscosity is unlikely to be as problematic from a theoretical standpoint, although we note that it was calibrated without waves and it is possible that including them might lead to different values of the constants, or a different form for the equations.On the other hand, given the controlled nature of the original calibrations, this may not make a large difference (Pacanowski and Philander, 1981).
Despite the issues we have illustrated, in all three studies big improvements were made in model performance.We note also that the B v model was applied with considerable success in modelling a controlled laboratory experiment for temperature diffusion in a wave tank (Dai et al., 2010).

Pleskachevsky et al.'s (2011) wave shear frequency approach
Pleskachevsky et al. ( 2011) derived their model quite differently.Their approach dealt simultaneously with both a modification to the turbulent viscosity of the flow and an augmentation to the turbulence production by wave-induced shear.Both effects were incorporated in their mixing model, using a two-equation turbulence model.For geophysical models where horizontal gradients vary slowly and vertical velocities are small, the production term for free-shear flows may be parameterised as (see Eq. A1 for the full equation) where M is the shear frequency.P thus represents the rate of production of turbulent kinetic energy by the mean current shear.They augmented the production by considering the addition to the shear from the waves, so P becomes P + P wave .P wave is the wave-induced shear production, parameterised as (5) and added to P .E AM is the ratio of wave energy lost to dissipation (or, equivalently, the production of turbulence) and the total wave energy.The wave mean shear frequency, M wave , is defined a bit differently to the usual shear frequency M, involving not just horizontal shear but vertical shear as well, in order to take into account the orbital motion of the water under the influence of the waves.It is also averaged over a wavelength in order to match the scale of the mean currents (the wave motion is assumed periodic in a timescale greater than the turbulence, but less than the mean currents).The wave shear frequency is expressed as The quantities U wave and W wave are component mean velocities of the wave motion, that is they represent the orbital motion with the turbulent fluctuations averaged out.The mean wave shear frequency is then the mean of M wave over a wavelength.
To augment the turbulent viscosity in these expressions, they defined a wave-induced turbulent viscosity which they added to the turbulent viscosity ν T ; it is proportional to the wave mean shear frequency and the square of a length scale, In a similar vein to the B v model, they assume that the length scale l wave is equal to the amplitude of the water particle orbits, a function of depth.However, in contrast to the B v model, Pleskachevsky et al. (2011) do not treat the full spectrum, and instead they work with bulk properties such as significant wave height, so they only deal with one value for l wave .
One of the problems with this approach is the same with that of the first B v implementation, which is that the waveinduced turbulent viscosity is added directly to the turbulent viscosity term in a two-equation model.However, unlike the B v model, this wave shear model also includes a turbulence production term, which both creates an extra problem and an opportunity.On the one hand, now that we are not only directly augmenting the viscosity but producing turbulence, the viscosity is effectively being counted twice.On the other hand, with an explicit production term we should not need to invoke any explicit addition to the turbulent viscosity, so we ought to be able to rectify the problem by simply removing the viscosity term.Note that this in itself is not to say that the formula proposed for turbulent viscosity is incorrect, but instead that it is not appropriate in the context of a two-equation turbulence model, other than perhaps as a diagnostic tool to check that the production is indeed leading to the right value for the turbulent viscosity.In other kinds of models, where the turbulent viscosity is not explicitly dependent on the turbulent kinetic energy, it could still be a useful formulation.
We applied this model for the wave-induced turbulence into the GOTM one-dimensional mixing model in, so far as we know, the same manner as that described by Pleskachevsky et al. (2011).This meant using the parameters pre-selected in the "wave breaking" test case supplied by the model developers, and using the k-model after Rodi.Using their formulas, however, it became clear on running the model that the turbulence production due to the waves, Eq. ( 5), was entirely negligible, meaning that the mixing was completely due to the direct addition of the wave-induced turbulent viscosity -the very term we had intended to remove on theoretical grounds.It would appear that the ratio E AM , which Pleskachevsky et al. (2011) set to 1.5 × 10 −4 , is too small.They indicated, however, that the variable could take a range of values, anywhere from 1.0 × 10 −3 to 1.0 × 10 −7 , depending on wave steepness.As a test we tried the higher value, and indeed the wave-induced mixing, now without any direct augmentation to the turbulent viscosity, became very strong.This suggests that a recalibration might fix the model, though we elected to defer this task.
What we are left with, then, is a direct parameterisation for the turbulent viscosity.It is easily shown that this is in fact identical to the monochromatic form for B v described by Qiao et al. (2010), provided that α = 1.In deep water, this is where A is the amplitude of the waves, and κ the wave number.This parameterisation was applied to a sediment suspension model and resulted in vastly improved modelling (Pleskachevsky et al., 2011); it has also been used to explain the under-prediction of mixed layer depth when waves are omitted from a mixing model (Toffoli et al., 2012).

Further comments
Despite some shortcomings we have discussed, the parameterisations explored in the previous section have led to improved model predictions, and it seems clear from previous studies that including some kind of extra turbulence due to the waves is preferable to leaving it out.Since the Pleskachevsky model as it stands is identical to B v for monochromatic waves, the work described above effectively leaves us with only one model.This describes the turbulent viscosity rather than the production of turbulence, so it is not really suited for use in a two-equation model.We therefore feel that there is room for another method, one that calculates turbulence production and which will work naturally with these models, and many other Reynolds-averaged turbulence models.Since the generation of wave-induced turbulent kinetic energy must by definition involve a transfer of energy from the waves, the same parameterisation should predict a component of the wave dissipation.Thus in Sect. 4 we test an alternative model, developed in terms of wave dissipation through turbulence generation.
4 One-dimensional modelling of the stratified ocean 4.1 An amplitude-scaled turbulence production parameterisation Babanin and Haus (2009) performed laboratory experiments with monochromatic waves and observed that the velocity spectra contained regions closely approximating the Kolmogorov interval.By measuring the spectrum at a depth just below that of the wave troughs, they determined a form for the dependence on amplitude of volumetric turbulent kinetic energy dissipation near the surface, where a 0 is the surface wave amplitude, and p = 3.0 ± 1.0 within the 95 % confidence limit, and b = 300.An earlier hypothesis (Babanin, 2006) on dimensional grounds supports the scaling with the cube of the amplitude, so p is set equal to 3. (This follows from a definition of a wave Reynolds number where the inertial forces scale with the square of amplitude, where V is the orbital velocity of the fluid particles, a is the amplitude of the motion (which scales exponentially with depth) and ν is the molecular viscosity of water.)Since the value of p = 3 is assumed to be constant, the quantities in expression Eq. ( 9) are dimensional.The dimensionless form was described by Babanin (2011), and in it b becomes where b 1 is a dimensionless parameter initially assumed to be constant.Turbulence was observed not more than once out of every ten wave periods for the steepest waves in the experiment, which were just below the breaking threshold.The value of b = 300 was determined for these isolated patches of turbulence, so over time the average value would be less than a tenth of this.(This value and the intermittent nature of the turbulence were found to be consistent with numerical experiments by Babanin and Chalikov, 2012.)Since the waves in the experiment were steep it was further suggested that for ocean swell the number could be reduced further.Ardhuin et al. (2009) had measured swell dissipation rates in the ocean and attributed them to a negative wind stress.By assuming instead that turbulence production in the ocean is the dominant sink for wave energy in swell, Babanin (2011) compared these swell dissipation rates with a depthintegrated form for the turbulence dissipation.The fit was good, except that it tended to over-estimate the swell dissipation, so a revised estimate of b was that it should be perhaps one twentieth of its laboratory-measured instantaneous value.Thus reckoned, b 1 was chosen to be 0.002.(We note that more recent measurements by Young et al., 2013 suggest a slightly smaller value of between 0.0014 and 0.0017, not too dissimilar from the results of Babanin, 2011, especially given that these values assume a constant value for b 1 .)Assuming a steady state where the turbulence dissipation and production are equal, the final form for the turbulence production due to the wave orbital motion is where H is the wave height, and z is measured upwards from the mean surface.On closer inspection the data appears to exhibit a quadratic dependence on wave steepness, thus b 1 becomes This form turns out not to match the lab experiments -it is far too high, in fact.There are several possible reasons for this: firstly, in the ocean there is always turbulence, whereas in the lab the water was left to rest.Since the rate of production of turbulence depends on the turbulent viscosity we would expect that the rate of production of oceanic turbulence would be much higher than the lab experiment for otherwise equivalent waves.Secondly, it is possible that the wave orbital motion is not the only, nor even the primary, source for swell dissipation.Thirdly, we have assumed that the ocean waves are monochromatic.The tank experiments were for monochromatic waves, but it is not entirely clear how a realistic ocean spectrum would behave.A fourth reason is, perhaps obviously, that the quadratic dependence is only really valid within the range of the data used to find it.All that really remains to be said is that the dimensional coefficients to b 1 were deduced independently of the experiments and are consistent with them, so we do not anticipate the problem to lie there.
The formula was adapted from experimental data for monochromatic waves, and in using it we assume that the wave spectrum is narrow enough such that we can use characteristic wave parameters such as significant wave height and peak period instead of specifying the production based on a complete wave spectrum.In general, particularly for the case of swell with little wind stress, this is a good approximation; however we note that, due to the cubic amplitude dependence, in a random wave field it may be that the largest waves, and not the significant waves, dominate the mixing.
The formulation of b 1 could well be enhanced with more swell data, but in its favour this model is based solely on the properties of the waves (there is no wind dependence, for example), it is scaled with the complete orbital motion (rather than the Stokes drift or some other property) and it models the turbulence production explicitly -and exclusivelyrather than the turbulent viscosity.We note that the model, with a different (and constant) b 1 based on much earlier estimations of swell dissipation, had been derived using different assumptions by Bowden (1950).We believe that, so far, this is the most likely formulation to succeed in a general ocean mixing model, especially, since it models turbulent kinetic energy production, when multiple sources of turbulence are introduced.There is a noteworthy deficiency, however, which is that the present model does not depend on turbulent viscosity.Indeed, in order to do so we would need either a more intimate knowledge of the physics which is thus far lacking, or a set of careful measurements showing how the rate of turbulence production changes with turbulent viscosity, also unavailable.One possible approach is to assume the validity of the wave mean shear approach (Eq.6) and attempt to find an empirical closure coefficient.On the other hand, the observational data from Young et al. (2013) suggest that the parameterisation does, on the mean, a good job of modelling the rate of swell dissipation in the ocean, so for ocean modelling applications it should be adequate.As the model time step is decreased, however, the model would probably benefit from this kind of refinement.

The mixing model and numerical experiments
To simulate the ocean mixing we use the one-dimensional mixing model GOTM (Umlauf and Burchard, 2005;Umlauf et al., 2011).This model is easy to set up, modify and use with a variety of different turbulence closure schemes and inputs, and can be coupled to other models, such as circulation models.It is not, however, a general turbulence model because the momentum transport equations are simplified so that they can efficiently model the oceanic boundary layer (Burchard, 2002) (if this were not the case the model would also have to be three-dimensional).This places some limits on the kind of approach needed to introduce a wave-induced turbulence production term into the mixing scheme, in particular with the method of Pleskachevsky et al. (2011) which uses the wave mean shear frequency, M wave , and therefore requires a horizontal shear as well as a vertical one.This means that we must add turbulence production due to M wave at a level higher, where it is completely parameterised based on mean wave properties rather than directly calculated from the flow.
While we do not anticipate that simply adding another source to the production term will be a problem, it is worth at least noting that this approach implicitly assumes that wave-induced turbulence is homogeneous on the scale of the waves.In their large eddy simulations Babanin and Chalikov (2012) showed that the turbulence is highly dependent on the gradients of the orbital velocities.Thus turbulence production occurs mostly in a narrow column immediately below the wave peaks.The turbulence may also be intermittent, as reported by Babanin and Haus (2009).This means that in order to include it in the GOTM model it must be presented in a space and time-averaged form, where the space and timescales are greater than the wavelength and period.In our application of the parameterisations this has been implicitly assumed, and is no different in that way to parameterising, say, turbulence from intermittent wave breaking.
The model setup was, just as with Pleskachevsky et al. (2011), based on the wave breaking scenario downloaded from the GOTM website.We removed the surface momentum flux which is meant to represent the effect of breaking waves and, in addition, changed the surface boundary condition from tke injection (intended for breaking waves) to logarithmic law of the wall.We tested two parameterisations: B v (the only one remaining after we discarded the others considered in Sect.3.3) and the new one described in Sect.4.1.
In the context of a turbulence model, we assume that the wave orbital shear behaves in a similar way to the mean shear, so for the present parameterisation the wave-induced turbulence production Eq. ( 12) was added directly to the shear production term, which is calculated in the GOTM production subroutine.Because we are using a twoequation turbulence model this is the only modification that should be made, apart from adjusting closure coefficients, especially the one relating to production of dissipation.how sensitive the model would be to this we varied this coefficient as far as the numerics would allow without generating errors and the differences, while measurable, were not so great as to overturn the main conclusions of this study; hence, without any field or laboratory measurements to guide us, we kept the default value.The B v model parameterises the turbulent viscosity directly, so is added to the viscosities calculated in the kolpran routine; we add the same value for B v , with α = 1, to the turbulent viscosities for momentum, heat and salinity.By eliminating all other sources of turbulence production this renders the turbulence transport equation irrelevant in the mixing model.More care would have to be taken if B v were used and another source of turbulence production were present, since the turbulent viscosities would not simply be additive.As we have noted above, our parameterisation lacks dependence on turbulent viscosity, but we do not anticipate that this will be of great importance in ocean applications except at the very highest model resolutions.
The two-equation turbulence scheme used is the Rodi kformulation.Tests with other two-equation models yielded results which differed a little, but compared with the difference between incorporating the wave-induced effect or omitting it entirely these differences were negligible.Apart from the waves there was no other forcing -no input or output of heat, no internal waves, surface wave breaking, biology, wind, currents or external sources of pressure.
The initial conditions consisted of prescribed salinity and temperature profiles which we sourced from the Global Temperature and Salinity Profile Program (GTSPP, 2006).We wanted to examine mixing due to non-breaking wave generated turbulence in isolation from other mixing mechanisms, so we sought a stable, stratified profile at mid-to-high latitudes in summer, so as to avoid the issues with convective mixing common in wintertime.It was also desirable to have something in deep water, well away from continents, as some of the wavelengths we tested were over 600 m.The profile we finally chose came from station 8541272, located some distance west of Chile at 115.9200 • W, 45.8840 • S, and measured at 03:24 UTC on 3 January 2010.The depth of the measured profile was 1813 m, and the model was set up to run with that depth of water (though for clarity we only show down to 250 m in the plots).
The results of the mixing on the temperature profile are shown in Figs.1-4 (a similar effect may be seen on the salinity profiles, but we omit the plots).In each the dashed line is the original temperature profile, and the solid lines are, from darkest to lightest, in two-hour intervals, the temperature profiles as the water column is mixed by turbulence induced by the waves.We ran the experiments for several different wave heights and wave numbers.In Fig. 1 the simulations were performed using the parameterisation Eq. ( 12) for a set of wave numbers and wave heights which sit within the range of data used to develop the parameterisation.The wave heights are 3, 4 and 5 m, and the wave numbers 0.014 and 0.02, corresponding to wave lengths of 449 and 314 m respectively.In Fig. 2 a wider range is used, where we have assumed that the parameterisation still holds.For this set of experiments the three wave heights are 2, 5 and 10 m and the two wave numbers 0.1 and 0.01 rad m −1 , corresponding respectively to wavelengths of 63 and 628 m.Note that tenmetre waves with a wave number of 0.1 are just above the breaking threshold, meaning that in any real situation there would be near-surface turbulence generated by wave breaking as well.Since we are investigating a source of turbulence in isolation (no currents, wind stress or wave breaking) these tests can be thought of as instructive, rather than a complete description of a column of ocean.
The results suggest that even for small H the wave motion may contribute substantially to the mixing.For wave heights of 5 m and a wave number of 0.02 it deepens the mixed layer by as much as a factor of four, and reduces the surface temperature by about three degrees after twelve hours.The largest part of the effect occurs in the first few hours, and significant cooling occurs after only minutes (not shown).Even the smaller, shorter waves have an impact on the mixed layer.Assuming an appropriately shallow mixed layer, an effect of this magnitude could strongly influence the intensity of tropical cyclones, especially where the sea surface temperature is close to the threshold for tropical cyclone generation.It may also be seen that the mixing extends much deeper than we would expect from wave breaking, which generates turbulence only to depths of the order of the wave height.
For the B v turbulent viscosity parameterisation Eq. ( 8), the results corresponding to the wave numbers and heights in Fig. 1 are shown in Fig. 3, with results for the wider range of wave parameters plotted in Fig. 4. The mixing is weaker, especially for the steeper waves, and qualitatively different: instead of a clear thermocline the temperature gradient is reduced.(This degradation of the thermocline was observed and discussed in more detail by Martin et al., 2011.)Even so, the surface temperature decrease is of the same magnitude as the new parameterisation with smaller waves, which are not far outside the range of the data of Ardhuin et al. (2009), and especially for shorter times; this is encouraging given that the two parameterisations tested here are for different quantities (one is the rate of production of turbulence, the other of turbulent viscosity), and were developed along different theoretical lines and validated with different kinds of observational data.For much larger waves, though, the results diverge greatly, and for the short waves of κ = 0.1 the B v turbulent viscosity is insufficient to affect the mixed layer depth, even for large wave heights -this is in no small part due to the initial depth of the mixed layer, and we would presumably have seen an effect if it were even shallower than it is.
The qualitative differences can be explored even further when we look at the plots for the turbulent viscosity in Fig. 5.For the B v parameterisation (the dotted curve in the plot) this is calculated directly from the wave parameters, which here are constant, and therefore do not change with time, while for the parameterisation of turbulence production the turbulent viscosity is calculated from the transport equations for turbulent kinetic energy and dissipation (see Eq. A4 in Appendix A).Its value depends not just on the wave-induced turbulence production but also on the stratification of the water column (hence the small bulge between 50 and 100 m depth in the plots, which corresponds to the end of the thermocline where the lower stability leads to increased mixing).
The viscosities generated by the two-equation model are a very different shape to those generated by the B v equation.While the B v model scales exponentially from a maximum at the surface, the two-equation model calculates the turbulent viscosity such that it is zero at the surface and there is a maximum at some depth.This is what we would expect since the size of the eddies near the surface should be constrained by the surface, leading to a reduced ability to mix.It should be noted that this characteristic is a result of using a logarithmic surface boundary condition in the mixing model, which was primarily intended for vertical shear generated turbulence, but the result is in any case more realistic than a viscosity that increases continuously to the surface.Perhaps when more is known about the physics of this source of turbulence a more appropriate boundary condition may be formulated for inclusion in the mixing model.
Extrapolating from these academic tests to more realistic ocean modelling applications is not necessarily straightforward, and we offer the following observations to explain our results more clearly in that context, particularly with regard to those situations where the two parameterisations produce disparate results.In the present idealised tests the parameterisations are used alone with no other sources of turbulence (except for a very small background turbulence required by the k-model).In a full ocean simulation they will be used in conjunction with other sources of turbulence; for the B v parameterisation this will mean that the value of B v will be added to whatever turbulent viscosity already exists, whereas our new parameterisation will add to the rate of turbulence production.The net result will be quite different for two reasons: firstly, adding turbulence sources does not result in an analogous addition of turbulent viscosities -put simply, doubling production does not double viscosity (an inspection of the turbulence transport equations in Appendix A will make The dotted line is the constant-in-time B v parameterisation, while the solid lines are calculated from the parameterisation based on wave dissipation (Babanin, 2011) using the k-model of Rodi (1987).The darkest line is that calculated two hours after the simulation started, with subsequent two-hour intervals indicated by progressively lighter lines.this clear).In more realistic applications, then, B v may be more effective than it seems based on some of the results shown in the figures.The second reason is that because the B v profile increases near the surface rather than decreasing to zero as is normally expected, it can help mix the very top layer of ocean into the deeper mixed layer even if mixing near the surface by currents alone is relatively inefficient.

Conclusions
While there are several different approaches available to modellers wishing to incorporate the effect of non-breaking wave-induced turbulence in their forecasting models, on closer inspection several of these cannot deal with situations of zero wind stress and here, in the context of wave-only mixing, had to be discarded from consideration.Leaving aside models invoking Langmuir circulation, which seems to represent only a part of the process, we identified two models which have been used successfully in modelling studies.One of these, the model proposed by Pleskachevsky et al. (2011), requires modification and at present represents a special case of the B v model.This latter model, first proposed by Qiao et al. (2004), models not the turbulence production but the turbulent viscosity.As an alternative to these we have proposed using the turbulence dissipation measurements of Babanin and Haus (2009) and the swell dissipation measurements of Ardhuin et al. (2009) as a source of turbulence production.Our tests suggest that, in an environment without any other sources of turbulence, the two models give broadly similar results in some situations, though with some important qualitative differences, especially with regard to the position and slope of the thermocline.However in other, more extreme situations their behaviour diverges considerably.What is needed now is more laboratory data in order to refine the theory, which is still not completely understood, and a good data set of simultaneous field measurements against which model results can be compared.
For the broader question of modelling the ocean and the atmosphere, our rather academic simulation was designed to test the efficacy of this source of turbulence in deepening a shallow mixed layer, and we have deliberately left out other sources of turbulence to get a measure of its effect in isolation.Of course in the ocean several turbulence sources coexist, including tangential wind stress, internal wave breaking, currents and the other surface wave-related sources: wave breaking and Langmuir circulation (though this latter source likely has been implicitly included in our model, at least in terms of the total energy), and a physical description of the ocean-atmosphere system is impossible without accounting for large-scale three-dimensional effects.Be that as it may, our simulation results, in the absence of surface fluxes and other sources of turbulence, suggest that the turbulence generated by the wave orbital motion, as reckoned from assuming that it is the primary mechanism for swell decay, can deepen an initially shallow mixed layer and reduce sea surface temperatures by, potentially, several degrees.Even a much smaller change would be significant in many modelling and forecasting applications, which indicates that it ought to be featured in mixing models -especially when modelling scenarios where there exists persistent or largeamplitude swell.This also implies that a wave model should be a component in any coupled model system where oceanatmosphere feedbacks are present.This finding aligns well with modelling studies by Jacobs (1978), Qiao et al. (2004Qiao et al. ( , 2010)), Wang et al. (2010), Huang and Qiao (2010), Huang et al. (2011) and Song et al. (2012).All showed improved correspondence between observations and models with the inclusion of turbulence attributed to the wave orbital motion, regardless of the details of the method chosen to parameterise it.We would anticipate that the model we have proposed, because of its more general approach, should aid future modelling studies.
As a final note, we would like to point out that from a modelling perspective the differences between the present model and other models such as B v , and even models based on Langmuir circulation such as that by Kantha and Clayson (2004), are more theoretical than practical at this stage.They all invoke waves as a source of turbulent kinetic energy and a way of distributing it through the water column, so in a sense they are all equivalent even if they are motivated by different physical assumptions.While the scaling in the present model is derived from experiments and theoretical assumptions which may not be always be applicable in the open ocean, the source of energy is, uniquely, explicitly quantified through measured swell dissipation, providing at least a cap on the available energy for mixing.Our 1-D tests when restricted to the regime for the available swell dissipation data (Figs. 1 and 3) indicate that the resultant mixing from the present model is comparable, though greater, than that produced by the B v model, which has been successfully applied in several different scenarios.Further research is needed to determine both the dominant physics of the mechanism and the partition between this and other possible sources of swell dissipation, such as that invoked by Ardhuin et al. (2009).

Two-equation models and artificially augmenting the turbulent viscosity
While there are many kinds of turbulence models, the socalled two-equation models occupy a special position in the modelling community.They are far more computationally efficient than large eddy simulations (which are too costly for many geophysical applications), and unlike algebraic and one-equation models they are "complete", which means that no prior knowledge of the flow need be had.They are the simplest such complete models, and their performance is not much hindered by their simplicity (Wilcox, 1998).There are a good number of formulations, but they all consist of two interdependent transport equations, one of which describes the evolution of the turbulent kinetic energy and the other that of a length-scale determining variable.The turbulent kinetic energy equation for an incompressible fluid takes the following form: where the rate-of-strain tensor is defined as Although some terms are more or less important depending on the application, this is the equation generally used (in simplified form) in geophysical turbulence modelling, and its essential features can be derived from the Reynolds stress equation.On the left-hand side we have the material derivative of the turbulent kinetic energy, k.On the right-hand side the first term is the production of turbulent kinetic energy due to shear stress in the fluid; the second term is the production due to buoyancy, where ρ is the potential density; the third term is a gradient diffusion term; and is the dissipation.ν T is the turbulent viscosity for momentum, and for scalars determining the buoyancy (salinity and temperature) the diffusivity is proportional to ν T , hence the scaling parameter σ T .The closure coefficient σ k is set empirically.Dimensional arguments may be used to determine that where l is a length scale, usually interpreted to be the characteristic size of the largest eddies.Further dimensional analysis provides an expression for the turbulent viscosity: where c µ is another empirical parameter.From substituting Eq. (A3) into Eq.(A4), it becomes clear that l and are interchangeable quantities, since ν T depends on k and only one of l or .It is also understood that they vary with the flow.
The motivation behind a two-equation model is to construct a second transport equation in order to determine the length scale l, which in turn enables us to fully determine the turbulent viscosity and the turbulent kinetic energy transport.The primary difference between the various two-equation models boils down to the choice of the variable used for the second transport equation, be it, for example, the dissipation rate as in the k-models, the specific dissipation rate, ω, as in k-ω models, or the product k l, as in the Mellor-Yamada model found in the Princeton Ocean Model used by Qiao et al. (2004).Any variable can in principle be used as long as it can unambiguously determine the length scale, although all of these models for the second equation are empirical, so they can differ substantially in both form and behaviour.
The dissipation transport equation in the Rodi k-model used in the numerical experiments in Sect.4.2 is (Rodi, 1987): where C 1 , C 2 , C 3 , σ are empirical constants.Note the existence of production terms here too; the right-hand side is analagous to that of Eq. (A1).Naturally any additional turbulence production due to waves must be present here just as for the equation for k.Now the main point we wish to make relates to Eq. (A4), which expresses the turbulent viscosity in terms of the turbulent kinetic energy k and the dissipation .Both of these quantities vary according to the rather complicated transport Eqs.(A1) and (A5), and both of the transport equations are functions of these quantities both explicitly and through the turbulent viscosity, so to simply add wave-induced turbulent viscosity to that obtained from other aspects of the flow (such as the current shear) not only amounts to effectively reformulating the equations, but it results in a kind of double counting of turbulent viscosity.Adding an extra source of turbulence production should be sufficient.This production term will appear in both the transport equations for turbulent kinetic energy and the dissipation (or another length-scale determining variable).Following Kantha and Clayson (2004), the only other change that one should wish to make would be to adjust the empirical coefficients.It is possible that this sort of approach might help to further improve the model presented in this paper, but it should be noted that in their parameterisation the augmentation to the length scale was necessary, since Stokes drift on its own adds relatively little to the current shear, and the primary production mechanism in their model, Langmuir circulation, involves length scales rather larger than the wave orbital motion.Hence we do not anticipate that in our model the adjustment will be as critical.

Appendix B
The experimental paper by Beyá et al. (2012) When Beyá et al. (2012) attempted to replicate the experiment of Babanin and Haus (2009) they did not find any evidence of wave-induced turbulence.Based on their observations and other arguments they concluded that no such effect exists.Since that paper has been gaining traction we felt that it was appropriate to explain why it does not invalidate the hypothesis of non-breaking wave-induced turbulence.Beyá et al. (2012) appear to have confused the papers by Babanin and Haus (2009) and Babanin (2006).They claim that Babanin and Haus (2009) established the Reynolds number threshold of 3000 for transition to fully developed turbulence, whereas it was actually done by Babanin (2006).Babanin and Haus (2009) investigated intermittent turbulence and could not confirm the threshold value either.Also, it was Babanin (2006) who conducted the laboratory experiments to illustrate the effect of fully developed turbulence on the dissolution of dye, but the experiments by Beyá et al. (2012) replicated the experimental conditions of Babanin and Haus (2009) where only intermittent turbulence was observed.Thus they did not and could not have seen the rapid dissolution of dye observed by Babanin (2006).
Secondly, Beyá et al. (2012) argue that the fact that Stokes waves are observed to high precision precludes the existence of vorticity, and hence turbulence.This is an apparent misreading of the basic theory of viscous wave motion where the viscosity, which leads to vorticity, reveals itself in the imaginary part of the solutions.This part leads to a slow decay of Stokes waves, but does not affect in any way the shape of such waves (see for example Kinsman, 1965).
Thirdly, the experiment conducted by Beyá et al. (2012) itself has serious shortcomings.In their experiment they observed dye in the presence of waves over 5 wave periods, for wave trains of 3 different steepnesses.The dye will disperse due to molecular diffusion even in the absence of turbulence, thus in order to determine the existence of turbulence a comparison between the rates of dissolution of dye with and without the presence of waves would be required.They made no such comparison, and indeed no measurement of molecular diffusion was attempted.Instead, the existence or nonexistence of turbulence was decided by a more subjective criterion, being whether or not the experimenters had judged Nonlin.Processes Geophys., 21, 325-338, 2014 www.nonlin-processes-geophys.net/21/325/2014/ the dye to have sufficiently dispersed.The weak, intermittent turbulence observed by Babanin and Haus would only be expected to create a visual difference after several hundred periods (see also Dai et al., 2010).Beyá et al. (2012) only observed their dye over 5-10 wave periods, so they could not have seen any significant effect visually.On the other hand, they did observe a reduction of thickness of their dye line, but their explanation contradicts their own data.In their Fig. 7, vertical profiles of dye line thickness are shown at zero time and after 5 wave periods, for three wave steepnesses: 0.17, 0.21 and 0.24, from the surface down to 7 cm below the surface.The dye line became twice as thin at the 6-7 cm level over the 5 periods, but notably not at the surface and not for the 0.21 steepness case.Beyá et al. (2012) explain the observed reduction of thickness by the effect of Stokes drift stretching the dye line.This is not supported by any estimates of expected thinning by this mechanism, but should it be the case then the thinning cannot be increasing away from the surface where the Stokes drift is decreasing exponentially.Furthermore, the Stokes drift could not have influenced the dye thickness at steepnesses 0.17 and 0.24 without affecting the dye at steepness 0.21.On the contrary, Fig. 7 of Beyá et al. (2012) much more likely reflects the impact of an intermittent turbulence in their experiment, and not the Stokes drift.Unless this figure is upside down, however, this turbulence could not have come from the waves and must have come from the bottom.
Finally, Beyá et al. (2012) comment that their observations of velocity time series do not "show any evidence of turbulence fluctuations", but they do not substantiate the claim.The original thesis (Beyá, 2011), however, provides more details on their observations of kinematics by acoustic doppler velocimeters (ADVs).The sampling frequency of the ADV was 50 Hz, which means that the Nyquist frequency of 25 Hz was below most of the frequency range of the Kolmogorov interval observed by Babanin and Haus (2009).More importantly, the Kolmogorov interval as described by Babanin and Haus only appeared over a fraction of the wave period and cannot be seen if the Fourier Transform is applied to the continuous time series.Beyá et al. (2012) did not attempt to locate segments of the intermittent turbulence, and therefore they could not have seen the Kolmogorov spectra and the turbulent oscillations even if their sampling frequency and resolution had been high enough.

Fig. 1 .
Fig. 1.Mixing in GOTM due to non-breaking wave-induced turbulence using the parameterisation for wave dissipation from Babanin (2011).The wave heights and wave numbers correspond to typical values from the data of Ardhuin et al. (2009).The dashed line is the initial temperature profile, and the solid lines for the profile at two-hour intervals, with lighter lines indicating more time having elapsed.

Fig. 2 .
Fig. 2. As for Fig. 1 but for a wider range of wave heights and wave numbers.

Fig. 3 .
Fig. 3. Mixing in GOTM due to non-breaking wave-induced turbulence using the B v parameterisation for monochromatic waves with α = 1 from Qiao et al. (2010).The range of wave heights and wave numbers is as for Fig. 1.

Fig. 4 .
Fig. 4. As for Fig. 3 but for a wider range of wave heights and wave numbers.

Fig. 5 .
Fig. 5. Turbulent viscosity due to non-breaking wave-induced turbulence for a wave height of 5 m and a wave number of 0.02 rad s −1 .The dotted line is the constant-in-time B v parameterisation, while the solid lines are calculated from the parameterisation based on wave dissipation(Babanin, 2011) using the k-model ofRodi (1987).The darkest line is that calculated two hours after the simulation started, with subsequent two-hour intervals indicated by progressively lighter lines.