## 1. Introduction

Estimation of the wind stress (drag coefficient *C*_{d} or roughness length *z*_{0}) over the sea surface is one of the most crucial issues in oceanic and atmospheric modeling, including tropical cyclone and storm surge modeling and forecasting. Although the wind stress has been a subject of study for over 50 yr, present parameterizations of the wind stress have still significant limitations, especially in high wind conditions (Jones and Toba 2001).

One of the main uncertainties regarding the roughness length estimation is the effect of the ocean surface wave field. There have been a number of studies that relate the equivalent roughness to wave parameters, such as the wave age (*c*_{p}/*u*∗), where *c*_{p} is the phase speed at peak frequency and *u*∗ is the friction velocity, representing the state of growth of wind waves relative to local wind forcing (e.g., Toba et al. 1990; Smith et al. 1992; Johnson et al. 1998; Drennan et al. 2003). But their results are far from conclusive.

Several approaches have been developed to predict the drag coefficient by explicitly calculating the wave-induced stress, that is, the stress supported by surface waves (Janssen 1989; Chalikov and Makin 1991; Makin et al. 1995; Makin and Mastenbroek 1996; Makin and Kudryavtsev 1999; Hara and Belcher 2004). These approaches are based on conservation of momentum over the ocean surface, which requires that the total stress is independent of height in the lower part of the atmospheric boundary layer (i.e., the wave boundary layer). The total stress inside the wave boundary layer is expressed as the sum of the turbulent stress and the wave-induced stress except inside the viscous sublayer. Since the wave-induced stress supports a significant part of the total stress near the water surface, the turbulent stress must be reduced. In most studies, eddy viscosity models are used to relate the reduced turbulent stress with the mean wind profile to estimate the drag coefficient (Janssen 1989; Chalikov and Makin 1991; Makin et al. 1995; Makin and Kudryavtsev 1999). Hara and Belcher (2004) introduced a new model based on the combination of momentum conservation and energy conservation inside the wave boundary layer. They used an analytical form of the equilibrium surface wave spectrum (Hara and Belcher 2002) to obtain an analytical expression for the equivalent surface roughness over mature seas, and investigated how the equivalent roughness depends on different external parameters. While most past studies treated surface waves as superposition of linear sinusoidal wave components, effectively neglecting the effect of surface breaking waves, Kudryavtsev and Makin (2001) and Makin and Kudryavtsev (2002) included the effect of wave breaking in a wave boundary layer model and predicted significant increase of the drag coefficient at high winds and over younger wave fields.

In order to predict the air–sea momentum flux in real wave fields using wave boundary layer models, prediction of realistic directional spectra is an important prerequisite. There have been considerable efforts made to predict the directional spectra of surface waves (e.g., WAMDI Group 1988; Komen et al. 1994). Recently, a new ocean surface wave model, WAVEWATCH III (Tolman 2002b), was developed at the National Oceanic and Atmospheric Administration/National Centers for Environmental Prediction (NOAA/NCEP) in the spirit of the wave model (WAM; WAMDI Group 1988). The WAVEWATCH III (WW3) has been validated over a global-scale wave forecast and a regional wave forecast (Tolman 1998, 2002a; Tolman et al. 2002; Wingeart et al. 2001) and consistently showed a good performance. Under hurricane conditions, where wave prediction is most difficult due to complicated and quickly varying wind forcing in time and space, Moon et al. (2003) have evaluated the performance of wave spectral simulation using the WW3. This study has shown that the model-simulated hurricane directional spectra are in very good agreement with the observed spectra obtained by National Aeronautics and Space Administration's (NASA's) scanning radar altimeter in Hurricane Bonnie (1998).

In this study we investigate the effect of surface waves on air–sea momentum flux over mature and growing seas. The atmospheric surface layer is assumed to be neutrally stable, and no stability corrections are made throughout this study. In the frequency range near the spectral peak, WW3 is used to estimate the wave spectra, and the spectra in the high-frequency range (equilibrium range) beyond the model resolution are parameterized by an analytical model developed by Hara and Belcher (2002). The full wave spectrum is then introduced to the wave boundary layer model of Hara and Belcher (2004) to calculate the drag coefficient and wind stress.

While Hara and Belcher (2004) treated the wind stress as a scalar quantity since the wave field was assumed to be symmetric relative to the mean wind direction, we treat the stress as a vector quantity. This is because, for the complex seas forced by nonuniform and nonstationary winds, the existence of dominant waves that propagate at a large angle to the local wind may influence both the magnitude of the wind stress and the angle between the wind stress and the mean wind (Smith 1980; Geernaert 1988; Rieder et al. 1994; Bourassa et al. 1999; Grachev et al. 2003). The method allows us to estimate the stress vectors and to investigate the alignment between the local wind and wind stress.

The present paper focuses on the effect of mature and growing seas forced by uniform winds on the momentum fluxes. For this purpose, a set of idealized experiments is designed for constant and spatially uniform winds from 10 to 45 m s^{−1}. The effects of complex seas forced by nonuniform wind, especially under tropical cyclones, are investigated in the companion paper (Moon et al. 2004).

A brief outline of the WW3 is given in section 2. The analytical model of wave spectra in equilibrium range is introduced in section 3. Section 4 describes the wave boundary layer model. Section 5 presents the detailed procedure to determine drag coefficient using the full wave spectrum and the wave boundary layer model. Sections 6 and 7 describe experiments with constant winds from 10 to 45 m s^{−1} to investigate drag behavior in mature and growing sea conditions, and these results are compared with available observational data. The summary and conclusions are given in the last section.

## 2. Wave spectrum near the peak

Surface wave spectra in the frequency range near the spectral peak are estimated by a well-tested ocean surface wave model, WAVEWATCH III (WW3). The WW3 explicitly accounts for wind input, wave–wave interaction, and dissipation due to whitecapping and wave–bottom interaction. It solves the spectral action density balance equation for directional wavenumber spectra. The implicit assumption of these equations is that the medium (depth and current) as well as the wave field varies in time and space scales that are much larger than the corresponding scales of a single wave. The source terms of the WW3 use nonlinear wave–wave interactions using a discrete interaction approximation (DIA) modified by Tolman and Chalikov (1996), input and dissipation from Tolman and Chalikov (1996), and bottom friction from the Joint North Sea Wave Project (JONSWAP) as used in the WAM (WAMDI Group 1988). A detailed description of the model is given by Tolman (2002b).

The parameters of the WW3 used in this study are 1800 s (time step and wind input interval), 24 directions (directional resolution), ^{1}_{12}° × ^{1}_{12}° (spatial grid resolution), and 2000 m (water depth). The grid is regularly spaced by a longitude–latitude grid. It extends 3000 km in the south–north direction and 1500 km in the east– west direction. The wave spectrum is discretized using 40 frequencies extending from 0.0285 to 1.1726 Hz (wavelength of 1.1–1920 m) with a logarithmic increment *f*_{n+1} = 1.1*f*_{n}, where *f*_{n} is the *n*th frequency. The WW3 provides output of two-dimensional directional wave spectra, as well as mean wave parameters such as significant wave heights (4*E**π**k*^{−1}*π**σ*^{−1}*f*_{p}), and peak direction. Here, *E* is *∫*^{2π}_{0}*∫*^{∞}_{0}*ψ*(*f,* *θ*) *df* *dθ,* *k* is the wavenumber, *σ* is the radian frequency, *ψ*(*f,* *θ*) is the frequency spectrum, and the overbar denotes *E*^{−1} *∫*^{2π}_{0}*∫*^{∞}_{0}*zψ*(*f,* *θ*) *df* *dθ* for a parameter *z* (Tolman 2002b). The peak (dominant) frequency is calculated from the one-dimensional frequency spectrum using a parabolic fit around the discrete peak. The peak (dominant) wavenumber (wavelength) is calculated from the peak frequency using the dispersion relation.

## 3. Tail of wave spectrum

*ψ*(

*k,*

*θ*) and the degree of saturation

*B*(

*k,*

*θ*) in the equilibrium range are expressed aswithwhere

*k*is the wavenumber,

*θ*is the wave direction relative to the mean wind direction,

*c*

_{β}is a coefficient in the wave growth rate parameterization, and

*h*(

*θ*) is the directionality of the wave growth rate. (See the appendix for more details of the model.) We set

*c*

_{β}= 40, following Plant (1982), and

*c*

_{θ}= 3

*π*/8 corresponding to

*h*(

*θ*) = cos

^{2}(

*θ*)(−

*π*/2 <

*θ*<

*π*/2),

*h*(

*θ*) = 0 (|

*θ*| ≥

*π*/2). The equilibrium spectrum depends on a single parameter

*k*

_{s}, called a sheltering wavenumber. In this study the sheltering wavenumber is determined empirically so that the spectral tail is smoothly connected to the WW3 spectrum as described in section 5.

## 4. Wave boundary layer model

*τ*_{tot}, which is the sum of the turbulent stress vector

*τ*_{t}(

*z*) and the wave induced stress vector

*τ*_{w}(

*z*), is constant inside the WBL:

*τ*_{tot}

*τ*_{t}

*z*

*τ*_{w}

*z*

*z*is the vertical coordinate measured (upward) from the instantaneous water surface. The energy conservation is expressed aswhere

**u**is the mean wind vector, Π is the vertical transport of the kinetic energy of the wave-induced motions, Π′ is the vertical transport of the turbulent kinetic energy,

*ρ*

_{a}is the air density, and ε is the viscous dissipation rate of the turbulent kinetic energy. The dissipation rate is parameterized in terms of the turbulent stress:where

*κ*is the von Kármán constant. Both the wave induced stress vector

*τ*_{w}(

*z*) and the kinetic energy transport Π(

*z*) of the wave-induced motions may be calculated explicitly if the wave spectrum is known. The gradient of the turbulent kinetic energy transport

*d*Π′/

*dz*is negligibly small relative to the other terms in (4) inside the WBL. Therefore, for a given wind stress vector

*τ*_{tot}and a given wave field, Eq. (4) can be integrated up to the top of the wave boundary layer to estimate the mean wind vector, the drag coefficient, and the equivalent surface roughness. The full details of the WBL model of Hara and Belcher (2004) are given in the appendix.

## 5. Procedure of stress vector calculation

As schematically shown in Fig. 1, the wind stress vector and the vertical wind profile are estimated as follows.

- The surface directional wave spectrum near the spectral peak is evaluated in WW3 with a given 10-m wind input.
- The WW3 output is truncated at a cutoff frequency
*f*_{c}. - For frequencies above
*f*_{c}, the tail part of the spectrum is specified by the equilibrium wave spectrum model. - The full wave spectrum is then introduced to the WBL model. The wind stress vector and the vertical wind profile are estimated in the following three steps.
- Starting with an initial estimate of the surface viscous stress vector, the wave-induced stress vector and the total wind stress vector are evaluated by integrating the product of the momentum input rate and the spectrum of the surface waves.
- The vertical mean wind profile is calculated based on the WBL model up to the top of the wave boundary layer.
- The 10-m wind speed vector is calculated and compared with the initial wind input of the WW3. If they do not agree, go back to step a with a modified estimate of the surface viscous stress vector. This iteration is repeated until the result converges using the Newton–Raphson method.

### a. Determination of cutoff frequency

In step 2 we define a cutoff frequency (wavenumber), where the WW3 spectrum is truncated and is connected to the equilibrium spectrum. This is because the WW3 uses its own parametric high-frequency tail in the equilibrium range, which is different from that used in this study. The WW3 defines its own cutoff frequency as 3.0*f*_{pi} to connect the prognostic spectrum and the parametric high-frequency tail. Here, *f*_{pi} is the peak input frequency of actively generated waves (Tolman 2002b). In the WW3, *f*_{pi} is estimated from the equivalent peak frequency of the positive part of the input source term to obtain consistent estimation of the *f*_{pi} even in complex multimodal spectra (Tolman and Chalikov 1996). In this study, we define the cutoff frequency *f*_{c} as 3.0*f*_{pi}, which is the same as that of the WW3. The corresponding cutoff wavenumber is defined by *k*_{c} = (2*πf*_{c})^{2}*g*^{−1}.

### b. Estimation of sheltering wavenumber

*k*

_{c}<

*k,*the sheltering wavenumber (

*k*

_{s}) is determined from the integration of Eq. (1) over all angles aswhere

*k*

_{c}is the cutoff wavenumber and

*ψ*(

*k*

_{c},

*θ*) is the WW3 directional wavenumber spectra at

*k*

_{c}.

More details of step 4 are given in the appendix. Important wavenumbers and stresses defined in this study and schematic picture of the wind stress calculation with height are presented in Fig. 2.

## 6. Results for mature seas

The effect of waves on air–sea momentum exchange for mature and growing seas is investigated by a series of numerical experiments, which are performed for spatially homogeneous winds from 10 to 45 m s^{−1}. We assume that northward winds blow over the whole model domain of 3000 km by 1500 km (in latitude and longitude directions) with 2000-m water depth. A central point at the northern part of the domain, where the effect of the model boundary is negligible, is selected to obtain the mean wave parameters and surface wave spectrum. The significant wave height (*H*_{s}) and mean wavelength (*L*) simulated with different wind speeds are presented in Fig. 3. They become constant within 72 h from the onset of the wind, suggesting that the wave field becomes fully developed. For the highest wind speed (45 m s^{−1}), the simulated values reach up to 47 and 960 m, respectively.

### a. Fully developed (mature) wave spectra

Mature (fully developed) surface wave spectra obtained at 72 h after the onset of wind are used for the present analysis. As explained in section 5, the spectral tail is attached to the WW3 spectra at *k* = *k*_{c}. Although this process ensures that the magnitude of *B*(*k*_{c}), which is the integrated degree of saturation over all angles at *k*_{c}, is matched between the WW3 output and the equilibrium range parameterization, the directionality of the spectra is not necessarily identical. Hara and Belcher (2004) set the directionality of the growth rate to be *h*(*θ*) = cos^{2}(*θ*), which yields the directional spreading of the equilibrium spectra being proportional to *h*^{1/2}(*θ*) = cos(*θ*).

Figure 4 shows a comparison between the directional-spreading function used in the equilibrium spectral model and the directional spreading simulated from the WW3 at different frequencies for various wind speeds. From the figure, the directional-spreading function of Hara and Belcher (2004) represents quite realistic directional spreading at the cutoff frequency of 3.0*f*_{pi}.

Figure 5 shows the degree of saturation *B*(*k*) over mature seas in the full range of wavenumbers at different wind speeds. The figure shows that *B*(*k*) in the equilibrium range increases, and the cutoff wavenumber decreases as wind speed increases.

The sheltering wavenumber *k*_{s} is a single dynamical variable to determine the equilibrium spectrum. For a fixed *k* the spectrum level increases as the sheltering wavenumber decreases. Figure 6 shows the sheltering wavenumber *k*_{s} estimated from the present approach versus the friction velocity *u*∗ for mature seas. The value is compared with two different estimates (solid and dash–dot lines) made by Hara and Belcher (2004), based on Banner and Peirson (1998) and Phillips (1985). The present results from the WW3 simulation are roughly in the middle of the lower and upper bounds of *k*_{s} estimated from the previous two approaches based on observational data. This agreement further confirms that the WW3 yields reliable estimates of the surface wave field.

### b. Momentum flux over mature seas

The nondimensional roughness length [or Charnock coefficient (Charnock 1955)] against the friction velocity for mature seas estimated from the present model is presented in Fig. 7. Our model results depend on one empirical coefficient *δ,* the normalized height of the inner layer (Belcher and Hunt 1993). Here, we have set *δ* = 0.01 so that our results of the Charnock coefficient best agree with previous observations of about 0.011 (e.g., Smith 1988; Fairall et al. 1996; Taylor and Yelland 2000). Hara and Belcher (2004) argue that *δ* should be of *O*(0.05) if it is interpreted as a blending height (Mason 1988). However, our estimates of the Charnock coefficient with *δ* = 0.05 are about 0.03 ∼ 0.04, which is somewhat larger than the commonly used values for mature seas.

There are several possible explanations why the estimates with *δ* = 0.05 suggested by Hara and Belcher (2004) overestimate the Charnock coefficient.

- The value of
*δ*is indeed smaller than 0.05; that is, the wave-induced stress decays faster than the prediction by Mason (1988). - The wave spectrum is overpredicted by WW3 or the equilibrium spectrum parameterization. If the true spectrum is lower than our model result, then
*δ*should be larger than 0.01 to obtain similar values of the Charnock coefficient. - The wave boundary model by Hara and Belcher (2004) overpredicts the surface wave effect; that is, the model underestimates the viscous dissipation rate of the turbulent kinetic energy (TKE) inside the wave boundary layer. If the true surface wave effect is smaller, then
*δ*should be larger than 0.01 to obtain similar values of the Charnock coefficient.

The present estimates in Fig. 7 show that the Charnock coefficient is mostly independent of wind speed except it is slightly larger at the lowest wind speed of 10 m s^{−1}. Our results therefore do not agree with those of Fairall et al. (2003), which suggest that the Charnock coefficient increases with wind speed from 10 to 18 m s^{−1}. However, as discussed in the same paper, the behavior of Charnock's coefficient at wind speeds above 10 m s^{−1} remains controversial because of the uncertainty of available field data.

Next, we present the estimates of drag coefficients *C*_{d} as a function of 10-m wind speed and compare them with empirical formulas based on various observations (Fig. 8). The *C*_{d} estimated from our model is consistent with most observation-based formulas for wind speeds below 26 m s^{−1}. There are few observations over 26 m s^{−1} wind speed for mature seas. Bulk formulas with a Charnock constant between 0.008 and 0.0185, extending up to 50 m s^{−1} wind, are most commonly used in the numerical modeling community (thick dash–dot– dot and and dash–dot lines, respectively). In high winds, *C*_{d} estimated from the present model lies in the middle of this range. The WW3 calculates drag coefficients internally, which are shown in Fig. 8 (+ symbol). These values are generally consistent with the bulk formula using a Charnock constant of 0.0185.

Figure 9 shows the mean wind profiles simulated with the wind speed of 15 and 45 m s^{−1} in mature sea conditions. Just below the top of the wave boundary layer (indicated by squares), the wind profile becomes slightly steeper because the loss of the kinetic energy of the mean flow is enhanced due to the energy flux into surface waves. However, the slope of the wind profile rapidly decreases toward the bottom of the wave boundary layer (indicated by diamonds) because the viscous dissipation of the TKE is reduced corresponding to the reduced turbulent stress. Above the top of the wave boundary layer, the wind profile is simply logarithmic corresponding to a constant turbulent stress. If this logarithmic profile is extended to zero wind speed, it gives the estimation of the equivalent roughness length *z*_{0}.

## 7. Results for growing seas

In the previous section, we have shown that the present model with the WW3 yields the Charnock coefficient that is consistent with the analytical results by Hara and Belcher (2004) for mature seas, except that the empirical coefficient *δ* must be reduced from 0.05 to 0.01. The model by Hara and Belcher (2004) is not applicable for growing seas since their model assumes that all wind-forced waves are in the equilibrium range; that is, the peak of the spectrum is outside the wind forcing range (*u*∗/*c*_{p} < 0.07). This condition is not satisfied for growing seas. With the WW3, however, we may explicitly calculate the spectrum near the peak and estimate the drag coefficient at any stages of the wave evolution.

### a. Wave spectrum of growing seas

As before, we first construct the full wave spectrum by connecting the WW3 directional spectra and the equilibrium spectra. Figure 10 shows the degree of saturation for constant winds of 15 and 45 m s^{−1} estimated at different wave development stages. The figure shows that for 15 m s^{−1} *B*(*k*) in the equilibrium range and *k*_{c} decreases as the wave field matures, while trends for 45 m s^{−1} are different and more complicated. In very strong wind conditions, *k*_{c} decreases with wave evolution (time), but *B*(*k*) in the equilibrium range increases with wave evolution in earlier stages and then reduces in later stages.

The estimated sheltering wavenumbers *k*_{s} for constant winds from 10 to 45 m s^{−1} are plotted against the wave age *c*_{p}/*u*∗ in Fig. 11. The figure shows that *k*_{s} varies with wave development from 10 min to 72 h after the onset of wind. At the lowest wind speed of 10 m s^{−1}, the wave age is already larger than 10 at 10 min. In contrast, at the highest wind speed of 45 m s^{−1}, the wave age is less than 3 at 10 min. It takes more than 9 h for the wave field to reach the wave age of 10. In general, the sheltering wavenumber *k*_{s} decreases with wave age for *c*_{p}/*u*∗ < 10 and increases with wave age for *c*_{p}/*u*∗ ≥ 10. Therefore, for a higher wind (≥20 m s^{−1}), *k*_{s} first decreases and then increases, while, for a lower wind (<20 m s^{−1}) *k*_{s} simply increases with wave evolution. This can be interpreted as the level of the equilibrium spectrum in very young waves is quite low at very high wind conditions.

### b. Drag coefficient for growing seas

Figure 12 shows the estimated nondimensional roughness length (Charnock coefficient) versus wave age *c*_{p}/*u*∗ at different wave development stages (same as in Fig. 11). With growing seas, present results for winds less than 30 m s^{−1} show that the Charnock coefficient is larger with younger seas except at very early wave stages, which is consistent with earlier studies. Dependence on the wave age is most pronounced at 10 m s^{−1} and gradually decreases as wind speed increases. For winds higher than 30 m s^{−1}, present results show a quite different trend; that is, very young waves yield a smaller Charnock coefficient and it increases with the wave age.

The Charnock coefficient is mainly determined by the two factors: the width and the level of the wind-forced part of the spectrum. The width monotonically increases as the wave age increases, while the level increases with the wave age for *c*_{p}/*u*∗ < 10 but decreases later. Therefore, for very young waves *c*_{p}/*u*∗ < 10 at higher winds, the Charnock coefficient must increase because both the width and the level increase. For *c*_{p}/*u*∗ > 10, the spectral width increases but the level decreases. Therefore, the trend of the Charnock coefficient is not trivial. Our model calculation suggests that at lower winds the decrease of the spectral level dominates and the Charnock coefficient decreases, while at higher winds the widening of the spectrum dominates and the Charnock coefficient increases.

This trend, indicating that younger waves produce lower *z*_{ch} at high wind speeds, is consistent with that of Toba et al. (1990), although their estimates are significantly larger. Makin and Kudryavtsev (2002) predicted that breaking waves might significantly enhance the wind stress (or Charnock coefficient), particularly over younger seas. This enhancement may provide an explanation of the difference between our results with Toba et al. (1990). We emphasize that this latter trend, that is, the Charnock coefficient increases with the wave age at an earlier stage, is observed only at very high winds because this earlier stage lasts less than 1 h and is simply not observable at lower winds.

Figure 13 shows drag coefficients against wind speeds at 10 m in different wave development stages (time). The estimated drag coefficients increase as wind speed increases and are between the two lines based on bulk formulas with the Charnock constants of 0.008 and 0.0185. Closer examination shows that young waves produce higher drag for winds less than 30 m s^{−1}, while young waves yield less drag for winds higher than 30 m s^{−1}. For young waves at very high winds, the estimated *C*_{d} shows a leveling-off with wind speed. This is because the wave-induced stress due to very young waves makes a relatively small contribution to the total wind stress in extremely high wind conditions as discussed earlier.

It is interesting that *C*_{d} estimated from the internal prediction of the WW3 shows a different trend with the present results; that is, *C*_{d} decreases significantly with the wave development, especially in high wind condition. In the WW3 *C*_{d} is estimated parametrically with the dependence of wave age following Janssen (1989) and Tolman and Chalikov (1996). This yields higher drag at young seas regardless of wind speed. The WW3 also shows a capping of *C*_{d}, especially in very early wave stages of high wind speed, but this is due to the limitation of a numerical range in internal calculation of *C*_{d}. Here it should be noted that the WW3 is a well-tuned surface wave prediction model and the dissipation term is empirically parameterized so that the results agree with observations. Since our model calculation suggests that the WW3 overestimates the wind stress and wind forcing, it is possible that the dissipation term is also overestimated in the model.

Recently, Powell et al. (2003) reported the drag coefficient derived from observed wind profiles using hundreds of GPS sondes launched from aircraft in tropical cyclones. It is among the first estimates made in a hurricane and in high wind speeds above 30 m s^{−1}. From this study, they found that the drag coefficient determined by the wind profile method above hurricane force shows a leveling off as the winds increase. This result is qualitatively consistent with the present trends showing a capping of *C*_{d} in very young waves under very high wind conditions.

## 8. Summary and conclusions

The effect of surface waves on air–sea momentum flux over mature and growing seas has been investigated by combining the WW3 ocean wave prediction model, the equilibrium spectrum model by Hara and Belcher (2002), and the wave boundary layer model by Hara and Belcher (2004). The combined model predicts the wind stress by explicitly calculating the wave-induced stress. This method allows us to estimate the wind stress vector for mature and growing seas as well as complex seas.

In numerical experiments performed for constant winds from 10 to 45 m s^{−1}, the effect of mature and growing seas on air–sea momentum flux is investigated. For mature seas, the Charnock coefficient is estimated to be about 0.01 ∼ 0.02 and the drag coefficient increases as wind speed increases, which are within the range of previous observational and theoretical studies. With growing seas, our results for winds less than 30 m s^{−1} show that the drag coefficient is larger for younger seas except at very early wave stages, as is consistent with earlier studies. For winds higher than 30 m s^{−1}, however, our results show a different trend; that is, very young waves yield less drag and it increases with wave ages. Our results, showing a capping of *C _{d}* in very young waves, are qualitatively consistent with recent estimates of Powell et al. (2003) at very high wind speed (over 40 m s

^{−1}) obtained from the measured wind profiles in tropical cyclones.

As suggested by Makin and Kudryavtsev (2002), breaking waves may significantly enhance the wind stress, particularly over younger seas. We expect that inclusion of breaking waves in the present model may significantly enhance the wind stress over younger seas, leading to more realistic estimation of the wind stress. However, even the most accurate measurement of breaking waves in the last few years indicates significant uncertainty. Our understanding of breaking processes is still very limited. Therefore, it is not feasible at present to include the breaking wave effect in our model and to quantify the breaking wave effect. Other factors, such as sea spray, may play an important role in the air–sea momentum exchange. Nevertheless, we emphasize that the trend of decreasing Charnock coefficient with very young seas found in this study is likely to be robust, simply because there are not as many waves (whether breaking or not), both in magnitude and in the wavenumber range, to support the momentum flux compared to the more developed stage.

## Acknowledgments

This work was supported by the U.S. National Science Foundation (Grant ATM 0001038). TH and SEB also thank the U.S. Office of Naval Research (CBLAST program, Grants N00014-011025 and N00014-0110133) and the U.S. National Science Foundation (Grant OCE0002314) for additional support.

## REFERENCES

Anderson, R. J., 1993: A study of wind stress and heat fluxes over the open ocean by the inertial dissipation method.

,*J. Phys. Oceanogr***23****,**2153–2161.Banner, M. L., , and W. L. Peirson, 1998: Tangential stress beneath wind-driven air–water interfaces.

,*J. Fluid Mech***364****,**115–145.Belcher, S. E., , and J. C. R. Hunt, 1993: Turbulent shear flow over slowly moving waves.

,*J. Fluid Mech***251****,**109–148.Bourassa, M. A., , D. G. Vincent, , and W. L. Wood, 1999: A flux parameterization including the effects of capillary waves and sea state.

,*J. Atmos. Sci***56****,**1123–1139.Chalikov, D. V., , and V. K. Makin, 1991: Models of the wave boundary layer.

,*Bound.-Layer Meteor***56****,**83–99.Charnock, H., 1955: Wind stress on a water surface.

,*Quart. J. Roy. Meteor. Soc***81****,**639–640.Donelan, M. A., 1990: Air–sea interaction.

*The Sea,*B. leMehaute and D. M. Hanes, Eds., Ocean Engineering Science, Vol. 9, Wiley and Sons, 239–292.Donelan, M. A., , and W. J. Pierson, 1987: Radar scattering and equilibrium ranges in wind-generated waves with applications in scatterometry.

,*J. Geophys. Res***92****,**4971–5029.Drennan, W. M., , G. C. Graber, , D. Hauser, , and C. Quentin, 2003: On the wave age dependence of wind stress over pure wind seas.

*J. Geophys. Res.,***108,**8062, doi:10.1029/2000JC000715.Fairall, C. W., , E. F. Bradley, , J. B. Edson, , and G. S. Young, 1996: Bulk parameterization of air–sea fluxes in TOGA COARE.

,*J. Geophys. Res***101****,**3747–3767.Fairall, C. W., , E. F. Bradley, , J. E. Hare, , A. A. Grachev, , and J. B. Edson, 2003: Bulk parameterization of air–sea fluxes: Updates and verification for the COARE algorithm.

,*J. Climate***16****,**571–591.Geernaert, G. L., 1988: Measurements of the angle between the wind vector and wind stress vector in the surface layer over the North Sea.

,*J. Geophys. Res***93****,**8215–8220.Grachev, A. A., , C. W. Fairall, , J. E. Hare, , J. B. Edson, , and S. D. Miller, 2003: Wind stress vector over ocean waves.

,*J. Phys. Oceanogr***33****,**2408–2429.Hara, T., , and S. E. Belcher, 2002: Wind forcing in the equilibrium range of wind-wave spectra.

,*J. Fluid Mech***470****,**223–245.Hara, T., , and S. E. Belcher, 2004: Wind profile and drag coefficient over mature ocean surface wave spectra.

*J. Phys. Oceanogr.,*in press.Janssen, P. A. E. M., 1989: Wave-induced stress and the drag of air flow over sea waves.

,*J. Phys. Oceanogr***19****,**745–754.Johnson, H. K., , J. Hfjstrup, , H. J. Vested, , and S. E. Larsen, 1998: On the dependence of sea surface roughness on wind waves.

,*J. Phys. Oceanogr***28****,**1702–1716.Jones, I. S. F., , and Y. Toba, 2001:

*Wind Stress over the Ocean*. Cambridge University Press, 303 pp.Komen, G. J., , L. Cavaleri, , M. Donelan, , K. Hasselmann, , S. Hasselmann, , and P. A. E. M. Janssen, 1994:

*Dynamics and Modeling of Ocean Waves*. Cambridge University Press, 520 pp.Kudryavtsev, V. N., , and V. K. Makin, 2001: The impact of air-flow separation on the drag of the sea surface.

,*Bound.-Layer Meteor***98****,**155–171.Large, W. G., , and S. Pond, 1981: Open ocean momentum flux measurements in moderate to strong wind.

,*J. Phys. Oceanogr***11****,**324–336.Makin, V. K., , and C. Mastenbroek, 1996: Impact of waves on air– sea exchange of sensible heat and momentum.

,*Bound.-Layer Meteor***79****,**279–300.Makin, V. K., , and V. N. Kudryavtsev, 1999: Coupled sea surface–atmosphere model. 1. Wind over waves coupling.

,*J. Geophys. Res***104****,**7613–7623.Makin, V. K., , and V. N. Kudryavtsev, 2002: Impact of dominant waves on sea drag.

,*Bound.-Layer Meteor***103****,**83–99.Makin, V. K., , V. N. Kudryavtsev, , and C. Mastenbroek, 1995: Drag of the sea surface.

,*Bound.-Layer Meteor***73****,**159–182.Mason, P. J., 1988: The formation of areally-averaged roughness lengths.

,*Quart. J. Roy. Meteor. Soc***114****,**399–420.Moon, I-J., , I. Ginis, , T. Hara, , H. Tolman, , C. W. Wright, , and E. J. Walsh, 2003: Numerical simulation of sea surface directional wave spectra under hurricane wind forcing.

,*J. Phys. Oceanogr***33****,**1680–1706.Moon, I-J., , I. Ginis, , and T. Hara, 2004: Effect of surface waves on air–sea momentum exchange. Part II: Behavior of drag coefficient under tropical cyclones.

,*J. Atmos. Sci***61****,**2334–2348.Phillips, O. M., 1977:

*The Dynamics of the Upper Ocean.*2d ed. Cambridge University Press, 336 pp.Phillips, O. M., 1985: Spectral and statistical properties of the equilibrium range in wind-generated gravity waves.

,*J. Fluid Mech***156****,**505–531.Plant, W. J., 1982: A relationship between wind stress and wave slope.

,*J. Geophys. Res***87****,**1961–1967.Powell, M. D., , P. J. Vickery, , and T. A. Reinhold, 2003: Reduced drag coefficient for high wind speeds in tropical cyclones.

,*Nature***422****,**279–283.Rieder, K. F., , J. A. Smith, , and R. A. Weller, 1994: Observed directional characteristics of the wind, wind stress and surface waves on the open ocean.

,*J. Geophys. Res***99****,**22596–22598.Smith, S. D., 1980: Wind stress and heat flux over the ocean in gale force winds.

,*J. Phys. Oceanogr***10****,**709–726.Smith, S. D., 1988: Coefficients for sea surface wind stress, heat flux, and wind profiles as a function of wind speed and temperature.

,*J. Geophys. Res***93****,**15467–15472.Smith, S. D., and Coauthors, 1992: Sea surface wind stress and drag coefficients: The HEXOS results.

,*Bound.-Layer Meteor***60****,**109–142.Taylor, P. K., , and M. A. Yelland, 2000: On the apparent “imbalance” term in the turbulent kinetic energy budget.

,*J. Atmos. Oceanic Technol***17****,**82–89.Taylor, P. K., , and M. A. Yelland, 2001: The dependence of sea surface roughness on the height and steepness of the waves.

,*J. Phys. Oceanogr***31****,**572–590.Toba, Y., , N. Iida, , H. Kawamura, , N. Ebuchi, , and I. S. F. Jones, 1990: The wave dependence of sea-surface wind stress.

,*J. Phys. Oceanogr***20****,**705–721.Tolman, H. L., 1998: Validation of a new global wave forecast system at NCEP.

*Ocean Wave Measurements and Analysis,*B. L. Edge and J. M. Helmsley, Eds., ASCE, 777–786.Tolman, H. L., 2002a: Validation of WAVEWATCH III version 1.15 for a global domain. NOAA/NWS/NCEP/OMB Tech. Note 213, 33 pp.

Tolman, H. L., 2002b: User manual and system documentation of WAVEWATCH-III version 2. 22. NOAA/NWS/NCEP/OMB Tech. Note 222, 133 pp.

Tolman, H. L., , and D. Chalikov, 1996: Source terms in a third-generation wind wave model.

,*J. Phys. Oceanogr***26****,**2497–2518.Tolman, H. L., , B. Balasubramaniyan, , L. D. Burroughs, , D. Chalikov, , Y. Y. Chao, , H. S. Chen, , and V. M. Gerald, 2002: Development and implementation of wind generated ocean surface wave models at NCEP.

,*Wea. Forecasting***17****,**311–333.WAMDI Group, 1988: The WAM model—A third generation ocean wave prediction model.

,*J. Phys. Oceanogr***18****,**1775–1810.Wingeart, K. M., , W. C. O'Reilly, , T. H. C. Herbers, , P. A. Wittmann, , R. E. Jenssen, , and H. L. Tolman, 2001: Validation of operational global wave prediction models with spectral buoy data.

*Ocean Wave Measurement and Analysis,*B. L. Edge and J. M. Hemsley, Eds., ASCE, 590–599.Yelland, M., , P. K. Taylor, , R. W. Pascal, , J. Hutchings, , and V. C. Cornell, 1998: Measurements of the open ocean drag coefficient corrected for airflow disturbance by the ship.

,*J. Phys. Oceanogr***28****,**1511–1526.

## APPENDIX

### Model of Wave Boundary Layer and Equilibrium Spectrum

#### Momentum conservation in the wave boundary layer

*τ*_{tot}

*τ*_{t}

*z*

*τ*_{w}

*z*

*β*

_{g}is the wave growth rate,

*ρ*

_{w}is the water density,

*σ*is the wave angular frequency,

*ψ*(

*k,*

*θ*) is the surface wave height spectrum,

*F*(

*z*;

*k*) is the vertical decay function of the wave induced stress due to waves at a wavenumber

*k,*and

*k*

_{m},

*k*

_{1}are the minimum and maximum wavenumbers of the wind-forced waves, respectively. If the decay function is approximated by a step function,

*F*(

*z*;

*k*) = 1 (

*z*≤

*δ*/

*k*) and

*F*(

*z*;

*k*) = 0 (

*z*>

*δ*/

*k*), where

*δ*is the normalized height of the inner layer (Belcher and Hunt 1993), the wave-induced stress at a height

*z*is simplythat is, it is equal to the integration of the momentum flux to surface waves in the wavenumber range of

*k*

_{m}<

*k*<

*δ*/

*z.*Above the wave boundary layer (

*δ*/

*k*

_{m}≤

*z*) the wave-induced stress is zero and the turbulent stress is equal to the total wind stress. Very near the water surface (

*z*<

*δ*/

*k*

_{1}) the turbulent stress is constant and is equal to the stress supported by viscosity inside the viscous sublayer, denoted by

*τ*_{ν}.

*β*

_{g}applied to waves at a particular wavenumber

*k*is assumed to be determined not by the total wind stress

*τ*_{tot}but by the turbulent stress evaluated at the top of the inner layer, which is called a

*local turbulent stress*and is denoted by

*τ*^{l}

_{t}

*k*). Using (A3) the local turbulent stress may be expressed asand the growth rate

