2005 Cambridge University Press doi:10.1017/S0022112005007044 Printed in the United Kingdom

*A numerical study on the evolution and*

*structure of a stress-driven free-surface*

*turbulent shear ﬂow*

*By W U -T I N G T S A I1,2 _{, S H I -M I N G C H E N}*2

*A N D* *C H I N -H O H M O E N G*3

1_{Institute of Hydrological Sciences, National Central University, Taoyuan 32001, Taiwan}
2_{Department of Civil Engineering, National Chiao Tung University, Hsinchu 30050, Taiwan}

3_{National Center for Atmospheric Research, Boulder, Colorado 80307, USA}

(Received7 October 2004 and in revised form 20 June 2005)

Turbulent shear ﬂow beneath a ﬂat free surface driven by a surface stress is simulated numerically to gain a better understanding of the hydrodynamic processes governing the scalar transfer across the air–water interface. The simulation is posed to mimic the subsequent development of a wind-driven shear layer as in a previous experiment except that the initiation of the surface waves is inhibited, thus focusing on the boundary eﬀect of the stress-imposed surface on the underlying turbulent ﬂow and vice versa. Despite the idealizations inherent in conducting the simulation, the computed ﬂow exhibits the major surface features, qualitatively similar to those that appear in the laboratory and ﬁeld experiments. Two distinct surface signatures, namely elongated high-speed streaks and localized low-speed spots, are observed in the simulated ﬂow. Including temperature as a passive tracer and describing an upward heat ﬂux at the surface, we obtain high-speed streaks that are colder and low-speed spots that are warmer than the surrounding regions. The high-low-speed streaks, arranged with somewhat equal cross-spacing of centimetres scale, are formed by an array of streamwise jets within the viscous sublayer immediately next to the surface. Beneath the streaks, counter-rotating streamwise vortex pairs are observed among other prevailing elongated vortices. However, they are signiﬁcantly shorter in length and more irregular than their corresponding high-speed streaks at the surface. Accompanying the more organized high-speed streaks, localized regions of low streamwise velocity emerge randomly on the surface. These low-speed spots are attributed to strong upwelling ﬂows which disrupt the viscous sublayer and also bring up the submerged ﬂuids of low streamwise velocity. The occasional interruptions of the streamwise high-speed jets by the upwelling ﬂows account for bifurcation or dislocation of the surface streaks. Statistics of the turbulence are presented and their implications for the formation of the ﬂow structures are discussed.

*1. Introduction*

The layer of water adjoining and including the air–sea interface controls the exchange of slightly soluble gases, such as carbon dioxide, between the atmosphere and the ocean. Such a liquid-phase controlled process is determined by the near-surface turbulence in the ocean side, of which the dominant forcing is the wind stress on the interface. The energetic, coherent vortical motions within the aqueous turbulent boundary layer, or known as the surface renewal eddies, therefore account for the

eﬀectiveness of transport across the air–water interface. Through interacting with the boundary, these coherent eddies also induced structural patterns on the water surface. A popular approach in recent years to estimate the air–water transfer rate is to infer the renewal events from remote-sensing images of water surface (e.g. Haussecker et al. 2002; Garbe, Schimpf & J¨ahne 2004; Schimpf, Garbe & J¨ahne 2004). The accuracy of the technique thus relies on relating quantitatively the small-scale surface signatures to the underlying turbulence processes.

The observation by Woodcock (1941) is among the pioneering experiments to study the small-scale motions on the water surface. He observed that when winds blow, high-speed surface water movements occur that subsequently form surface streaming in lines and streaks roughly parallel to the wind direction. The water in the streams moves down-wind much faster than that between the streams. The streaks are approximately a centimetre across and with a distance of the order of centimetres between the streaks. The subsequent visualization experiments of Gemmrich & Hasse (1992) and Kenny (1993) all observed similar small-scale surface streaks on natural water surfaces. These high-speed streaks, which are the most frequently appearing patterns of surface water movement at wind speeds of 2 or 3 cm s−1 or higher, were taken by Gemmrich & Hasse (1992) as indication of surface renewals that might enhance air–sea gas exchange. They interpreted the appearance of these surface streaks with a cross-spacing of the order of centimetres either as a result of wave breaking or as an indication of the existence of counter-rotating helical vortex rolls similar to those in Langmuir circulations, which have a typical spacing of metres to kilometres (Langmuir 1938; or see the reviews by Leibovich 1983 and Thorpe 2004). The observation of Kenney (1993), nevertheless, indicated the contrary. No evidence was found of ﬁxed circulations beneath the streaky surface. In addition, no signiﬁcant surface waves were observed and the streaks could appear on a calm water surface at low wind condition, indicating the absence of wave–current interaction for generating Langmuir circulations.

In a laboratory experiment by Melville, Shear & Veron (1998) (and also Veron & Melville 2001), predominant streaks aligned in the wind direction were also observed at a wind-driven water surface. The streaky surface undergoes several stages in its evolution which begins with a quiescent water surface at the onset of the wind. The appearance of the surface streaks was shown to follow closely the initial growth of the surface waves and the inception of the Langmuir circulations. (In another laboratory experiment by Caulliez, Ricci & Dupont 1998, transition of the initial laminar ﬂow to a turbulent ﬁeld was observed, which then leads to an explosive growth of the initial wind-generated ripples and the formation of surface streaks.) In addition to the surface streaks, local upwellings and divergence of the velocity ﬁeld at the surface were also observed by Melville et al. (1998), which caused bifurcation or dislocation of the streaks.

The acting of wind stress on the water surface induces both surface waves and turbulent shear ﬂow. It is possible that the coherent surface signatures observed in various experiments might result from several diﬀerent mechanisms, including the interaction among surface waves, the interaction between waves and shear current, and also the development of the turbulent boundary layer itself. None of these mechanisms has been demonstrated to the exclusion of the others. The formation of a wind-driven shear current always accompanies the inception of the surface waves. Accordingly, it is diﬃcult, if not impossible, to distinguish various mechanisms in the experiments, in particular those pertinent to surface waves. Attempts to separate the generation of surface waves from the wind-driven shear current by supplied

surfactants to inhibit the wave initiation (Faller & Perini 1984; Veron & Melville 2001), however, would also modify the sub-surface ﬂow structures and, consequently, aﬀect the induced surface signatures (McKenna & McGillis 2004; Tsai & Yue 1995; Tsai 1996).

In this study, the diﬃculty encountered in the experiments is resolved by numerical simulation of a canonical problem. Speciﬁcally, turbulent shear layer driven by a shear stress acting at a water-surface boundary is simulated numerically without the formation of the surface waves. Such a model problem has previously been studied for stress-driven oceanic mixed layers by large-eddy simulations (e.g. Skyllingstad & Denbo 1995; McWilliams, Sullivan & Moeng 1997; Sander et al. 2000), which require parameterization of the near-surface ﬂow using the law-of-the-wall theory. The approach adopted here is a direct numerical simulation of the wind-driven turbulent shear ﬂow using neither any turbulence closure to model the subgrid ﬂow processes nor a ‘vortex force’ to initiate the mechanism of wave–current interaction as in these large-eddy simulations. Therefore, with direct numerical simulation we can resolve the ﬂow scales down to that of the viscous sublayer and examine the detailed ﬂow ﬁeld immediately next to the surface. The present study thus attempts to deduce the possible correlation between the observed centimetre-scale surface signatures and the coherent ﬂow processes induced by shear turbulence. The simulation also focuses on providing further insight into the detailed hydrodynamic mechanism, particularly the coherent ﬂow structures that govern the process of scalar exchange across the air–sea interface.

