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 N2
A N D C H I N -H O H M O E N G3
1Institute of Hydrological Sciences, National Central University, Taoyuan 32001, Taiwan 2Department of Civil Engineering, National Chiao Tung University, Hsinchu 30050, Taiwan
3National 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.
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: ∇ · v = 0, (2.1) ∂v ∂t +v · ∇v = − 1 ρ∇p + ν∇ 2v, (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)
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 , 06 ζ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 t0 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) = At0 (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 t0. 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
t0= 45, 30 and 20 s. The characteristic depths Ls= 2(νt0)1/2 and velocities us= At0 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-ergyu2+v2z= 0, and the volume-averaged turbulent kinetic energyu2+v2+w2V,
tUs/Ls 0 100 2000 0 100 200 t (s) f ug z=0 /us f KE gz=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 gz=0 /us 2 f KE g /us 2 0 0.05 0.10 f KE gz=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 u2+ v2z= 0/u2s (dashed curves), and the
volume-averaged turbulent kinetic energy u2+ v2V/u2
s (dash-dotted curves) for the wind
speeds ua= (a) 5 m s−1, (b) 4 m s−1and (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 z0 , (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 z0 by z0= (ν/u∗) exp(−κψ+). The greater the constant ψ+, the smaller the surface roughness length z0. 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 u2 ∗ . (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(−tr/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∗ z0 δ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
the increase in equivalent surface roughness z0. 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 z0has 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 z0 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, ui
∗, 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 w2/u2
∗ 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, u2(z)/u2∗
∗ (dash-dotted curves), w2(z)/u2∗ (dashed curves) and the averaged
turbulent kinetic energy,u2+ v2+ w2(z)/u2
∗ (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), ω2y(z) and ωz2(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,ωx2(z)
(solid curve), ﬂuctuating spanwise vorticity,ω2y(z) (dash-dotted curve), and vertical vorticity,
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
∂u2 ∂t =−2u w∂u ∂z P11 −∂u2w ∂z T11 +2 ρ p∂u ∂x Π11 +ν∂ 2u2 ∂z2 D11 −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: P11, 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 u3 s/Ls. ∂v2 ∂t =− ∂v2w ∂z T22 +2 ρ p∂v ∂y Π22 +ν∂ 2v2 ∂y2 D22 −2ν ∂v ∂xk ∂v ∂xk 22 , (5.10) ∂w2 ∂t =− ∂w3 ∂z T33 −2 ρ ∂pw ∂z + 2 ρ p∂w ∂z Π33 +ν∂ 2w2 ∂z2 D33 −2ν ∂w ∂xk ∂w ∂xk 33 , (5.11)
where P11 is the shear production rate, Tii the turbulent advection rate, Πii the velocity–pressure correlation term, Diithe viscous diﬀusion rate, and iithe 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 P11term, 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 u2 component to v2 and w2 and hence is a consumption for u2 and a major source for the other two components, except near the viscous sublayer for w2. 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, T11, 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 P11 in the submerged region. Such a signiﬁcant increase in near-surface T11 is balanced by pressure transport (Π11) and an enhancement in viscous diﬀusion (D11). 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 u2 and v2, and is negligibly small for vertical component w2. However, for turbulence bounded by a sheared surface, both the streamwise and spanwise diﬀusion rates, D11 and D22, 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, D33, 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. 1282uniformly 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−1and (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 = λ+(1 + ψ2
λ)−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, us(y) =
u(x, y, 0) dx. The distributions of the normalized average spectra, us(λ)/u2, 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
us(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, us(λ)/u2, 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 us(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.