*β*

_{g}is written aswhere

*ρ*

_{a}is the air density,

*c*is the wave phase speed, and

*θ*

_{τ}is the direction of the local turbulent stress. Following Plant (1982), we set

*c*

_{β}= 40 and

*k*

_{m}= 0.07

^{2}

*g*/

*u*

^{2}

_{∗}

*τ*_{tot}| =

*ρ*

_{a}

*u*

^{2}

_{∗}

*h*(

*θ*) of the growth rate is set to be

*h*(

*θ*) = cos

^{2}(

*θ*)(−

*π*/2 <

*θ*<

*π*/2), and

*h*(

*θ*) = 0 (|

*θ*| ≥

*π*/2). In Hara and Belcher (2004)

*δ*is set to be 0.05 following Mason (1988). However, in this study we determine

*δ*empirically so that the drag coefficient at lower wind speeds is consistent with the previous observations and parameterizations. This is the only tuning parameter in the model since the other parameters are all well constrained.

Equations (A4) and (A5) show that if the surface viscous stress *τ*_{ν} and the surface wave spectrum are given, the local turbulent stress *τ*^{l}_{t}*k*), the turbulent stress *τ*^{l}_{t}*z*), and the total wind stress *τ*_{tot} can be calculated by integrating (A4) from *k* = *k*_{1} to *k* = *k*_{m}. In this study the upper bound of the wavenumber range is set as *k*_{1} = 400 rad m^{−1} because Hara and Belcher (2004) have shown that the results of the drag coefficient are insensitive to the choice of *k*_{1} provided it is chosen to be above 400 rad m^{−1}. Important wavenumbers defined in this study and a schematic picture of the wind stress calculation are presented in Fig. 2.