The numerical aspects of the simulation, including the model formulation, the numerical method and the simulation implementation are described in § § 2, 3 and 4, respectively. Statistics of the turbulence are ﬁrst discussed in § 5. The numerical results of the ﬁrst-order mean proﬁles are compared with those derived from the surface renewal theory (Liu & Businger 1975) and the experimental measurements (Wu 1984). Turbulent energetics budgets of the present shear turbulence bounded by a stress-driven boundary, compared to those next to either a shear-free boundary or a no-slip wall, are discussed. Two-dimensional surface signatures of the velocity and scalar ﬁelds arising from the development of the turbulent boundary layer are then presented in § 6, and compared with those observed in the experiment of Melville et al. (1998). The underlying three-dimensional ﬂow structures, which induce the surface signatures, are ﬁnally examined and discussed in§ 7. Some of the preliminary results from this investigation were reported in letters (Tsai 2001; Tsai et al. 2003).

*2. Mathematical formulation*

*We consider a turbulent shear ﬂow beneath a ﬂat free surface at z = 0. The*
* coordinate axes x = (x, y, z) are in the streamwise, spanwise and upwards vertical*
direction, respectively, with the corresponding velocities

**v = (u, v, w). The turbulent***shear ﬂow is driven by a shear stress τs*acting at the boundary surface. The ﬂuid is assumed to be incompressible and Newtonian, such that the velocity

**v and the***dynamic pressure p are governed by the solenoidal condition and the Navier–Stokes*equations:

*(2.1)*

**∇ · v = 0,***∂*

**v***∂t*+

*1*

**v · ∇v = −***ρ*2

**∇p + ν∇**

_{v,}_{(2.2)}

With the focus of this study being on the boundary eﬀect of the stress-imposed
surface on the underlying shear turbulence and vice versa, an idealized boundary of
a ﬂat free surface is considered (e.g. Walker, Leighton & Garza-Rios 1996; Calmet
& Magnaudet 2003) in this study. By a ﬂat free surface, what is meant is that the
exerted stress is balanced by local surface tangential stress and that the surface is
free to ﬂow horizontally, but not vertically. A previous study (Tsai 2001) using more
realistic free-surface boundary conditions (Borue, Orszag & Staroselsky 1995; Tsai
1998) indicates only a minor impact of the small-amplitude surface motions on both
the surface and subsurface turbulent structures. At the surface, the tangential surface
stress is balanced by the mean exerting shear stress, whereas the normal surface stress
vanishes. Accordingly, the corresponding kinematic and dynamic conditions on the
*mean surface, z = 0, are*

*w= 0,* *∂u*
*∂z* =
*τs*
*ρν,*
*∂v*
*∂z* *= 0,*
*∂p*
*∂z* *= 0.* (2.3)

Note that in contrast to the Dirichlet pressure condition on a free surface with
deformation (e.g. Tsai 1998), the normal-stress condition on a ﬂat free surface results
*in a homogeneous Neumann pressure condition. With w = 0, the surface gravity waves*
are excluded from our simulations. Although the idealized free surface diﬀers from
that of the experiment of Melville et al. (1998) and also from the real ocean surface,
the exclusion of surface waves allows us to isolate the eﬀect of shear turbulence on
small-scale signatures at the air–sea interface without the complicating eﬀects of the
mean-ﬂow/surface-wave interaction. Our simulations also diﬀer from that of a
wall-bounded shear turbulent boundary layer (e.g. Mansour, Kim & Moin 1988; Moin
& Mahesh 1998). In contrast to shear turbulence over a no-slip wall, our simulated
ﬂows allow for tangential motions on the surface in response to the driven surface
stress and the underlying turbulence.

To study the eﬀects of the near-surface ﬂow processes on the transport of passive
scalars, advection–diﬀusion equations governing the evolutions of the temperature
* ﬁeld θ (x, t) and the gas concentration c(x, t) are also integrated in time with the ﬂow*
simulation. The boundary condition at the surface, however, depends on the transport

*process considered. A given uniform heat ﬂux Q gives rise to a Neumann condition*for the temperature ﬁeld at the surface as

*∂θ(x, y, 0, t)*

*∂z* =−

*Q*

*ρcpνθ*

*,* (2.4)

*where cp* *and νθ* are speciﬁc heat and thermal diﬀusivity, respectively. In our
simulations, the temperature ﬁeld is treated as a passive tracer; the buoyancy eﬀect
due to temperature ﬂuctuations and the imposed surface heat ﬂux does not modify
the vertical momentum equation in (2.2).

In the case of transfer of sparingly soluble gas, the resistance is dominated by the subsurface aqueous ﬂow and the atmosphere can be regarded as an inﬁnite reservoir with a constant concentration. This results in a Dirichlet condition for the concentration ﬁeld at the surface as

*c(x, y, 0, t) = cs.* (2.5)

*3. Numerical method*

The numerical method used in the present simulation resembles that in Tsai (1998) but with modiﬁcations to improve the near-surface resolution and to impose the

solvability condition for the pressure equation. A mixed scheme of the
pseudo-spectral method in the two horizontal directions and second-order ﬁnite diﬀerencing
in the vertical is adopted for approximation of the spatial-diﬀerential operators in the
equations and boundary conditions. A stretched grid system is employed such that
the discretization grids cluster when approaching the surface so as to allow proper
resolution of the viscous sublayer adjacent to the upper surface. The vertical domain
of physical coordinate *−h 6 z 6 0 is discretized according to the transformation*
function
*zk*=*−h*
1− *tanh aζk*
*tanh a*
*,* 0*6 ζk* *6 1,* (3.1)

*where the grids of the numerical coordinate ζ are uniformly distributed. The degree*
*of cluster is determined by the parameter a. A typical value used in the present*
*computations is a = 1.8417 (Gavrilakis 1992).*

A low-storage Runge–Kutta scheme (Spalart, Moser & Rogers 1991) is adopted for temporal integration of the Navier–Stokes equation (2.2) for the velocity ﬁeld and the transport equations for the scalar ﬁelds. The pressure ﬁeld is ﬁrst solved at each interval step in the Runge–Kutta time integration before updating towards the next time step for the new velocities. The Poisson equation for the pressure ﬁeld is obtained by taking the divergence of the temporal-discretized Navier–Stokes equation and applying the solenoidal constraint to the velocity ﬁeld at the new time interval. Since the pressure conditions on both the ﬂat free surface and the free-slip bottom are of the Neumann type, an additional solvability condition is required to make the solution unique when solving the pressure Poisson equation. In the present numerical implementation, a convenient value for pressure is assigned to a chosen grid point. For spectral representation of the pressure ﬁeld, this corresponds to explicitly specifying the constant mode. Such a solvability condition, however, would not aﬀect the solution as only the pressure gradients are required in integrating the velocity equations.

*4. Posing the simulation*

We pose the simulations to mimic the subsequent development of the ﬂow of water
which has been driven by a wind ﬁeld accelerating constantly from rest to a reference
*speed of ua* *in a time interval of t*0 as in the experiment of Melville et al. (1998).
Accordingly, numerical simulation of the ﬂow is initiated with a mean streamwise
velocity distribution (in the unit of cm s−1):

*U(z) = At*0
*(1 + 2η*2*)erfc(η)*− 2
π*1/2η*exp(−η
2_{)}
*,* (4.1)

*where η =* *− z/2(νt)1/2*_{. Equation (4.1) corresponds to the velocity proﬁle of a plane}
laminar shear layer driven by a surface stress, which induces a linear increase in the
*surface velocity with a constant acceleration A for a period of t*0. Assuming continuity
in the evolution of the air ﬂow, the shear stress acting at the surface,

*τs* *= ρν*
*dU*
*dz*
*z*=0
= 2
π*1/2ρA(νt)*
*1/2 _{,}*

_{(4.2)}

*therefore increases with time in the rate proportional to t1/2* _{from the start of}
the simulation. The surface velocity develops freely and the maximum turbulence
*production occurs near the water surface. The constant acceleration A depends on*

*the reference wind speed ua*. For the three reference wind speeds considered in the
*experiments of Melville et al. (1998), ua*= 3, 4 and 5 m s−1, the estimated acceleration
*constants are A = 0.27, 0.5 and 1 cm s*−2 with the corresponding acceleration periods

