Wave Boundary Layer Turbulence over Surface Waves in a Strongly Forced Condition

Accurate predictions of the sea state–dependent air–sea momentum ﬂux require a thorough understanding of the wave boundary layer turbulence over surface waves. A set of momentum and energy equations is derived to formulate and analyze wave boundary layer turbulence. The equations are written in wave-followingcoordinates,andallvariablesaredecomposedintohorizontalmean,waveﬂuctuation,andturbulentﬂuctuation.Theformulationdeﬁnesthewave-inducedstressasasumofthewaveﬂuctuationstress(because of the ﬂuctuating velocity components) and a pressure stress (pressure acting on a tilted surface). The for- mulations can be constructed with different choices of mapping. Next, a large-eddy simulation result for wind over a sinusoidal wave train under a strongly forced condition is analyzed using the proposed formulation. The result clariﬁes how surface waves increase the effective roughness length and the drag coefﬁcient. Spe-ciﬁcally, the enhanced wave-induced stress close to the water surface reduces the turbulent stress (satisfying themomentumbudget).Thereducedturbulentstressiscorrelatedwiththereducedviscousdissipationrateoftheturbulentkineticenergy.Thelatterisbalancedbythereducedmeanwindshear(satisfyingtheenergy budget), which causes the equivalent surface roughness to increase. Interestingly, there is a small region farther above where the turbulent stress, dissipation rate, and mean wind shear are all enhanced. The observed strong correlation between the turbulent stress and the dissipation rate suggests that existing turbu- lence closure models that parameterize the latter based on the former are reasonably accurate.


Introduction
The wind stress (or the drag coefficient) at the ocean surface is an important parameter needed for ocean, atmosphere, and surface wave models. When a surface ocean wave field is fully developed, that is, is in equilibrium with local wind forcing, the wind stress is a function of local neutral wind speed (corrected for stability) and can be parameterized using a bulk formula. However, if the wave field is not in equilibrium, which is the norm rather than the exception, the wind stress may deviate significantly from the bulk parameterization and may require sea state-dependent parameterization with concurrent predictions of surface wave fields.
Many previous modeling studies have investigated how the wind stress and drag coefficient are modified by different sea states, including growing seas (e.g., Makin and Kudryavtsev 2002;Moon et al. 2004b; Kukulka and Hara 2008;Mueller and Veron 2009) and complex seas (e.g., Moon et al. 2004a;Donelan et al. 2012;Reichl et al. 2014). They all start with the momentum conservation constraint that the wind stress is equal to a sum of the momentum flux into surfaces waves (form drag of surface waves) and the momentum flux directly into the subsurface currents through viscous stress. The momentum flux into waves is normally evaluated by integrating the fluxes to all wave spectral components and may include explicitly the enhanced form drag due to breaking waves (e.g., Kudryavtsev and Makin 2001;Makin and Kudryavtsev 2002;Kukulka and Hara 2008;Mueller and Veron 2009;Banner and Morison 2010). The next step of the drag coefficient estimation is to model the feedback of the wave form drag on the mean wind profile. This step is needed to establish a relationship between the wind stress and the wind speed (normally at 10-m height). The wind profile in some studies is simply approximated using log-layer vertical wind profiles (e.g., Kudryavtsev and Makin 2001;Mueller and Veron 2009;Donelan et al. 2012). In this case, the wind profile is dependent only on the surface roughness parameter z 0 , that is, the feedback appears only in the parameterization of the sea state-dependent z 0 . Other studies explicitly account for the feedback of the modified turbulent stress due to wave form drag on the mean wind shear in the wave boundary layer using various turbulence closure models (e.g., Makin and Kudryavtsev 1999;Hara and Belcher 2004;Kukulka and Hara 2008). In such studies it is often assumed that the turbulent stress is reduced because of the increased wave-induced stress inside the wave boundary layer, where the total wind stress remains constant, and that this reduced turbulent stress is responsible for the reduction of the mean wind shear and the increase of the equivalent surface roughness (or the drag coefficient).
Although the wave boundary layer turbulence model is an essential component of the estimation of sea statedependent drag coefficients, its validity has not been thoroughly investigated either numerically or experimentally. This is because the wave modulation of turbulence mainly occurs very close to the water surface, often below the level of wave crests. Turbulence observations are extremely difficult to carry out very close to moving water surfaces. While numerical studies, such as direct numerical simulation (DNS) and large-eddy simulation (LES), can be carried out over wavy surfaces, interpretations of such results are not trivial. For example, the traditional definition of the wave-induced stress in terms of the wave-correlated velocity components (e.g., Makin and Kudryavtsev 1999;Hara and Belcher 2004) breaks down below the level of the wave crests. The more recent modeling studies of Sullivan et al. (2000), Chalikov and Rainchik (2011), and others, formulated in a wavefollowing coordinate system, clarify how the momentum flux is partitioned into the contribution of the wavecorrelated fluctuating velocities and the contribution of the wave-correlated pressure acting on a sloped surface and that the latter may become increasingly important very close to the water surface.
The main objectives of this study are 1) to develop a robust theoretical framework to describe and interpret wave boundary layer turbulence, which is applicable in areas very close to the water surface, and 2) to investigate the wave boundary layer turbulence and its impact on the mean wind profile and the drag coefficient using LES results.