Once the vertical profile of *τ*_{t}(*z*) is determined, the mean wind profile is obtained from the energy conservation constraint as described next.

#### Energy conservation in the wave boundary layer

*F*(

*k,*

*z*) has been applied to the pressure-transport as to the wave-induced stress. We further assume that the divergence of the turbulent transport term is smaller than the other terms in the energy balance, as in a homogeneous rough wall boundary layer.

*z*is parameterized in terms of the turbulent stress at the same height:where

*κ*is the von Kármán constant.

*z*

_{ν}of the viscous sublayer,where

*ν*

_{a}is the air viscosity. On integrating (A10)– (A12) in

*z*from the lower boundary to a 10-m height with

**u**being continuous at

*k*=

*k*

_{1}/

*δ*and

*k*=

*k*

_{m}/

*δ,*we may estimate the mean wind vector at 10 m.

In this study it is necessary to calculate the wind stress vector for a given 10-m wind speed vector. We therefore start with an initial estimate of the surface viscous stress vector *τ*_{ν} and calculate *τ*_{t}(*z*) and *τ*_{tot} using (A4) and (A5). The results are then introduced to (A10)–(A12) to estimate the 10-m wind speed vector. The calculated 10-m wind speed vector is compared with that used for the WW3 input wind vector. If they do not agree, we go back to the calculation of the total and turbulent stress with a modified estimate of the surface viscous stress vector.