*t*0*= 45, 30 and 20 s. The characteristic depths Ls= 2(νt*0)*1/2* *and velocities us= At*0 of
the shear layers at the beginning of the simulation are 1.4, 1.1 and 0.9 cm, and 12.15,
15 and 20 cm s−1, respectively. These result in the corresponding Reynolds numbers
*of the initial driven shear ﬂows, Re = usLs/ν*= 1630, 1650 and 1800, for the three
*reference wind speeds ua*= 3, 4 and 5 m s−1, respectively. From (4.2), the initial surface
*friction velocities are determined as u*_{∗}*= (τs/ρ*)*1/2= 0.45, 0.56 and 0.71 cm s*−1, which
*result in the wall Reynolds numbers Re*_{∗}*= u*_{∗}*Ls/ν= 60.4, 61.9 and 63.6 for ua*= 3, 4
and 5 m s−1, respectively.

In the experiment of Melville et al. (1998), predominant streamwise-elongated
high-speed streaks (see also the present simulation results in§ 6) were observed. The streaks
reveal a nearly equal spacing, which varies with the wind speeds. Accordingly, the
length and width of the computational domain are chosen to be at least six times
*that of the observed characteristic steak spacing; the domain size in both x and*

*y* is therefore 45, 33.6 and 23 cm, respectively, for the three reference wind speeds

*ua*= 3, 4 and 5 m s−1. All the simulations are carried out with 1283 grid points. This
means that the lateral region between two streaks is resolved by at least twenty grids.
The corresponding non-dimensional streamwise and spanwise grid sizes in wall unit,

*x*+ * _{and y}*+

_{, are 15.89, 14.54 and 12.87 for u}*a*= 3, 4 and 5 m s−1, respectively. In
the vertical direction near the upper surface, the stretched grid distribution allows
suﬃcient resolution of the viscous layer. For the three simulations considered, there
*are at least ten grids in the near-surface region (z*+_{6 10).}

In the depth direction, the computational domain is closed by imposing a free-slip bottom condition at a depth of a quarter of the length (and the width) of the domain to emulate the inﬁnite depth. Thus, the depth of the turbulent layer will grow in time through entrainment in proportion to the applied wind stress. The thickness of the shear layer develops to about a quarter of the depth of the computational domain at the end of the simulation. For all the simulations conducted in the present study, the numerical depth is suﬃciently large not to aﬀect the near-surface ﬂow.

In order to generate a realistic turbulent velocity ﬁeld, a ﬂuctuating but solenoidal
vector ﬁeld is ﬁrst generated by taking the curl of another random vector ﬁeld.
The simulation is then conducted using the ﬂuctuating solenoidal vector ﬁeld as the
initial velocity condition without the mean streamwise velocity. As the integration
of the momentum equations proceeds, the simulated ﬂow ﬁeld adjusts itself from an
initial artiﬁcial random vector ﬁeld to a realistic turbulence. The ‘spin-up’ simulation
is performed until the simulated ﬂow reaches an equilibrium state in which the
inertial and dissipating ranges of the energy spectra converge. During this stage,
the turbulence ﬁeld near the surface become eﬀectively anisotropic, with a rapid
attenuation of the vertical ﬂuctuating velocities and an accompanying increase in
the horizontal turbulence intensities. The turbulent velocity ﬂuctuations, which are
*homogeneous in the horizontal plane ((x, y)-plane), is then added to the mean ﬂow*
(4.1) to form the initial velocity ﬁeld for the simulation.

The abrupt superposition of the turbulent velocity ﬁeld onto the mean streamwise
velocity would require additional spin-up time before meaningful ﬂow statistics and
structures can be obtained from the simulations. Evolutions of three averaged ﬂow
*properties, the surface-averaged mean velocity us*=*uz*= 0 and turbulent kinetic
en-ergy*u*2*+v*2*z*= 0, and the volume-averaged turbulent kinetic energy*u*2*+v*2*+w*2V,

*tUs/Ls*
0 100 2000
0 100 200
*t (s)*
f
*ug*
*z=0*
*/us*
f
KE
g*z=0*
*/us*
2
f
KE
g
*/us*
2
0 200 300
0 10
0.4
0.6
0.8
1.0
f
*ug*
*z=0*
*/us*
0.4
0.6
0.8
1.0
f
*ug*
*z=0*
*/us*
0.4
0.6
0.8
1.0
0
0.05
0.10
f
KE
g*z=0*
*/us*
2
f
KE
g
*/us*
2
0
0.05
0.10
f
KE
g*z=0*
*/us*
2
f
KE
g
*/us*
2
0.05
0.10
*(a)*
*(b)*
*(c)*
5
0 5 10
0 5 10 15 20
100

Figure 1. Evolutions of the surface-averaged streamwise velocity *uz*_{= 0}*/us* (solid curves),

the surface-averaged turbulent kinetic energy *u*2*+ v*2z= 0*/u*2*s* (dashed curves), and the

volume-averaged turbulent kinetic energy *u*2*+ v*2_{V}*/u*2

*s* (dash-dotted curves) for the wind

*speeds ua= (a) 5 m s*−1*, (b) 4 m s*−1*and (c) 3 m s*−1.

*are shown in ﬁgure 1, where the turbulent ﬂuctuating quantity q**= q− q, and ·*
and ·_{V} represent the averaging over the horizontal plane and the entire volume,
respectively. Immediately after the start of the simulation, both the surface mean
velocity and ﬂuctuation intensity decrease drastically showing a rapid adjustment of
the imposed initial ﬂow to turbulence. The subsequent evolution of the ﬂow after the
surface mean velocity and turbulent intensity reaching the minima can be divided into
*two stages. During the ﬁrst stage (non-dimensional time tus/Ls*≈ 30 to 120), the mean
surface velocity rises in a higher rate while the turbulent intensity evolves gradually.
*After the developing stage (tus/Ls* & 120), the ﬂow can be considered as that of fully

developed turbulence in which the mean surface velocity increases in a reduced rate
while the surface turbulent intensity approaches a constant level. To determine if the
simulated ﬂow has reached an equilibrium state and has become representative of
turbulence, the evolutions of the one-dimensional turbulent energy spectra at various
depths are examined. (The details are the same as those in Tsai 1998.) Once the
surface turbulent intensities reach the fully-developed stage, the distributions of the
inertial and dissipation ranges of the energy spectra also reach the steady states at
both the shallow and submerged depths, indicative of the equilibrium of turbulence.
The turbulence statistics, which are discussed in§ 5, are space and time averages (over
*the whole horizontal plane and the large-scale time tus/Ls*= 120∼ 200) of the ﬂow
ﬁelds in this fully developed stage.

*5. Flow statistics*

5.1. Mean proﬁles

For the temporal mean current beneath a wind-blown air–water interface, a two-layer velocity proﬁle, analogous to that of the turbulent wall layer, has been observed in the laboratory (Wu 1975, 1984; McLeish & Putland 1975) and ﬁeld (Churchill & Csanady 1983) experiments. The mean current near, but not immediately below, the interface is found to be characterized by a logarithmic velocity proﬁle:

*us− u(z)*
*u*_{∗} =
1
*κ* *ln z*
+* _{+ ψ}*+

_{=}1

*κ*ln

*z*

*z*0

*,*(5.1)

*where κ is the von K ´arm ´an constant, the friction velocity u*_{∗}*= (τs/ρ*)*1/2*, the wall
*coordinate z*+_{= zu}

∗*/ν, and the constant of integration ψ*+ which is related to surface
*roughness length z*0 *by z*0*= (ν/u*∗) exp(*−κψ*+*). The greater the constant ψ*+, the
*smaller the surface roughness length z*0. The universal function of (5.1) is also
appli-cable to the proﬁles of the mean temperature*θ(z) and dissolved gas concentration*
*c(z) as*
*θ(z) − θs*
*θ*_{∗} =
1
*κθ*
*ln z*+* _{+ ψ}*+

*θ*(5.2) and

*c(z) − cs*

*c*

_{∗}= 1

*κc*

*ln z*+

*+*