Governing equations for wave-induced motions
Let us consider air with a constant density r a and a constant kinematic viscosity n a . Since this study focuses on processes inside the thin wave boundary layer just above the water surface, density stratification of the air and the Coriolis effect are ignored. We start with a rectangular coordinate (x 1 , x 2 , x 3 ) 5 (x, y, z), where x and y are horizontal and z is vertically upward, with z 5 0 at the mean water surface. The air velocities in (x 1 , x 2 , x 3 ) 5 (x, y, z) directions are denoted by (u 1 , u 2 , u 3 ) 5 (u, y, w). The continuity, momentum, and energy equations are written as where S ij 5 1 2 is the strain rate tensor, s ij 5 22nS ij is the viscous stress, p 5 p total /r a 1 gz is the dynamic pressure divided by r a , and p total is the total pressure. Next, we introduce the Reynolds decomposition of a variable a: where the overbar denotes Reynolds (ensemble) average, and the prime denotes a turbulent fluctuation. The Reynolds averaged equations of continuity, momentum, mean energy, and turbulence kinetic energy (TKE) become where t ij 5 u 0 i u 0 j , t ij 5 u 0 i u 0 j is the Reynolds stress, E 5 (1/2)u i u i is the mean kinetic energy, e 5 (1/2)t ii 5 (1/2)u 0 i u 0 i is the TKE, and is the viscous dissipation of the TKE. Since all the viscous terms except « are negligible outside the viscous sublayer, they have been omitted for simplicity. In this study, a simple periodic surface wave train is considered. The wave train is treated as a deterministic motion, that is, the wave motion is retained after ensemble averaging is taken. We introduce a second averaging, denoted by brackets h i applied in a horizontal x 2 y plane to filter out the wave-induced motions. (In general, this second averaging can be taken in time t instead of in x and y. However, when a coordinate system moving with a wave train is introduced later, the time averaging does not filter out the wave-induced motions.) The triple decomposition of a variable a is defined as a 5 a 1 a 0 5 hai 1ã 1 a 0 , where the Reynolds (ensemble) average a is split into the horizontal mean hai and the wave fluctuationã. The former is a function of z only. Then the continuity and momentum equations of the wave fluctuation are and where t w ij 5ũ iũj is the momentum flux due to wave fluctuations, andt w ij 5 t w ij 2 ht w ij i. The continuity equation of the horizontal mean requires that hwi 5 0, and the momentum equation governing the horizontal mean is The vertical (i 5 3) component of Eq. (14) yields ht w 33 i 1 h pi 1 ht 33 i 5 0, while the horizontal components yield that is, the wind stress is constant in z and is equal to a sum of the horizontally averaged turbulent stress ht i3 i and the horizontally averaged stress due to wave fluctuations ht w i3 i; the latter is often called ''wave-induced stress.'' Let us define the kinetic energy of the horizontal mean E m 5 (1/2)hu i ihu i i and the kinetic energy of the wave fluctuation E w 5 (1/2)t w ii 5 (1/2)ũ iũi . The energy equations for E m and E w are obtained by multiplying Eqs. (14) and (13) by hu i i andũ i , respectively: and the turbulence kinetic energy Eq. (9) becomes Here, 2ht w i3 i(›hu i i/›z), 2ht i3 i(›hu i i/›z), and 2ht ij (›ũ i /›x j )i denote energy transfers from horizontal mean to wave fluctuation, from horizontal mean to turbulence, and from wave fluctuation to turbulence, respectively. Vertical energy fluxes of horizontal mean, wave fluctuation, and turbulence are denoted by F m , F w , and F t , respectively.
If Eqs. (17) and (18) are added together, we obtain the energy equation for the sum of E w and e: Here, the first term is the total shear production (i.e., total loss of the horizontal mean energy), the second term is the total transport term, and the third term is the viscous dissipation. If we integrate Eq. (19) from z 5 0 to a reference height z 5 z r that is located above the wave boundary layer (where F w and F t are negligibly small), and set the surface current to be zero, we obtain where (F w 1F t ) z50 is equal to the energy flux into surface waves. Therefore, the relationship between wind stress and wind speed at z 5 z r can be obtained if the energy flux into the surface waves and the TKE viscous dissipation « below z 5 z r are known (Hara and Belcher 2004). In previous studies (e.g., Makin and Kudryavtsev 1999;Hara and Belcher 2004), Eqs. (12) to (18) are the basis of a wave boundary layer model used to estimate how the mean wind profile and drag coefficient are modified by surface waves. However, in rectangular coordinates these equations are valid only above the highest wave crest because the horizontally averaged variables, such as the mean wind speed hu i i, cannot be defined below. Therefore, the validity of these wave boundary layer turbulence models over real waves (whose amplitude is not infinitesimal) has not been investigated either experimentally or numerically (i.e., against direct numerical simulations or large-eddy simulations). This is particularly problematic in strongly forced conditions (when the wind speed is much larger than the wave phase speed) because the mean wind profile is modified mostly in a very thin layer whose height is often smaller than the wave amplitude.