In the present model, the height of the WBL is always much less than 10 m [for the highest wind speed (45 m s^{−1}), the height of the WBL is about 1.4 m]. Note that the WBL model uses a wave-following coordinate system. Our 10-m wind speed is a quantity calculated from the wind stress and the equivalent roughness *z*_{0}. Only if the significant wave height is much less the 10 m, our 10-m wind speed actually corresponds to the wind speed measured at 10 m above the mean water level. In addition, the wave boundary layer is so thin that the stability correction is not necessary inside. The WBL model basically determines the equivalent roughness length (*z*_{0}). In this study, we assume that the atmospheric surface layer is neutrally stable. However, if the atmospheric surface layer is stable or unstable, our results can be used to estimate the roughness length (or the neutral drag coefficient) and the stability correction may be made in a straightforward manner.

#### Equilibrium spectrum

In the above calculation of the wind stress vector and wind speed vector, we need the surface wave spectrum in the entire wavenumber range (*k*_{m} < *k* ≤ *k*_{1}). While the spectrum near the peak is estimated by the WW3, the tail part is parameterized using the analytical model of the equilibrium spectrum by Hara and Belcher (2002).

*θ*= 0 as the mean wind direction and consider the stress component in the wind direction only. Following Phillips (1985), it is assumed that the wind input

