## Abstract

Under high-wind conditions, breaking surface waves likely play an important role in the air–sea momentum flux. A coupled wind–wave model is developed based on the assumption that in the equilibrium range of surface wave spectra the wind stress is dominated by the form drag of breaking waves. By conserving both momentum and energy in the air and also imposing the wave energy balance, coupled equations are derived governing the turbulent stress, wind speed, and the breaking-wave distribution (total breaking crest length per unit surface area as a function of wavenumber). It is assumed that smaller-scale breaking waves are sheltered from wind forcing if they are in airflow separation regions of longer breaking waves (spatial sheltering effect). Without this spatial sheltering, exact analytic solutions are obtained; with spatial sheltering asymptotic solutions for small- and large-scale breakers are derived. In both cases, the breaking-wave distribution approaches a constant value for large wavenumbers (small-scale breakers). For low wavenumbers, the breaking-wave distribution strongly depends on wind forcing. If the equilibrium range model is extended to the spectral peak, the model yields the normalized roughness length (Charnock coefficient) of growing seas, which increases with wave age and is roughly consistent with earlier laboratory observations. However, the model does not yield physical solutions beyond a critical wave age, implying that the wind input to the wave field cannot be dominated by breaking waves at all wavenumbers for developed seas (including field conditions).

## 1. Introduction

Understanding the role of breaking surface waves in air–sea exchange processes is integral to improving parameterization schemes for coupled ocean–atmosphere models, which are commonly used for weather and climate predictions (Komen et al. 1996). Although there exists evidence that wind-generated breaking surface waves enhance air–sea fluxes of heat, gases, and momentum (Melville 1996), the mechanisms that lead to this enhancement are poorly understood. Wind–wave interaction is a highly complex phenomenon. Part of the complexity arises from interactions of waves with the turbulent wind, with subsurface currents, and with other waves of multiple scales. Each of these processes not only influences the dynamical balances of waves and wind, but can also lead to different routes to breaking (Melville 1996). In addition, it is observationally very challenging to quantify wave breaking (Banner and Peregrine 1993) or turbulent properties in the lower part of the atmospheric boundary layer influenced by waves (in the “wave boundary layer”).

The goal here is to develop a simple model for strongly wind forced surface waves. The length of breaking crests per unit surface area strongly increases with wind speed, as casual observations or detailed field investigations (Melville and Matusov 2002) suggest. Furthermore, laboratory studies have shown that the form drag of waves increases significantly over breakers (Banner 1990). Therefore, we postulate that breaking waves dominate the wave form drag under high-wind conditions. Our model represents a limiting configuration with the maximal possible influence of breaking waves. Such a model shall provide a guidance to investigate the effects of breaking waves on the air–sea momentum flux and wave field.

An important concept for modeling wind waves is the “equilibrium range” of surface wave spectra, in which the wave field is nearly stationary and therefore in a dynamical balance between wind input, nonlinear interactions, and dissipation. Phillips (1985) determined the wave spectrum as function of total wind stress by assuming that all three contributions, represented by simple parameterizations, are comparable and proportional to each other. In the model of Belcher and Vassilicos (1997), sharp-crested breaking waves control the equilibrium range, which is assumed scale invariant, so that the shape of breaking waves is self-similar and the distribution of breaking waves varies as power law in the wavenumber *k*. Based on these assumptions the wave spectrum must also be governed by a power law in *k*, whose exponent was determined by imposing a wave energy balance between nonlinear interactions and dissipation.

Wave dynamics are coupled with the air dynamics, because waves extract energy and momentum from the wind and thereby reduce the wind forcing on the waves (see, e.g., Makin and Mastenbroek 1996; Makin and Kudryavtsev 1999). Within the wave boundary layer and outside the viscous sublayer, the wind stress partitions into turbulent and wave-induced components. Because shorter waves are forced by a reduced turbulent stress due to the momentum uptake of longer waves, the short wave spectrum deviates significantly from the spectrum proposed by Phillips (Hara and Belcher 2002). The mean wind profile in the wave boundary layer can be determined from turbulence closure schemes, based on the energy balances of mean, wave-induced, and turbulent motions. Makin and Mastenbroek (1996) employ the full balance equations of turbulent kinetic energy and its dissipation to calculate drag coefficients over fully developed seas. Balancing locally turbulent production and dissipation, Makin and Kudryavtsev (1999) estimate an eddy viscosity, which is also used to parameterize the turbulent dissipation rate. Such a closure model, however, generally does not conserve the total energy within the wave boundary layer (Hara and Belcher 2004). To satisfy energy conservation, Hara and Belcher (2004) instead relate the turbulent dissipation rate to the reduced turbulent stress in the wave boundary layer.

In most of the previous wind–wave models, only the wind input to nonbreaking waves has been considered (see, e.g., discussion of Komen et al. 1996). A notable exception is the investigation by Kudryavtsev and Makin (2001), who estimated that a considerable fraction (up to 50%) of the wind stress is supported by breaking waves. The total wind stress at the surface has contributions from the viscous stress and the wave form drag, which can be split into two parts due to breaking and nonbreaking waves. In the study by Kudryavtsev and Makin (2001), the parameterization of the momentum flux from wind to nonbreaking waves depends on the total wind stress, rather than on a reduced turbulent stress. Such a parameterization, however, generally violates momentum conservation (see, e.g., Hara and Belcher 2002). Also, the dependency on the reduced turbulent stress is crucial to accurately predict the stress partitioning between viscous stress and the form drag of nonbreaking waves (e.g., Kukulka and Hara 2005).