_{+ ψ}*c,*(5.3)

*where θs*=*θ(0), cs*=*c(0), θ*∗= *− Q/(ρcpu*∗*), c*∗= *− J/u*∗*, and Q and J are the*
mean heat and gas ﬂuxes, respectively.

Immediately below the interface, the vertical component of the turbulence is largely suppressed so that a thin viscous sublayer exists where molecular viscosity dominates the momentum transport. This viscous sublayer was found to be thinner than the analogous sublayer adjacent to a solid no-slip wall (Csanady 1984; Wu 1984). Adopting the surface renewal model that the viscous sublayer undergoes cyclic growth and disruption, Liu & Businger (1975) proposed an exponential proﬁle for the mean velocity distribution within the sublayer as

*us− u(z)*
*us− ub*
= 1− exp
−*z*
*δu*
*,* (5.4)

*where ub* *denotes the bulk mean velocity. δu* is a scaling depth equivalent to the
thickness of the viscous sublayer, and relates to the average residence duration of

*renewal ﬂuid elements tr*∗ through
*δu= (νtr*∗)
*1/2 _{= ν}us− ub*

*u*2 ∗

*.*(5.5)

In developing (5.4), Liu & Businger (1975) assumed that the probability of renewal of
any ﬂuid element is independent of the duration for which it has resided in contact
with the surface (Danckwerts 1951), and the probability distribution of the residence
*time tr* is described by an exponential function, exp(−t*r/tr*∗*)/tr*∗. Equation (5.4) can be
re-written in terms of the wall coordinate as

*us− u(z)*
*u*_{∗} *= ξu*
1− exp
−*z*+
*ξu*
*,* (5.6)

*where ξu= δuu*∗*/νis the sublayer parameter. The asymptotic limit of (5.6) as (z*+*/ξu*)→
*0 becomes the familiar linear proﬁle (us− u)/u*∗*= z*+.

For heat and gas transport across the conduction and diﬀusion sublayers, similar exponential proﬁles can be obtained for temperature and concentration distributions within the sublayers as

*θ(z) − θs*
*θ*_{∗} *= ξθ*
1− exp
*−P r* *z*+
*ξθ*
*,* (5.7)
*c(z) − cs*
*c*_{∗} *= ξc*
1− exp
*−Scz*+
*ξc*
*,* (5.8)

*where the Prandtl number P r = ν/νθ* *and the Schmidt number Sc = ν/νc, νθ* and

*νc* are the diﬀusivities of heat and dissolved gas, respectively, the non-dimensional
*parameters ξθ= δθu*_{∗}*/νθand ξc= δcu*_{∗}*/νc*, the thicknesses of the thermal and molecular
*sublayers δθ= (νθtθ*)*1/2* *and δc= (νctc*)*1/2*.

Comparisons between the theoretical and simulated vertical mean distributions of
the streamwise velocity, the temperature and the gas concentration are shown in
*ﬁgure 2 for the three wind speeds ua*= 5, 4 and 3 m s−1. As shown in ﬁgure 2, the
simulated mean velocity, temperature and concentration distributions are all well
represented by the two-layer proﬁles. Our simulated ﬂows reveal an exponential
*viscous sublayer above z*+ * _{≈ 10 and a logarithmic inertial layer between z}*+

_{≈ 20}and 100. For the velocity comparison, the available laboratory measurements by

*Wu (1984) at the wind speeds of 4.8, 4.0 and 3.0 m s*−1 are also plotted. The experimental measurements, the numerical simulations and the parameterized proﬁles are comparative for the whole sublayer range extending to the buﬀer layer.

*In determining the non-dimensional constants κ and ψ*+ _{in (5.1), least-squares}
*ﬁtting to the simulated mean streamwise velocities for the depths z*+_{= 20 to 100 is}
employed. The exponential and logarithmic proﬁles are then matched by equating
the right-hand sides and their ﬁrst derivatives of (5.1) and (5.4) at the matching
*depth z*+* _{= ζ}*+

_{(Liu, Katsaros & Businger 1979; Kraus & Businger 1994). This yields}

*two nonlinear equations for the sublayer parameter ξu*and the matching coordinate

*ζ*+_{. The same procedure is employed to reach the matched temperature and gas}
concentration proﬁles.

The corresponding parameters of the matched mean velocity proﬁles in ﬁgure 2
for the three wind speeds are listed in table 1. For logarithmic ﬁtting of the averaged
*velocity distributions, the typical computed values of κ are between 0.35 and 0.40,*
which are close to the universal K ´arm ´an constant, 0.4. The typical values of the
*computed ψ*+_{, however, are much less than 5.0 of a no-slip smooth wall, indicating}

0 10 15
100
101
102
*(a)* *(d)* *(g)*
*(c)* *( f )* *(i)*
*(b)* *(e)* *(h)*
*z+*
5 0 10 15
100
101
102
5 10 15
100
101
102
5
0 10 15
100
101
102
*z+*
5 0 10 15
100
101
102
5 10 15
100
101
102
5
0 10 15
100
101
102
*z+*
5 0 10 15
100
101
102
5 10 15
100
101
102
5

Figure 2.*(a–c) Mean proﬁles of the streamwise velocity (us−u(z))/u*_{∗}*, (d–f ) the temperature*

(*θ(z)−θs)/θ*_{∗}, (g–i) the gas concentration (c(z)−cs*)/c*_{∗}*for various wind speeds: ua*= 5 m s−1

*(a, d, g), 4 m s*−1 *(b, e, h), and 3 m s*−1 *(c, f, i). The computed results are represented by solid*

curves. The dash-dotted curves are the exponential proﬁles (5.6) of the viscous sublayer developed based on the surface renewal model, and the dash-dot-dotted curves are the linear

approximation. The dashed curves denote the matched logarithmic proﬁles (5.1).䊊, laboratory

*measurements by Wu (1984) at wind speeds of 4.8, 4.0 and 3.0 m s*−1.

*ua* *us* *u*∗ *z*0 *δu* *tr*∗
(m s−1) (cm s−1) (cm s−1) *κ* *ψ*+ _{(mm)} * _{ζ}*+

_{ξ}*u*(cm) (s) 5.0 13.0 0.74 0.39 1.1 0.089 20.9 9.8 0.13 1.8 4.0 10.7 0.61 0.34 1.5 0.092 23.3 10.5 0.17 2.9 3.0 8.2 0.47 0.36 1.9 0.098 26.2 11.2 0.24 5.6

Table 1. Computed parameters in the mean velocity proﬁles (5.1) and (5.4) for the three

wind speeds.

*the increase in equivalent surface roughness z*0. For the ﬂow beneath a wind-blown
sea surface, this is caused by direct energy supply to the near-surface turbulence
through inﬂuences such as breaking wavelets (e.g. Csanady 1984). In the present
case of stress-driven shear ﬂow, the increase infers the enhancement of near-surface

horizontal turbulence intensity by the presence of the upper surface. (This will be
discussed in*§ 5.2.) The increase in surface roughness z*0has also been noted in a recent
direct numerical simulation of an oceanic boundary layer driven by wave breaking
(Sullivan, McWilliams & Melville 2004). For the limiting case with no breaking
*surface forcing, their simulation also indicates a higher eﬀective roughness z*0 than
that of wall-bounded turbulence. Their simulations with breaking surface forcing
further reveal that wave breaking would eﬀectively increase the surface roughness.

*The computed non-dimensional viscous sublayer parameters ξu* are around 10,
smaller than that immediately next to a smooth no-slip wall (≈16), and increase
*slightly (within the range of O(10*−1)) with decreasing wind speeds. This is
quanti-tatively similar to that observed in the experiments (Csanady 1984; Wu 1984). The
*increase in ξu* with decreasing wind speed also implies the thickening of the viscous
*sublayer δu* and the growth of resistance against the vertical momentum transport.

5.2. Turbulence intensities and production

The vertical distributions of the turbulent velocity variance normalized by the surface
friction velocity, *u**i*

2* _{(z)/u}*2