*S*

_{ω}to the wave action is proportional to the cube of the local wave spectrum,

*S*

_{ω}

*β*

_{g}

*k,*

*θ*

*g*

^{1/2}

*k*

^{−9/2}

*B*

*k,*

*θ*

*αgk*

^{−4}

*B*

^{3}

*k,*

*θ*

*B*(

*k,*

*θ*) =

*k*

^{4}

*ψ*(

*k,*

*θ*) is the degree of saturation and

*α*is a nondimensional proportionality constant. From (A14), the degree of saturation is expressed in terms of the growth rate as

*B*

*k,*

*θ*

*α*

^{−1}

*β*

_{g}

*k,*

*θ*

*g*

^{−1/2}

*k*

^{−1/2}

^{1/2}

*τ*

^{l}

_{t}

*k*), which can be solved asand the corresponding degree of saturation is obtained aswithHere,

*k*

_{s}is called the

*sheltering*wavenumber and represents the wavenumber at which the local turbulent stress begins to be affected by sheltering by the longer wavelength waves. At low wavenumbers (

*k*≪

*k*

_{s}),

*τ*

^{l}

_{t}

*k*) becomes constant and

*B*(

*k,*

*θ*) is proportional to

*k*

^{1/2}, while

*τ*

^{l}

_{t}

*k*) decreases like

*k*

^{−1}and

*B*(

*k,*

*θ*) becomes independent of

*k*at high wavenumbers (

*k*≫

*k*

_{s}).

In summary, the equilibrium spectrum of Hara and Belcher (2002) is determined by a single dynamical variable, *k*_{s}. For a fixed wavenumber *k,* the degree of saturation *B*(*k*) monotonically increases as the sheltering wavenumber *k*_{s} decreases. This spectrum is used to represent the tail part unresolved by the WW3 and is connected with the WW3 spectrum to construct the full wave spectrum, as described in section 5.