Parameterization of Wave Boundary Layer

Abstract

It is known that drag coefficient varies in broad limits depending on wind velocity and wave age as well as on wave spectrum and some other parameters. All those effects produce large scatter of the drag coefficient, so, the data is plotted as a function of wind velocity forming a cloud of points with no distinct regularities. Such uncertainty can be overcome by the implementation of the WBL model instead of the calculations of drag with different formulas. The paper is devoted to the formulation of the Wave Boundary Layer (WBL) model for the parameterization of the ocean-atmosphere interactions in coupled ocean-atmosphere models and wave prediction models. The equations explicitly take into account the vertical flux of momentum generated by the wave-produced fluctuations of pressure, velocity and stresses (WPMF). Their surface values are calculated with the use of the spectral beta-functions whose expression was obtained by means of the 2-D simulation of the WBL. Hence, the model directly connects the properties of the WBL with an arbitrary wave spectrum. The spectral and direct wave modeling should also take into account the momentum flux to a subgrid part of the spectrum. The parameterization of this effect in the present paper is formulated in terms of wind and cut-off frequency of the spectrum.

1. Introduction

A theory of the surface layer for the case of the flat underlying surface is well developed and confirmed by experimental data. The Monin-Obukhov similarity theory for the stratified surface layer is used in all atmospheric models. The lowest layer of such numerical model is usually identified as the surface layer. It is supposed that the surface layer is close to the state of local horizontal homogeneity and stationarity and, consequently, the vertical divergence of momentum, energy and vapor is equal to zero. It means that the vertical fluxes remain constant by height. As a rule, in the lowest layer all other terms of equations are also taken into account, but their influence is small, so, the vertical structure of layer can be described with the Monin-Obukhov similarity theory. Above the surface layer up to the height of the order of one kilometer, extends the atmospheric boundary layer. Attempts were also made to parameterize the effects of this layer (which is often called the ‘Ekman layer’), but generally with no success. This was because in contrast to the surface layer, all terms of the equations, including the Corioles force, are of the same order of value. Hence, the structure of the Ekman layer cannot be in the least degree universal. The so called Ekman layer turns out to be more complicated than ‘free atmosphere’ because of the strong influence of the vertical turbulent fluxes produced by velocity shear and convective motions.

The surface layer above the sea is generally similar to the surface layer above a solid surface, so for the calculation of fluxes, the same methods are usually accepted. However, there is a significant difference in the presence of a moving and curvilinear lower boundary. The most popular way to take this property into account is the introduction of the Charnock relation for the roughness parameter of sea surface $z_0$:

$$z_0=m{v}_{*}^2/g$$

(1)

where $v_{*}={\left|\tau \right|}^{1/2}$ is the friction velocity ($\tau$ is a vector of the vertical kinematic momentum flux at the upper level of surface layer); $m$ is an empirical coefficient of the order of 0.01, g is the gravity acceleration. The length scale $Ch=v_{*}^2/g$ was introduced in a famous one-page paper by Charnock [1]. This scale always appears in a theory of wind-wave interaction, so it is reasonable to call it the ‘Charnock scale’. Formula (1) can be used for estimation of the drag coefficient $C_h$ at height $z=H:$

$$C_h=\frac{{\kappa }^2}{{\ln }^2\left(H/z_0\right)}$$

(2)

where $\kappa =0.41$ is von Karman constant.

Formula (1) does not include any characteristics of the wave field. It is implicitly assumed that waves are fully developed and defined by wind velocity through $v_{*}.$ Such a suggestion does not seem senseless because on average stronger wind produce higher waves. This is why the Charnock’s scaling slightly reduces the scatter of dependences on the sea surface roughness parameter. Such property is a likely reason for the popularity of the Charnock formula (1). Uncertainty of such an approach, however, is expressed in a large range of values for coefficient $m$, which according to different sources varies within the interval 0.01–0.05. It can be expected that the value of this coefficient can influence a wind forecast above sea. It looks like the Charnock scale can be replaced by the significant wave height $H_s$:

$$H_s=4{\left({\displaystyle \underset{0}{\overset{{\omega }_r}{\int }}{\displaystyle \underset{0}{\overset{2\pi }{\int }}S\left(\omega ,\theta \right)d\theta d\omega }}\right)}^{1/2}$$

(3)

(where $S\left(\omega ,\theta \right)$ $\left(m^2{s}^{-1}\right)$ is the wave spectrum in polar coordinates: $\theta$ (radians) is a direction counted counterclockwise from east and $\omega$ ($rad{s}^{-1}$) is frequency). Value of $H_s$ characterizes wave energy explicitly; however, such scaling turns out to be worse than the Charnock scale. As shown in this paper, such an effect can be explained by the strong dependence of friction on high-frequency waves which produce most of the drag but contribute little to $H_s$, while big waves define $H_s$ but contribute a much smaller part of the drag.

The results of this paper show that, for the same wind speed, the coefficient $m$ increases with the growth of wave age $c_p/U_{10}$ ($c_p$ is a phase velocity in a peak of the spectrum, $U_{10}$ is a wind velocity at height 10 m), which means that the connection between wind and roughness is not direct. It is evident that the roughness of the sea surface depends on a specific shape of the wave spectrum $S\left(\omega ,\theta \right)$ ($\omega$ is wave frequency, $\theta$ is wave direction). Many other factors can influence the value the of drag coefficient.

At present, many branches of geophysical fluid mechanics need a more detailed description of air-sea interactions. The most important are the modeling of coupled circulation of ocean and atmosphere (including climate models and the medium-range weather forecast models), and the models for surface wave forecast. The coupled models require accurate formulations of the heat, moisture and momentum exchanges between ocean and atmosphere. The models of wave forecast are based on the equation of evolution of wave spectrum $S(x,y,t,\omega ,\theta )$ (x and y are coordinates, t is time). Their task is calculations of the evolution of spectrum at space grid under the action of energy input, dissipation and advection of wave energy. Energy input depends essentially on the structure of boundary layer above waves.

Investigation of the wind and wave coupling was initiated by Ursell [2] who made the only conclusion that this process was not understood. His review stimulated two theories of the generation of wind waves by Miles [3] and Phillips [4]. The first of those theories explained the role of the pressure disturbances produced by waves in the air, while the second one explained the initial generation of waves by turbulent fluctuations of air pressure in the wind.

All further models designed for the description of interaction of waves and wind were based on Miles’ linear theory in which the finite amplitude waves were actually absent. Such models considered the flow structure above the single mode while for calculation of the integral effect a linear superposition was used. Hence, the model is constructed in the coordinates linked with the mean level. However, the mean level is not a physical object contrary to the real physical interface where a kinematic surface condition is valid. This condition assumes the absence of exchange by finite volumes of fluid through the surface. The assumptions used in linear theory oversimplify the process. In reality the fluctuating boundary layer rises above the moving interface which becomes a viscous sub-layer close to the surface. Formally, the averaging of the Navier-Stokes equations in Cartesian coordinates (which is possible to perform above sea crests) gives nothing except the turbulent fluxes of momentum, i.e., the fluxes produced by wave fluctuations of velocity, pressure and stresses are absent.

2. Wave Boundary Layer (WBL)

The energy is transferred to waves by non-uniform distribution of the non-static pressure on the surface. The formation of spatial pressure distribution occurs over the entire Wave Boundary Layer (WBL). The WBL is defined as the lowest part of the atmospheric boundary layer where the fluctuations produced by waves are clearly manifested [5,6,7,8,9]. The height of WBL is equal to several significant wave heights $H_s$. A typical height of the WBL for medium waves is several meters, while for big waves it can reach 20 meters [10]. At the bottom, the WBL contacts with the wave surface, while on top it merges with the Monin-Obukhov stratified boundary layer [11]. In most cases the WBL height does not exceed the dynamic sub-layer height; hence, a direct influence of stratification on wind-wave interaction is small. Inside WBL the movement is influenced by surface waves. As waves produce fluctuations in the WBL, its structure depends on the characteristics of a wave field.

The fully nonlinear modeling of wind and wave interactions, as well as the concept of WBL was presented in papers [5,6,12]. A 2-D model of the WBL was used for the investigation of some problems of wind-wave interaction (see review [8]) and for formulation of a first version of the input energy in the WAVEWATCH model [13,14]. Later the approach was improved by the use of the conformal coordinates and the Fourier-transform method and by the coupling of WBL with the conformal wave model [15,16]. This approach is based on several important principles: (1) the model is formulated in a curvilinear surface-following coordinate system; (2) the waves are the object of modeling, i.e., the nonlinear equations of potential motion are solved simultaneously with the equation for the WBL and with matching of the solutions on interface; (3) due to the periodicity condition the Fourier-transform method is used; (4) the problem is formulated as a typical problem of statistical fluid mechanics, i.e., the model generates extended ensembles of 2-D fields of main variables which then are processed with the many-way statistical processing. A coupled model was used for a detailed investigation of the 2-D structure of the WBL and the development of a new scheme for the calculation of the energy input ([17], see also book [18]). Currently, a new version of the 2-D coupled model is underway.