∗, are shown in ﬁgure 3 for the three reference wind speeds

*ua*= 5, 4 and 3 m s−1. Away from the viscous layer, the magnitudes of the turbulent
velocity variances are close to those next to a no-slip solid wall (e.g. Kim, Moin
& Moser 1987). In particular, the upper bound of *w*2*/u*2

∗ is about 1; it does not exhibit a large increase like that in the large-eddy simulations of oceanic mixed layers (Skyllingstad & Denbo 1995; McWilliams et al. 1997), in which the Craik– Leibovich ‘vortex force’ attributed to the wave–current interaction mechanism (Craik & Leibovich 1976) is incorporated into the simulations.

The turbulence ﬁeld near the interface is eﬀectively anisotropic and two-dimensional,
with a rapid attenuation of the vertical velocity ﬂuctuation and an accompanying
increase in the streamwise component on approaching the interface. In comparison
with the streamwise and vertical components, the spanwise turbulence intensity shows
relatively mild variation. This is contrary to the observations of the turbulent shear
ﬂows bounded by shear-free surfaces (such as the free-surface jet ﬂow in Walker, Chen,
& Willmarth 1995; the free-surface wake ﬂow in Tsai 1998; and the open-channel
ﬂow in Calmet & Magnaudet 2003), in which the transfer of the turbulence energy
to the streamwise component is much more signiﬁcant than that to the spanwise one.
*Immediately beneath the upper surface and within the viscous sublayer (z*+_{. 10), both}
streamwise and spanwise turbulence intensities exhibit salient features, in particular for
*the case of ua*= 3 m s−1*(ﬁgure 3c). (The depths of the viscous sublayers, z*+*= ξ*+≈ 10,
are marked in the ﬁgures.)

For turbulence bounded by a shear-free surface, the normal turbulent ﬂuctuations
are forced to decay towards the boundary. This is a manifestation of the blocking eﬀect
brought about by the kinematic condition of vanished normal velocity at the surface.
The shear-free boundary conditions, nevertheless, give rise to much weaker constraints
on the tangential velocities and lead to the dissipation rate close to the surface actually
becoming smaller than in the bulk of the ﬂow (Teixeira & Belcher 2000). Accordingly,
the horizontal velocity ﬂuctuations increase upon approaching a shear-free surface
(e.g. Perot & Moin 1995; Walker et al. 1996; Tsai 1998; Calmet & Magnaudet 2003).
On the present stress-imposed surface, the constraints of the horizontal ﬂuctuating
*velocities, ∂u**/∂z= ∂v**/∂z*= 0, are identical to those of a shear-free surface although
*the mean velocity ﬁeld is subjected to a constant vertical gradient, ∂u/∂z = τs/µ*.
As ﬁgure 3 shows, the intensities of the horizontal turbulence are also ampliﬁed as
they approach a stress-imposed surface. This is also because of the decrease in the

0 5 10
0
50
100
*(a)*
*(b)*
*(c)*
*z*+
0 5 10
0
50
100
*z*+
0 5 10
0
50
100
*z*+
— — — —
— —
— — — —
— —

Figure 3.Vertical distributions of the normalized turbulent velocity variances, *u*2*(z)/u*2_{∗}

(solid curves),*v*2*(z)/u*2

∗ (dash-dotted curves), *w*2*(z)/u*2∗ (dashed curves) and the averaged

turbulent kinetic energy,*u*2*+ v*2*+ w*2*(z)/u*2

∗ (thick solid curves) for various reference wind

*speeds (a) ua*= 5 m s−1*, (b) 4 m s*−1 *and (c) 3 m s*−1. The depths of the viscous sublayers,

*z*+*= ξ*+_{≈ 10, are marked by dashes in the ﬁgures.}

near-surface dissipation proﬁles, which will be discussed later. Such an enhancement of the horizontal turbulence intensities has also been observed in the recent simulation by Sullivan et al. (2004).

The drastic changes within the viscous sublayer are also discernible in vertical
variations of the vorticity variances,*ω*2

*x(z), ω*2*y(z) and ωz*2*(z), as well as the vertical*
strain rate (horizontal divergence), *(∂w/∂z)*2* _{ = (∂u/∂x + ∂v/∂y)}*2

_{, as shown in}ﬁgure 4. Both the vertical and horizontal vorticity variances, in particular the spanwise component, reach the maxima within the viscous sublayer. At the interface, the streamwise and spanwise vorticities decrease to negligible levels, whereas the vertical

*component remains ﬁnite. Away from the interface (z*+

_{& 40; the characteristic depth}

*of the shear layer is approximately at z*+

_{≈ 70), the horizontal vortical ﬁeld becomes}essentially isotropic and the vorticity variances continue reducing. No signiﬁcant

0 0.5 1.0
0
50
100
*(a)* *(d )*
*(b)* *(e)*
*(c)* *( f )*
*z*+
—
—
—
—
—
—
—
——
—
—
—
0.1 0.2
0
50
100
0 0.5 1.0
0
50
100
*z*+
0.1 0.2
0
50
100
0 0.5 1.0
0
50
100
*z*+
0.1 0.2
0
50
100

Figure 4. *Vertical distributions of the vorticity variances (a–c): streamwise vorticity,ω _{x}*2

*(z)*

(solid curve), ﬂuctuating spanwise vorticity,*ω*2_{y}(z) (dash-dotted curve), and vertical vorticity,

*ω*2

*z(z) (dashed curve); and the vertical velocity–strain variance, (∂w/∂z)*2*(z) (d–f ) for the*

*wind speeds of ua*= 5 m s−1*(a, d), 4 m s*−1*(b, e) and 3 m s*−1*(c, f ). All the terms are normalized*

*by us/Ls. The depths of the viscous sublayers, z*+* _{= ξ}*+

_{≈ 10, are marked by dashes.}

change in the streamwise vorticity variance is observed in these submerged depths for all the wind speeds considered throughout the simulations, implying the absence of predominant streamwise vortices.

To further elucidate the impact of the sheared surface on the transport of
near-surface turbulence energy, the energetics budgets are examined next. Assuming
equilibrium, the mean budget equations for turbulent ﬂow with two-dimensional
*mean shear ∂u(z)/∂z are*

*∂u*2
*∂t* =*−2u*
_{w}_{}*∂u*
*∂z*
*P*11
−*∂u*2*w*
*∂z*
*T*11
+2
*ρ*
*p**∂u*
*∂x*
*Π11*
*+ν∂*
2* _{u}*2

_{}

*∂z*2

*D*11

*−2ν*

*∂u*

*∂xk*

*∂u*

*∂xk*

*11*

*,*(5.9)

–0.001 0 0.001
0
20
40
*(a)*
*(b)*
*(c)*
*z*+
–0.0005 0 0.0005
–0.0002 0 0.0002
0
20
40
*z*+
0
20
40
*z*+

Figure 5.*Vertical distributions of the various production/consumption terms of (a) *

*stream-wise, (b) spanwise and (c) vertical turbulent intensities in the Reynolds-stress equations (5.9),*

*(5.10) and (5.11) for the reference wind speed ua*= 5 m s−1*: P*11*, thick solid curves; Tii, *

*dash-dot-dotted curves; Πii, solid curves; Dii, dash-dash-dot-dotted curves; ii, dashed curves. All the terms are*

*normalized by u*3
*s/Ls.*
*∂v*2
*∂t* =−
*∂v*2*w*
*∂z*
*T*22
+2
*ρ*
*p**∂v*
*∂y*
*Π*22
*+ν∂*
2* _{v}*2

_{}

*∂y*2

*D*22

*−2ν*

*∂v*

*∂xk*

*∂v*

*∂xk*

*22*

*,*(5.10)

*∂w*2

*∂t*=−

*∂w*3

*∂z*

*T33*−2

*ρ*

*∂p*

*w*

*∂z*+ 2

*ρ*

*p*

*∂w*

*∂z*

*Π33*

*+ν∂*2

*2*

_{w}_{}

*∂z*2