Governing equations in wave-following coordinates
To investigate the wave-induced motions near or below the wave crest, we need to introduce a wave-following coordinate system. In this study we focus on a simple problem of a single periodic wave train with a fixed wavelength and a fixed phase speed. The wave shape is assumed unchanged as the waves propagate. We first introduce a frame of reference moving with the wave phase speed so that the wave motion becomes steady in time. Next, we introduce a coordinate mapping (without time dependence) from a wavy surface to a flat surface. (Although it is straightforward to introduce a time-dependent coordinate mapping following a time-dependent surface elevation field with multiwave components, such an approach will be discussed in a future study.) Let us introduce a coordinate system (j 1 , j 2 , j 3 ) 5 (j, h, z), j5j(x, y, z), h 5 h(x, y, z), z5z(x, y, z) , (21) such that z 5 0 is at the water surface and z 5 z as z/'. At this stage we do not need to specify the functional form of Eq. (21). It is expected that the horizontal coordinates (j, h) are either the same as or slightly depart from (x, y) and that constant z planes smoothly transition from the actual wavy water surface to a flat surface as z increases. We will later demonstrate that our results are relatively insensitive to different choices of mapping. The continuity and momentum equations are now written as (Anderson et al. 1984) where is the contravariant flux velocity perpendicular to a constant j i surface, and p, u i and s ij are the same variables as in rectangular (Cartesian) coordinates except they are now functions of j i and t. Notice that the momentum Eq. (24) is written for the Cartesian momentum component r a u i but now varies with the mapped coordinates j i . This expression is more convenient when the wind stress (vertical flux of x and y momentum) is considered later. Physically, the first term in the bracket u i U j represents a flux of x i momentum (momentum in the x i direction) across a constant j j plane due to an advective velocity U j , and the second term (1/J)p(›j j /›x i ) is a flux of x i momentum in the j j direction due to the pressure force applied on a constant j j plane. When the constant j j plane is not parallel to the x i axis, this pressure term introduces tangential stress. Finally, the energy equation can be derived by multiplying Eq. (24) by u i : Let us introduce the Reynolds decomposition as before. Since the wave field is independent of t in the mapped coordinate, all Reynolds averaged variables are independent of t. The Reynolds averaged equations of continuity, momentum, mean energy, and TKE become where t ij 5 u 0 i U 0 j , t ij 5 u 0 i U 0 j is equivalent to the Reynolds stress t ij in the rectangular coordinate, and t p ij 5 (1/J)p(›j j /›x i ) is the stress due to the Reynolds averaged pressure.
Next horizontal averaging, denoted by brackets h i, is introduced with the averaging performed in j and h (at a fixed z) to filter out wave-induced motions. Then, the continuity and momentum equations of the wave fluctuation become where t w ij 5ũ iŨj is the momentum flux due to wave fluctuations,t w ij 5 t w ij 2 ht w ij i, andt p ij 5 t p ij 2 ht p ij i. The horizontally averaged continuity equation requires that hWi 5 0. Note that in general hwi 6 ¼ 0. (For example, if we introduce a surface drift current that is larger in the windward side than in the lee side of the crest, hwi . 0 at z 5 0.) The horizontally averaged momentum equation becomes where the quantity inside the bracket is a flux of x i momentum in the z direction. The vertical component of Eq. (33) yields ht w 33 i 1 ht p 33 i 1 ht 33 i 5 0, and the horizontal components of Eq. (33) yield Again, this equation shows how the wind stress (horizontally averaged flux of x and y momentum in the z direction) is realized in the wave boundary layer. In contrast to Eq. (15) in rectangular coordinates, where the wind stress is a sum of the turbulent stress and stress due to wave fluctuations (wave-induced stress), in Eq. (34) the wind stress is a sum of the turbulent stress ht i3 i and the two wave-induced terms ht w i3 i and ht p i3 i. The first term ht w i3 i represents a flux due to wave fluctuations and is equivalent to ht w i3 i in the rectangular coordinates defined in the previous section; we call this stress ''wave fluctuation stress.'' The second term ht p i3 i is a flux due to pressure applied on a tilted constant z plane (tilted because of the waves); we call this stress ''pressure stress.'' (Note that this pressure stress is defined in the entire domain above the water surface, not just at the water surface.) Therefore, it is natural to define a sum of these two wave-induced terms ht w i3 i 1 ht p i3 i as wave-induced stress. Very close to the water surface (z/0) the wave fluctuation stress approaches zero (because W/0, i.e., the velocity normal to the stationary water surface must approach zero; see Sullivan et al. 2014), and the pressure stress dominates the wave-induced stress as expected. Far from the water surface (z/'), where a constant z plane becomes flat (z approaches z), the pressure stress becomes zero; hence, the wave fluctuation stress alone determines the wave-induced stress. Between these two limits both wave-induced terms may contribute to the wave-induced stress.
The energy equations for horizontal mean and wave fluctuation are obtained by multiplying Eqs. (33) and (32) by hu i i andũ i , respectively, and by using the following identity The resulting equations are The TKE equation is obtained from Eq. (30): If we compare Eqs. (36) to (38) with Eqs. (16) to (18), it is clear that the wave-induced stress ht w i3 i in the rectangular coordinate is replaced by the wave-induced stress (ht w i3 i 1 ht p i3 i) in the mapped coordinate. Adding Eqs. (37) and (38) yields the equation for the sum of the TKE and wave fluctuation energy: In summary, if a wave-following coordinate is introduced, it is possible to decompose a variable into the horizontal mean, wave fluctuation, and turbulence components and to derive the continuity, momentum, and energy equations everywhere, including regions below the wave crest level. In particular, the wave-induced stress is naturally defined such that it is a sum of the wave fluctuation stress (i.e., Reynolds-like stress) and the pressure stress. A wave-following coordinate also allows us to examine the energy budget (including the TKE dissipation rate) and the mean wind profile very close to the water surface and to clarify how surface waves modify the equivalent roughness length and the drag coefficient.
Note that the formulation in a wave-following coordinate is not new (e.g., Sullivan et al. 2000), and the same momentum Eq. (33) has been obtained by Chalikov and Rainchik (2011) with a particular choice of wave-following coordinates. The derivation of the energy equation in a wave-following coordinate system was made in earlier studies as well (e.g., Hsu et al. 1981).
Here, a similar but more general approach has been applied (with different choices of mapping) and has been extended to include energy Eqs. (36) to (39). Naturally, the definitions of the horizontal mean and wave fluctuations depend on a particular choice of mapping. We will therefore employ different mapping approaches and investigate the sensitivity of the results.

LES of wind over a periodic wave train
Next, we analyze an LES of wind over a periodic wave train using the formulation in wave-following coordinates with the triple decomposition of the variables outlined previously. Both mean wind and waves are assumed to be in the x direction, that is, they are aligned. The location of the water surface is specified as and the velocities at the water surface are set, using the linear deep-water wave solutions, as where a is the wave amplitude, k is the wavenumber, v 5 ffiffiffiffiffi ffi gk p is the angular frequency, and c 5 v/k is the phase speed. Therefore, the surface boundary condition accounts for the wave orbital velocity but no additional mean surface currents (drift velocities) are added. Both for the simulations and for the data analysis we introduce the following coordinate mapping: or inversely x 5 j, y 5 h, z 5 z(j, z) .
This mapping does not change the horizontal coordinates and only stretches or shrinks the vertical axis according to the water surface elevation. Such a mapping is preferable because it is straightforward to extend it to a wavy surface with multiwave components. Other commonly used coordinates over a single wave train include the area-conserving mapping (e.g., Belcher and Hunt 1993) and the conformal mapping (e.g., Benjamin 1959). These approaches modify the horizontal coordinate in a manner related to the wavelength, and the formulation becomes complex (i.e., not as convenient for application) if two or more wave trains of different wave lengths coexist. We will not examine such coordinates in this study. In the LES, z 5 z(x, z) is determined numerically, and in the analysis different choices of z 5 z(x, z) are examined.
The actual LES calculations are performed in a fixed frame of reference, that is, the waves vary in time and continually propagate through the computational box. The size of the LES computational domain is l j 3 l h 3 l z , with l j 5 l h 5 5l, and l z 5 l, where l 5 2p/k is the wavelength. Doubly periodic boundary conditions are imposed in the horizontal directions. At the top of the box, a slip (no tangential stress) condition is imposed, while at the water surface the roughness of the smaller unresolved waves is parameterized by setting the equivalent roughness length z 0 . The discretization employs (N j , N h , N z ) 5 (256, 256, 128) grid points. The vertical distribution of points in computational space is nonuniform. The spacing ratio between neighboring cells is held constant at 1.0028, with the first point off the water surface located at z 1 5 0. 0065l. The mapping between physical and computational space vertical coordinates is where the shape of the underlying wave is The wind forcing is applied by an external pressure gradient ›P/›x that yields a surface stress t s 5 (›P/›x)l and a surface friction velocity u * s 5 jt s j 1/2 . The simulation is carried forward for 50 nondimensional large-scale turnover times (t 5 50), approximately 130 000 time steps. The surface stress becomes statistically steady at nondimensional timet 5 25, and flow statistics are then computed over the intervalt 5 [25,50]. The details of the LES algorithm and numerical methods used to solve the governing equations are fully described in Sullivan et al. (2014).
If we normalize all the variables using a length scale 1/k, a velocity scale u * s , and a time scale 1/ku * s , the problem depends on three nondimensional parameters, that is, the wave steepness ka, the normalized roughness length kz 0 , and the normalized wind forcing u * s /c. In this study we choose a strongly forced condition where u * s /c 5 0. 632 with a finite wave steepness ka 5 0.226, which is typical for dominant wind waves in laboratory conditions and is also applicable to small-scale waves in the open ocean. The roughness length of the unresolved waves is set kz 0 5 0.002 70 such that the resulting wind field is consistent with typical observed conditions.
To analyze the LES results we do not use the mapping Eq. (44) of the LES calculations but employ the following vertical mapping: with varying s, such that z 5 0 is exactly at the water surface and z 5 z if kz 1. We start with s 5 1: The Jacobian of this transformation is calculated to be J 5 ›z ›z 5 1 1 2 ka cos(kj)e 2kz .
Since the flow is forced by a constant horizontal pressure gradient ›P/›x, the stress is not constant in z or in z. The horizontally averaged x momentum equation is modified to When the mapping Eq. (47) is used, h1/Ji 5 1, and the wind stress t wind 13 varies linearly in z: if normalized by the wind stress.
The horizontally averaged energy equations are also modified to If we add Eqs. (53) and (54), we obtain the energy equation for the sum of wave fluctuation and turbulence: where the first term is the energy input from the external pressure force, the second term is the shear production (conversion from the mean energy), the third term is the transport term, and the last term is the TKE viscous dissipation. If Eqs. (53) to (55) are multiplied by 2kz/u 3 * , with u 2 * 5 jt wind 13 j, the normalized energy equations become and In the constant stress layer over a flat surface, where the pressure forcing is zero and the mean wind profile is logarithmic, it is known that the normalized wind shear (›hui/›z)(kz/u * ) and the normalized dissipation h(1/J)«i(kz/u 3 * ) are both close to 1 and the transport term is small. Even if the stress is not strictly constant, the normalized wind shear is expected to be close to 1 over a flat surface (according to the mixing length scaling), provided the normalization is done using the height-dependent friction velocity u * instead of the constant surface friction velocity u * s .

a. Overview of LES results
In Fig. 1 the results of the LES are first presented in physical space x-z coordinates moving at the phase speed c. The streamlines of the (ensemble averaged) wind are shown in Fig. 1a. The figure clearly shows a cat's eye pattern (closed stream lines) on the lee side of the wave crest. This pattern arises because the wind blows to the right, but air velocity along the wave surface is always negative (to the left) in the coordinate system moving with the wave. Above the cat's eye the streamline is significantly modulated relative to the wave shape, although the flow is not separated in a sense that both the velocity and tangential stress along the wave surface are always negative. The pressure plot in Fig. 1b shows that the location of maximum pressure has significantly migrated downwind from the wave trough, that is, the high pressure acts on a positive surface slope and pushes the wave to the right (i.e., contributes to the air-water momentum flux). Figure 1c shows that the TKE dissipation is large near the surface as expected. (It varies like 1/z over a flat surface.) However, the high dissipation rate (strong turbulence) region appears to be advected by the mean flow (indicated by the streamline) and is detached from the surface at the location of the cat's eye. Below the cat's eye the dissipation rate is significantly reduced near the surface, suggesting that turbulence is very weak there. The TKE plot in Fig. 1d shows that the TKE is nearly constant (which is expected in a constant stress layer) but is also significantly reduced near the surface below the cat's eye. We will next show that this weakening of turbulence near the surface is related to the modification of the mean wind profile and the increase of the equivalent surface roughness.

b. Mean wind profile
Next, we introduce the mapping Eq. (47). This mapping allows us to define the horizontal mean and the wave fluctuations everywhere above the wave surface. In Fig. 2, the computed normalized mean wind speed hui/u * s and the computed normalized mean wind shear (›hui/›z)(kz/u * ) are plotted as a function of normalized height kz. From here on kz is always plotted in a log scale since we focus on the processes very close to the wave surface. The computed normalized mean velocity hui/u * s matches the theoretical surface value 2c/u * s 5 21. 58 at the roughness height kz 5 kz 0 5 0.002 70 as expected. It is seen that the computed normalized shear is close to 1 when kz is roughly between 0.7 and 3, suggesting that the mean wind profile is similar to that over a flat surface above kz 5 0.7. (Above kz 5 3, the results are affected by the top boundary and should be ignored.) Since there is a well-established region where the mean wind profile and the mean shear behave like those over a flat surface, The normalized wind shear has been defined such that the area integral of the normalized shear in Fig. 2b is approximately proportional to the increase of the normalized wind speed in Fig. 2a, that is, the area shaded in blue is approximately equal to the area with green hatches. (If the wind stress is constant in height, these two areas are exactly the same.) It is apparent that the increase of the equivalent surface roughness is mainly caused by the reduction of the wind shear (relative to that over a flat surface) very close to the surface (kz , 0.15). It is interesting that there is a small region where the wind shear is enhanced (kz between 0.15 and 0.7). However, the equivalent roughness length increases because the decrease of the wind shear below kz 5 0.15 is more significant than its increase above kz 5 0.15.

c. Momentum flux budget
Next we examine the momentum flux budget. In Fig. 3a the computed normalized turbulent stress ht 13 i/t wind 13 and computed normalized wave-induced stress (ht w 13 i 1 ht p 13 i)/t wind 13 , as defined in Eq. (51), are plotted against normalized height. Here, the turbulent stress is a sum of the resolved stress and the SGS stress.
The computed normalized total stress (a sum of the turbulent stress and the wave-induced stress), shown with a red line, is very close to 1 everywhere. This assures that the LES conserves momentum accurately.
It is evident that the wave-induced stress increases and the turbulent stress decreases very close to the surface (kz , 0.2). The wave-induced stress is roughly 40% of the total stress at/near the surface, that is, the form drag of the surface waves supports about 40% of the total wind stress. Interestingly, there is a small region where the normalized wave-induced stress is negative (momentum flux is upward) and the turbulent stress is enhanced (larger than the total wind stress) around 0.2 , kz , 0.7. This range of kz is similar to the range where the mean wind shear is enhanced in Fig. 2b. Let us examine the two components of the wave-induced stress separately. As discussed earlier, the wave-induced stress is a sum of the wave fluctuation stress ht w 13 i/t wind 13 and the pressure stress ht p 13 i/t wind 13 . Very close to the surface the pressure stress dominates as expected. It is always positive (downward momentum flux) and monotonically decays with height. In contrast, the wave fluctuation stress is always negative (upward momentum flux). Its magnitude is the largest around kz 5 0.35 and approaches zero near the surface and far from the surface. Figures 3c-f show the spatial distribution of the stress components with their horizontal averages given in Fig. 3a. Note that Fig. 3d shows the excess normalized turbulent stress (t 13 /t wind 13 2 1) rather than the total turbulent stress. The streamlines in the mapped coordinate are shown in Fig. 3b   shows that the pressure stress is the largest at the wave surface where the high pressure acts on a positive surface slope (roughly p , kj , 13p/8). It is also positive where the negative pressure acts on a negative surface slope (roughly 0 , kj , 7p/8). It monotonically decays with height. Figure 3e shows that the wave fluctuation stress is significant only at midheights (around 0. 05 , kz , 1). The two strongly negative regions appear where the streamline leaves the surface (with positiveũ and positiveW) and where the streamline approaches the surface (with negativeũ and negativeW ). It is interesting that the location of the enhanced turbulent stress in Fig. 3d is not correlated with the locations of the negative wave-induced stress in Fig. 3c. However, they exactly compensate each other when they are averaged horizontally in Fig. 3a.

d. Energy budget
The energy budgets of the wave fluctuation energy E w , TKE e, and the sum of the two (E w 1 e) are examined in Fig. 4. All the terms in Eqs. (56) to (58) Figure 4a shows that the wave fluctuation energy E w is generated by shear production (conversion of mean energy) near the surface, is transferred upward, and then is converted back to the mean energy (negative shear production). In addition, some E w is converted to TKE near the surface, but some TKE is converted back to the E w farther above. The effect of the pressure forcing (dark green) is negligibly small, and the overall energy conservation (red) is reasonably accurate. All the terms are negligible above around kz 5 1.
The TKE budget is shown in Fig. 4b. The normalized shear production term and the normalized viscous dissipation term are both close to 1, and the transfer term is small above around kz 5 0.7, suggesting that the TKE budget is similar to that over a flat surface. The shear production is significantly enhanced around kz 5 0.35 because both the mean wind shear and the turbulent stress are enhanced compared to those over a flat surface. The TKE dissipation is also enhanced in this region, although some of the TKE is transferred above and below instead of being dissipated. In contrast, below about kz 5 0.15 both the shear production and the viscous dissipation are significantly reduced.
Finally, the budget of the wave fluctuation energy plus TKE (E w 1 e) is examined in Fig. 4c. Note that the shear production term (blue) is now identical to the normalized mean wind shear examined in Fig. 2b. Therefore, this plot helps us understand how the mean wind profile is modified by the surface waves. It is clear that the shear production term (blue) is mostly balanced by the normalized viscous dissipation term (magenta). The contribution of the flux term (cyan) is not negligible but is relatively small. This suggests that the reduction/ enhancement of the mean wind shear is closely related to the reduction/enhancement of the TKE viscous dissipation. Figure 4d shows the spatial distribution of the magnitude of the excess normalized TKE viscous dissipation («kz/Ju 3 * 2 1), such that its horizontal average plus 1 is identical to the negative of the magenta line in Fig. 4a. This figure shows significant reduction of the dissipation rate near the surface, particularly at the location of the cat's eye pattern, and enhancement of the dissipation rate just downstream of the top of the cat's eye pattern. This spatial distribution of the excess normalized TKE viscous dissipation is quite similar to that of the normalized excess turbulent stress (Fig. 3d), suggesting that the reduction/enhancement of the TKE viscous dissipation is correlated with the reduction/enhancement of the turbulent stress.
Based on the above analyses we can summarize the relationship between the modification of the mean profile and the modification of the turbulence due to surface waves. When the wave-induced stress increases and the turbulent stress decreases (from the momentum flux budget) very close to the surface (roughly kz , 0.15), the TKE dissipation rate also decreases. The reduction of the TKE dissipation rate is balanced by the reduction of the mean wind shear (from the energy budget). This reduction of the mean wind shear makes the equivalent surface roughness increase. Interestingly, exactly opposite trends (decrease of the wave-induced stress, increase of the turbulent stress, increase of the TKE dissipation, and increase of the mean wind shear) appear around kz 5 0.35. However, the effect of the enhanced wave-induced stress (roughly kz , 0.15) is stronger than the effect of the reduced wave-induced stress (around kz 5 0.35), and the equivalent roughness length increases because of the surface waves.

e. Discussion of different mappings
In this subsection, we investigate how the above analyses of the LES results changes if a different mapping is introduced. It is expected that the most significant changes occur if the wavy surface at the air-water interface transitions to a flat surface at a different rate as z increases. In particular, it is possible to introduce a vertical coordinate that is not stretched or compressed but simply translates up and down as the surface moves. With such a mapping the wave fluctuation terms appear at all heights, even at a large height where true wave effects are negligible. (For example, a simple uniform horizontal wind velocity u would introduce a wave fluctuationW far away from the surface with this mapping.) Let us examine two more cases of s 5 2 and s 5 0.1 in the mapping Eq. (46). The former transitions from a wavy surface to a flat surface twice as fast, while the latter transitions 10 times slower. The latter case is very similar to the translating vertical coordinate (no stretching or compressing) discussed above. The results of the momentum budget and the energy budget with different mappings (different s) are presented in Fig. 5. Figure 5a shows that the energy budget terms are quite insensitive to different mappings. Both the normalized mean wind shear profile (blue) and the profile of the normalized dissipation (magenta) are hardly affected. Figure 5b shows that both the turbulent stress (blue) and waveinduced stress (magenta) are quite insensitive to different mapping as well. They are enhanced and reduced in almost identical manners regardless of the mapping. The profiles of the pressure stress (cyan) and wave fluctuation stress (green) are significantly modified by different mapping. As s decreases (as the waviness persists to higher elevations), the magnitude of both stress components increases at midheights, and its decay is much slower with height. This is not surprising since smaller s tends to introduce artificially large wavy fluctuations as explained above. Nevertheless, once the two stress components are added, the effects of different mapping mostly cancel out and the resulting wave-induced stress is quite independent of mapping. From the momentum conservation, the turbulent stress profile is quite independent of mapping as well.
In summary, most of the results presented in this study are quite robust since they are not significantly affected by different choices of mapping, provided the wavy airwater interface smoothly transitions to a flat surface. The enhanced wave-induced stress, the reduced turbulent stress, the reduced TKE dissipation rate, and the reduced mean wind shear are all robust features inside the wave boundary layer very close to the air-water interface, and they explain how the equivalent surface roughness is increased by surface waves.

f. Implications for observations made from a moving platform
Field measurements of the wind stress are sometimes performed using anemometer measurements from moving platforms, such as ships and buoys. Although measured velocities are carefully motion corrected before the stress calculations are made, such estimates may still be different from those from a fixed platform if the wind measurement is performed at different elevations depending on the phase of the surface wave. For example, if the platform is wave following, the wind measurement is effectively performed at a constant z with s 1 instead of at a constant z. Therefore, it is interesting to simulate such observations using the LES results. Assuming that the stress estimation is made from measured u and w (by subtracting their time mean, taking their product, and taking its time mean), the resulting stress estimate corresponds to hu 0 w 0 1ũwi. In Fig. 5c, this simulated wind stress estimate hu 0 w 0 1ũwi/u 2 * [as well as the sum of the residual terms (t wind 13 2 hu 0 w 0 1ũwi)/u 2 * ] is shown for different mappings. Above kz 5 2 the simulated wind stress is very accurate. Between kz 5 0.5 and 2 the simulated wind stress remains quite accurate. The error is larger with a smaller s but does not exceed 7%. In contrast, the stress below kz 5 0.4 is significantly underestimated by this simulation. In particular, the simulation almost entirely misses the flux around 0.05 , kz , 0.1, perhaps because the simulated wind stress misses the important contribution of the pressure stress very close to the surface.
It should be emphasized here that this is a single result with a particular wind and wave condition. More LES simulations with different wind forcing conditions must be performed before any conclusions are drawn regarding the accuracy of wind stress measurements from a moving platform.

g. Turbulence closure inside the wave boundary layer
As discussed in section 1, the predictions of the sea state-dependent drag coefficient are often carried out by first estimating the total wave-induced stress by integrating the contributions from all spectral components and then imposing a turbulence closure model to relate the reduced turbulent stress and the modified mean wind profile. The first step assumes that the waveinduced stress from different wave components can be simply summed up, that is, the wave-induced stress is mainly determined by the first harmonics of velocities and pressure, which are linearly correlated with the wave elevation. We may test the validity of this assumption using the LES results. In Fig. 6a, the calculated wave-induced stress components (wave fluctuation stress, pressure stress, and total) using the first harmonics only (dashed lines) are compared with the original calculations (solid lines) that include the higher harmonics. The two results are almost identical. This is surprising since the spatial distribution of these stress components, shown in Figs. 3c-f, look quite nonlinear. Nevertheless, this result suggests that estimations of the wave-induced stress and the reduced turbulent stress over multiwave components (as routinely done in many modeling studies) are reasonably accurate even if the steepness of each wave component is not very small.
We next test the validity of some of the existing turbulence closure models. While some studies employ rather complex schemes [e.g., the higher-order turbulence closure scheme by Chalikov and Rainchik (2011)] that are difficult to test, others use simple parameterizations of the eddy viscosity and/or the viscous dissipation rate « in terms of the reduced turbulence stress, which are easily tested using the LES results. For example, Makin and Kudryavtsev (1999) parameterize the dissipation rate «(z) as proportional to [t(z)/t] 3/4 , while Hara and Belcher (2004) parameterize «(z) as proportional to [t(z)/t] 3/2 , where t is the total wind stress and t(z) is the reduced turbulent stress at a height z. If these parameterizations are introduced in the present analysis in mapped coordinates, the normalized dissipation h(1/J)«i(kz/u 3 * ) is parameterized as These two parameterizations are tested in Fig. 6b. It is seen that the second parameterization works quite well, while the first parameterization underestimates the wave effect. In general, the observed strong correlation between the reduced turbulent stress and the reduced TKE dissipation rate suggests that the parameterization of the latter based on the former is a reasonable approach. Of course, more LES simulations with different wind forcing conditions are needed to obtain more conclusive results.
h. Surface stress and wave growth rate Since the energy flux EF into surface waves is entirely because of the normal stress t n and the tangential stress t t applied on the tilted wave surface, it can be expressed as EF 5 EF n 1 EF t , EF n 5 u n t n cosu , EF t 5 u t t t cosu , where u is the angle of the surface tilt, u n and u t are the normal and tangential components of the wave orbital velocity, and EF n and EF t are the energy fluxes due to normal stress and tangential stress, respectively. Note that the factor 1/cosu is needed to account for the increase of the surface area due to the surface tilt. If the wave surface is smooth, there is a viscous sublayer along the wave surface, and the tangential stress is determined by the tangential viscous stress. The normal viscous stress is zero, and the normal stress is equal to the surface pressure. In this study, however, we parameterize the impact of unresolved small scale waves using a prescribed roughness length. Therefore, the surface tangential stress is determined by the roughness length and the horizontal velocity at the first grid level (using the law of the wall), while the total normal stress is a sum of the pressure and the turbulent normal stress, which is evaluated at the first grid level rather than at the wave surface. (It is expected that the total normal stress is almost constant in the local normal direction above the wave surface over a distance that is much smaller than the wavelength. We have ascertained that the LES result indeed satisfies this expectation over a first few grid levels above the water surface. Therefore, it can be assumed that the total normal stress evaluated at the first grid level is almost identical to that at the true water surface.) In Fig. 7, we plot the tangential stress along the wave surface as well as the normal stress components (pressure, turbulent normal, and total) evaluated along the first grid level. As expected, the total normal stress variation is significantly larger than the tangential stress variations. Below the cat's eye both stresses are close to zero. Using these surface stress results, the normalized energy flux (EF/u 3 * s ) into the surface waves is calculated to be EF 5 0.782, EF n 5 0.676, and EF t 5 0.106. Therefore, the tangential stress accounts for 14% of the total energy flux. However, this number likely varies if the equivalent surface roughness is allowed to vary along the wavy surface (reflecting the modulation of small-scale waves). Interestingly, the DNS by Yang and Shen (2010) of turbulent wind over waves with a similar wave age shows a comparable contribution of the tangential stress to the energy flux into waves, even though the physical meaning of the tangential stress is very different. (The tangential stress in the DNS is the viscous stress applied on a smooth wavy surface, while our tangential stress is mostly the form drag due to unresolved smaller waves.) It is noteworthy that the contribution of the total normal stress (EF n 5 0.676) consists of the pressure contribution (0.603) and the turbulent normal stress contribution (0.073). Therefore, the pressure contribution alone accounts for 77% of the total energy flux, missing 14% by the turbulent tangential stress contribution and 9% by the turbulent normal stress contribution.
The wave growth rate b (in a dimensional form) is often expressed as where c b is a nondimensional coefficient, and r w is water density. Since the (dimensional) energy flux is equal to a product of the wave energy (1/2)r w ga 2 and the growth rate b, the coefficient c b can be expressed as The LES result then yields c b 5 19.4. If the tangential stress contribution is ignored, c b 5 16.7. If the contribution of pressure alone is used, c b 5 14.9. These numbers are near the lower end of observational results and are quite consistent with previous theoretical and numerical results of the linear wave growth rate in the strongly forced conditions (see Belcher 1999), even though our LES calculation has been performed with a relatively large wave steepness.

Concluding remarks
We have derived the momentum and energy equations inside the wave boundary layer by introducing a wave-following coordinate and triple decomposition (horizontal mean, wave fluctuation, and turbulent fluctuation) of variables. The formulations are valid very close to the water surface, even below the wave crest level, and can be derived with different choices of mapping. The formulations naturally define the waveinduced stress as a sum of the wave fluctuation stress and the pressure stress and show that the sum of the waveinduced and turbulent stresses remains constant with height. They also describe the energy balance occurring inside the wave boundary layer.
Next, an LES result of wind over the sinusoidal finiteamplitude wave train (in a strongly forced condition) has been analyzed using the proposed formulations with three different coordinate-mapping choices. The results show that the enhanced wave-induced stress very close to the water surface (around kz , 0.15) reduces the turbulent stress (from the momentum budget). The reduced turbulent stress is correlated with the reduced TKE viscous dissipation rate. The reduced dissipation rate is then balanced by the reduced mean wind shear (from the energy budget), which causes the equivalent roughness length to increase. Interestingly, the exactly opposite trends (increased turbulent stress, increased dissipation rate, and increased mean wind shear) occur around 0.15 , kz , 0.7 and reduces the overall increase of the roughness length and drag coefficient. These results are quite insensitive to different choices of mapping. The observed strong correlation between the dissipation rate and the turbulent stress suggests that the existing parameterization of the former in terms of the latter is a reasonable approach.
There are many remaining questions to be answered. So far, only one case of wind forcing u * s /c, normalized roughness length kz 0 , and wave steepness ka has been studied. Clearly, more LES experiments varying all these parameters are needed to fully understand wave boundary layer turbulence. Furthermore, a surface drift velocity can be added, and both the roughness length and the drift velocity can be made variable along the wavy surface. Conditions of misaligned wind and waves are of significant interest as well.
Acknowledgments. TH thanks the National Center for Atmospheric Research (NCAR) and University of Rhode Island for supporting his collaboration with PPS at NCAR in 2012. TH also acknowledges support of the National Science Foundation (Physical Oceanography). FIG. 7. Surface stress components (normalized by the surface friction velocity squared): turbulent tangential stress (magenta), total normal stress (solid blue), pressure (dotted blue), and turbulent normal stress (dashed blue). The tangential stress is evaluated along the water surface, while the normal stress components are evaluated along the first grid level.