The problem is significantly simplified if the momentum flux to nonbreaking waves is neglected. We will closely follow the approach by Kudryavtsev and Makin (2001) to parameterize the form drag due to breaking waves. Unlike Kudryavtsev and Makin (2001), however, we will not prescribe the mean wind by a logarithmic profile, which generally does not satisfy energy conservation in the wave boundary layer. Instead, we will derive coupled equations governing the wind speed, turbulent stress, and the wave-breaking distribution, based on the conservation of wave energy, as well as air-side momentum and energy (section 2). Furthermore, in our approach smaller-scale breaking waves are sheltered from wind forcing, if they are in airflow separation regions of longer breaking waves (spatial sheltering effect). After deriving solutions without the spatial sheltering effect (section 3), we will investigate how spatial sheltering affects the dynamics (section 4). Before concluding, we will apply our model to describe strongly forced, developing wave fields (section 5).

## 2. Theory

To investigate the effects of breaking waves under high-wind conditions, we introduce a hypothetical “high-wind equilibrium range” of the surface wavenumber spectrum, in which the wind input is dominated by breaking-wave dynamics for wavenumbers between *k*_{0} and *k*_{1}. For mathematical convenience we let *k*_{1} → ∞, as high *k* asymptotes are derived in the next sections. Physical solutions then need to be truncated at a finite *k*_{1}, where the effects of surface tension or viscosity become important.

Airflow separation ahead of the breaking-wave crests leads to a pressure drop, causing the transfer of momentum from wind to waves. In sections 2a–c, we will develop conservation equations for the air-side momentum and energy, as well as the wave energy. Throughout the development, we assume 1) the breaking-wave crests are orthogonal to the wind direction, 2) the shape of breaking waves is self-similar, and 3) the mean wind flow *u* in the horizontal *x* direction is stationary and horizontally homogeneous. The first assumption is not critical for deriving the main conclusions of this paper (see appendix B). Laboratory studies of short gravity wind waves support the validity of the second assumption (e.g., Caulliez 2002). The third assumption is a common simplification for idealized wind–wave models and is based on a local approximation. For simplicity, we will at first neglect the effect of airflow separation due to longer breaking waves on shorter waves. This “spatial sheltering effect” will later be introduced in section 2d.

### a. Air-side momentum

The velocity vector (𝒰, 𝒲) is decomposed into mean and turbulent fluctuation (denoted by primes):

The momentum balance of a volume of air at height between *z* − *dz*/2 and *z* + *dz*/2 (see Fig. 1) may be expressed as

that is, the divergence of the turbulent momentum flux *τ _{t}*(

*z*) = −

*ρ*is balanced by the resistance force per unit volume

_{a}*f*due to breaking waves. This resistance force is a result of airflow separation ahead of the breaking-wave crest, leading to a pressure drop at the forward face of the breaking-wave front. Our approach resembles the urban canopy model by Coceal and Belcher (2004), in which individual canopy elements give rise to a localized drag.

_{b}Let us first consider a monochromatic wave field with a wavenumber *k*. The breaking waves at wavenumber *k* are obstacles to the wind with a frontal area 2*a*(*k*)*l*(*k*), where *a* is the amplitude of the breaking waves, and *l* is the total length of breaking-wave crests per unit surface area. Because breaking waves are assumed self-similar, the slope of breaking waves *ε* = *ak*, is fixed, that is, *a*(*k*) = *ε*/*k*. The wind force *τ _{b}* (per unit surface area) that acts on the breaking-wave crests is

where Δ*p* denotes the pressure drop due to airflow separation ahead of the wave crest. Scaling arguments as well as previous experiments suggest that the pressure drop can be parameterized by

where *c* is the phase speed, *u*(*k*) = *u*[*z* = *a*(*k*)] denotes the wind speed evaluated at the height of the breaking waves (see Fig. 1), and *c _{d}* is the drag coefficient with

*c*= 0 for

_{d}*u*(

*k*) <

*c*.

For a spectrum of waves, the average length of breaking-wave crests per unit surface area for waves with wavenumbers between *k* − *dk*/2 and *k* + *dk*/2 is given by Λ(*k*)*dk*, where Λ is the distribution function of breaking-wave crests (“breaking-wave distribution”; Phillips 1985). Then the momentum flux into wavenumbers between *k* − *dk*/2 and *k* + *dk*/2, denoted by *dτ _{b}*(

*k*), is

Most of *dτ _{b}*(

*k*) is transferred to the wave crest at

*z*=

*a*(

*k*) =

*ε*/

*k*, since the pressure drop is proportional to the squared wind speed, which rapidly decreases toward the surface. Therefore,

*dτ*(

_{b}*k*) is approximately the momentum flux

*dτ*(

_{b}*z*) at heights between

*z*−

*dz*/2 and

*z*+

*dz*/2, that is,

where the height *z* is uniquely related to the wavenumber *k* through

Then, the momentum equation in (2) can be rewritten

or

if we introduce a turbulent stress *τ _{t}*(

*k*) =

*τ*[

_{t}*z*=

*a*(

*k*)] as a function of

*k*.

where *β*_{Λ} = 2*c _{d}ε*.

### b. Air-side energy

The Reynolds-averaged total energy equation of the airflow over breaking waves includes an additional term representing energy loss to breaking waves:

The first three terms are familiar (see, e.g., Cohen and Kundu 2002), with *uτ _{t}* describing the energy flux due to the turbulent stress acting on the mean flow, Π = [ +

*ρ*(1/2)], the mean energy flux associated with turbulent fluctuations, and the third term ɛ denoting the viscous dissipation rate of turbulent eddies. Here, we neglect the divergence of the flux Π, similar to traditional turbulence models close to a solid wall. The simulation of airflow over nonbreaking waves by Makin and Mastenbroek (1996) suggests that the divergence of Π is not the dominant term in the balance equation of turbulent kinetic energy.

_{a}The loss of wind energy to breaking waves is to first order given by *cf _{b}*, since the wind force

*f*mostly acts on breaking-wave crests, where the water parcels move approximately with the phase speed