*D33*

*−2ν*

*∂w*

*∂xk*

*∂w*

*∂xk*

*33*

*,*(5.11)

*where P*11 *is the shear production rate, Tii* *the turbulent advection rate, Πii* the
*velocity–pressure correlation term, Diithe viscous diﬀusion rate, and ii*the dissipation
rate. Figure 5 shows the vertical variations of the various terms in (5.9), (5.10) and
*(5.11) for the reference wind speed ua*= 5 m s−1. The distributions for other wind
speeds resemble those in ﬁgure 5.

The distributions of the energy budgets shown in ﬁgure 5 reveal some general
features known for shear turbulence bounded by a no-slip wall (e.g. Mansour
*et al. 1988): the shear production rate, the P*11term, is the major turbulence production
for the streamwise budget, which reaches its maximum level within the viscous
*sublayer before diminishing at the surface. The pressure–strain term, Π*11, converts
*the u*2 *component to v*2 *and w*2 *and hence is a consumption for u*2 and a major
*source for the other two components, except near the viscous sublayer for w*2. The
*dissipation terms, ii*, account for the major loss for all three energy components. The
*contributions of turbulent advection, Tii* *terms, and viscous diﬀusion, Dii* terms, are
in general negligibly small except within the viscous sublayer. These general features
of the turbulent kinetic energy budget in the region of submerged water have also
been observed in the simulation of a stress-driven oceanic boundary layer by Sullivan
et al. (2004).

The main diﬀerences in the energetics budgets between the present
stress-driven turbulence and the wall-bounded turbulence are within the viscous sublayer
attributed to the diﬀerent boundary conditions. For a free surface (sheared or
*shear-free), the vertical ﬂuctuating strain rate, ∂w**/∂z*= *− (∂u**/∂x+ ∂v**/∂y*), increases
rapidly approaching the boundary (ﬁgure 4). Concurring with the ampliﬁcation of
near-surface streamwise turbulent velocities (ﬁgure 3), the convective transport of
*streamwise turbulence, T*11, increases from consumption in the submerged water to
production within the sublayer, and reaches its maximum at the upper surface.
The near-surface streamwise turbulent transport is of a magnitude comparable to
*that of shear production P*11 in the submerged region. Such a signiﬁcant increase
*in near-surface T*11 *is balanced by pressure transport (Π*11) and an enhancement in
*viscous diﬀusion (D*11). The present near-surface turbulent energy budget is somewhat
diﬀerent from that in a similar simulation by Sullivan et al. (2004), but with a higher
Reynolds number, in which negligible contributions from turbulent transport and
viscous diﬀusion have been observed.

The variations in the dissipation rates within the sublayer next to a stress-imposed
surface are completely opposite to those next to a wall boundary. The dissipation
*rates of both the streamwise and spanwise turbulent energies, *11 *and *22, decrease
drastically when approaching the upper surface. This results in a corresponding
increase in the horizontal ﬂuctuating velocities as explained by the theoretical model
of Teixeira & Belcher (2000). The dissipation rate of the vertical turbulent energy,

33, increases rapidly and reaches its maximum level at the surface, indicative of an abrupt reduction in vertical velocity towards the surface.

For shear turbulence near a no-slip wall, viscous diﬀusion dominates the productive
*transport of u*2 *and v*2*, and is negligibly small for vertical component w*2. However,
for turbulence bounded by a sheared surface, both the streamwise and spanwise
*diﬀusion rates, D*11 *and D*22, become the primary consumption of the horizontal
components within the viscous sublayer and reach their negative minima at the
*surface. On the other hand, viscous diﬀusion, D*33, is the major term for vertical
transport immediately beneath the surface. This conﬁrms that the vertical exchange
across the interface is controlled by viscous diﬀusion within the viscous sublayer.

*6. Surface characteristics*

6.1. Streaky appearance

The ﬁrst clue of the existence of organized structures within the stress-driven surface turbulent boundary layer is the emergence of elongated ﬁlaments of high-speed streaks

with somewhat equal transverse spacing on the water surface as have been observed
in both ﬁeld and laboratory experiments. In the experiment of Melville et al. (1998)
the typical spanwise spacing between the streaks ranges from∼3 cm to ∼7 cm as the
reference wind speed decreases from 5 m s−1 to 3 m s−1. These streamwise elongated
*streaks with transverse spacing of O(1 cm) to O(10 cm) have also been observed*
in the ocean (Gemmrich & Hasse 1992) and lake (Kenney 1993; Woodcock 1941).
Such a surface signature arises from the wind-driven shear layer and is geometrically
similar to the low-speed streaks near a stationary wall, which develop within the
viscous sublayer and extend into the logarithmic region (e.g. Kline et al. 1967; Smith
& Metzler 1983).

In the present simulation of surface-driven turbulent shear ﬂow, high-speed surface streaming, similar to that observed in the laboratory and ﬁeld experiments, emerges on the surface boundary. Instantaneous distributions of the surface streamwise velocity and the temperature contours for the three wind speeds are shown in ﬁgure 6. The distribution patterns clearly demonstrate the existence of elongated surface streaks carrying cooler and faster-moving ﬂuids. The major streaks are arranged with somewhat equal transverse spacing, and the spacing increases with decreasing wind speed.

To further visualize the structures and also the inception and evolution of the surface
streaming, 1282 _{uniformly distributed ﬂoating Lagrangian particles are released on}
the water surface and their trajectories are tracked. Figure 7 depicts the evolution
*of the Lagrangian particles for the wind speed of ua*= 5 m s−1. Immediately after
the particles are released, they aggregate and form a ‘ﬁsh-skin’ pattern as shown in
*ﬁgure 7(a). The merged particles further accumulate into several elongated groups*
*and align along the fast-moving velocity streaks (ﬁgures 7b and 7c). The rapidity of*
the formation of this surface streaky pattern is the same as the observation of Melville
et al. (1998) and as that described by Woodcock (1941): ‘When the powder is dusted
over the water it moves into lines so quickly that, even in moderate breeze, the eye
can hardly detect the transition from random scattering to the linear pattern.’ It takes
less than two seconds for the particles to evolve from the initial uniform distribution
*to the streaky pattern in ﬁgure 7(c). The aligned particles continue travelling along*
the high-speed streaky regions although some particles of the individual streets may
diverge and join other ones.

As time proceeds, some segments of the particle streets break, and the particles
*further cluster in a streamwise fashion and form shorter streets, as shown in ﬁgures 7(e)*
*and 7(f ). These shorter streets, however, still travel along the same velocity streaks. To*
conﬁrm that the elongated velocity streaks persist though the streets of Lagrangian
particles seem to have dislocated, a new set of uniformly distributed particles are
*released on the surface immediately after the time of ﬁgure 7(d) (5.364 s from the*
release of the ﬁrst set of particles). The distributions of the newly released particles at
the corresponding escaping times in ﬁgure 7 are shown in ﬁgure 8. Similar evolution
patterns are observed for the two sets of surface ﬂoating particles released at diﬀerent
times.

Accompanying the more organized fast-moving cooler streaks are slowly-moving
warmer localized spots appearing intermittently at the surface, as shown in ﬁgure 6.
The appearance of these random surface spots can be educed from the initial
*formation of the ‘ﬁsh-skin’ pattern in ﬁgure 7(a). At the early stage after the particles*
are released, the ﬂoating particles bypass the localized regions of lower streamwise
velocities and begin to accumulate around the upstream edges of the slowly moving
*spots. The process is focused on by showing in ﬁgure 9(a) the distribution of the*

10
*(a)*
5
0
–5
–10
–10 –5 0 5 10
10
5
0
–5
–10
–10 –5 0 5 10
*y (cm)*
10
*(b)*
5
0
–5
–10
–10 –5 0 5
10
5
0
–5
–10
–10 –5 0 5
*y (cm)*
10
*(c)*
*(d)*
*(e)*
*( f )*
5
0
–5
–10
–10 –5 0 5
10
5
0
–5
–10
–10 –5 0 5
*y (cm)*
*x (cm)* *x (cm)*