The WBL model uses turbulent closure schemes based on the equation of the turbulent energy balance and additional hypothesis. It means that such a model does not include a direct interaction of turbulence with waves (Phillips mechanism). It does not mean that the energy of parameterized turbulence is not influenced by waves since it is generated and transferred by the full velocity field which includes the wave-produced components of velocity. Hence, it is assumed that the flow is affected by strong viscosity, the characteristics of which depend on the mean velocity field. A much more general approach is represented by three-dimensional models. Actually, the 3-D modeling of the wind-wave interaction was launched with so-called direct (or turbulence-resolving) models (see [19,20,21,22,23,24]. This type of model does not require any scheme for parameterization of turbulence, since it is assumed that the fluctuations of all scales are simulated explicitly, while stabilization of the statistically stationary regime occurs under the action of molecular viscosity. Naturally, this type of simulation can be applied to very small volumes only. Such models can be useful for the investigation of capillary waves with air flow, but mostly they are used for the development of a numerical scheme for the 3-D Euler equation in the curvilinear coordinate system. The LES modeling assumes that motion can be separated into two parts: the explicitly simulated and parameterized fluctuations. The better the resolution, the more detailed the description of the statistical structure of flow. In practice the boundaries between the simulated and subgrid ranges are set for reasons of the computational realism. It is desirable that the energy of individually described turbulence would be much larger than the parameterized energy; however, it is not always possible. The presence of multi-scale waves imposes additional troubles, since the resolution of the LES model should be in compliance with the wave spectrum.

Unfortunately, the coupled LES modeling requires extremely large computer resources, so the usual practice of geophysical computer modeling consisting of multiple runs with variation and tuning of parameters is impossible. A most advanced LES model [21] is designed for the simulation of turbulent flow in volume $3000\times 3000\times 500m^3$. The effects of Coriolis force and stratification are taken into account. The drag produced by wave breaking is described with a simple speculative scheme. Waves are assigned as a superposition of linear waves with random phases and a prescribed spectrum. The statistical properties of such wave field can be taken into account by joining with a 3-D phase-resolving wave model (for example, the HOS model, or the model described in [18]). The aim of this work is the creation of a ’piece of nature’ that can be investigated in its full complexity. It is a very promising direction, but it should be combined with a detailed investigation of different sides of the process using more specific models.

3. Wind and Wave Dynamic Interaction

Energy and momentum are transferred from the air to the water by the surface pressure field and by tangent stress. The tangent stress cannot be assimilated in a potential wave model. Even in a non-potential wave model the flux of energy to waves through the tangent stress is very small, since the local stress is weakly correlated with the surface orbital velocities. Hence, such momentum is transferred not to waves but to surface currents. When considering the formation of the tangent stress it is necessary to take into account its spectral properties. If the resolution of wave spectrum is very high, we can suggest that the momentum is transferred through molecular viscosity at very large local gradients of a tangent component of velocity.

According to the linear theory [3], the Fourier components of surface pressure $p$ are connected with those of surface elevation through the following expression:

$$p_{k,l}+\mathrm{i}{p}_{-k,-l}=\frac{{\rho }_a}{{\rho }_w}\left({\beta }_{k,l}+\mathrm{i}{\beta }_{-k,-l}\right)\left({\eta }_{k,l}+\mathrm{i}{\eta }_{-k,-l}\right),$$

(4)

where $h_{k,l},h_{-k,-l},{\beta }_{k,l},{\beta }_{-k,-l}$, are real and imaginary parts of elevation $\eta$ and the so-called $\beta$-function (k and l denote COSINE Fourier mode and −k and −l denote SINE mode; ${\rho }_a/{\rho }_w$ is a ratio of the air and water densities). Equation (4) is a standard presentation of pressure in a wave number space above a multi-mode surface. It means that every wave mode with amplitude ${\left(h_{k,l}^2+h_{-k,-l}^2\right)}^{1/2}$ initiates the pressure mode with amplitude ${\left(p_{k,l}^2+p_{-k,-l}^2\right)}^{1/2}$ shifted off the phase of wave mode by angle $\alpha =\mathrm{atan}\frac{{\beta }_{-k,-l}}{{\beta }_{k,l}}$. Both $\beta$ coefficients are the functions of the nondimensional frequency:

$${\mathsf{\Omega }}_{k,l}={\omega }_{k,l}U/g$$

(5)

(${\omega }_{k,l}$ is a dimensional frequency, $s^{-1}$ which characterizes a ratio of wind velocity at reference level to phase velocity of $c_{k,l}=g/{\omega }_{k,l}$:

$${\mathsf{\Omega }}_{k,l}=U/c_{k,l}$$

(6)

Since the supply of wave with energy and momentum occurs in the layer, the height of which is proportional to the wave length, it is reasonable to suggest that a reference height for wind velocity should be different for a different virtual wave length (distance ${\lambda }_{k,l}/\cos {\theta }_i$ between wave peaks in wind direction, index $i$ denotes a direction of mode). Wind velocity can be found by interpolation or extrapolation to the reference level:

$$z_{i,k,l}=0.5{\lambda }_{k,l}/\cos {\theta }_i$$

(7)

The definition of $\mathsf{\Omega }$ should take into account an angle ${\theta }_i$ between the wind vector $U$ and the direction of wave mode. Finally, the virtual nondimensional frequency takes the form:

$${\mathsf{\Omega }}_{i,k,l}={\omega }_{i,k,l}\cos {\theta }_{i}U(z_{i,k,l})/g=\cos {\theta }_{i}U\left(z_{i,k,l}\right)/c_{k,l}$$

(8)

The distribution of surface pressure and the spectral structure of momentum and energy exchange depend on the shape of $\beta$-function. Spectral models use an imaginary part of $\beta$-function only, which describes the pressure mode shifted by $\pi /2$ relatively wave mode. The nonlinear phase-resolving models use both parts of $\beta$-function.

The concept of $\beta$-function was advanced by Miles [3,25,26]; its shape was investigated theoretically in [27,28] and others. Such schemes qualitatively explained the mechanism of the wind-wave interaction, but they were not able to provide a solid basis for the calculation of the momentum and energy exchange between wind and waves. For the experimental derivation of the $\beta$-function shape it is necessary to simultaneously measure the wave surface elevation and non-static pressure on the surface. The measurement of surface pressure is a very difficult problem since it should be done very close to a moving surface, preferably with a surface-following sensor. Such measurements are done quite seldom, especially in the field. The measurements were carried out for the first time by a team of authors both in the laboratory and field [29,30,31,32,33]. The data obtained in this way allowed them to construct an imaginary part of the beta-function used in some versions of wave forecasting models [34,35]. Such measurements and further data processing are quite complicated since the wave-produced pressure fluctuations are masked by the turbulent pressure fluctuations.

Different versions of the experimental $\beta$-function are given in Figure 1. As seen, the values differ by several times. The function calculated in [17] which is discussed below is very close to a relatively old version [30]. Naturally, such an agreement does not confirm the trueness of both approximations. Figure 1 demonstrates that the problem of wind-wave interaction is very far from being resolved. However, it is not surprising that the wave prediction models using such different recommendations generally give satisfactory results. The evolution of waves occurs under the influence of two competing processes: input energy and its dissipation. Dissipation is even more complicated than the input energy. Therefore, the parameters regulating dissipation as well as the entire scheme can be modified in broad limits. Finally, the comparison of wave forecast with the available experimental data gives the possibility to compensate inaccuracies in both schemes for obtaining an acceptable result.

The second way of the $\beta$-function evaluation is based on the results of the numerical investigations of the statistical structure of the boundary layer above waves with use of Reynolds equations and an appropriate closure scheme. In general, this method works so well that many problems in the technical fluid mechanics are often solved using numerical models, not experimentally [12,38]. This method was being developed beginning from [6,8], followed by [9,39,40,41]. The results were implemented in the WAVEWATCH model, i.e., a third-generation wave forecast model [13], and thoroughly validated against the experimental data in the course of developing WAVEWATCH-III [14]. This method was later improved on the basis of the more advanced coupled modeling of waves and boundary layer [17], while the beta-function used in WAVEWATCH-III was corrected and extended up to high frequencies.

In our earlier works (see review [8]) the surface was assigned as a superposition of linear modes with a prescribed spectrum. Actually, the shape of the real wave is influenced by nonlinearity: the wave crests are high and sharp, and the probability distribution of elevation deviates considerably from Raleigh distribution. Also, the waves demonstrate a tendency for nonlinear sharpening [42] and breaking. All these features strongly modify the wave drag. A more realistic description was used in a model [17] where the model of WBL was coupled with a nonlinear phase-resolving wave model. Both models were written in the surface -following conformal coordinates. The 2-D WBL model was based on Reynolds equations and the energy/dissipation closure scheme; the solutions were matched at the interface. The initial conditions for waves were assigned as a random superposition of Stokes waves which amplitudes corresponded to JONSWAP spectrum. The description of the models and the numerical scheme can be found in [18].

One of the aims of the above investigation was the evaluation of function $\beta$ in a wide range of nondimensional frequency $\mathsf{\Omega }$. For such purposes, 47 long-term (up to several hundred peak wave periods) numerical runs with the different wind velocity and spectral resolution were performed using the coupled wind-wave model. Wind velocity was assigned through the different values of wind stress $\tau =v_{*}\left|v_{*}\right|$ at the upper level directed along and against the general direction of waves. The non-dimensional peak wave number ${\mathsf{\Omega }}_p=U/c_p$ was in a range (−28, 16), while the non-dimension wave numbers were in a range (−50, 50). The total number of Fourier amplitudes for the surface pressure and surface elevation falling in this interval was equal to 1,390,100. For the estimation of $\beta$ values, the above interval was separated into 100 bins with width $\Delta \mathsf{\Omega }=0.1$. The number of points falling within each interval was roughly the same, i.e., around 1400. Generation of such large ensemble was necessary because together with the wave-produced pressure fluctuations the signal included a wide range of the pressure fluctuation caused by nonlinearity.

The structure of flow in WBL is illustrated by tracers that were initially placed in the cells of the regular conformal grid. The trajectories of tracers are not affected directly by turbulence, otherwise they would be less regular, and so, they are carried by the averaged flow. It is seen that separation of flow is a frequent phenomenon; it often happens behind more or less steep waves, but not necessarily after the breaking. We believe that the overflow of breaking waves is very similar to that of steep non-breaking waves. The events of breaking essentially modify the distribution of pressure and all these cases are, in fact, taken into account in the parameterization of the surface pressure distribution described below. Not that a picture in Figure 2 is remarkably similar to the pictures with tracers obtained in a laboratory wind-wave channel (described in [43]).

We found that the processing of the signal done in a paper [17] was not accurate enough and was repeated here with some improvements. Note that the data described above were obtained with an undirected model. It is why in the equations below, the index $l$ is omitted. Before calculating of $\beta$ values we performed filtration of the data excluding the points with an excessive deviation of a mean value. In one-dimensional wave field the one-point values of real ${\beta }_k^i$ and imaginary ${\beta }_{-k}^i$ parts of $\beta$-function can be calculated through the Fourier coefficients for $p_0$ and $\eta$ for each bin $\Delta {\mathsf{\Omega }}_i$:

$${\beta }_{-k}^l={\left|{\eta }_k^i\right|}^2\left(p_{-k}^i{\eta }_k^i-p_k^i{\eta }_{-k}^i\right),{\beta }_k^i={\left|{\eta }_k^i\right|}^{-2}\left(p_k^i{\eta }_k^i+p_{-k}^i{\eta }_{-k}^i\right),$$

(9)

while the averaged values of ${\beta }_k$ and ${\beta }_{-k}$ were calculated in each bin with i number using rms method:

$${\beta }_k={\displaystyle \sum _{i=1}^{N_i}\left({\left|{\eta }_k^i\right|}^{-2}\left(p_k^i{\eta }_k^i-p_{-k}^i{\eta }_{-k}^i\right)\right)}{\left({\displaystyle \sum _{i=1}^{N_i}{\left|{\eta }_k^i\right|}^4}\right)}^{-1}$$

(10)

$${\beta }_{-k}={\displaystyle \sum _{i=1}^{N_i}\left({\left|{\eta }_k^i\right|}^{-2}\left(p_{-k}^i{\eta }_k^i-p_k^i{\eta }_{-k}^i\right)\right)}{\left({\displaystyle \sum _{i=1}^{N_i}{\left|{\eta }_k^i\right|}^4}\right)}^{-1}$$

(11)

where $N_i$ is the number of points falling in each $\Delta {\mathsf{\Omega }}_i$ bin.

The data marked by dots and by the + sign show the imaginary and real part of function $\beta$. A fragment in the middle of Figure 3 shows an imaginary part of function $\beta$.

The data obtained in [17] fall in the interval $\left(-50<\mathsf{\Omega }<50\right)$. However, the scatter of the data outside this interval $\left(-20<\mathsf{\Omega }<20\right)$ is large, hence, the approximation of $\beta$-function may not be accurate (though, in our opinion, it is still acceptable). Usually, such high values seldom occur. The imaginary part of function $\beta$ is negative in the interval $\left(0<\mathsf{\Omega }<0.7\right)$ where the phase velocity of waves is larger than the wind velocity; hence, the momentum and energy are transferred from waves to wind. It means that the pressure at the upwind slope of waves is less than the pressure at the downwind slope; hence, waves accelerate wind and lose their momentum and energy in a vicinity of $\mathsf{\Omega }=0$.

In all other intervals of $\mathsf{\Omega }$ the flux of energy and momentum is directed from wind to waves, i.e., the pressure at the upwind slope of a wave is larger than at the downwind slope what decreases the absolute value of wind velocity. The flux of momentum is directed to waves, but the wave momentum is negative, hence, the waves decelerate, and the energy dissipates. In the cases when wind blows against waves, the function is larger than in the opposite case with the same value of $\left|\mathsf{\Omega }\right|$, so the function $\beta$ is not symmetric regarding $\mathsf{\Omega }=0$. The processing of the large volume of data generated in [17] based on Equations (9)–(11) gave the approximation:

$${\beta }_i=\{\begin{array}{cc}{\beta }_0+a_0\left(\mathsf{\Omega }-{\mathsf{\Omega }}_0\right)+a_1{\left(\mathsf{\Omega }-{\mathsf{\Omega }}_0\right)}^2& \mathsf{\Omega }>{\mathsf{\Omega }}_0\\ {\beta }_0-a_0\left(\mathsf{\Omega }-{\mathsf{\Omega }}_0\right)+a_1{\left(\mathsf{\Omega }-{\mathsf{\Omega }}_0\right)}^2& \mathsf{\Omega }<{\mathsf{\Omega }}_0\end{array}$$

(12)

$${\beta }_r=\{\begin{array}{cc}{\beta }_1+a_3\left(\mathsf{\Omega }-{\mathsf{\Omega }}_2\right)& \mathsf{\Omega }<{\mathsf{\Omega }}_2\\ a_2{\left(\mathsf{\Omega }-{\mathsf{\Omega }}_1\right)}^2& {\mathsf{\Omega }}_2<\mathsf{\Omega }<{\mathsf{\Omega }}_3\\ {\beta }_1-a_3\left(\mathsf{\Omega }-{\mathsf{\Omega }}_3\right)& \mathsf{\Omega }>{\mathsf{\Omega }}_3\end{array}$$

(13)

where ${\mathsf{\Omega }}_0=0.355$, ${\mathsf{\Omega }}_1=1.2,$ ${\mathsf{\Omega }}_2=-18.8$, $a_0=0.0228$, $a_1=0.0948$, $a_2=-0.0372,$ $a_3=14.8$, ${\beta }_0=-0.02,$ ${\beta }_1=-148.0$.

The $\beta$-function can be represented as a dependence of time scale for energy input $T_{inp}$ for the modes with nondimensional wave number $\mathsf{\Omega }$, normalized by period of mode $T_k$:

$$\tilde{T}=\frac{T_{inp}}{T_k}=\frac{{\rho }_w}{2{\rho }_a\pi {\beta }_i}$$

(14)

(here the imaginary part of $\beta$ is included; ${\rho }_a$ and ${\rho }_w$ are densities of air and water). The Equation (14) gives an estimation of the timescale for the rate of wave growth in terms of wave period $T_p$ as a function of $\mathsf{\Omega }$. This function varies by orders and has both positive and negative values. Therefore, it is convenient to use the following transformation:

$$R=\{\begin{array}{cc}{\lg }_{10}\left(1+\tilde{T}\right)& \tilde{T}>0\\ -{\lg }_{10}\left(1-\tilde{T}\right)& \tilde{T}\le 0\end{array}$$

(15)

For estimations, formula $\left|\tilde{T}\right|\approx {10}^R$ can be used. Vertical lines at $\mathsf{\Omega }=0$ and $\mathsf{\Omega }=0.7$ correspond to the zero values of $\beta$ (see inset in Figure 3) and infinite time scale. Timescale is everywhere positive, which means that energy goes to waves. It is negative (an inverse flux of energy) in a narrow interval $\left(0<\mathsf{\Omega }<0.7\right)$ where waves give their energy to air at the timescale of the order of ${10}^4{T}_k$.

The same order of value has a positive timescale for small frequencies. For $\left|\mathsf{\Omega }\right|\approx 10$ the timescale for input is of the order of 10, while for $\left|\mathsf{\Omega }\right|\approx 20$ the timescale is as small as few periods of mode. It means that only low-frequency waves, typically in the peak of the spectrum, are conservative, but high-frequency modes in the tail of the wave spectrum are short-living. Because the shape of the tail is more or less steady, the rate of dissipation in the spectral tail is also high. Consequently, the suggestion that the spectral tail is mostly formed of the nonlinear cascade of energy cannot be confirmed. Note that the approximation suggested in [33], (Figure 1) gives unrealistic small-time scales, i.e., by orders smaller than in Figure 4.

The non-stationarity of wave modes as well as the non-linearity of wave field as a whole (e.g., the tendency for being a superposition of Stokes-like modes), variability of phase velocities (see [44]) and the presence of the broad spectrum of pressure fluctuations of non-wave nature all indicate that the connection between $\eta$ and $p$ is not deterministic. This is why the calculations, (10), (11) based on the results of many runs with different initial conditions and wind velocities give an effective shape of $\beta$-function, which, in fact, implicitly takes into account most of the complicating factors mentioned above, including the effects of random flow separation. The time scales of these processes are much smaller than those for energy growth [44]. This is why we tend to treat Equation (4) not as a connection between instantaneous amplitudes of pressure and elevation but rather as a connection between the quantities averaged over time and over some frequency intervals.

4. Wave-Produced Momentum Flux

The existence of interface $\eta \left(x,y\right)$ dictates to use the surface-following coordinate system where the vertical coordinate $\zeta =0$ at $z=\eta \left(x,y\right)$. Above the flat surface $\eta =0$ the natural boundary condition is a zero vertical velocity, i.e., $w=0$. In the conformal coordinates the role of vertical velocity is played by the contravariant component of velocity $W$ ([15,18]). The kinematic condition on surface $W=0$ suggests the absence of exchange through the surface by finite volumes.

An equation of evolution of horizontal momentum can be obtained by averaging of equation for $u$ component of velocity along the coordinate lines $\zeta =const$. Brackets $\langle \rangle$ denote the averaging (see details in [18])

$$\begin{array}{cccc}\frac{\partial \langle u\rangle }{\partial t}& =\frac{\partial }{\partial z}(-\langle uW\rangle & +\langle {\eta }_{x}p\rangle & +\langle -{\eta }_x\overline{u^{\prime }{u}^{\prime }}+\overline{u^{\prime }{w}^{\prime }}\rangle )\\ \langle 1\rangle & \langle 2\rangle & \langle 3\rangle & \langle 4\rangle \end{array}$$

(16)

Since $\langle \varsigma \rangle =z$, and $\langle \partial /\partial \varsigma \rangle =d/dz$ $\langle \partial /\partial \zeta \rangle =\partial /\partial z$, Equation (16) can be considered as written in the Cartesian coordinate system.

According to (16) the rate of evolution of horizontal momentum $\langle 1\rangle$ depends on vertical divergence of vertical momentum flux generated by wave-produced components of velocity $\langle 2\rangle$, by pressure $\langle 3\rangle$ and by the averaged turbulent flux $\langle 4\rangle$.

Equation (16) is analogous to the equation of momentum balance above flat surface

$$\frac{\partial u}{\partial \tau }=\frac{\partial \langle u^{\prime }{w}^{\prime }\rangle }{\partial z},$$

(17)

which, in fact, is used in linear theories and can be obtained by averaging in the Cartesian coordinate system above sea crests. Equation (17) is useless since it does not include wave effects. Note that the averaging in a surface-following coordinate system is physically performable both in the laboratory [32] and the field [33].

Let us represent Equation (16) in the form

$$\frac{\partial U}{\partial t}=\frac{\partial }{\partial z}\left(\tau +{\tau }_w\right)$$

(18)

where $U=\langle u\rangle$ is a wind profile, $\tau =\langle \overline{u^{\prime }{w}^{\prime }}\rangle$ is a turbulent momentum flux, and ${\tau }_w=\langle {\eta }_{x}p+-{\eta }_x\overline{u^{\prime }{u}^{\prime }}+\overline{u^{\prime }{w}^{\prime }}\rangle$ is the Wave Produced Momentum Flux (WPMF). Vectors $U,\tau$ and ${\tau }_w$ are directed along x axis.

The vertical profiles of WPMF were investigated in detail in [17], where the numerical experiments with the coupled wind/wave model were described. A wave boundary layer (WBL) model is based on Reynolds equations with the $K-\epsilon$ closure scheme. The solutions for air and water matched through the interface. The structure of the WBL and vertical profiles of the wave-produced momentum flux (WPMF) in a long-term simulation of the coupled dynamics are investigated and parameterized. The averaged profiles of spectral components of WPMF calculated in [17] are shown in the top panel of Figure 5 reproduced from [17]. As seen, the shape of those profiles depends on the wavenumber. These data can be arranged by normalizing each profile by its surface ${\tau }_{w0}^k$ value and representing them as a function of nondimensional height. After that transformation the profiles converge. The convergence can be even closer when by introducing the following relation:

$${\tau }_w^k={\tau }_{w0}^k\exp \left(Gkz\right),$$

(19)

where G is a function of $\tilde{\omega }=\mathsf{\Omega }/{\mathsf{\Omega }}_p=\omega /{\omega }_p:$

$$G=0.985+0.4{\tilde{\omega }}^{0.81}$$

(20)

Note that at the moderate spectral resolution the introduction of function G in (19) changes the result insignificantly, hence, in such cases it is possible to put $G=1$. The profiles of spectral components of ${\left({\tau }_{\omega }/{\tau }_{w0}\right)}_k$ are given in the bottom panel of Figure 5 as a function of $Gkz$(20). As seen, Formula (19) demonstrates a good approximation of vertical profiles of WPMF.

5. One-Dimensional Model of the Boundary Layer above Waves

The WBL is defined as the lowest part of the atmospheric boundary layer where a considerable part of momentum and energy is transferred by fluctuations of velocity and pressure produced by waves. The WBL height is of the order of several significant wave heights, typically of the order of 10 m. In the bottom, WBL adjoins the wave surface and in top it merges with the stratified Monin-Obukhov boundary layer. In the majority of cases the WBL height does not exceed the height of the dynamic sub-layer, so it is not necessary to take into account the effects of stratification, though these effects can be easily included in the model below. The equations of non-dimensional WBL were derived in [17] with an averaging of two-dimensional equations in the conformal coordinates along the coordinate lines. Here we present a more general case considering velocity as a vector velocity $u$ with components $\left(u,v\right)$:

$$\begin{array}{l}\frac{\partial u}{\partial t}=\frac{\partial }{\partial z}\left(K\frac{\partial u}{\partial z}+{\tau }_{wx}\right)\\ \frac{\partial v}{\partial t}=\frac{\partial }{\partial z}\left(K\frac{\partial v}{\partial z}+{\tau }_{wy}\right)\end{array}$$

(21)

$$\frac{\partial e}{\partial t}=\frac{\partial }{\partial z}{K}_e\frac{\partial e}{\partial z}+P-\epsilon ,$$

(22)

where $\left({\tau }_{wx},{\tau }_{wy}\right)$ is a vector of vertical WPMF; e is kinetic energy of turbulence; $K=\kappa z{e}^{1/2}$ is the coefficient of turbulent viscosity ($\kappa$ is Karman constant); τ_w—wave-produced momentum flux; P—a rate of turbulent energy generation:

$$P=\frac{\partial u}{\partial z}\left(K\frac{\partial u}{\partial z}+{\tau }_{wx}\right)+\frac{\partial v}{\partial z}\left(K\frac{\partial v}{\partial z}+{\tau }_{wy}\right)$$

(23)

In a previous version of the model the system of equations included an equation of evolution of the dissipation rate $\epsilon$ together with the formula written with dimension consideration. Fortunately enough for the current problem this equation is not required as it does not change the results noticeably. In fact, this equation does not seem to be physical, and in complicated problems it often gives unrealistic values removed with the use of the so-called ‘limiters’.

In the current work the averaged over horizontal domain structure of the wave boundary layer above wave surface assigned by a two-dimensional wave spectrum $S({\omega }_i,{\theta }_j)$($m^{2}s$) is under consideration. Axis $x$ is directed to East, while axis $y$—to North. The frequency ${\omega }_i$ varies in the interval $\left(0,{\omega }_M\right)$, the wave direction $\theta$ changes counterclockwise in a range $\left(0,{\theta }_M\right)$. The frequency grid is set with the three parameters: number of intervals $M$, top value of frequency ${\omega }_M$ and stretching coefficient ${\gamma }_{\omega }.$

At the upper bound of domain $z=H$ the components of wind velocity are assigned:

$$u_{z=H}=U_h\cos \psi ,v_{z=H}=U_h\sin \psi ,$$

(24)

where $U_h$ is wind velocity; $\psi$ is direction. In the description of wind-wave interaction in Section 2 it was accepted that direction of wind coincides with k-axis. It is assumed that WPMF attenuates at $z=H$ and the value of turbulent energy at $z=H$ is proportional to the turbulent stress ${\tau }_{tH}:$

$$e_{z=H}=c_1{\left({\tau }_{tx}^2+{\tau }_{ty}^2\right)}_{z=H},$$

(25)

where ${\tau }_{tx}$ and ${\tau }_{ty}$ are the components of horizontal stress; $c_1$ is a constant. The same rule is used for calculation of the energy of turbulence at the lowest level $z=z_1$:

$$e_{z=z_1}=c_1{\left({\tau }_{tx}^2+{\tau }_{ty}^2\right)}_{z=z_1}.$$

(26)

The equations of the WBL model are approximated at the irregular stretched grid $\Delta z_{j+1}={\gamma }_z{\Delta }_j$ where ${\gamma }_z$ is a stretching coefficient. In a research version of the model considered here, the height of the domain was 10 m, the number of levels $L$ in the air was equal to 50 and the stretching coefficient was taken as ${\gamma }_z=1.10$. The variables $u,$ $v$ and $e$ $\mathrm{were}$ placed to the middle of the layers; their fluxes and wave-produced momentum flux were placed in accordance to the intermediate levels. A high value of ${\gamma }_z$ was used for setting a high resolution close to the surface so that the thickness of the lowest level was of the order of height of the viscous sub-layer $z_{\nu }=60\nu /v_{*0}$ ($\nu =0.15\cdot {10}^{-4}{m}^2/s$ is a molecular air viscosity and $v_{*0}$ is the local friction velocity). In this case an effective roughness parameter can be calculated with Formula [11].

$$z_{00}=0.11\nu /v_{*0}$$

(27)

The local surface drag coefficient $C_0$, surface tangent stress ${\tau }_{t0}$ and local friction $v_{*0}$ are calculated by a formula:

$$C_0={\left(\frac{\kappa }{\ln \raisebox{1ex}{$z_1$}\!\left/ \!\raisebox{-1ex}{$z_{00}$}\right.}\right)}^2,{\tau }_{t0}=C_0{u}_1^2,v_{*0}={\tau }_{t0}^{1/2};$$

(28)

The method of WPMF profile calculation requires more detailed consideration.

Typically, the wave spectrum $S\left(\omega ,\theta \right)$ contains thousands of values of wave energy. For calculation of the energy input to each interval $\Delta {\omega }_i\Delta {\theta }_j$ it is necessary to calculate the wind velocity value at level ${\lambda }_i/2$. Since the vertical grid is irregular, this simple operation procedure turns into a big problem, especially, if the WBL model is used for the coupling of the atmosphere with the wave model or ocean model. A simple way to overcome this problem is the approximation of wind profile $U={(u^2+v^2)}^{1/2}$, for example by a relation:

$$U(z)=W-A\ln \left(\raisebox{1ex}{$H$}\!\left/ \!\raisebox{-1ex}{$2z$}\right.\right),$$

(29)

where

$$W=U_h\frac{\ln \left(\raisebox{1ex}{$z$}\!\left/ \!\raisebox{-1ex}{$z_0$}\right.\right)}{\ln \left(\raisebox{1ex}{$H$}\!\left/ \!\raisebox{-1ex}{$z_0$}\right.\right)},U_h={\left(u_h^2+v_h^2\right)}^{1/2},$$

(30)

where coefficient $A$ is calculated by the root mean square method:

$$A=\frac{{\displaystyle \sum _1^L{W}_i\ln \left(H/2./z_i\right)}}{{\displaystyle \sum _1^L{\left(\ln \left(H/2./z_i\right)\right)}^2}}$$

(31)

where the summation is performed over the vertical coordinate. Formula (29) with coefficient A is used for calculation of the interpolated and extrapolated values of wind velocity $U_k$:

$$U_k=U_h\frac{\ln \left(z_k/z_0\right)}{\ln \left(H/z_0\right)}-A\ln \left(H/2/z_k\right),$$

(32)

at height:

$$z_k={\lambda }_k\left|\sec (\theta -\psi )\right|/2,$$

(33)

where $\psi =\arctan \left(u_h/v_h\right)$, and $\theta -\psi$ is an angle between the wind direction and the direction of kth mode. The function $\sec (\theta -\psi )$ is included in a reference level $z_k$ because waves moving at angle $\theta -\psi$ to the wind direction have a corresponding longer virtual length. Note that Equation (32) can give unnaturally high values of $z_k$, especially for the low-frequency waves. We restricted this quantity by limit $z_k=100,$ which did not influence the results.

The nondimensional frequency $\mathsf{\Omega }$ is calculated by a formula:

$${\mathsf{\Omega }}_{k,j}={\omega }_j{U}_k\cos \left(\theta -\psi \right)/g$$

(34)

These values of ${\mathsf{\Omega }}_{k,j}$ are used for the calculation of an imaginary part of $\beta$-function (12). The components of surface momentum fluxes from air to waves $m^2{s}^{-1}$ are calculated with the formulae:

$${\tau }_{wx0}(k,j)={\omega }_k^2{S}_{k,j}{\beta }_{k,j}\cos \psi ,$$

(35)

$${\tau }_{wy0}(k,j)={\omega }_k^2{S}_{k,j}{\beta }_{k,j}\sin \psi$$

(36)

The flux of energy to waves in the spectral cell $\Delta {\omega }_k\Delta {\theta }_j$ $\left(m^3{s}^{-2}\right)$ is calculated by formula:

$$F_{k,j}=\frac{{\rho }_w}{{\rho }_a}g{\omega }_k{S}_{k,j}{\beta }_{k,j}$$

(37)

Here, in Equations (35)–(37) ${\beta }_{k,j}$ is the imaginary part of $\beta$-function of ${\mathsf{\Omega }}_{k,j}$ (34) calculated by (12). The 1-D spectra of the surface wave flux of momentum ${\tau }_{wx0}(k)$ and ${\tau }_{wy0}(k)$ can be calculated by summation over angle $\theta$

$${\tau }_{wx0}\left(k\right)={\displaystyle \sum _{j=1}^M{\tau }_{wx0}\left(k,j\right)},{\tau }_{wy0}\left(k\right)={\displaystyle \sum _{j=1}^M{\tau }_{wy0}\left(k,j\right)}$$

(38)

while the components of total wave flux—by summation over frequencies:

$${\tau }_{wx0}={\displaystyle \sum _{k=1}^N{\tau }_{wx0}\left(k\right)}\Delta {\omega }_k,{\tau }_{wy0}={\displaystyle \sum _{k=1}^N{\tau }_{wy0}\left(k\right)}\Delta {\omega }_k.$$

(39)

The total flux of energy to waves is calculated by summation over angles $\theta$ and frequencies $\omega$

$$F={\displaystyle \sum _{k=1}^N{\displaystyle \sum _{j=1}^M{F}_{k,j}\Delta {\omega }_k}}.$$

(40)

The vertical profiles of WPMF components are calculated taking into account the exponential attenuation (19), (20):

$${\tau }_{wx}\left(z_i\right)={\displaystyle \sum _{k=1}^N{\displaystyle \sum _{j=1}^M{\tau }_{wx0}(k,j)\Delta {\omega }_k\exp \left(-K_{k,j}{G}_k{z}_i\right)}},$$

(41)

$${\tau }_{wx}\left(z_i\right)={\displaystyle \sum _{k=1}^N{\displaystyle \sum _{j=1}^M{\tau }_{wy0}(k,j)\Delta {\omega }_k\exp \left(-K_{k,j}{G}_k{z}_i\right)}}.$$

(42)

where $K_{k,j}$ is the virtual wavenumber (i.e., taking into account the direction $\psi$):

$$K_{k,j}={\omega }_k^2{g}^{-1}\left|\sec (\psi )\right|$$

(43)

and $G_k$ is the function (19).

Note that in the above equations, a variable with two indexes denotes matrix; with one index it denotes the same variable but vector obtained by summation over angle; and with no indexes–the variable denotes scalar obtained by summation over angles and frequencies.

6. Examples of the Calculations for Idealized Spectrum

A model formulated above is designed for the calculation of the drag coefficient that can be used for estimation of fluxes specific feature of the model is that it links the structure of WBL with the upper boundary condition and any arbitrary wave spectrum.

Below are the results of the calculations performed for the cases when surface waves are assigned by JONSWAP spectrum $S\left(\omega ,\theta \right)$ [45] for different wind velocities and inverse wave age ${\mathsf{\Omega }}_p=U_{10}/c_p$ ($U_{10}$ is a wind velocity at height $H=10\mathrm{m}$, $c_p$ is a phase velocity in spectral peak).

$$S\left(\omega ,\theta \right)=S_1\left(\omega \right)D\left(\tilde{\omega },\theta \right).$$

(44)

Here $S_1\left(\omega \right)$ is a one-dimensional spectrum and $D\left(\theta \right)$ is a function describing the angle distribution of the potential energy density.

$$S_J\left(\omega \right)=\alpha g^2{\omega }^{-5}\exp \left(-1.25{\tilde{\omega }}^{-4}\right){\gamma }^{\Gamma }$$

(45)

where $\tilde{\omega }=\omega /{\omega }_p$ (${\omega }_p$ is frequency in spectral peak), $\alpha$ is a coefficient depending on $\mathsf{\Omega }$ [46]

$$\alpha =0.01{\mathsf{\Omega }}^{0.66},$$

(46)

and $\mathsf{\Gamma }$ is a function of $\tilde{\omega }$:

$$\mathsf{\Gamma }=\exp \left({\left(\frac{\tilde{\omega }-1}{\sqrt{2}\sigma }\right)}^2\right),\mathrm{where}\sigma =\{\begin{array}{cc}0.7& \tilde{\omega }<1\\ 0.9& \tilde{\omega }\ge 1\end{array}$$

(47)

The angle spreading is taken in a form suggested by Donelan (1999).

The shape of spectrum depends on ${\mathsf{\Omega }}_p$ and on the frequency in spectral peak

$${\omega }_p={\mathsf{\Omega }}_{0}g/U_{10}.$$

(48)

A wave spectrum was assigned at grid $\left({\omega }_k,{\theta }_j\right),$ $k=1,2,3\dots N,N=100$ at stretched grid with a stretching coefficient ${\gamma }_{\omega }=1.03$; $j=1,2,3\dots M,M=90$ in the interval $\left(0<\theta <2\pi \right)$ with a constant step over angle $\Delta \theta =2\pi /M$. The WBL model was approximated at grid with a number of levels $L=100$ at $H=10m$ with a stretching coefficient ${\gamma }_z=1.10$. Such a high number of level and value ${\gamma }_s$ provided a high resolution close to the surface which justified a use of the roughness parameter for smooth surface (27).

A steady numerical solution of Equations (21) and (22) was obtained with a second order scheme over $z$ and an explicit time integration scheme. Such a scheme requires a very small time step $\Delta t$ estimated by condition:

$$\Delta t<0.25\min \left({\left(\Delta z\right)}^2/K\right),$$

(49)

For steady conditions the total flux of momentum $\tau$ should be constant over height

$$\tau =K\frac{\partial u}{\partial z}+{\tau }_w.$$

(50)

Note that a term $K\partial u/\partial z$ in Equation (50) at $z=0$ should be replaced by the term ${\tau }_{t0}=C_0{u}_1^2$ (see Equation (28)). The criterion of steadiness was chosen is a form:

$$\left(\max \left(\tau \right)-\min \left(\tau \right)\right)/{\tau }_H<0.001$$

where ${\tau }_H$ is the total flux of momentum at $z=H$. An explicit scheme was used here because the solution was very fast. However, when the model is used for massive calculations together with another model (for example, with a wave forecast model), it is reasonable to use the TDMA method [47].

The structure of WBL is illustrated in the cases shown in Figure 6. The area of unrealistic value of the significant wave height $H_s=11m$ corresponding to the wave height about 33 m is separated by a dotted line. This area will be shown for other characteristics of WBL.

The wind profiles shown in Figure 7 demonstrate a considerable deviation of the logarithmic behavior (described by the second term in Equation (29)). For a smaller wave age (larger ${\mathsf{\Omega }}_p$) the deviation is larger and can be as big as $10\mathrm{m}/\mathrm{s}$ for $U_{10}=30\mathrm{m}/\mathrm{s}$.

The specific property of WBL is that the mechanism of momentum transfer to wave develops very close to the surface. Figure 7 shows that deviations of the logarithmic wave profile are concentrated below the height of 1 m. It proves that the main part of momentum is transferred by small waves at high frequencies. The main advantage of a one-dimensional approach is that it allows a broad range of frequencies to be taken into account and most part of the wave produced flux. Contrary to the turbulent stress, the wave produced flux of momentum is generated directly by an external factor, i.e., by curvilinearity of a moving surface. In a steady boundary layer the vertical flux of momentum should be constant over height; hence, the turbulent momentum flux should decrease when approaching the surface, and the connection between the stress and wind profile is modified (Figure 8).

As can be seen, a decrease of turbulent momentum flux is completely compensated by the increase of WPMF that when close to the surface, can reach values as high as 0.7${\tau }_H$. For the wind that is less than 20 m/s, the WPMF turns into 0 at the upper level. For higher winds, the WPMF can reach a value of $0.1{\tau }_H$, which means that the total height of WBL is higher than $z=10\mathrm{m}$. The calculations with the increased height of domain up $H=20\mathrm{m}$ showed that the results in the layer $0-10\mathrm{m}$ obtained at complete attenuation of WPMF to the upper level are very close to the current results. Since $z=10\mathrm{m}$ is a standard reference level, we keep to this value for all variants of calculations.

The role of surface tangent stress ${\tau }_{t0}$ calculated with Equation (28) is illustrated in Figure 9. The contribution of tangent stress in the total surface stress is large for small wind and increases with growth of inverse wave ${\mathsf{\Omega }}_p$.

The dependence of these quantities on nondimensional frequency $\omega /{\omega }_p$ is shown in Figure 10 and Figure 11. As seen, both pictures demonstrate a great contribution of high-frequency modes to the fluxes. The spectral composition of the surface momentum ${\tau }_{w0}$ and energy fluxes $F_{w0}$ is convenient to consider their integrated values ${\tau }_{w0}^k$ and $F_{w0}^k$ as a function of the upper limit $\omega /{\omega }_p$:

$${\tau }_{w0}^k={\tau }_H^{-1}{\displaystyle \sum _{k=1}^k{\tau }_{w0}^k}\Delta {\omega }_k,$$

(51)

$$F_{w0}^k=F_H^{-1}{\displaystyle \sum _{k=1}^k{F}_{w0}^k}\Delta {\omega }_k,$$

(52)

where $F_H=u_h{\tau }_h$ is the total vertical flux of energy at $z=H$.

The flux of energy comes to the saturation level at $\omega /{\omega }_p\approx 100$ while the flux of momentum remains growing up to $\omega /{\omega }_p\sim 100$. It can be explained just by dimensional consideration: Equations (35) and (36) contain ${\omega }^2$, while an expression (37) contains $\omega$. The tendency for the tangent stress to drop at high frequencies in Figure 9 and the tendency of fluxes to increase in Figure 10 and Figure 11 can be explained by an insufficient resolution. This effect could be reduced by introducing a higher resolution in frequency, but a shape of the $\beta$-function for high frequencies is not well known. However, the total stress as a sum of and ${\tau }_{t0}$ ${\tau }_{w0}$ is calculated correctly.

The stationary solution for WBL can be obtained without separation of subgrid stress into a tangent and WPMF part by assuming that a total stress on surface is equal to the outer stress ${\tau }_H$.

$${\tau }_{t0}+{\tau }_{w0}={\tau }_H$$

(53)

Such formulation of the WBL problem does not require the use of high resolution in $z$ and $\omega$. However, this scheme cannot be recommended for the precise calculation of wind-wave interaction, since WBL is not always fully adjusted to waves, and the sum of turbulent momentum flux and WPMF is not constant over height (see Section 5). In general cases, the WBL can be unsteady. The nondimensional roughness parameter $z_{0n}=z_{0}g/v_{*}^2$ is calculated by relation

$$z_0=H\exp \left(-\frac{\kappa U_H}{v_{*}}\right),$$

(54)

here $v_{*}={\tau }_H^{1/2}$ is plotted in Figure 12. As seen, $z_0$ is not a constant but significantly depends on wind velocity and inverse wave age ${\mathsf{\Omega }}_p$.

The JONSWAP wave spectrum depends on the two parameters, i.e., wind velocity and inverse wave age ${\mathsf{\Omega }}_p$. Hence, the drag coefficient also depends on these parameters: the dependence on $U_{10}$ is direct, while on inverse wave age ${\mathsf{\Omega }}_p$ it depends because the larger ${\mathsf{\Omega }}_p$, the larger the steepness of the wave field which can be defined as follows:

$$St={\displaystyle \underset{0}{\overset{{\omega }_N}{\int }}{\displaystyle \underset{0}{\overset{2\pi }{\int }}{k}^{2}S(\omega ,\theta )d\omega d\theta }},$$

(55)

where $k={\omega }^2/g$ is a wavenumber. Naturally, $St$ is nondimensional, but for $S\propto {\omega }^{-5}$ its value depends logarithmically on cut-frequency ${\omega }_N$. This is why $St$ cannot serve as an absolute characteristic of wave field, but it is useful for the comparison of the wave fields with the same spectral resolution. In our case for ${\mathsf{\Omega }}_p=1$ it equals 0.25 and for ${\mathsf{\Omega }}_p=0.22$ it is 0.3. This difference is large enough to consider that the dependence on ${\mathsf{\Omega }}_p$ might be significant. The dependence of $C_{10}$ on $U_{10}$ turns out to be quadratic (Figure 13). For the developed sea (${\mathsf{\Omega }}_p=1$) it does not exceed $2\cdot {10}^{-3}$, while for the young sea it reaches the value $3.5\cdot {10}^{-3}$. A comparison of these results with the experimental data is easy, as the existing data due to a large scatter can confirm any reasonable results.

The dependence of drag coefficient on wind was studied in laboratory ([48,49]) and field ([50,51,52,53]). In general, these data show that at least up to $U_{10}=30m/s$ the drag coefficient grows (see the data collected in [54]). The values $C_{10}$ at wind velocity $25-30м/c$ fall in the interval $\left(1.5-3\right){10}^{-3}$. Above the value $U_{10}=30m$ the drag coefficient in compliance with the majority of recent data starts to drop. According to [54,55] the drag coefficient can be influenced by swell. Different functions $C_{10}\left(U_{10}\right)$ were tested being implemented in coupled models [56,57,58,59,60,61]. Those studies showed that introduction of a drag coefficient decreased at strong wind improves the forecast results.

The data from different sources being plotted together (see [54]) look chaotic. The data can equally prove and disprove any reasonable theoretical scheme. A reason for this is an absence of ideal conditions when the theory is applicable. There are many reasons why the vertical flux of momentum is not a constant over height. For example, in a specially selected area for conducting a study of the Monin-Obukhov boundary layer the conditions during the vertical flux of momentum changing twice over the distance of 10 m, was not unusual. The cause of such a phenomenon was a forest at a distance of 2 km from the point of observations, which violated the condition of horizontal homogeneity (A.M. Obukhov, private communication). Finally, it looked as if the drag coefficient depended on wind direction. Sea conditions are usually closer to homogeneity, but the deviation can be created by meso-scale wind and sea surface temperature disturbances as well as by secondary circulations in the planetary boundary layer and upper ocean, and for many other reasons. If the wind profile is not logarithmic, the value of drag coefficient can change unpredictably in a broad range.

If the sea surface is flat at the conditions of stationarity and homogeneity, the structure of the boundary layer is very simple and can be described by the well-known formulae. Hence, even at ideal conditions in the presence of waves the drag coefficient is a product of waves. Wave spectrum is seldom described by the JONSWAP approximation [45]; it often contains more than one peak. Also, a general direction of waves can deviate from the wind direction. The angle distribution even for the JONSWAP approximation is not well known. The drag coefficient depends on many factors:

1. A specific shape of the spectrum. The spectrum can have more than one maximum; it contains swell and high-frequency tail. The spectrum can be asymmetric about wind direction, and the general direction of energy-containing waves can deviate from the wind direction. The data on time scales for the energy input (Figure 4) show that high-frequency waves should be adjusted to the local wind conditions and are most likely generated by the local in time and space wind.

2. The local non-stationarity and non-homogeneity are caused by change of wave spectra and wind velocity. In this case wind profile is not adjusted to the local wave spectra and can considerably deviate from the structure shown in Figure 7. At the time scales of a governing model the WBL can be also nonstationary, hence, we cannot rely on stationary solution.

Great variations of drag coefficient can be observed in the complicated processes such as tropical storms where waves are not directly connected with local wind, and the conditions are very far from being homogeneous and stationary. In this case, the characteristics of ocean-atmosphere interaction are the product of local wind and waves as well as of the prehistory of the process. Definitely, drag coefficient should be different in the different phases of storm. It means that this process should be described not by the selected formulae and rules, but by a system of the equations describing a dynamic interaction between the atmospheric boundary layer and upper ocean taking into account sea waves.

7. Pragmatic Model of the WBL

A one-dimensional WBL model described above is not convenient for regular use in a frame of the ocean-atmosphere-sea wave models because it is based on direct calculations of high-frequency waves and tangent stress in a viscous sub-layer. Hence, the formulation of such a model suggests a high resolution for wave spectrum and for vertical coordinate. Such a model is well suited for research. The spectral wave models usually do not include high-frequency waves. Besides, the high resolution WBL model being included in a large model takes considerable time, which can slow down the entire process. The straight way to avoid such a problem is to assume the stationarity of WBL. In this case WPMF ${\tau }_w\left(z\right)$ for the wave modes included in wave spectrum are calculated directly with Equations (41) and (42), while a surface boundary condition for turbulent stress is taken in a form:

$${\tau }_{tx0}={\tau }_{Hx}-{\tau }_{wx0},{\tau }_{ty0}={\tau }_{Hy}-{\tau }_{wy0}$$

(56)

following from the assumption of the homogeneity and stationarity of the WBL. With such formulation the WBL model can use a poor resolution by vertical (about 10–20 levels). This approach is too simple to be universal, because the real process is seldom close to the conditions assumed, hence Equation (56) is, in general, inapplicable. However, the ocean-atmosphere and wave modeling require the formulation of a simplified approach, what we call here a ‘pragmatic’ model.

Let us consider a model of the WBL (21)–(23) for the calculation of the WBL structure above waves assigned by an arbitrary wave spectrum $S(\omega ,\theta )$. The upper boundary conditions (24) and (25) remain unchanged. The principal difference between the research model and pragmatic model is that wave spectrum is assigned up to the relatively low frequencies, of the order of 1 Hz, i.e., the smallest wave has the length of the order of $10m$. We should admit that such spectral resolution is too poor for a good description of the WBL structure, because even the smallest wave with the frequency of 1Hz, is supplied by energy from a layer with the thickness of the order of WBL height. In fact, wave prediction models are the models of big waves. However, we hope that the resolution of wave models will be gradually improving, at least for research purposes.

Let us suppose that the wave model uses a grid with the maximum frequency ${\omega }_m$. It means that the influence of waves with smaller frequencies is described individually in terms of the spectral components of momentum flux calculated with a use of $\beta$-function. The influence of the subgrid modes with $\omega >{\omega }_m$ should be taken into account parametrically. Let us suggest that above frequency ${\omega }_m$ waves are generated by current wind; the spectrum is with respect to wind direction and the energy is falling down to ${\omega }^{-5}$. Such assumptions are not voluntary, since according to Figure 4, high-frequency waves have a small timescale, i.e., it would be reasonable to suggest that they are generated by local wind. Such situation can happen, for example, in a tropical storm ${\tau }_{0w}$ when the wind quickly changes direction and starts to generate new waves. The problem is describing the total momentum flux ${\tau }_{0m}$ to waves in the interval $\omega >{\omega }_m$:

$${\tau }_{0m}={\tau }_{0w}+{\tau }_{0t},$$

(57)

where

$${\tau }_{0w}={\displaystyle \sum _{k_m}^N{\tau }_{0w}^k}\Delta {\omega }_k,$$

(58)

$k_m$ is an index of frequency ${\omega }_m$ and $N$ is a number of frequencies in the research model. The data on ${\tau }_{0w}$ and ${\tau }_{0t}$ were generated in a series of calculations with a high-resolution model for the JONSWAP spectrum. The parameterization means that the values of ${\tau }_{0m}$ and ${\tau }_{0t}$ should be represented in terms of the individually calculated variables of pragmatic model. There are several ways of obtaining such an approximation, and we chose a version that gives the smallest errors. We found that the stress produced by subgrid waves can be approximated by relation:

$${\tau }_{wa}=0.4v_{*}^2\frac{A_0\tilde{\omega }+A_1}{\tilde{\omega }+A_2},$$

(59)

where $\tilde{\omega }=\omega /{\omega }_p$ (${\omega }_p$ is peak frequency), $v_{*}$ is the external friction velocity, while $A_0,A_1$ and $A_2$ are the functions of $v_{*}\left(m/s\right)$:

$$\begin{array}{l}{A}_0=-0.0440+0.0175v_{*}^{0.760}\\ A_1=-11.203+14.056v_{*}^{0.0543}\\ A_3=-0.319+2.649v_{*}^{0.2389}\end{array}$$

(60)

The tangent stress ${\tau }_{t0}$ is approximated by relation:

$${\tau }_{ta}=v_s^2\left(B_0+B_1{L}_{*}+B_2{L}_{*}^2\right),$$

(61)

where $L_{*}=\ln \left(v_{*}\right)$ and

$$\begin{array}{l}{B}_0=0.3953-0.1706{\mathsf{\Omega }}_p+0.0301{\mathsf{\Omega }}_p^2\\ B_1=-0.1043+0.0133{\mathsf{\Omega }}_p-0.0026{\mathsf{\Omega }}_p^2\\ B_2=0.0392-0.0037{\mathsf{\Omega }}_p+0.0007{\mathsf{\Omega }}_p^2\end{array}$$

(62)

The approximations (59)–(62) are valid in a range $\left(0.9<{\mathsf{\Omega }}_p<3.\right)$ and for the wind velocity $\left(3<U_{10}<35\right)$. The accuracy of the approximations (59)–(62) is demonstrated in Figure 14 where the values calculated with (55)–(58) are compared with the model calculated with high resolution. It should be underlined that parameterization (59)–(62) is based on the assumption that parameterized high-frequency stress is produced by local wind.

8. Dependence of Drag Coefficient on Wind Direction

The model described above was used for illustration of dependence of drag coefficient on wind direction. It seldom occurs that wave spectrum is symmetrical and is produced by constant wind. In reality, wind often rotates and starts to generate its own wave system and provide by energy the existing waves. In a most clear form it happens in tropical storms where waves are a product of a rapidly changing wind system.

For the calculations a simplified model with modest resolution was used: the number of frequencies $N=50$; maximum frequency is ${\omega }_m=6{\omega }_p$; the stretching coefficient ${\gamma }_{\omega }=1.03$; the angle resolution $\Delta \theta =2\pi /M,$ $M=90$. The WBL model was approximated at a grid with $L=20$ levels; the time step was estimated with Equation (48). The surface boundary conditions for the subgrid wave-produced momentum flux ${\tau }_{w0}$ and ${\tau }_{to}$ were calculated by Equations (58)–(62). The wave spectrum in a range $0<\omega <{\omega }_m$ was assigned by the JONSWAP approximation. It was assumed that at higher frequencies $\omega >{\omega }_m$ the wave spectrum is also described by the JONSWAP asymptotic, but a general direction of such high-frequency waves coincides with the wind direction. The calculations were made for three values of inverse wave age: ${\mathsf{\Omega }}_p=1,\mathsf{\Omega }=1.5$ and ${\mathsf{\Omega }}_p=2$; for three directions of wind: $\mathsf{\Psi }=0$ (the directions of wind and waves are the same, curves 1 in Figure 15); $\mathsf{\Psi }=\pi /2$ (vector of wind is perpendicular to waves which are treated as a multi-mode swell); merging curves 2; ${\mathsf{\Omega }}_p=\pi$ (wind is directed against swell, curves 3). Curves from Group 1 are close to the same curves in Figure 13. Case 3 corresponds to the waves moving mostly perpendicular to wind. In this case the drag coefficient is supported by subgrid stress only. This is why the drag coefficient is small and it is the same for different wave ages. The curves in Group 3 relate to the case of opposite wind (high-frequency waves are still directed along the wind). In this case, in calculation of $\mathsf{\Omega }$ the overflow is a sum of wind and phase velocities (the values of an imaginary part of $\beta$ for $\mathsf{\Omega }<0$ are larger than $\beta$ for $\left|\mathsf{\Omega }\right|$). As seen, the drag coefficient for opposite wind can reach very high values.

Figure 15 confirms that the instantaneous value of drag coefficient for any wind velocity is a product of wave spectrum. Another example of the strong influence of wave spectrum on drag coefficient is given in Figure 16. The calculations were done with a set of parameters used for the calculations shown in Figure 15.

As seen, drag coefficient can change significantly above a variable wave spectrum. All of the calculations in this Section were done up to the state of momentum balance (Equation (50)). In general case, it is incorrect, since the boundary layer can remain in nonstationary conditions. In our opinion, the time step should be comparable with the time steps for the wave and atmospheric model. In our calculations the explicit time scheme was used, since the rate of calculations did not matter. When a model is intended as a link between the atmosphere and ocean, it is reasonable to use an implicit time step. In such a mode the time spent for the WBL modeling is very small compared with the time for the calculation of atmospheric physics. Note that in the nonstationary WBL above complicated wavefield the connection between wind velocity and drag coefficient can vary in broad ranges.

9. Thoughts on Drag Coefficient at Strong Wind

As demonstrated in Figure 12, the current theory predicts a monotonic growth of drag coefficient with an increase of wind speed. However, there exist some data revealing that for a wind speed exceeding 25–30 ${\mathrm{m}\mathrm{s}}^{-1}$, the drag coefficient reaches an upper limit (Powell et al. 2003, Donelan et al. 2004), while as the wind speed increases further, the drag coefficient can decrease. Earlier such an effect was also noted when analyzing the tropical cyclone development (Emanuel, 1995). Currently, some attempts are made to explain such behavior of drag coefficient on the basis of a ‘droplet theory’ [43,62,63,64]. According to this theory, the drops generated as a result of the splitting of falling water volumes intensify the dissipation of turbulence, which causes reduction of drag coefficient. Such scheme looks reasonable, although a model of drop generation (as well as the mechanism of drops-turbulence interaction) uses too many arbitrary assumptions and ‘constants’, leaving a great possibility for fitting the theory to the observation data. In fact, this mechanics was suggested by Barenblatt and Golitsyn [65] for explanation of a dust storm phenomenon (and for liquid fuel droplets in rocket engines before that). According to their theory, sand particles increase the dissipation of turbulence and decrease friction, which accelerates the flow. The difference between the sea condition and desert is that the specific weight of sand particles is 3—4 times larger than that of water, and the concentration of sand particles is very high (visibility in dust storms is usually zero). It is easy to estimate that the density of sand fraction in a dust storm is by orders of magnitude higher than the density of water in a flow above sea waves.

According to our data, an effect of the wave-drag reduction at high wind speeds can be easily explained through the influence of high-frequency waves. As seen in Figure 10, about 50–70% of WPMF is transferred to waves at frequencies $\omega /{\omega }_p>10$. Hence, the variation of stress in a high-frequency part of the spectrum can modify the drag coefficient. It would be appropriate to suggest that the energy of short waves at high winds decreases due to two factors: the presence of bubble foam that suppresses short waves, and the flow deceleration in troughs due to the flow separation. The latter can be investigated using the coupled wind-wave model directly, though the model must take into account a broad range of wave spectrum from peak waves to capillary waves. A high wind speed and the necessity to use a high vertical resolution in the WBL make this problem quite time-consuming; however, such calculations are possible from a technical point of view. Note that the reduction of drag can also be caused by ‘blowing away’ sharp crests: high winds can smooth the surface and remove the elements responsible for stress. This effect was directly observed in the wind-wave tunnel (J. Troitskaya—private communication).

To prove that the reduction of drag coefficient is caused by the suppression of small waves, the additional calculations of drag coefficient were performed with a use of a non simplified WBL model described in Section 3 for the modified JONSWAP spectrum. A somewhat arbitrary assumption suggests that waves, whose frequency exceeds some value of ${\omega }_f$, are absent. In other words, we just cut the spectral tail in small pieces starting from $\omega ={\omega }_N$ and calculated the stationary structure of the WBL and drag coefficient. The calculations were started from $U_{10}=25m/s$ when a drag coefficient was equal to $2.6\cdot {10}^{-3}$. Then the calculations were repeated for the next wind value, with cutting a new portion of wave tail. The sizes of cuts were chosen for good qualitative approximation of data obtained on [51,55,66]. The final dependence $C_{10}\left(U_{10}\right)$ is shown by dots in the lower panel of Figure 17. In the upper panel a relative decline of total wave energy $\Delta E_p/E_p(\%)$ is shown. It is seen that a considerable decrease of drag can be reached by a very small reduction of energy in a high-frequency part of the spectrum. Note that the form of reduction does not influence general results: essentially the same dependence $C_{10}\left(U_{10}\right)$ can be obtained by increase of the rate of energy drop in the spectral tail.

It should be emphasized that the result presented in Figure 17 is purely qualitative, since the exact shape of spectrum at high wind speeds is unknown. Figure 17 just illustrates a simple explanation of the drag coefficient reduction at high wind speeds on the basis of the high-frequency wave spectrum modification. To take into account this effect, it is necessary to study the exact shape of spectrum tail, i.e., to find a connection between the wind velocity and energy in a high-frequency part of spectrum. This problem seems to be simpler than defining a set of ‘constants’ in the equations of the ‘droplet’ theory.

10. Conclusions

Contrary to the turbulent transfer, the wave-produced momentum flux (WPMF) is not an internal factor of the wave boundary layer (WBL), since the structure and magnitude of the WPMF are controlled by the wave field. Only the stationary and horizontally homogeneous flow above smooth or rough surface has the universal structure. Even above land such conditions are seldom, because in reality the surface is inhomogeneous, and the external flow is nonstationary. Above sea the variable wind and advection of wave energy produce infinite diversity of wave spectra. The dynamic structure of wave boundary layer considered in this paper is completely defined by upper boundary conditions (i.e., wind velocity or turbulent stress) and wave field properties. This is true if the wind or wave spectrum does not change; otherwise the structure of the WBL depends on a relatively short prehistory. All these processes can be well simulated by coupling the WBL 2-D model or 3-D model with turbulence closure or LES schemes [67]—with the phase resolving model of waves (like the HOS model [68] or the finite difference model [69] or model [70].

Spectral wave models [14], as well as 2-D [17] and 3-D LES [21] modeling [67], or modeling with a turbulent closure scheme [17] require the parameterization of tangent stress and subgrid waves effects. The scheme (59–62) was developed specifically for these purposes. A more practical approach, such as wave prediction modeling [14] or coupled ocean-waves-atmosphere modeling [71] is based on a spectral presentation of wave fields. Parameterizations of the dynamic WBL should solve two main tasks: (1) calculate the flux of momentum lost by the atmosphere model and obtained by the ocean; (2) calculate the flux of energy to waves. Both of these problems are far from being resolved. Until recently numerous attempts have been made to find the best approximation for dependence of surface stress on wind velocity (see [72]). However, it became clear many years ago that such a unique dependence cannot exist (e.g., [73]). The drag coefficient at any given wind velocity is a product of the wave spectrum which can be quite different. For example, let us consider an ideal wave field consisting of one single mode. If wind direction coincides with the direction of wave propagation, the wave will produce some wave drag which is superimposed upon tangent drag. If the wave front is parallel to the wind direction, the wave-produced drag is equal to zero, and friction is supported by the tangent stress only. If the direction of the wind is the opposite to wave propagation, the drag will be larger than in the first case, because the rate of overflow is a sum of wind and phase velocities. In a multimode wave field, waves of different directions are present, but the effect of deviation of wave direction from the wind direction remains important (see Figure 16). Naturally, the magnitude and direction of wave drag strongly depends on thea shape of the two-dimensional wave spectrum $S\left(\omega ,\theta \right)$.

Unfortunately, the data on wind vector and spectral shape cannot be converted directly into the magnitude of drag coefficient. For such purposes it is necessary to use a model of WBL that calculates the structure of the WBL and the rate of exchange by momentum and energy between wind and waves. The construction of the WBL model became possible following the investigation of $\beta$ function performed in [17] and the recent new post-processing of their results. For approximation (12), (13) of $\beta$-function 1,400,000 individual values in the range of $-50<\mathsf{\Omega }<50$ were used. The data on $\beta$ at $\left|\mathsf{\Omega }\right|>20$ had a large scatter and were considered as not quite reliable. Such values of $\mathsf{\Omega }$ are used infrequently. $\beta$-function was defined by the processing of the simultaneous local values of elevation and surface pressure. The pressure was calculated by a standard method based on the Poisson equation, and to some degree can be considered as exact. However, the results depend on the vertical and horizontal resolution of WBL, and the significance of this dependence is not clear enough. We were not able check the effect of this resolution because the calculations similar to those described in [17] normally take several months. Note that the entire approach demonstrated in the current paper requires further investigation of $\beta$-function calculated with a high-resolution coupled model. This problem is far from being solved; the current paper is just a first step in this direction.

The simplest way of derivation of the 1-D WBL model (Equations (21) and (23)) should be based on the complete Euler equations extended by an appropriate closure scheme. The averaging of the equation for horizontal momentum gives a 1-D equation that contains the components of WPMF. The vertical profiles of these components can be represented in a convenient form as a sum of fluxes produced by separate wave modes (see Figure 5). Finally, we obtain the equations which directly connect the structure of the WBL with an arbitrary wave spectrum. To avoid uncertainty with a lower boundary condition we used a very high resolution for frequency and height, which allowed us placing the lowest layer in a viscous sub-layer and calculate the tangent stress in the near-surface layer with the usage of molecular viscosity. The calculations showed that the results are not too sensitive to resolution. As an example, the solution of equations was done for the symmetrical JONSWAP spectrum. The structure of the WBL was illustrated by profiles of wind velocity and WPMF, as well as by energy of turbulence and some other characteristics. It is shown that the Charnock’s scaling decreases the dependence on wind velocity, but the coefficient in this formula still remains dependent on wind velocity and wave age. Finally, the dependence of drag coefficient on wind velocity and wave age was presented. Naturally, this result is valid only for the situations which are similar to the idealized conditions of the JONSWAP experiment. In reality, the general direction of the wind wave spectrum deviates from wind direction, and the spectrum can contain the waves not generated by local wind. Each spectrum at given velocity generates its own drag coefficient. If the conditions change quickly, the structure of WBL is not stationary, and the drag coefficient is not fully defined by the local spectrum and wind.

A detailed model that uses high resolution over wave frequency and elevation above the water surface cannot be used as an element of a complicated system similar to the wave prediction model or coupled ocean-atmosphere model. However, this model can be used for the development of a simplified approach which does not require considerable computational resources. The simplification of the model can be done by the sharp reduction of frequency range (engineering models like WAM or WAVEWATCH do not consider high-frequency waves). Parameterization is based on the evident properties of small waves: they are mostly generated by current wind; hence, their structure is more or less universal. We assume that an asymptotic form of wave spectrum at high frequencies is similar to that of the JONSWAP spectrum. It was shown before that the friction generated by subgrid waves and tangent stress can be parameterized in terms of cut-frequency, external friction velocity and wave age. The calculations of the WBL structure with a simplified model proved that such parameterization is acceptable.

The 1-D model of WBL is useful for the investigation of a broad range of problems. The main advantage of the model is a direct connection of the WBL structure with an arbitrary wave spectrum. In general cases, the WBL problem is nonstationary: the boundary layer is all the time adjusting to the changing wave field and wind velocity. If the model is used as a link in the atmosphere-waves-upper ocean system, the stationary solution is not always reachable, since the boundary conditions are typically unsteady.

This paper is devoted to the purely-dynamic interaction. The model will be generalized by including stratification effects and the equations for temperature and vapor transfer.