_{b}*c.*Therefore, the power per unit volume transmitted by the wind to the breaking waves is

*cf*. Following Hara and Belcher (2004), ɛ is related to the turbulent stress as

_{b}where *κ* = 0.4 is the von Kármán constant.

As we have done for the momentum equation, we rewrite the energy equation in *k* using the same relationship *z* = *ε*/*k*, yielding

or

is the spectral energy flux density into breaking waves at *k*.

### c. Wave energy

The pressure asymmetry on the waves leads to energy input to the waves. For a stationary and homogeneous wave field, the energy input rate *I _{b}*(

*k*) from wind to waves balances the rate of wave energy loss, and the redistribution of energy due to nonlinear interactions, the sum of the latter two terms is denoted by

*L*(

*k*), that is

Frequently, the parameterization of the nonlinear interaction term is based on weak resonant interactions of four gravity waves (Komen et al. 1996, sections II.3 and III.3.2). In the regime, however, where breaking-wave dynamics dominates the wind input, this weak nonlinear interaction mechanism might be secondary. Therefore, in the equilibrium range we assume that the loss term in *L* due to wave breaking is dominant. Following Phillips (1985), the energy loss due to wave breaking is related to the breaking distribution:

where *b* denotes a coefficient that is independent of *k*. This parameterization is also consistent with the study by Thorpe (1993), who investigated the energy loss due to breaking waves as a function of wind speed. Note that *b* could be adjusted to account for the possible imbalance that results from neglecting the nonlinear interaction term, if the dissipation and nonlinear interaction terms were approximately proportional to each other. Substitution for *I _{b}* and

*L*in (16) yields the wave energy balance

Notice that (18) does not depend explicitly on Λ, but is determined by the coupled system of equations, summarized at the end of this section. To understand physically the complex coupling between the air-side and wave dynamics, a time-dependent problem will be investigated in section 3c.

### d. Spatial sheltering

Ahead of breaking-wave crests the airflow separates. If a shorter breaking wave is just ahead of a longer breaking wave and therefore is inside the airflow separation region (or short “sheltered area”) due to the longer breaking wave, the wind forcing on the shorter breaking wave must be reduced (“spatial sheltering effect”). We include this effect in our model by simply setting the wind forcing on a shorter breaking wave to zero if it is located in a sheltered area of a longer breaking wave. The vertical cross section of an airflow separation region is modeled as a rectangular area, which extends with length *D* ahead of the breaking-wave crests (see Fig. 2). Since the breaking waves are assumed self-similar, *D* scales as *D* = *ε̂k*^{−1}, where ε̂ is independent of *k* and is called the spatial sheltering coefficient. With a spectrum of breaking waves we define the fractional area (area per unit surface area) free of airflow separation *α*(*z*). At the top of the wave boundary layer, where the airflow is not influenced by waves, there exists no sheltered area, so that *α* = 1. Below, *α*(*z*) monotonically decreases as *z* decreases. Let us introduce *α*(*k*) = *α*(*z* = *ε*/*k*), which is the fractional area free of airflow separation due to breaking waves with wavenumbers less than *k*. This function monotonically decreases with increasing *k*.

A simple way of modeling *α*(*k*) is to assume that breaking-wave crests with a wavenumber *k* generate new sheltered area, as long as the new area does not overlap with already existing sheltered area (overlapping area in medium gray in Fig. 2) due to longer breaking waves. Then the expected value that the area *D*Λ*dk* falls outside already existing airflow separation regions is given by *α*(*k*)*D*Λ*dk*. Hence, breaking waves at *k* change the fractional area free of sheltered area by

Or, with *D* = *ε̂k*^{−1}, the last equation becomes

Note that based on this *α* model, breaking crests could erroneously create new sheltered area, if they break within sheltered area (hatched area in Fig. 2 is included as sheltered area). This, however, likely is only a small source of error and is neglected in the following calculation.

The air-side dynamics discussed in previous subsections holds only in regions that are free of airflow separation. Therefore, the force on breaking waves, the viscous dissipation of turbulent kinetic energy, and the momentum flux due to turbulent eddies must all be reduced by the fraction *α*. Hence, with the spatial sheltering effect included, the air-side momentum and energy, and wave energy equations are, respectively,

and

where *u* is the wind outside sheltered area; that is, the volume flux due to the wind is *αu*.

### e. Set of governing equations

In this section we summarize the governing equations for *u*, *τ _{t}*, Λ, as well as

*α*, and present them in a convenient form for further analysis. After substitution of (

*u*−

*c*)

^{2}from (18) in (10) and (

*u*−

*c*)

^{2}Λ from (10) in (14), the set of governing equations with spatial sheltering effect, that is, the air-side momentum and energy, the wave energy, and the spatial sheltering equations, can be written respectively as

and

where

denotes a coefficient that is related to the ratio of energy input to dissipation.

We will reduce the four governing equations to two equations for *ατ _{t}* and

*u*only. First the expression for

*α*from (26) is substituted into (27) to obtain Λ as function of

*u*and

*k*only:

The phase speed is estimated by *c* = *n*(*g*/*k*)^{1/2}, where *g* is the gravitational constant. The enhanced phase speed of nonlinear waves is modeled by assuming that the nonlinear phase speed is augmented by a fixed fraction, *n*. Note that *n* = 1 yields the linear phase speed of surface gravity waves.

where

is proportional to the spatial sheltering coefficient. Replacing *du*/*dk*, *α*, and Λ with expressions from (30), (26), and (29), respectively, in the air-side energy equation in (25) yields the second equation that determines *ατ _{t}* and

*u*, namely,

The governing equations in (30) and (32) are more conveniently expressed after introducing the nondimensional variables

The dependent variables *S* and *U* can be interpreted as normalized wind forcing and wind speed, respectively. Note that *U* has the following properties: 0 < *U* < 1 for forced waves, *U* = const. for *u* ∝ *c*, and *U* ≈ 1 for *u* ≫ *c*. The new independent variable can be expressed as *K* = ln(*k*/*k*_{ref}), where *k*_{ref} is some reference wavenumber. For convenience, we set *k*_{ref} = *k*_{0}, so that the solutions coincide for fixed boundary conditions *S*(*k*_{0}) = *S*_{0} and *U*(*k*_{0}) = *U*_{0}. Notice that the governing equations are independent of the phase speed factor *n*. Once *S* is known, the breaking distribution can be determined through (24), which can be expressed with *S* as

where *b*′ = *bρ _{w}*/

*ρ*. The solutions to the nonlinear autonomous system (36) and (37) depend only on the boundary conditions, as well as on the parameters

_{a}*γ*and

*δ*.

Before deriving solutions, it is useful to provide an approximate range of the experimental coefficients *ε*, *b* (or *b*′), *c _{d}*, and

*ε̂*, as well as the model coefficients

*γ*and

*δ*. The slope of a breaking wave depends on the breaking mechanism, but is generally confined between

*ε*= 0.1 and 0.5 (see discussion by Melville 1996). The wave energy dissipation coefficient has been determined in the range of

*b*= 0.003–0.07 for breaking laboratory waves, with lower values found for unsteady breakers (Melville 1996). The study from Phillips et al. (2001) suggests that

*b*could be as small as

*b*= (0.7–1.3) × 10

^{−3}for open-ocean conditions, but we use the range suggested by Melville since it is based on direct observations of the dissipation rate. We adopt

*c*≈ 0.5 from the study of Kudryavtsev and Makin (2001), which is based on previous laboratory studies of reattaching flow over a backward facing step. An upper bound of

_{d}*ε̂*can be estimated by assuming that

*D*(

*k*) cannot exceed one wavelength, so that

*ε̂*

_{max}< 2

*π*. Laboratory investigations by Reul et al. (1999) suggest that the recirculation cell extends approximately from wave crest to trough, so that

*ε̂*≈

*π*. To confine

*ε̂*, we will arbitrarily set

*ε̂*

_{min}=

*π*/4.

From the range of experimental coefficients, one can estimate a rough upper and lower bound for *γ*, *δ*, and *b*′, as *γ* = 0.04–0.5, *δ* = 0.01–3, and *b*′ = 2–60. Note that *γ* = 0.27 for the values of the parameters used in the study by Kudryavtsev and Makin (2001). In the later presentation we will assign by default the following reference values, *ε* = 0.3, *b* = 0.01, *ε̂* = 1, corresponding to *γ* = 0.2, *b*′ = 8, *δ* = 0.1. We will also carefully investigate the dependency of the solutions on the model coefficients, *γ* and *δ*. Although *γ* and *δ* depend on *b* and cannot be chosen independently, we will assume for mathematical convenience that these model coefficients are independent of one another, which effectively increases the model uncertainty.

## 3. Model without spatial sheltering

It is insightful to first analyze the behavior of *u*, Λ, and *τ _{t}* without spatial sheltering, that is,

*δ*= 0. In this case, the governing equations in (36) and (37) simplify to

For this system, exact analytic solutions can be easily obtained.

### a. Solutions

From (35) with *α* = 1, the nondimensional wind speed is constant at

with the solution

where *A* is a coefficient of integration that is determined by the boundary condition. The breaking distribution is obtained from (38):

In dimensional variables and with the local friction velocity defined as

the solution can be rewritten as

and

The dimensional wind speed is proportional to the phase speed

hence, the wind profile is proportional to *z*^{1/2} rather than logarithmic.

### b. Discussion

Because of the simple wave energy balance between input and dissipation, with both terms being proportional to Λ, the wind speed must be proportional to *c*. According to (48), *u* always exceeds *c* by a factor [1 + (1/*γ*)], so that all waves in the equilibrium range are forced by the wind. The coefficient *γ* (ratio of input to dissipation) can be interpreted as the “efficiency” of the wind speed in maintaining the wave energy balance. Note that as a consequence of the wave energy conservation (18), the pressure drop must adjust to balance the dissipation, which is proportional to *b*. For larger *b*, the pressure drop could only be elevated by a larger input coefficient, *β*_{Λ}, or an increase in *u*. Therefore, for *γ* ≫ 1 (large *β*_{Λ}, small *b*), a relatively small wind close to *c* can maintain the wave energy balance. On the other hand, if the input is weak relative to dissipation, *γ* → 0, *u* must be very large (*u* → ∞), so that sufficient energy is fluxed into waves.

According to (43) and (44), *S* and Λ monotonically increase (decrease) for *A* > 0 (*A* < 0) from the boundary values at *k* = *k*_{0}, where *S* = *S*_{0} and

to the asymptotic values for *k* → ∞,

and

For *A* = 0, that is, *S*_{0} = *S*_{∞}, *S* and Λ remain constant (see Fig. 3). Note that if *S* is constant, the turbulent stress decays as 1/*k*. If *S*_{0} < *S*_{∞}, Λ_{0} must be larger than Λ_{∞}, according to (49). The elevated Λ increases the momentum flux into waves, so that the turbulent stress in (24) decreases more strongly relative to the case of *S*_{0} = *S*_{∞}. Since *S* is proportional to the turbulent stress, *S* must decrease until the solution converges to *S*_{∞}. The opposite is true for *S*_{0} < *S*_{∞}, where *S* must increase to approach *S*_{∞}.

For *k* → ∞, both *S*_{∞} and Λ_{∞} decrease with increasing *γ* for fixed *b*′ (black solid line in Fig. 4). The constant Λ at high *k* agrees with the model by Belcher and Vassilicos (1997). The power-law distribution (*k*^{0}) infers scale invariance, which is one of the basic assumptions made by Belcher and Vassilicos (1997). The zero exponent found in their study, however, follows from a dynamical balance between wave dissipation rate and the divergence of nonlinear wave interaction flux. Interestingly, the model by Hara and Belcher (2002), which is based on a three-way wave energy balance between dissipation, nonlinear interactions, and wind input to nonbreaking waves only, also predicts a constant Λ for high *k*.

According to (49), *S*_{0} needs to exceed a threshold, so that Λ cannot take unphysical negative values. One can show that Λ is positive for any *A* < 0 or *A* > 0 and

Recall that

so that *S*_{0} is the normalized wind forcing applied at the height *z* = *ε*/*k*_{0}. The inequality (52) suggests that the applied wind forcing, *S*_{0}, must exceed a threshold for admissible solutions for a given *γ* (see Fig. 5). This is consistent with the assumption that our model is only applicable for strongly wind forced surface waves. Notice also that *S*^{1/2}_{0} defines an inverse wave age. Therefore, our theory is only applicable for “young,” that is, strongly forced wind waves. This will be discussed in more detail in section 5.

### c. Time-dependent problem

To gain better physical insight, we will next examine how the solutions Λ, *u*, and *τ _{t}* respond in time if Λ is slightly perturbed. Let us denote the perturbed breaking-wave distribution by Λ′ and the resulting perturbed wind speed and turbulent stress by

*u*′ and

*τ′*, respectively. If the adjustment time scale of the wind is much shorter than the adjustment time scale of gravity waves, the air-side momentum and energy balances adjust more rapidly in comparison with the wave energy balance. Hence, the time derivative terms can be neglected in the air-side momentum and energy equations, so that (10) and (14) approximately hold for the perturbed wind. The wave energy, on the other hand, changes with time; that is, (18) becomes

_{t}where *E* denotes the spectral wave energy density. If the magnitude of the first term on the right-hand side exceeds the magnitude of the second term, that is, if *u*′(*k*) > *u*(*k*), the wave field will grow (∂*E*/∂*t* > 0) and the breaking distribution will increase (∂Λ/∂*t* > 0), provided Λ monotonically increases with *E*. Therefore, the solutions are stable if the wind speed decreases (increases), that is, *u*′ < *u* (*u*′ > *u*), in response to increased (decreased) breaking distribution, that is, Λ′ > Λ (Λ′ < Λ) for all *k*.

In the following analysis, the wavenumber of the longest wave, *k*_{0}, and the boundary value of the turbulent stress, *τ _{t}*(

*k*=

*k*

_{0}), are held constant. Then, the perturbation in Λ and the resulting wind perturbation can be simulated by considering a second solution of the steady problem, that is, (10), (14), and (18), with a slightly different value for

*b*, namely,

*b*+ Δ

*b*and the same boundary condition

*S*

_{0}. This is possible because the parameter

*b*appears only in the wave equation but not in the air-side equations in (10) and (14). The second solution therefore satisfies the air-side momentum and energy equations exactly but does not satisfy the wave energy equation with the original choice of

*b*.

One can show that for Δ*b* > 0, Λ′(*k*) < Λ(*k*) and *u*′ > *u*, while for Δ*b* < 0, Λ′(*k*) > Λ(*k*) and *u*′ < *u* for all *k*. Therefore, in light of (53), the perturbed Λ′ will always go back to Λ, so that the solution is stable. Furthermore, one can show that for Λ′ < Λ, the turbulent stress is enhanced, that is, *τ′ _{t}* >

*τ*, below the boundary height

_{t}*z*=

*ε*/

*k*

_{0}. However, it is not straightforward to physically explain why the wind speed and wind stress increase in response to the reduced breaking distribution, because the coupling between the wave field and the wind is very complex.

## 4. Model with spatial sheltering

### a. Solutions

As before, we will set the boundary conditions at *k*_{0}. In general, breaking waves with *k* < *k*_{0} will introduce some sheltered area, that is, *α*(*k*_{0}) < 1. However, we only need to focus on solutions with *α*(*k*_{0}) = 1, since the family of solutions obtained here also contains solutions for *α*(*k*_{0}) < 1, if the lower boundary of *k* is redefined as *k*′_{0} with *k*′_{0} > *k*_{0}. Before solving numerically the full system (36) and (37), we will examine the behavior of the solutions for *k* close to *k*_{0} and for large *k* (*k* → ∞).

#### 1) Asymptotes for *k* close to *k*_{0} with *α*(*k*_{0}) = 1

For *k* close to *k*_{0}, only a small amount of sheltered area covers the surface area, that is, *α* ≈ 1, so that *U* can be approximated by

Furthermore, we may linearize (37) in the small perturbation *S*′ = *S* − *S*_{0}:

where the constant coefficients *q*_{1} and *q*_{2} are defined by

The asymptotic solution for *S* is

With *S*, Λ can be obtained from (38). Note that there are now two conditions required for Λ > 0, namely, either

or

The former condition (58) suggests that the normalized wind forcing must be bound not only by the lower limit found in the nonsheltering case, but also by an upper limit. This will be discussed in more detail in section 5. The latter condition (59) generally does not yield a physical solution, which can be shown by a critical point analysis as discussed in appendix A and the next section.

#### 2) Asymptotes for *k* → ∞

To determine the solution for large *k* we calculate the critical points (*U*_{∞}, *S*_{∞}) of the system, where the solution is stationary. For *dU*/*dK* = *dS*/*dK* = 0, (36) and (37) become

We reject *S*_{∞} = 0 as unphysical solution. Then *U*_{∞} = 1 and

Notice that we chose the positive sign in front of the square root, such that *S*^{1/2}_{∞} > 0. By means of critical point analysis (see appendix A), one can identify this solution as a stable node to which physical solutions converge, as *k* → ∞. Last, the asymptotic physical variables can be summarized as follows:

#### 3) Numeric solutions

The system (36) and (37) can be integrated with a standard numerical ordinary equation solver. An example of the family of solutions for a particular boundary condition, a fixed *γ*, and several values of *δ* is shown in Figs. 6 and 7. Note that asymptotic solutions coincide with the numerical ones in the valid range. The peak in Λ, however, could only be revealed through the full numeric solutions.

### b. Discussion

Before comparing the solutions with previous findings, we will here highlight the differences between sheltering and nonsheltering models. First, we will investigate if the two models are consistent, in the sense that they converge to the same solutions for *δ* → 0 and *α*(*k*_{0}) = 1. In general, the two models are not consistent. For simplicity, we will only compare the two boundary values Λ_{0} and Λ_{∞} = *b*′^{−1}*S*_{∞}. In virtue of (38), it is easy to show that Λ_{0} with sheltering converges to the value of Λ_{0} without sheltering as *δ* → 0. For *k* → ∞, on the other hand, *S*_{∞} without sheltering, denoted by *S*_{∞ns}, is given by (50), while one can show [e.g., with (62)] that with spatial sheltering

Therefore, *S*_{∞} converges to *S*_{∞ns} as *δ* → 0 only if *γ* → 0. This suggests that both models are only consistent for *γ* → 0. Notice that this condition implies *u* ≫ *c* for the nonsheltering model, which is required for the sheltering asymptote *U*_{∞} = 1. Practically, as Fig. 4 indicates, for *γ* < 1, the sheltering asymptote Λ_{∞} closely approaches the nonsheltering value as *δ* → 0.

According to (38), an enhanced *δ* increases Λ_{0} for fixed *S*_{0} (see also Fig. 7). This increase might be surprising, because at *k* = *k*_{0} there is no sheltered area, since we set *α*(*k*_{0}) = 1. Notice, however, that spatial sheltering modifies the divergence of flux terms and therefore plays a role in the dynamics of the system (24)–(27), even if *α* = 1.

At the high *k* limit, the “effective” airspace is reduced by spatial sheltering, which decreases Λ_{∞} relative to the solution without spatial sheltering (see Figs. 4 and 7). The reduction of the effective airspace causes increases in *u* to maintain the wave energy balance (26) (Fig. 6). One can show that *u* decreases monotonically toward the surface only for *δ* < (*γ*/*κ*)^{2} (see also Fig. 6). However, the average wind speed, that is, the volume flux *αu*, always decreases toward the surface. Therefore, an increase in *u* toward the surface could occur only locally.

An interesting feature of the sheltering solution is the peak of Λ (Fig. 7). The peak is due to a transition from nonsheltering to sheltering dynamics. For low *k*, the sheltering effect is first weak, and therefore Λ is closely governed by the asymptotic solution with *α* ≈ 1, discussed above. As *k* increases, more and more area is covered by sheltered area. After a critical transition point the wind speed follows closely the asymptotic solution for high *k *(68).

## 5. Developing young wave fields

Our model has been formulated such that it is valid only for the equilibrium range, in which the input from wind to waves is dominated by the form drag of breaking waves and is balanced by the dissipation of breaking waves. In principle, to model wind and waves accurately, wind interactions with the whole spectrum of waves, including the spectral peak, need to be considered. The application of our equilibrium range model requires knowledge of the boundary condition

For a given total stress, the local friction velocity *u _{l}*

_{*}(

*k*

_{0}) can be determined by integrating the momentum fluxes into all waves with

*k*<

*k*

_{0}, which includes fluxes to the spectral peak. The dynamics of the spectral peak may differ from the dynamics of the equilibrium range discussed in this paper, because nonlinear wave interactions likely play a role in shaping the spectral peak (see, e.g., Komen et al. 1996).

Nevertheless, as a first step to investigate developing young wave-breaking fields, we will simply extend the equilibrium model to the spectral peak and set the spectrum to zero for waves with *k* < *k _{p}* =

*k*

_{0}, where

*k*is the peak wavenumber. This simplification is based on the assumption that shorter wind waves support a dominant fraction of the overall wave-induced momentum flux and that the detailed shape of the breaking distribution near the spectral peak is not important. Consequently, the wave age,

_{p}*σ*, is expressed as

where *u*_{*} = *τ*_{0}/*ρ _{a}* denotes the total friction velocity and

*τ*

_{0}is the total wind stress.

Since *S*_{0} must be bound by (58), our model applies only for a limited range of wave ages between *σ*_{min} = *δ*/2 and *σ*_{max} = 2*γ*/*κ*(1 + *γ*). A possible interpretation of the upper limit for *S* (lower limit for *σ*) is that our model cannot be applied for very young initial wavelets generated just after the onset of the wind. For *δ* < 0.1, most of the observed wave conditions have wave ages exceeding the lower bound (see also discussion below).

Figures 8 and 9 illustrates how the wave field develops according to our theory for a fixed total stress. As the wave field develops, the peak phase speed increases. Since faster waves feel a reduced wind forcing (i.e., smaller *S*_{0} and larger *σ*), Λ_{0} decreases. The wave age *σ*_{max} defines an upper limit, where Λ_{0} approaches zero. The value of *σ*_{max} is estimated to be 1.7 or less. Therefore, our theory suggests that the dynamics of more developed seas (*σ* > *σ*_{max}) cannot be dominated by breaking waves, at least not at smaller wavenumbers close to *k*_{0}. Likely, nonbreaking waves play a significant role in the momentum transfer from wind to waves for older seas. A simple physical interpretation of this result could be that longer waves of older seas propagate too fast, so that the wind speed is inefficient in forcing breaking waves.

For younger wave fields, the enhanced Λ causes a more rapid coverage with airflow separation regions (Fig. 8, left bottom panel). The area is quickly covered by separation regions as the surface is approached (i.e., for increasing *k*). Notice that the high-*k* asymptote Λ_{∞} is smaller with spatial sheltering than without (Fig. 9). Furthermore, the breaking distribution with spatial sheltering differs from the nonsheltering solution in that the response of longer waves to the wind forcing is more sensitive as the wave age increases.

It would be interesting to compare the breaking distribution determined from our simple model with previous observations. Although limited observations exist for open-ocean conditions (e.g., Gemmrich 2005; Melville and Matusov 2002; Phillips et al. 2001), to date there are no quantitative data for young, strongly forced wind waves.

To demonstrate a potential application of our model, we will predict the total length of breaking-wave crests per unit surface area as well as the breaking rate for a realistic high-wind-condition laboratory study. The experiment from Donelan et al. (2004), who measured the drag coefficient in a wind–wave tank, was conducted for *u*_{*} up to 2.5 m s^{−1} with *σ* ≈ 0.5. As the wave field develops from wave ages of *σ* = 0.3, 0.5, 0.8 (see Figs. 8 and 9), we find that the length of total breaking-wave fronts close to the spectral peak decreases from (44, 9, 0.6) m m^{−2}, respectively. For this calculation, Λ has been integrated from *k _{p}* = 22, 10, 4 m

^{−1}(for

*σ*= 0.3, 0.5, 0.8, respectively) to 3

*k*and

_{p}*u*

_{*}is set to

*u*

_{*}= 2 m s

^{−1}. The breaking rate is here defined by the average number of breaking crests with wavenumbers between

*k*and 3

_{p}*k*that pass through a vertical line per unit time and can be calculated by integrating

_{p}*c*Λ over this wavenumber range. The breaking rate is (23, 6, 0.6) s

^{−1}for

*σ*= (0.3, 0.5, 0.8), respectively. These estimates could be an order of magnitude smaller or larger depending on the particular choice of parameters. Interestingly, for young, strongly forced wind waves a small change in wave age has a large impact on the number of breaking waves. For such a high coverage of breakers, the water just below the air–sea interface likely is well mixed and turbulent as a result of the overturning of breaking waves, which has important implications for air–sea gas transfer (Melville 1996). Enhanced turbulence and kinetic energy dissipation under breaking waves have been previously observed (see, e.g., discussion by Thorpe 2005, chapter 9.2).

One important observable, which can be compared with our model results, is the nondimensional roughness length (or Charnock coefficient),

where *z*_{0} is the roughness length. The Charnock coefficient is often used to parameterize the air–sea momentum flux and has been determined for a wide range of wind and wave conditions. A simple expression for *r* can be derived here, by assuming that the total stress on the sea surface is solely determined by the form drag of breaking gravity waves. The turbulent wind profile for a neutrally stable atmosphere above the wave boundary layer, namely,

makes a continuous transition to the wind profile of the wave boundary layer, so that

where z* _{T}* =

*ε*/

*k*

_{0}denotes the height of the wave boundary layer. Solving for

*r*yields

Interestingly, provided that *α*(*k*_{0}) = 1, spatial sheltering has no net effect on the Charnock coefficient. Since the total energy flux into the wave boundary layer, *u*(*z _{T}*)

*τ*

_{0}, is the same with and without sheltering, the total air energy losses to waves and viscous dissipation must be the same also. Therefore, spatial sheltering changes the air dynamics only internally inside the wave boundary layer. Note that the enhanced phase speed of nonlinear gravity waves decreases the Charnock coefficient.

As Fig. 10 indicates, our *r*(*σ*) roughly agrees with previous laboratory observations of very young wave fields. The relation from Toba et al. (1990) results from dimensional analysis and a best fit to composite data from laboratory and field studies, including storms with wind speeds up to 25 m s^{−1}. The composite data from Toba and Ebuchi (1991) also stem from field and laboratory experiments for *σ* > 1, and only from laboratory experiments for *σ* < 1. As the wind speed increases from 30 to 50 m s^{−1}, Donelan et al. (2004) found that the dimensional roughness length approaches a constant value of 3.35 mm. For completeness, we also show data from the recent field experiment by Drennan et al. (2003) for more developed seas.

From (74) we can also derive a simple expression relating the roughness length to the amplitude of the longest wave,

Therefore, according to our theory, the roughness length is a fraction of the peak wave amplitude between exp(−2) ≈ 0.14 and one with *n* = 1. The expression (75) is qualitatively consistent with previous expressions that relate the roughness length to the significant wave height (see, e.g., review in Jones and Toba 2001).

Since the complete wind speed profile is determined with the assumptions discussed above, we can examine what happens if the steady solution is perturbed by adding a single wave component at _{i}_{0} that is slightly below the lowest wavenumber *k*_{0}, that is, _{i}_{0} < *k*_{0} and _{i}_{0} → *k*_{0}. As long as the disturbance is small enough, the wind speed at _{i}_{0} does not change significantly, so that

If the energy input to the wave exceeds the dissipation term in (18) the disturbance may grow (see discussion in section 3c). The temporal change in wave energy is given by

After introducing (72) into the last equation and linearizing in (_{i}_{0} − *k*_{0}), one finds that

with *σ̃* = *c*(_{i}_{0})/*u*_{*}. Therefore, disturbances with wave ages smaller than *σ*_{max} will grow and the wave field may develop. Disturbances with wave ages greater than *σ*_{max} on the other hand will decay, so that the wave field cannot develop beyond the wave age of *σ*_{max}.

## 6. Conclusions

We have developed a coupled wind and wave model for strongly forced wind waves, based on the assumption that in the equilibrium range of surface wave spectra the wind stress is dominated by the form drag of breaking waves. On the foundation of wave energy conservation, as well as conservation of momentum and energy in the air, expressions for the turbulent stress, wind speed, and the distribution of breaking-wave crests have been derived and their solutions have been obtained. We have also considered the spatial sheltering of shorter waves, which cannot be forced in airflow separation regions of longer waves.

Next, the equilibrium range model has been extended to the spectral peak and growing (fetch dependent) breaking-wave fields have been examined. Our model results indicate that the breaking-wave distribution near the spectral peak decreases rapidly as the wave field develops, but they approach a constant at high wavenumbers, being consistent with earlier theoretical studies. The model also yields a simple analytic formula for the Charnock coefficient, which roughly agrees with previous laboratory observations for young wave fields. For intermediate to fully developed seas, our model does not yield physical solutions. Therefore, for older seas the wind input to the wave field cannot be dominated by breaking waves at all wavenumbers; likely, the input to nonbreaking waves also needs to be considered.

As a next step, we plan to construct a model that includes the wind input to both breaking and nonbreaking waves. Starting from the previous approaches by Hara and Belcher (2002, 2004) that include wind input to nonbreaking waves only, we could add the breaking-wave effect, which is parameterized as in this study. In this case, the present model represents an important guidance as a limiting configuration, revealing the maximum possible influence of breaking waves.

The present model depends on a set of parameters that are not well constrained. Therefore, the estimated numerical values of the breaking-wave distribution and the drag coefficient include large uncertainties. It is critically important to constrain these parameters using a more complete set of observations including the direct estimates of the breaking-wave distribution. Potentially, the problem needs to be revisited to include the effects of sea sprays (Andreas 2004) and resulting ill-defined air–sea interfaces under high-wind conditions.

## Acknowledgments

We thank two anonymous reviewers for carefully reading our manuscript and providing helpful comments. This work was supported by the U.S. National Science Foundation (Grant OCE-0526177) and the U.S. Office of Naval Research (Grant N00014-06-10729). Author TK was in part supported by a University of Rhode Island Graduate Fellowship.

## REFERENCES

**,**

**,**

**,**

**,**

**,**

**,**

**.**

**.**

**,**

**,**

**,**

**.**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

### APPENDIX A

#### Critical Point Analysis

To examine the nature of the critical point (*U*_{∞}, *S*_{∞}) found in section 4a(2), we introduce the small perturbations *U* ′ and *S*′, so that (*U*, *S*) = (*U*_{∞} + *U* ′, *S*_{∞} + *S*′). Substituting this into (36) and (37) and linearizing results in

The first equation determines *U* ′. Because (*δ*/2)*S*_{∞} > 0, *U* ′ converges to zero for *K* → ∞. The solution to *S*′, determined from the second equation, can be decomposed into homogeneous and inhomogeneous parts. Since *U* ′ converges to zero, we only need to consider the homogeneous part, governed by

Because the second bracket is always positive, as one can show as follows:

*S*′ converges only to zero if the first bracket is also positive, that is, for

Now the condition (58) leads to the inequality

### APPENDIX B

#### Directional Effects

In this section we will investigate how breaking waves traveling in oblique directions relative to the wind modify the dynamics. Let the two-dimensional breaking distribution be

where *θ* denotes the angle of propagation relative to the wind direction and *h* is the directionality of the breaking distribution with

We assume that the angle between the wind and the propagation direction of the waves does not exceed 90°. For simplicity we will focus on the model without spatial sheltering and, in the case of spatial sheltering, on the high *k* limit. If *h* is independent of *k*, the solutions obtained above are still valid in both cases, after redefining the model coefficients. The directional solutions converge to the unidirectional solutions discussed above, as the directionality approaches the unidirectional distribution.

First we examine the case without spatial sheltering, in which the air momentum and energy and wave energy equations are expressed, respectively, as

and

where

If *h* is independent of *k*, *u* = (1 + *γ*′^{−1})*c* is proportional to *c*, as in the unidirectional case; compare with (48). The constant coefficient *γ*′ is determined by the wave energy equation

To combine the governing equations to a single equation for *S*, in analogy to (42), we approximate from the wave energy equation

where the coefficient *d _{n}* is defined by

Following similar steps as in section 2e, we find that *S* is governed by

where *γ*″ = *γ*′[(*d*_{1} − 1)*γ*′ + *d*_{1}]^{−1} and *κ*′ = *κ*/*d*_{1}. Note that the last equation is identical to (42), after redefining the coefficients. For example, if we set *h* ∝ cos^{3}*θ* and *γ* = 0.3, *d*_{1} = (9/32)*π* ≈ 0.88, *γ*′ = 0.26, *κ*′ = 0.45, and *γ*″ = 0.31. The error introduced by (B8) is roughly 5%. The coefficients in (B7) and (B10) change at most by about 13%.

Next we investigate the effect of spatial sheltering, which is more pronounced at high *k*, where a significant aerial fraction is covered by separation regions. Therefore, we will focus on the case *α* ≪ 1 for *k* → ∞. Including the directionality, the wave energy equation in (26) is

For *α* ≪ 1 this balance can only be established for *u* ≫ *c*, and hence arcos(*c*/*u*) ≈ *π*/2. Neglecting small terms, the wave energy equation becomes

Following similar steps as in section 2e and section 4a(2), the asymptotic solutions for *U* and *S* are *U*_{∞} = 1 and

with *δ*′ = *δd*_{1} and *γ*‴= *γ**d*_{2}. The last equation is identical to the unidirectional case (62), after redefining the model coefficients. For *h*∝ cos^{3}*θ*, *d*_{2} = 4/5 ≈ 0.89, so that the difference in the coefficients does not exceed 12%.

## Footnotes

*Corresponding author address:* Tobias Kukulka, Graduate School of Oceanography, University of Rhode Island, Narragansett, RI 02882. Email: tkukulka@gso.uri.edu