Figure 6. *Representative contour distributions of the streamwise velocity (a–c) and *

*tempera-ture (d–f ) at the free surface for the three reference wind speeds: (a, b) 5 m s*−1*; (c, d) 4 m s*−1;

*(e, f ) 3 m s*−1. The ﬂow travels from left to right. The warm/cold-colour areas represent regions

*(a)* *(b)*
10
5
0
–5
–5 0 5 10
–10
10
5
0
–5
–10
–10
*(c)*
10
5
0
–5
–5 0 5 10
–10
–10
*(d)*
10
5
0
–5
–5 0 5 10
–10
–10
*( f )*
10
5
0
–5
–5 0 5 10
–10
–10
*(e)*
10
5
0
–5
–5 0 5 10
–10
–10
–5 0 5 10
–10

Figure 7.Instantaneous distribution of the surface ﬂoating Lagrangian particles at times

*(a) 0.224 s, (b) 0.671 s, (c) 0.894 s, (d) 1.788 s (e) 3.129 s and (f ) 5.364 s after the release of the*

particles. 1282_{uniformly distributed particles are released.}

Lagrangian particles superposed on the contours of the streamwise velocity shortly
after the release of the particles. A high-speed streak bifurcates into branches when it
encounters a localized low-speed spot. This can be seen from the streamwise velocity
*contours (ﬁgure 6) as well as the distribution of Lagrangian particles (ﬁgure 7c, d).*
A close-up superposed distribution of streamwise velocity and particles revealing
*the bifurcation events is shown in ﬁgure 9(b). Such a feature is identical to the*
streak dislocation or bifurcation observed by Melville et al. (1998). The merging and
bifurcation of Lagrangian particles also look similar to those shown in McWilliams
et al. (1997, ﬁgure 20), which they claimed are due to Langmuir cells.

*(a)* *(b)*
10
5
0
–5
–5 0 5 10
–10
10
5
0
–5
–10
–10
*(c)*
10
5
0
–5
–5 0 5 10
–10
–10
*(d)*
10
5
0
–5
–5 0 5 10
–10
–10
*( f )*
10
5
0
–5
–5 0 5 10
–10
–10
*(e)*
10
5
0
–5
–5 0 5 10
–10
–10
–5 0 5 10
–10

Figure 8. Instantaneous distribution of the surface ﬂoating Lagrangian particles at times

*(a) 0.224 s, (b) 0.671 s, (c) 0.894 s, (d) 1.788 s (e) 3.129 s and (f ) 5.364 s after the release of the*
second set of uniformly distributed particles. The particles are released immediately after the
*time of ﬁgure 7(f ).*

6.2. Spanwise distribution of streaks

For wall-bounded turbulent ﬂows, low-speed streaks similar in geometry to the present
high-speed counterparts on a stress-imposed surface, have been identiﬁed within the
viscous sublayer in close proximity to the no-slip wall. The non-dimensional mean
*spanwise spacing, λ*+_{= λu}

∗*/ν*, between these low-speed streaks next to a no-slip wall
is generally accepted to be around the value of 100 (e.g. Kline et al. 1967), where

*λ is the dimensional mean spanwise spacing. Further study by Smith & Metzler*

(1983) indicated that such a non-dimensional mean value is essentially invariant with Reynolds number, and the streak spacing exhibits a log-normal distribution as suggested previously by Nakagawa & Nezu (1981).

*(a)* *(b)*
0
–5
–5
0
–10 –5 0

Figure 9.Close-up snapshots of the Lagrangian particles superposed on the streamwise

*velocity contours. Immediately after the particles are released (a), the ﬂoating particles bypass*
the localized regions of lower streamwise velocities (blue areas). The particles eventually
*accumulate and travel along the fast-moving streaks (red areas) as shown in (b), where the*
fast-moving streak bifurcates into branches when it encounters a localized slowly moving spot.

*-y+*
10 20 30 40 50 60 70 80 90 100
0
500
1000
1500
Frame number

Figure 10.*Sequence of the streak marks for the reference wind speed ua*= 5 m s−1. Each

vertical line corresponds to one image frame for streak identiﬁcation. The time interval between two subsequent frames is about 0.0224 s. The dashes indicate the transverse locations of the identiﬁed streaks.

Following the study of Smith & Metzler (1983), we determine the streak-spacing
values from visual identiﬁcation using the image sequence of the contours of
streamwise velocity and the distributions of Lagrangian particles, as presented in
ﬁgures 6 and 7. The images of particle distributions are superposed on the velocity
contours to allow a cross-examination when the streak positions are determined. For
each frame of the images, a reference spanwise axis at a ﬁxed streamwise coordinate
is chosen where the streak positions are to be marked. A streak intersecting the
reference spanwise axis is identiﬁed when the Lagrangian particles aggregate along an
elongated well-deﬁned region of high streamwise velocity with a streamwise extent of
*at least x*+_{= 300. The locations of the streak intersections along the transverse axis}
are then marked using a graphic digitizer software, and the transverse coordinates
are stored. An example of the sequence of streak marks for the reference wind speed

*P(λ*+_{)}
*P(λ*+_{)}
*P(λ*+_{)}
200 400 600
0
0.002
0.004
0.006
0.008
*(a)*
*(b)*
*(c)*
*λ*+_{ = 183}*σ*+
*λ *= 90
200 400 600
0
0.002
0.004
0.006
0.008
*λ*+
0 200 400 600
0.002
0.004
0.006
*λ*+_{ = 161}*σ*+
*λ *= 63
*λ*+_{ = 201}*σ*+
*λ *= 77

Figure 11. Probability density histograms of the spanwise streak spacing and the log–normal

*distributions with the corresponding mean value λ*+ * _{and standard deviation σ}*+

*λ* for the wind

*speed (a) ua*= 5 m s−1*, (b) 4 m s*−1*and (c) 3 m s*−1.

frame for streak identiﬁcation, and the dashes indicate the transverse locations of the
*identiﬁed streaks. The transverse distances between the neighbouring streaks, λ, are*
then measured from the sequence of streak marks.

Distribution histograms of the spanwise spacing between two streamwise streaks
at the water surface determined from the described procedure of streak identiﬁcation
and counting are shown in ﬁgure 11 for the three reference wind speeds considered.
*The streak-spacing values are presented in non-dimensional form λ*+_{= λu}

∗*/ν*. The
*three distributions all skew toward values lower than the means, λ*+_{, similar to that}
of low-speed streaks within a wall boundary layer. The determined non-dimensional
*values of the mean spacing λ*+ _{are 183, 161 and 201, with the corresponding standard}
*deviations σ*+

*λ* *= 90, 63 and 77, respectively, for the three wind speeds ua*= 5, 4 and
3 m s−1*, as plotted in ﬁgure 12(a). These results are about twice the values of λ*+_{≈ 100}
*and σ*+

*λ* ≈ 40 for the mean spacing of low-speed streaks near a no-slip wall (Smith &
Metzler 1983).

*ua *(cm s–1)
*λ (cm)*
2
4
6
8
10
3 4 5
3 4 5
0
100
200
300
*(a)*
*(b)*
*λ*+

Figure 12.*(a) Variation of the non-dimensional mean streak spacing, λ*+, determined by the

*visual procedure with the wind speed. (b) Comparison of the dimensional streak spacings*

determined by both the visual procedure () and the spectral analysis () in the present

simulations with those from the experiments of Melville et al. (1998) () for the three wind

*speeds. The corresponding spacing of the second peaks for ua*= 5 m s−1 from the spectral

analysis is beyond the upper range and is not shown in the ﬁgure.

As suggested by Nakagawa & Nezu (1981) and also Smith & Metzler (1983), the
skewed distribution of the streak spacing can be characterized statistically using a
*log–normal probability density function determined from the values of λ*+ * _{and σ}*+

*λ* in
the form:
*P(λ*+) = √ 1
2π*λ*+* _{ψ}*
0
exp
−1
2
1

*ψ*0 ln

*λ*+

*λ*+

_{0}2

*,*(6.1)

*where λ*+

_{0}

*= λ*+

*2*

_{(1 + ψ}*λ*)*−1/2* *is the median value of λ*+*, ψ*0=

*ln(1 + ψ*2
*λ*)

*1/2*
is the
*coeﬃcient of variation of ln λ*+_{, and ψ}

*λ= σλ*+*/λ*+. The ﬁtted log–normal probability
density functions are superposed over the corresponding histograms in ﬁgure 11.
The comparison between the histograms and the log–normal distribution functions
appear to be quite good. These results also suggest the possible similarity between
the small-scale high-speed streaks on a wind-blown water surface and the low-speed
streaks next to a wall boundary.

The characteristics of the surface streak spacing can also be estimated by examining
the spanwise spectrum of the streamwise-averaged surface ﬂuctuating streamwise
*velocity, u**s(y) =*

*u**(x, y, 0) dx. The distributions of the normalized average spectra,*
*u**s(λ)/u*2, for the three reference wind speeds are shown in ﬁgure 13. The power
spectra are obtained from the ensemble average of the streamwise-averaged velocity

*u**s(y) at various times within the interval when the surface turbulence reaches the*
stationary state. Despite the initial broadband spectra without any predominant
span-wise wavelengths, the streamspan-wise ﬂuctuating velocities are dominated by components

200 1000 1800
0.01
0.03
0.05
*(a)*
*(b)*
*(c)*
200 1000 1800
0.01
0.03
0.05
200 1000 1800
0.01
0.03
*λ+*

Figure 13. Distributions of the normalized average spectra, _{u}* _{s}(λ)/u*2, for the three

*reference wind speeds of (a) 5 m s*−1*, (b) 4 m s*−1 *and (c) 3 m s*−1. The power spectra are

*obtained from ensemble averages of the streamwise-averaged velocity u*_{s}(y) at various times

within the intervals when the surface turbulence reaches the stationary state.

with preferential cross-spacing wavelengths when surface streaks appear, as shown in
ﬁgure 13. The spectra all display two characteristic peaks for the three wind speeds
*considered. The non-dimensional average streak spacing λ*+_{, which can be considered}
as the shorter wavelength of the ﬁrst spectrum peak, are 219, 249 and 270, respectively,
*for the wind speeds ua*= 5, 4 and 3 ms−1. Consistent with the ﬁnding for low-speed
streaks next to a no-slip wall (Kline et al. 1967), the spacing values determined in
this more objective manner are greater than those determined by the visual counting
scheme. The longer wavelength implies possible spanwise modulation of the streaks
as shown in ﬁgure 6.

While the present numerical simulations are not completely compatible with
the conditions in the experiments of Melville et al. (1998), the simulated surface
signatures, however, quantitatively resemble the observations. It is therefore worth
comparing the characteristic cross-spacing of the surface streaks arising from
the simulated boundary layer with that from a wind-driven surface shear layer.
Summary of the surface streak spacing determined by both the visual procedure
and the spectral analysis in the present study, and from the experiment of Melville
*et al. (1998) is shown in ﬁgure 12(b). The spacing values observed in the experiment*

of Melville et al. (1998) are in general slightly larger than those determined in the
*present numerical simulations. For the wind speed of ua*= 5 m s−1, the numerical
predictions are close to the measurement of Melville et al. (1998). Both numerical
and experimental results clearly show the same trend: the streak spacing decreases as
the wind speed increases.

*7. Coherent ﬂow structures*

The surface features discussed in§ 6 are strongly linked to the coherent structures
in the subsurface ﬂow. To reveal the underlying structures, three-dimensional
*cut-away views of instantaneous streamwise velocity u and temperature θ distributions*
*are shown in ﬁgures 14 and 15 for the wind speed ua*= 5 m s−1. The corresponding
*distributions of vertical velocity w and gas concentration c on the facing cross-stream*
planes are also depicted in ﬁgures 14 and 15, respectively.

It is discernible that the longitudinal streaks are the surface signatures of an array of high-speed jets strongly contiguous to the surface. The vertical depths of the jets are limited within the viscous sublayer. Looking into the distribution of the vertical velocity, downdrafts occur beneath the longitudinal streaks (marked with downward arrows in the ﬁgure). The cold water collects along the high-speed jets which are characterized by cold streaks at the surface. The well-known cool-skin eﬀect (e.g. Katsaros 1980), therefore, is attributed to these coherent thermal structures within the surface sublayer. Despite the diﬀerent surface boundary conditions between the temperature and the gas concentration ﬁelds, the two scalar distributions are virtually identical except in the region immediately next to the water surface.

Beneath the sublayer, distinct ‘tongues’ of warmer water with lower streamwise velocities stick vertically out of the well-mixed submerged region. In contrast to the more organized cold high-speed jets, these warm low-speed tongues appear randomly with no spatial pattern. Comparing the distributions of the streamwise and vertical velocities on the cross-stream vertical plane, it can be seen that the coherent tongues are formed by upward advection of the warm and low-speed water underneath. When the upwelling is strong enough to penetrate the viscous sublayer (one particular event is highlighted by a circle in ﬁgure 14), it results in a localized low-speed spot at the surface. The bifurcation of the elongated surface streak occurs (ﬁgure 9) when a high-speed jet encounters such a localized low-speed updraft.

*An explanation of this local upwelling ﬂow is the emergence of Ω-shaped horseshoe*
vortices in the turbulent shear layer by turning and stretching of the spanwise vortices
(Tsai 1998). These coherent vortices move upwards and impinge on the surface. The
*upper heads of the Ω-shaped vortices thereby induce the local upwelling ﬂows and*
bring up the submerged slowly moving ﬂuids. Such an impinging process can also
be deduced from the enhancement of horizontal divergence of the velocity ﬁeld
immediately beneath the surface (ﬁgure 4).

To assess the presence of circulatory ﬂow that forms the surface streaks, we plot
*in ﬁgure 16 the instantaneous iso-surfaces of streamwise vorticity ωx*, the contours of

*ωx* *on various vertical cross-sections (y, z-plane) along with the streamwise velocity*
contour on the surface. In sharp contrast to the more organized high-speed streaks
at the surface, the streamwise vortical ﬁeld display is irregularly spaced and is much
shorter in length. The spanwise spacing between the vortex tubes is also narrower
than that between the streaks. Furthermore, the predominant streamwise vortices are
strongly conﬁned to near the surface, highly consistent with the vertical distribution
of streamwise vorticity (ﬁgure 4). While the pairing of the counter-rotating vortices

–15
–10
–5
0
5
10
15 15
10
*y/L*
*x/L*
5
0
–5
–10
–15
0
20
*40 z*+
60

Figure 14. *Three-dimensional cut-away view of streamwise velocity u for the reference wind*

*speed ua*= 5 m s−1*. The corresponding contours of vertical velocity w on the facing cross-stream*

plane are also depicted on an additional vertical section. Note that the non-dimensional friction

*coordinate z*+ _{is used for the vertical axis to demonstrate better the vertical scale of the ﬂow}

structures. The downward arrow indicates the occurrence of downdraft beneath the high-speed streaks. The circle highlights a particular region where the upwelling tongue penetrates the viscous sublayer, and the upward advection carries the submerged slowly-moving water and induces a low-speed spot on the free surface.

–15
–10
–5
0
5
10
15 15
10
*y/L*
*x/L* _{5} 0
–5
–10
–15
0
20
*40 z*+
60

Figure 15. *Three-dimensional cut-away view of temperature θ for the reference wind speed*

*ua*= 5 m s−1*. The corresponding contours of gas concentration c on the facing cross-stream*

plane are also depicted on an additional vertical section. Note that the non-dimensional friction

*coordinate z*+ _{is used for the vertical axis to demonstrate better the vertical scale of the ﬂow}

structures. The circle highlights the region where the upwelling penetrates the conductive sublayer and induces a warm spot on the free surface.