# “Choppy wave” model for nonlinear gravity waves

Frédéric Nouguier, Charles-Antoine Guérin, Bertrand Chapron

► To cite this version:

Frédéric Nouguier, Charles-Antoine Guérin, Bertrand Chapron. “Choppy wave” model for non-linear gravity waves. Journal of Geophysical Research. Oceans, 2009, 114 (C9), pp. C09012. 10.1029/2008JC004984. [hal-01345520]

HAL Id: hal-01345520

https://hal.science/hal-01345520

Submitted on 13 Jul 2016

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# “Choppy wave” model for nonlinear gravity waves

Frédéric Nouguier,<fnref n="1" /> Charles-Antoine Guérin,<fnref n="2" /> and Bertrand Chapron<fnref n="3" />

[1] We investigate the statistical properties of a three-dimensional simple and versatile model for weakly nonlinear gravity waves in infinite depth, referred to as the “choppy wave model” (CWM). This model is analytically tractable, numerically efficient, and robust to the inclusion of high frequencies. It is based on horizontal rather than vertical local displacement of a linear surface and is a priori not restricted to large wavelengths. Under the assumption of space and time stationarity, we establish the complete first- and second-order statistical properties of surface random elevations and slopes for long-crested as well as fully two-dimensional surfaces, and we provide some characteristics of the surface variation rate and frequency spectrum. We establish a relationship between the so-called “dressed spectrum,” which is the enriched wave number spectrum of the nonlinear surface, and the “undressed” one, which is the spectrum of the underlying linear surface. The obtained results compare favorably with other classical analytical nonlinear theories. The slope statistics are further found to exhibit non-Gaussian peakedness characteristics. Compared to observations, the measured non-Gaussian omnidirectional slope statistics can only be explained by non-Gaussian effects and are consistently approached by the CWM.

## 1. Introduction

[2] The development of fully consistent inversions of sea surface short wave characteristics via the ever increasing capabilities (radiometric precision, spatial resolution) of remote sensing measurements has considerably advanced. Yet, difficulties remain, mostly associated to stringent requirements to have adequate understandings and means to describe very precisely the sea surface statistical properties in relation to surface wave dynamics. The simplest linear superposition and Gaussian models remain in common use. Such models provide insight and are often accurate enough for many practical purposes. Yet, common visual inspections of natural ocean surface waves often reveal geometrical asymmetries. Namely, when the steepness of a wave locally increases, its crest becomes sharper and its trough flatter. Harmonic phase couplings occur, and an ocean surface wave field can become rapidly a non-Gaussian random process. For remote sensing applications and model developments, the statistical description of random nonlinear gravity waves is then certainly not straightforward, but must be taken into account to improve uses and interpretation of measurements, e.g., to correct for the sea state bias in altimetry, to explain the upwind/downwind asymmetry of

the radar cross section or to interpret the role of fast scatterer in Doppler spectra.

[3] As usually described, nonlinear surface gravity waves are generally prescribed in the context of the potential flow of an ideal fluid. For small wave steepness, the resulting nonlinear evolution equations can first been solved by means of a perturbation expansion [Tick, 1959]. This approach consists in finding iteratively a perturbative solution of the equations of motion for both the surface elevation and the velocity potential, by matching the boundary conditions at the bottom and at the free surface [Hasselmann, 1962; Longuet-Higgins, 1963; Weber and Barrick, 1977]. Following an other approach, Zakharov [1968] showed that the wave height and velocity potential evaluated on the free surface are canonically conjugate variables. This helps to uniquely formulate the water wave equations as a Hamiltonian system. For water waves, the Hamiltonian is the total energy E of the fluid. The Hamiltonian approach is based on operators expansions technique [Zakharov, 1968; Creamer et al., 1989; Watson and West, 1975; West et al., 1987; Fructus et al., 2005], albeit using truncated Hamiltonian. We refer to Elfouhaily [2000] for a comparison and discussion between the two approaches. For two-dimensional water waves, where the free surface evolves as a function of one variable in space, effective methods have been improved and include conformal mapping variables [Zakharov et al., 2002; Ruban, 2005; Chalikov and Sheinin, 2005]. A recent review on numerical methods for irrotational waves can be found in the paper by Dias and Bridges [2006]. For the three-dimensional problem, one loses the possibility to employ complex analysis, except to still consider a quasi-planar approximation, i.e., very long

<footnote n="1">¹ Institut Fresnel, UMR 6133, Faculté de Saint-Jérôme, Université Paul Cézanne, CNRS, Marseille, France.</footnote>

<footnote n="2">²LSEET, UMR 6017, Université du Sud-Toulon-Var, CNRS, La Garde, France.</footnote>

<footnote n="3">³Laboratoire d'Océanographie Spatiale, IFREMER, Plouzané, France.</footnote>
crested waves. Consequently, for the general problem, the perturbative technique has the advantage of simplicity, but remains essentially a low-frequency expansion and produces some nonphysical effects at higher frequencies, such as the divergence of the second-order spectrum. The Hamiltonian approach will be capable of handling stronger nonlinearities but is more tedious, remains essentially numerical and does not provide explicit statistical formulas. Finally, a Lagrangian description of surface wave motion may be more appropriate to describe steep waves [Chalikov and Sheinin, 2005]. In such a context, the Gerstner wave [Gerstner, 1809] is a first well-known exact solution for rotational waves in deep water, and Stokes [1847] derived a second-order Lagrangian approximation for irrotational waves leading to a well-known and observed net mass transport, the Stokes drift phenomenon, in the direction of the wave propagation.

[4] The aim of this paper is to build on this latter simplified phase perturbation methodology to propose a simple, versatile model, that can reproduce the lowest-order nonlinearity of the perturbative expansion but does not suffer from its related shortcomings. This analytical model is certainly not properly new, as it is widely used by the computer graphics community [Fournier and Reeves, 1986; J. Tessendorf, unpublished data, 2004] to produce real-time realistic looking sea surfaces. The terminology choppy wave model (henceforth abbreviated to “CWM”) originates from the visual effect imposed by the transformation compared to linear waves. In addition to gravity waves nonlinear interaction, the model can incorporate further physical features such as the horizontal skewness induced by wind action over the waves, an effect that we will not consider in this paper and which will be left for subsequent work.

[5] On the mathematical level, the model identifies completely with the perturbative expansion in Lagrangian coordinates as proposed four decades ago by Pierson [1962, 1961]. In the case of a single wave, it coincides with the Gerstner solution and is consistent with the Stokes expansion [Stokes, 1880] at third order in slope. Our present contribution is to provide a complete, nontrivial statistical study of this model and a comparison with the classical approaches. As understood, the CWM does not claim to compete with Hamiltonian-based methods and is in fact limited to the lowest-order nonlinearity. Its main strength is to provide a good compromise between simplicity, stability and accuracy. More precisely, it is (1) numerically efficient, as time evolving sample surfaces can be generated by FFT; (2) analytically tractable, as it provides explicit formulas for the first- and second-order point statistics; and (3) robust to the frequency regime, as it is found to be equivalent to the canonical approach [Creamer et al., 1989] at low frequencies while remaining stable at higher frequencies.

[6] In the following we have studied the two- and three-dimensional case pertaining to long-crested or truly two-dimensional sea surfaces, which from now on we will refer to as the 2-D and 3-D case. Since the methodology remains the same in both instances, we have chosen to give a complete exposure of the technique in the 2-D case which is considerably simpler. All the analytical results of the 2-D case (section 2) have their counterpart in the 3-D case (section 3). In the subsequent study, the emphasis will be put on the spatial properties of a “frozen” surface, even

through some temporal properties will also be discussed. Using a phase perturbation in the Fourier domain, the nonlinear local transformation simply consists in shifting the horizontal surface coordinates. Starting with a linear, reference surface, assumed to be a second-order Gaussian stationary process in space and time with given power spectrum, the complete first- and second-order properties of the resulting, non-Gaussian, random process is derived and related to the statistics of the reference surface. In particular, the resulting spectrum, which we refer to as dressed, has been related to the reference spectrum, termed undressed, in a way which is found to be very similar to Weber and Barrick’s [1977] and Creamer et al.’s [1989], but corrects the former and extends the latter to the 3-D case. As well, the sea surface slope statistical description is modified to exhibit a non-Gaussian behavior with a measurable peakedness effect, i.e., an excess of zero and steep slopes. A comparison with recent airborne laser measurements which allows to discriminate the slope statistics of gravity waves from smaller, short gravity, and capillary waves, is presented in section 4. As found, the CWM brings the excess kurtosis of omnidirectional slopes significantly closer to the data.

## 2. The 2-D Model

### 2.1. Definition

[7] As mentioned above, our goal is to assess the statistical properties of a nonlinear random process resulting from shifting horizontal coordinates. For a Gerstner wave in deep water only, the coordinates $(x, z)$ of particles at the free surface have the following parameterization in time $t$:

$$ 
\begin{aligned}
x &= x_0 - a \sin(kx_0 - \omega t) \\
z &= a \cos(kx_0 - \omega t),
\end{aligned}
$$

where the points $(x_0, 0)$ labels the undisturbed surface and $\omega = \sqrt{g|k|}$ satisfies the gravity-waves dispersion relation ($g = 9.81 \text{ m s}^{-2}$ is the acceleration due to gravity). This model introduces the horizontal displacement $D(x, t) = -a\sin(kx_0 - \omega t)$. At a given time, the locus of the points on the free surface describes a trochoïd. Strictly, this solution is not physical since it produces rotational motion. However, the vorticity is of order $(ka)^2$, and the solution is expected to be accurate for small slope parameter $ka \ll 1$. The obvious generalization to multiple waves writes,

$$ 
\begin{aligned}
x &= x_0 - \sum_j a_j \sin(k_j x_0 - \omega_j t + \phi_j) \\
z &= \sum_j a_j \cos(k_j x_0 - \omega_j t + \phi_j),
\end{aligned}
$$

where $\omega_j = \sqrt{g|k_j|}$ and $\phi_j$ are random phases. Such superpositions are known to be the solutions of the linearized equations of motion in Lagrangian coordinates [Pierson, 1961, 1962; Gjosund, 2003] but include effects that are nonlinear in the Eulerian formulation.

[8] Note that the horizontal displacement of the particles can be achieved through the Hilbert transform of the vertical coordinate, since this operation turns sine into cosine functions. Accordingly, the proposed nonlinear superposition also compares to the improved linear representation derived by Creamer et al. [1989]. A discrete or continuous
superposition can thus be realized though the following nonlinear transformation:

$$ (x, h(x, t)) \mapsto (x + D(x, t), h(x, t)) \quad (1) $$

where the horizontal displacement

$$ D(x, t) = \int_{-\infty}^{+\infty} dk \text{ i sign } (k) e^{ikx} \hat{h}(k, t) \quad (2) $$

is the Hilbert transform of the Gaussian elevation profile $\hat{h}(x, t)$. Here the function $\hat{h}$ is the spatial Fourier transform of the surface elevation:

$$ \hat{h}(k, t) = \frac{1}{2\pi} \int_{-\infty}^{+\infty} dx \, e^{-ikx} \hat{h}(x, t). \quad (3) $$

The relation

$$ \tilde{h}(x + D(x, t), t) = h(x, t) \quad (4) $$

implicitly defines a function $\tilde{h}$ for the displaced surface, provided the map $x \mapsto x - D(x, t)$ is one-to-one, an assumption that will be made systematically in the following. This is the case if the space derivative $D'$ remains smaller than one in magnitude.

## 2.2. Statistical Properties of the Space Process

[9] We will now study the spatial statistical properties of the displaced surface $\tilde{h}$ at a given time, say $t = 0$. The time dependence will from now on be omitted. The underlying reference surface $h(x)$ is assumed to be a stationary centered Gaussian process, which results from the summation of a sufficient number of free waves. Under this assumption, the process $D'(x)$ is again stationary and Gaussian, with the same variance as the slope process $h'$. Hence, the model is expected to hold for moderate slopes, for which the threshold $|D'| = 1$ is attained with exponentially small probability.

[10] We will denote $C$ and $\Gamma$ the spatial correlation function and spectrum of $h$, respectively:

$$ C(x) = \langle h(x)h(0) \rangle, \quad \Gamma(k) = \frac{1}{2\pi} \int_{-\infty}^{+\infty} e^{-ikx} C(x)dx \quad (5) $$

where the bracket denotes the ensemble average. Even though the function $\tilde{h}$ is not explicit, its first- and second-order statistical properties can be established analytically. In order not to go too far off the leading path of the paper, we have chosen to restrict the study to the statistical quantities which are truly needed for the scattering problem, namely the first- and second-order properties of the surface process. We will thus derive the distribution of elevations and slopes, together with the wave number spectrum. In the long-crested case, we will also provide the distribution of variation rate of elevation at a given location, together with the frequency spectrum. These quantities will be later compared with those derived from classical theoretical approaches. There is also an abundant literature of nonlinear wave amplitudes (crests, troughs or crest-to-trough amplitudes). Many studies have dealt with intercomparison of wave height distributions after approximate solutions, exact numerical models or experi-

mental measurements. A recent survey are given by Tayfun and Fedele [2007]. We will, however, not discuss here the height-amplitude distribution after the CWM, a study which is left for further investigation.

### 2.2.1. First-Order Properties

[11] The one-point characteristic function of the nonlinear surface is given by:

$$ \Phi(v) = \langle e^{iv\tilde{h}} \rangle. \quad (6) $$

Since the process $\tilde{h}$ is stationary, we can rewrite:

$$ \Phi(v) = \lim_{L \to \infty} \frac{1}{2L} \int_{-L}^{+L} \langle e^{iv\tilde{h}(x)} \rangle dx \quad (7) $$

Now, operating the change of variable $x \mapsto x + D(x)$ makes it possible to eliminate the implicit function $\tilde{h}$:

$$ \Phi(v) = \lim_{L \to \infty} \frac{1}{2L} \int_{-L}^{L} \langle e^{iv\tilde{h}(x)}(1 + D'(x)) \rangle dx \quad (8) $$

As $h$ and $D'$ are Gaussian stationary processes, the term in bracket can be easily evaluated yielding to:

$$ \Phi(v) = (1 - iv\sigma_1^2) \exp\left(-\frac{1}{2}v^2\sigma_0^2\right) \quad (9) $$

Here, we have introduced the absolute moments of the spectrum:

$$ \sigma_n^2 = \int_{-\infty}^{+\infty} |k^n| \Gamma(k) dk \quad (10) $$

Note that $\sigma_0^2$ and $\sigma_2^2$ are the mean squared height (msh) and mean squared slope (mss) parameters of the surface, respectively. A Fourier inversion of (9) provides the probability distribution function (pdf) of elevations:

$$ \tilde{P}_0(z) = P_0(z) \left(1 - \frac{\sigma_1^2}{\sigma_0^2} z\right) \quad (11) $$

where

$$ P_n(z) = \frac{1}{\sqrt{2\pi\sigma_n^2}} \exp\left(-\frac{z^2}{2\sigma_n^2}\right) \quad (12) $$

is the centered normal law with variance $\sigma_n^2$. The evaluation of the characteristic function together with its successive derivatives at the origin provides the first few moments and cumulants ($\tilde{\kappa}_n$) of the transformed process:

$$ \begin{aligned} \langle \tilde{h} \rangle &= -\sigma_1^2, & \langle \tilde{h}^2 \rangle &= \sigma_0^2, & \langle \tilde{h}^3 \rangle &= -3\sigma_0^2\sigma_1^2, & \langle \tilde{h}^4 \rangle &= 3\sigma_0^4 \\ \tilde{\kappa}_1 &= -\sigma_1^2, & \tilde{\kappa}_2 &= \sigma_0^2 - \sigma_1^4, & \tilde{\kappa}_3 &= -2\sigma_1^6, & \tilde{\kappa}_4 &= -6\sigma_1^8 + 3\kappa_2^2 \end{aligned} \quad (13) $$

[12] Starting with a zero-mean linear surface, the resulting nonlinear surface becomes a nonzero mean random non-
<figure><img src="image_2.png" /><figcaption>Figure 1. Comparison of the distribution of elevations for the CWM and the Tayfun distribution, for typical parameters $\nu = 0.1$ and $\sigma_0 = 0.5$ m. The Gaussian distribution is given as reference.</figcaption></figure>

Gaussian process. The corresponding skewness is slightly negative,

$$ \tilde{\lambda}_3 = \frac{\tilde{\kappa}_3}{\tilde{\kappa}_2^{3/2}} = \frac{-2\sigma_1^6}{(\sigma_0^2 - \sigma_1^4)^{3/2}}, \quad (14) $$

and the msh is slightly diminished. There is no, however, significant creation of kurtosis with respect to the Gaussian case:

$$ \tilde{\lambda}_4 = \frac{\tilde{\kappa}_4}{\tilde{\kappa}_2^2} = \frac{-6\sigma_1^8}{(\sigma_0^2 - \sigma_1^4)^2} + 3 \quad (15) $$

Hence the transformed surface is shifted toward negative values and skewed. This is natural since the transformation tends to sharpen the crests and enlarge the troughs, unbalancing the contribution of top and bottom points in favor of the latter. The obtained distribution of elevation (11) can be compared with the well-known Tayfun distribution for narrow spectra [Tayfun, 1980], rewritten for the pdf with our notations:

$$ P_{\text{Tayfun}}(z) = \frac{1}{\pi \sigma_0 \nu} \int_0^{+\infty} e^{-\frac{\xi^2}{2\nu^2}} \left( e^{-\frac{(1-C(\xi))^2}{2\nu^2}} + e^{-\frac{(1+C(\xi))^2}{2\nu^2}} \right) \frac{d\xi}{C(\xi)} \quad (16) $$

where $\nu = k_0\sigma_0$ is a small dimensionless parameter, $k_0$ is the central wave number of the narrow spectrum, and $C(\xi) = (1 + 2k_0z + \xi^2)^{1/2}$. For narrow spectra, note that $\sigma_1^2 \simeq k_0\sigma_0^2$, to that the CWM distribution (11) is also a function of the parameters $\nu$ and $k_0$. Figure 1 gives a comparison of the Gaussian reference distribution, the distribution arising (11) from the CWM and the Tayfun distribution for typical parameters $\nu = 0.1$ and $\sigma_0 = 0.5$ m. The CWM is extremely

close to the Tayfun distribution in the first standard deviation interval.

[13] Differentiating equation (4) provides an implicit definition of the slopes of transform process:

$$ \frac{d\tilde{h}}{dx}(x + D(x)) = \frac{h'(x)}{1 + D'(x)} \quad (17) $$

We have not been able to calculate explicitly the characteristic function of the slopes process. However, $h'(x)$ and $D'(x)$ are two independent random variables, and we can use a formula for the distribution of quotient to derive the pdf. of slopes as:

$$ \begin{aligned} \tilde{P}_2(z) &= \int_{-\infty}^{+\infty} dx \, x|x| P_2(zx) P_2(x-1) \\ &= \frac{e^{-1/2\sigma_2^2}}{\pi(1+z^2)^2} + \frac{1}{\sqrt{2\pi\sigma_2^2}} \frac{\sigma_2^2(1+z^2) + 1}{(1+z^2)^{5/2}} \\ & \quad \times \text{Erf}\left(\frac{1}{\sqrt{2\sigma_2^2(1+z^2)}}\right) \exp\left(-\frac{1}{2\sigma_2^2}\left(\frac{z^2}{1+z^2}\right)\right) \end{aligned} \quad (18) $$

where Erf is the error function. Note that this distribution is even. The transformed slopes are thus centered and unskewed. The fourth moment of this distribution is unbounded, making the tail of the distribution unrealistic. The distribution should thus be truncated beyond a given threshold value. It can be checked that this truncation has a negligible impact on the normalization of the distribution, since the steepest events are very rare. It is interesting to note that a very resembling expression was recently obtained for the distribution of slope at a level upcrossing in the framework of a similar Lagrangian model, the main ingredient of the proof being Rice's formula for level crossings [Aberg, 2007].

[14] To test the shape of the distribution, we can compare it with the classical Gram-Charlier expansion as used by Cox and Munk [1954] to analyze ocean glitter distribution,

$$ \tilde{P}_2(z) = \frac{1}{\sqrt{2\pi\sigma_2^2}} e^{-z^2/2\sigma_2^2} \times \left(1 + \frac{c_4}{24} \times \left(\frac{z^4}{\sigma_2^4} - 6\frac{z^2}{\sigma_2^2} + 3\right)\right) \quad (19) $$

Figure 2 shows the different distributions for a typical mss slope value $\tilde{\sigma}_2^2 = 0.03$. (Second-order moment of CWM and Cox and Munk slopes pdf's are set equals.) The agreement between a Gram-Charlier expansion and CWM is found excellent with a clear departure from the Gaussian distribution. For the chosen mss, the agreement is found with $c_4 = \lambda_4 - 3 \simeq 0.27$. This excess of kurtosis is comparable to Cox and Munk reported values. Since the CWM slopes excess kurtosis is unbounded, we truncated the CWM distribution at a realistic maximum slope value, here $z_{max} = 0.7$. As derived, this excess of kurtosis is a consequence of the geometrical wave profile asymmetries, but also on the implicit modulation of the shorter waves by much longer waves [e.g., Creamer et al., 1989, paragraph 5]. These
<figure><img src="image_3.png" /><figcaption>**Figure 2.** Slopes' distributions after CWM and Gram-Charlier expansion. The Gaussian distribution is given as reference.</figcaption></figure>

interactions can then further lead to an excess of kurtosis [Chapron et al., 2000].

### 2.2.2. Second-Order Properties

[15] The second-order statistical properties are completely characterized by the two-point characteristic function $\langle \exp(\mathrm{iv}_1h(x_1) + \mathrm{iv}_2h(x_2)) \rangle$. We did not find it possible to obtain the latter analytically. However, we can derive a related function, namely its one-dimensional Fourier transform on the diagonal:

$$ \Psi(u; v) = \int_{-\infty}^{+\infty} e^{iux} \left( \left\langle e^{iv(\tilde{h}(x) - \tilde{h}(0))} \right\rangle - \left\langle e^{iv\tilde{h}(x)} \right\rangle \left\langle e^{-iv\tilde{h}(0)} \right\rangle \right) dx. \quad (20) $$

Introducing the structure function:

$$ S_0(x) = \langle |h(x) - h(0)|^2 \rangle = 2[\sigma_0^2 - C(x)], \quad (21) $$

applying the change of variable as in equation (8) and using standard properties of Gaussian processes [Papoulis, 1965], we obtain:

$$ \Psi(u; v) = \int_{-\infty}^{+\infty} e^{iux} \left[ \exp\left(-\left(u^2 + v^2\right) \frac{S_0}{2}\right) \left[1 - 2iuC' - C'' - u^2C'^2 + \frac{1}{4}v^2S_1^2\right] - \exp\left(-\left(u^2 + v^2\right)\sigma_0^2\right) \left(1 + v^2\sigma_1^4\right) \right] dx \quad (22) $$

Here we have introduced the first and second derivative of the correlation function ($C'$ and $C''$, respectively) and the structure function:

$$ S_1(x) = 2[\sigma_1^2 - C_1(x)] \quad (23) $$

where $C_1$ is the so-called Gilbert transform of the correlation function.

$$ C_1(x) = \int_{-\infty}^{+\infty} dk |k| \Gamma(k) e^{ikx}. \quad (24) $$

[16] Now denote $\tilde{C}$ and $\tilde{\Gamma}$ the centered correlation function and spectrum of the nonlinear process,

$$ \begin{aligned} \tilde{C}(x) &= \langle \tilde{h}(x) \tilde{h}(0) \rangle - \langle \tilde{h}(x) \rangle \langle \tilde{h}(0) \rangle, \\ \tilde{\Gamma}(k) &= \frac{1}{2\pi} \int_{-\infty}^{+\infty} e^{-ikx} \tilde{C}(x) dx \end{aligned} \quad (25) $$

We will make use of the terminology introduced by Elfouhaily et al. [1999] and Soriano et al. [2006] to designate the quantities pertaining to the linear or transformed processes. The “measured,” “output,” or dressed spectrum denotes the spectrum which is actually measured experimentally on the true ocean surface, including nonlinearities ($\tilde{h}$). The “bare,” “input,” or undressed spectrum refers to the linear surface that underlies the nonlinear process ($\tilde{h}$). To be able to generate realistic nonlinear surfaces, it is important to have a relationship between dressed and undressed quantities. For this we observe that:

$$ \tilde{\Gamma}(k) = \frac{1}{4\pi} \left[ \frac{\partial^2 \Psi(k, v)}{\partial v^2} \right]_{v=0} \quad (26) $$

resulting in the following expression of the dressed spectrum:

$$ \tilde{\Gamma}(k) = \frac{1}{2\pi} \int_{-\infty}^{+\infty} dx e^{ikx} \left\{ e^{-k^2 \sigma_0^2} (\sigma_0^2 - \sigma_1^4) - e^{-\frac{1}{2}k^2 S_0} \left[ \frac{1}{2} S_0 (1 - 2ikC' - C'' - k^2 C'^2) - \frac{1}{4} S_1^2 \right] \right\} \quad (27) $$

The dependence in the space variable is implicit in the involved functions. The oscillating nature and the slow decay of the correlation function for sea spectra makes the numerical evaluation of the above integral challenging. However, the formula can be simplified considerably by investigating the different frequency regimes.

### 2.2.3. Low-Frequency Asymptotics

[17] For small values of $k\sigma_0$ the real exponentials arising in formula (27) can be linearized. Reordering the different terms in the expansion according to powers of $k\sigma_0$, and using the identity:

$$ k^2(\Gamma*\Gamma) = 2(k\Gamma)*(k\Gamma) + 2(\Gamma)*(k^2\Gamma), \quad (28) $$

we obtain at the lowest corrective order the following expression for the dressed spectrum:

$$ \begin{aligned} \tilde{\Gamma}(k) = \Gamma(k) &+ \frac{1}{2} \int_{-\infty}^{+\infty} dk' \Gamma(k') \{k^2 \Gamma(k-k') - 2k^2 \Gamma(k)\} \\ &+ \int_{-\infty}^{+\infty} dk' \Gamma(k') \Gamma(k-k') [|k'| |k-k'| - k'(k-k')] \\ &- 2 \int_{-\infty}^{+\infty} dk' |k'| |k| \Gamma(k) \Gamma(k') \end{aligned} \quad (29) $$

A consistency check can be performed by integrating this modified spectrum and comparing it with the variance of elevation found previously. Even though the above expansion holds at low frequency only, the variance of elevations is imposed by large scales and thus the comparison is
<figure><img src="image_4.png" /><figcaption>Figure 3. Dressed and undressed curvature for a Phillips spectrum at 7 m s⁻¹ versus wave number *k* (in decimal log-log scale). The different approximations of the dressed spectrum (low- and high-frequency regime) are shown in their respective range of validity. The high-frequency regime depends on the chosen cutoff *l* for the smaller gravity waves. HF, high frequency; LF, low frequency.</figcaption></figure>

meaningful. The integration leads after simple manipulations to:

$$ \int_{-\infty}^{+\infty} dk \, \tilde{\Gamma}(k) = \int_{-\infty}^{+\infty} dk \, \Gamma(k) - \left( \int_{-\infty}^{+\infty} dk \, |k|\Gamma(k) \right)^2, \quad (30) $$

which coincides with the variance of elevation predicted by (13). The low-frequency formula will be compared in section 4 with classical expansions of the literature.

### 2.2.4. High-Frequency Asymptotics

[18] For $k\sigma_0 \sim 1$, the integrand in (27) contribute mainly through their behavior around the origin. Since the functions $C'$ and $S_1$ vanish at zero and $\sigma_2^2 \ll 1$, $\sigma_1^4 \ll \sigma_0^2$, the dressed spectrum may be approximated by:

$$ \tilde{\Gamma}(k) = \frac{1}{2\pi} \int_{-\infty}^{+\infty} dx \, e^{ikx} \left[ e^{-k^2 \sigma_0^2} \sigma_0^2 - e^{-\frac{1}{2}k^2 S_0} \frac{1}{2} S_0 \right] \quad (31) $$

[19] Figure 3 shows an example of dressed spectrum for a Phillips undressed spectrum by a wind of 7 m s⁻¹ ($\Gamma(k) = 0.0025 |k|^{-3}$, $k_p < k < k_u$), where $k_p = 0.14$ rad m⁻¹ is the peak frequency and $k_u$ is a high-frequency cutoff for gravity waves. The transition between gravity and capillarity waves is not sharp and involves surface wavelengths *l* between 1 and 5 cm, so that the value of $k_u = 2\pi/l$ lies in the range from 125 to 630 rad m⁻¹.

[20] To highlight the difference with the undressed spectrum, it is the curvature $k^3\Gamma(k)$ which is plotted. The straight horizontal line corresponds to the constant curvature 0.025 of the linear surface. The different approximations of the dressed spectrum are shown in their respective range of validity. The dressed spectrum has an enhanced curvature of

1.5–2 dB at higher frequencies, with a peak depending on the chosen value of $k_u$. Note that the dressed spectrum is nonvanishing at $k_u$, which means that nonlinearities have added high-frequency components. The low-frequency formula is consistent with the high-frequency expansion but starts diverging around $k = 10$ rad m⁻¹.

## 2.3. Statistical Properties of the Time Process

### 2.3.1. First-Order Properties

[21] The technique which has been used to derive spatial first-order distribution functions can also be employed to obtain first-order statistical properties in the time domain, assuming the process is stationary in time. As an example we will derive the probability density function of the variation rate $\partial_t \tilde{h}$ at a given location (the superscript “t” refers to time-dependent quantities):

$$ \tilde{P}_2^t(\tau) = \frac{1}{2\pi} \int_{-\infty}^{+\infty} \langle e^{i\nu\partial_t\tilde{h}(x,t)} \rangle e^{-i\nu\tau} d\nu \quad (32) $$

For one-sided time spectra, that is for waves traveling in one single direction, say to the right, we might write the time process in the form:

$$ h(x, t) = \int_{-\infty}^{+\infty} dk \, e^{i(kx - \operatorname{sign}(k)\omega t)} \hat{h}(k), \quad (33) $$

where as usual $\omega = \sqrt{g|k|}$. Using the same spatial averaging as in (8), we obtain in a similar way as previously the following distribution of the variation rate:

$$ \tilde{P}_2^t(\tau) = \int_{\mathbb{R}^2} |x_2|(1+x_1)P_2^t(\tau x_2)G(x_1, x_2)dx_1dx_2, \quad (34) $$
<figure><img src="image_5.png" /><figcaption>Figure 4. Distribution of variation rate $\partial_t \tilde{h}$ after CWM transformation. The Gaussian distribution of the underlying linear surface is given as reference (only the positive part of the symmetric distribution is shown).</figcaption></figure>

where $P_2^t$ is the Gaussian distribution of $\partial_t h$ (mss: $g\sigma_1^2$) and $G$ is the bivariate Gaussian distribution with covariance matrix:

$$N = \begin{bmatrix} \sigma_2^2 & -\sqrt{g}\sigma_{3/2}^2 \\ -\sqrt{g}\sigma_{3/2}^2 & g\sigma_1^2 \end{bmatrix}$$

We have not been able to push the analytical calculation further, but this double integral can be easily implemented numerically. It is parameterized by the absolute moments $\sigma_n$ of the k spectrum. Figure 4 shows the time-slope probability density function for a 10 m s⁻¹ wind speed, omnidirectional fully developed Elfouhaily spectrum. A comparison is given with the corresponding Gaussian distribution $P_2^t$. The slope distribution $\tilde{P}_2^t$ of the CWM can be fitted with a striking accuracy with the Cauchy distribution, $p(\tau) = 1/\pi(1 + \tau^2)$. Hence, it has a slow $\tau^{-2}$ decay at large arguments, thereby rendering quite probable the occurrence of very large slopes. A recent experimental study by Joelson and Néel [2008] has shown that the distribution of variation rate measured in a tank can actually be well fitted by heavy tail distributions such as alpha-stable laws.

### 2.3.2. Second-Order Properties

[22] Even though the surface evolution is essentially governed by the linear dispersion relationship, the occurrence of nonlinear interactions alters the latter. The main contribution to the time spectrum energy at frequency $\omega$ is due to the waves of length $\lambda = 2\pi/k = 2\pi g/\omega^2$, but different scales are also involved through nonlinear effects. The simple example of the Gerstner wave is illuminating in that respect, as it contains different spatial scales (the bound waves) traveling with the same velocity. Hence, nonlinearities render the definition of the dispersion relationship ambiguous. However, experimental set-ups very often record surface elevation variation with time at a prescribed location and therefore provide estimation of time domain spectra. As mentioned, the presence of nonlinearities makes the link with space domain spectra difficult. Hence, incorporating time domain spectra in the model is highly desirable.

We have not, however, been able to fully mimic the procedure that was adopted for space spectra. However, the time spectrum can be easily estimated under the assumption of small displacements, an hypothesis which is valid if the surface has only low-frequency components. Denote $C^t$ and $\Gamma^t$, respectively, the temporal correlation function and frequency spectrum,

$$C^t(\tau) = \langle h(0, \tau)h(0, 0) \rangle, \quad \Gamma^t(\Omega) = \frac{1}{2\pi} \int_{-\infty}^{+\infty} e^{-i\Omega\tau} C^t(\tau)d\tau, \qquad (35)$$

as well as their nonlinear counterparts $\tilde{C}^t$ and $\tilde{\Gamma}^t$:

$$\begin{aligned} \tilde{C}^t(\tau) &= \langle \tilde{h}(0, \tau) \tilde{h}(0, 0) \rangle - (\langle \tilde{h}(0, 0) \rangle)^2, \\ \tilde{\Gamma}^t(\Omega) &= \frac{1}{2\pi} \int_{-\infty}^{+\infty} e^{-i\Omega\tau} \tilde{C}^t(\tau) d\tau \end{aligned} \qquad (35')$$

[23] For small displacements $D$, we may approximate

$$\tilde{h}(x, t) \approx h(x, t) - D(x, t) \partial_x h(x, t), \qquad (36)$$

which entails:

$$\tilde{C}^t = C^t \left( 1 + \frac{1}{g^2} \partial_{\tau}^4 C^t \right) + \frac{1}{g^2} \left( \partial_{\tau}^2 C^t \right)^2, \qquad (37)$$

or, equivalently, in the frequency domain ($*$ is the convolution in time):

$$\tilde{\Gamma}^t = \Gamma^t + \frac{1}{g^2} [\Gamma^t * (\Omega^4 \Gamma^t) + (\Omega^2 \Gamma^t) * (\Omega^2 \Gamma^t)]. \qquad (38)$$

The frequency dressed spectrum thus enjoys a similar relationship as the wave number spectrum in the low-frequency regime (29).

## 3. The 3-D Model

[24] In the 3-D case, Pierson [1961] has provided the solution of the linearized equations of motion for an inviscid irrotational fluid in Lagrangian coordinates. In deep water, the particle positions at the free surface have following parameterization:

$$\begin{aligned} x &= x_0 - \sum_j a_j \hat{\mathbf{k}}_j \cdot \mathbf{x}_0 \sin(\mathbf{k}_j \cdot \mathbf{r}_0 - \omega_j t + \phi_j) \\ y &= y_0 - \sum_j a_j \hat{\mathbf{k}}_j \cdot \mathbf{y}_0 \sin(\mathbf{k}_j \cdot \mathbf{r}_0 - \omega_j t + \phi_j) \\ z &= \sum_j a_j \cos(\mathbf{k}_j \cdot \mathbf{r}_0 - \omega_j t + \phi_j), \end{aligned}$$

where $\hat{\mathbf{k}}_j$ is a two-dimensional vector, $\hat{\mathbf{k}}_j = \mathbf{k}_j/|\mathbf{k}_j|$ and $\mathbf{r}_0 = (x_0, y_0)$ labels the particles at rest on the flat surface. Similarly to the 2-D case, the corresponding surface can be realized through horizontal displacements of a reference, linear, surface:

$$(\mathbf{r}, h(\mathbf{r}, t)) \mapsto (\mathbf{r} + \mathbf{D}(\mathbf{r}, t), h(\mathbf{r}, t)) \qquad (39)$$

where $\mathbf{r} = (x, y)$ is the horizontal coordinate. The function

$$\mathbf{D}(\mathbf{r}, t) = i \int_{\mathbb{R}^2} e^{i\mathbf{k} \cdot \mathbf{r}} \hat{h}(\mathbf{k}, t) \hat{\mathbf{k}} \, d\mathbf{k} \qquad (40)$$
<figure><img src="image_6.png" /><figcaption>**Figure 5.** The pdf of elevations of a 3-D linear (h) and CWM (h̃) surfaces for a wind of 10 m s⁻¹ (Elfouhaily spectrum). The nonlinear height distribution is shifted toward negative values.</figcaption></figure>

is the so-called Riesz transform of the function h, and

$$ \hat{h}(\mathbf{k}, t) = \frac{1}{(2\pi)^2} \int_{\mathbb{R}^2} \mathrm{d}\mathbf{r} \, e^{i\mathbf{k} \cdot \mathbf{r}} h(\mathbf{r}, t) \quad (41) $$

is its two-dimensional spatial Fourier transform.

### 3.1. First-Order Properties of the Space Process

[25] The calculations are similar to the 2-D case, although more involved.

[26] Let us introduce the partial and total absolute moments of the spectrum:

$$ \sigma_{\alpha\beta}^2 = \int_{\mathbb{R}^2} \frac{|k_x|^\alpha |k_y|^\beta}{|\mathbf{k}|^\gamma} \Gamma(\mathbf{k}) \mathrm{d}\mathbf{k}, \quad \sigma_n^2 = \int_{\mathbb{R}^2} |\mathbf{k}|^n \Gamma(\mathbf{k}) \mathrm{d}\mathbf{k}, \quad (42) $$

[27] Standard calculations lead to the following expression for the characteristic function (43) of elevations:

$$ \Phi(v) = (1 - iv\sigma_1^2 + v^2\Sigma_1) \exp\left(-\frac{1}{2}v^2\sigma_0^2\right) \quad (43) $$

with $\Sigma_1 = \sigma_{111}^4 - \sigma_{201}^2\sigma_{021}^2$. An example of deviation from the normal distribution is shown on Figure 5 for an input linear surface with directional Elfouhaily spectrum [Elfouhaily et al., 1997] at 10 m s⁻¹ wind speed.

[28] From the characteristic function, the following moments are easily obtained:

$$ \begin{aligned} \langle \tilde{h} \rangle &= -\sigma_1^2, & \langle \tilde{h}^2 \rangle &= \sigma_0^2 - 2\Sigma_1 \\ \langle \tilde{h}^3 \rangle &= -3\sigma_0^2\sigma_1^2, & \langle \tilde{h}^4 \rangle &= 3\sigma_0^4 \left(1 - 4\frac{\Sigma_1}{\sigma_0^2}\right), \end{aligned} \quad (44) $$

as well as the pdf of elevations:

$$ \tilde{P}_0(z) = P_0(z) \left( 1 + \frac{\Sigma_1}{\sigma_0^2} - \frac{\sigma_1^2}{\sigma_0^2}z - \frac{\Sigma_1}{\sigma_0^4}z^2 \right) \quad (45) $$

where as before $P_0$ is the Gaussian pdf of the linear surface.

[29] We can also derive the skewness ($\tilde{\lambda}_3$) and the kurtosis ($\tilde{\lambda}_4$) of elevation. The respective values for isotropic spectra are given in parenthesis:

$$ \begin{aligned} \tilde{\lambda}_3 &= \frac{-2\sigma_1^2[\sigma_1^4 + 3\Sigma_1]}{[\sigma_0^2 - (\sigma_1^4 + 2\Sigma_1)]^{3/2}} \quad \left( = -\frac{\sigma_1^6}{2} \left( \sigma_0^2 - \frac{\sigma_1^4}{2} \right)^{-3/2} \right) \\ \tilde{\lambda}_4 &= 3 \frac{\sigma_0^4 \left( 1 - 4\frac{\Sigma_1}{\sigma_0^2} \right) - \sigma_1^4 (\sigma_1^4 + 4\Sigma_1 + 2\sigma_0^2)}{(\sigma_0^2 - (\sigma_1^4 + 2\Sigma_1))^2} \quad (46) \\ & \quad \cdot \left( = 3\sigma_0^2(\sigma_0^2 - \sigma_1^4) \left( \sigma_0^2 - \frac{\sigma_1^4}{2} \right)^{-2} \right) \end{aligned} $$

[30] Again, there is a negative skewness and a positive excess of kurtosis, and the msh is diminished by a negligible amount.

[31] We have not been able to calculate explicitly the pdf of slopes $\tilde{P}_2(\mathbf{z})$. However, we could establish the following integral representation, which can be estimated numerically:

$$ \tilde{P}_2(\mathbf{z}) = \int_{\mathbb{R}^3} \frac{\operatorname{sign}(|\mathbf{J}|)|\mathbf{J}|^2}{(2\pi)^{5/2} \sqrt{|S_1|} \sqrt{|S_2|}} \exp\left\{-\frac{1}{2}X^T S_1^{-1} X\right\} \cdot \exp\left\{-\frac{1}{2}\mathbf{z}^T J^T S_2^{-1} J \mathbf{z}\right\} dX \quad (47) $$

with

$$ \begin{aligned} S_1 &= \begin{bmatrix} \sigma_{402}^2 & \sigma_{312}^2 & \sigma_{222}^2 \\ \sigma_{312}^2 & \sigma_{222}^2 & \sigma_{132}^2 \\ \sigma_{222}^2 & \sigma_{132}^2 & \sigma_{042}^2 \end{bmatrix}, & S_2 &= \begin{bmatrix} \sigma_{200}^2 & \sigma_{110}^2 \\ \sigma_{110}^2 & \sigma_{020}^2 \end{bmatrix}, \\ X &= \begin{bmatrix} x_1 \\ x_2 \\ x_3 \end{bmatrix}, & J &= \begin{bmatrix} 1 + x_1 & x_2 \\ x_2 & 1 + x_3 \end{bmatrix} \end{aligned} $$

where $|M|$ denote the determinant of the matrix M. Figure 6 displays the pdf of slopes in the upwind and crosswind direction for a directional Elfouhaily spectrum [Elfouhaily et al., 1997]. A comparison is given with the associated Gaussian distribution. The tail of the distribution decreases slower for the CWM and is significantly higher than the Gaussian tail for slope magnitudes beyond 0.5. Again, the slopes larger than some threshold (about 0.7) are not physical and the distribution must be truncated beyond this value.

### 3.2. Second-Order Properties of the Space Process

[32] As in the 2-D case, we can derive the two-dimensional Fourier transform of the two-point characteristic function on the diagonal, namely:

$$ \Psi(\mathbf{u}; v) = \int_{\mathbb{R}^2} e^{i\mathbf{u} \cdot \mathbf{r}} \left( \left\langle e^{iv(\tilde{h}(\mathbf{r}) - \tilde{h}(0))} \right\rangle - \left\langle e^{iv\tilde{h}(\mathbf{r})} \right\rangle \left\langle e^{-iv\tilde{h}(0)} \right\rangle \right) \mathrm{d}\mathbf{r}. \quad (48) $$
<figure><img src="image_7.png" /><figcaption>**Figure 6.** Distribution of slopes for linear and CWM surfaces at a wind of 12 m s⁻¹ (Elfouhaily spectrum). The logarithm is plotted to highlight the difference at large arguments.</figcaption></figure>

Operating the change of variable **r** $\mapsto$ **r** + **D**(**r**) we obtain:

$$ \Psi(\mathbf{u}; v) = \int_{\mathbb{R}^2} \mathbf{dr} e^{i\mathbf{u} \cdot \mathbf{r}} \left[ - \left| \left\langle e^{iv(h(\mathbf{r})+i\mathbf{u} \cdot \mathbf{D}(\mathbf{r}))} | J(\mathbf{r}) \rangle \right\rangle \right|^2 \right. \\ \left. + \left\langle e^{iv(h(\mathbf{r})-h(\mathbf{0}))+i\mathbf{u} \cdot (\mathbf{D}(\mathbf{r})-\mathbf{D}(\mathbf{0}))} | J(\mathbf{r}) | J(\mathbf{0}) \rangle \right\rangle \right] \quad (49) $$

[33] Here, J is the Jacobian matrix:

$$ J(\mathbf{r}) = \begin{bmatrix} 1 + \partial_x D_x(\mathbf{r}) & \partial_x D_y(\mathbf{r}) \\ \partial_y D_x(\mathbf{r}) & 1 + \partial_y D_y(\mathbf{r}) \end{bmatrix} $$

Discarding the quadratic terms in the Jacobian,

$$ |J(\mathbf{r})| |J(\mathbf{0})| \approx (1 + \nabla \cdot \mathbf{D}(\mathbf{r}))(1 + \nabla \cdot \mathbf{D}(\mathbf{0})) \quad (50) $$

and using standard properties of Gaussian processes [e.g., Papoulis, 1965] we obtain after tedious but straightforward calculations the following expression for the functional Ψ:

$$ \Psi(\mathbf{u}; v) = \int_{\mathbb{R}^2} \mathbf{dr} e^{i\mathbf{u} \cdot \mathbf{r}} \left[ \exp \left( -\frac{v^2}{2} S_0 - \frac{|\mathbf{u}|^2}{2} \Phi_{\mathbf{u}} \right) \right. \\ \left. \times \left( 1 - 2i\mathbf{u} \cdot \nabla C - (\mathbf{u} \cdot \nabla C)^2 - \Delta C + \frac{v^2}{4} S_1^2 \right) \right. \\ \left. - (1 + v^2 \sigma_1^4) e^{-v^2 \sigma_0^2} e^{-|\mathbf{u}|^2 \sigma_u^2} \right] \quad (51) $$

Here $\nabla C$ and $\Delta C$ are the gradient and the Laplacian of the correlation function, respectively. The dependence in the space variable is implicit. The auxiliary functions $\Phi_u$ and $S_1$ are defined by:

$$ S_1(\mathbf{r}) = 2 \int_{\mathbb{R}^2} \mathbf{dk}' |\mathbf{k}'| \Gamma(\mathbf{k}') [1 - e^{i\mathbf{k}' \cdot \mathbf{r}}] \quad (52) $$

$$ \Phi_{\mathbf{u}}(\mathbf{r}) = 2 \int_{\mathbb{R}^2} \mathbf{dk}' (\hat{\mathbf{u}} \cdot \hat{\mathbf{k}}')^2 \Gamma(\mathbf{k}') [1 - e^{i\mathbf{k}' \cdot \mathbf{r}}] \quad (53) $$

$$ \sigma_{\mathbf{u}}^2 = \int_{\mathbb{R}^2} \mathbf{dk}' (\hat{\mathbf{u}} \cdot \hat{\mathbf{k}}')^2 \Gamma(\mathbf{k}') \quad (54) $$

[34] Using the same technique as in 2-D we obtain for the dressed spectrum:

$$ \begin{aligned} \bar{\Gamma}(\mathbf{k}) = & \frac{1}{(2\pi)^2} \int_{\mathbb{R}^2} \mathbf{dr} e^{i\mathbf{k} \cdot \mathbf{r}} \left\{ -(\sigma_1^4 - \sigma_0^2) e^{-|\mathbf{k}|^2 \sigma_{\mathbf{k}}^2} \right. \\ & \left. + \left[ \frac{-S_0}{2} (1 - \Delta C - 2i\mathbf{k} \cdot \nabla C - (\mathbf{k} \cdot \nabla C)^2) + \frac{S_1^2}{4} \right] e^{-\frac{|\mathbf{k}|^2}{2} \phi_{\mathbf{k}}} \right\} \end{aligned} \quad (55) $$

The calculation of the low-frequency expansion of the dressed spectrum is similar to the 2-D case, leading to:

$$ \begin{aligned} \bar{\Gamma}(\mathbf{k}) = & \Gamma(\mathbf{k}) \\ & + \frac{1}{2} |\mathbf{k}|^2 \int \mathbf{dk}' \Gamma(\mathbf{k}') \left\{ \Gamma(\mathbf{k}'') - 2 (\hat{\mathbf{k}} \cdot \hat{\mathbf{k}}')^2 \Gamma(\mathbf{k}) \right\} \\ & + \int \mathbf{dk}' |\mathbf{k}'| \Gamma(\mathbf{k}') \Gamma(\mathbf{k}'') |\mathbf{k}''| \left[ (\hat{\mathbf{k}}' \cdot \hat{\mathbf{k}}'')^2 - \hat{\mathbf{k}}' \cdot \hat{\mathbf{k}}'' \right] \\ & - 2 \int \mathbf{dk}' |\mathbf{k}'| |\mathbf{k}'| (\hat{\mathbf{k}} \cdot \hat{\mathbf{k}}')^2 \Gamma(\mathbf{k}') \Gamma(\mathbf{k}), \end{aligned} \quad (56) $$

with $\mathbf{k}'' = \mathbf{k} - \mathbf{k}'$.

### 3.3. Undressing the Spectrum

[35] As can be seen on Figure 3, the dressed spectrum has an enhanced curvature with respect to the undressed one. This is natural since the inclusion of bound waves enriches the high-frequency content of the spectrum. In the 3-D case, we might also expect an enhancement of the spreading function at high frequencies through the nonlinear interaction of strongly directive long waves and weakly directive short waves. Now it is the dressed spectrum which is measured experimentally. To generate a nonlinear surface with a preassigned spectrum, it is thus necessary to go through an undressing procedure of the latter. The CWM transformation of the linear, fictitious, surface with undressed spectrum will eventually produce a nonlinear surface with suitable dressed spectrum.

[36] Soriano et al. [2006] introduced a simple undressing method assuming a power law form of the high-frequency part of the undressed spectrum. The parameters were fitted in such a way that the dressed spectrum leads to the correct values of the mean square height and slope after the nonlinear transformation proposed by Creamer et al. [1989]. Elfouhaily et al. [1999] used another method to retrieve the lowest-order cumulants of the nonlinear surface.

[37] The equation (55) can be incorporated in a simple iterative procedure to undress a spectrum with prescribed curvature and spreading function $\tilde{B}_{\text{target}}$, $\tilde{\Delta}_{\text{target}}$. Assuming a
<figure><img src="image_8.png" /><figcaption>Figure 7. Undressed and dressed curvature and spreading functions after two iterations (Elfouhaily spectrum).</figcaption></figure>

second harmonic azimuthal expansion of the dressed and undressed spectra:

$$2\pi k^4 \Gamma(\mathbf{k}) = B(\mathbf{k})(1 + \Delta(k) \cos(2(\phi_{\mathbf{k}} - \phi_{\text{wind}}))), \quad (57)$$

the iterative procedure to find undressed curvature (B) and spreading (Δ) functions runs as follows:

$$B^{(n+1)} = B^{(n)} - dB^{(n)}, \quad dB^{(n)} = \tilde{B}^{(n)} - \tilde{B}_{\text{target}} \quad (58)$$

$$\Delta^{(n+1)} = \Delta^{(n)} - d\Delta^{(n)}, \quad d\Delta^{(n)} = \tilde{\Delta}^{(n)} - \tilde{\Delta}_{\text{target}} \quad (59)$$

with $B^{(0)} = \tilde{B}_{\text{target}}$ and $\Delta^{(0)} = \tilde{\Delta}_{\text{target}}$. As an example, Figure 7 shows the first few iterates for a fully developed Elfouhaily dressed spectrum by a $U_{10} = 11 \text{ m s}^{-1}$ wind.

## 3.4. Numerical Surface Generation

### 3.4.1. Frozen Surface

[38] Sample nonlinear surfaces at a given time can be generated efficiently at the cost of three successive two-dimensional fast Fourier transforms: one for the spectral representation of the linear surface $h(\mathbf{r})$ and the other two for its Riesz transform $\mathbf{D}(\mathbf{r})$:

$$h(\mathbf{r}_{mn}) = \Re \sum_{ij} e^{i\mathbf{k}_{ij} \cdot \mathbf{r}_{mn}} \sqrt{\Gamma(\mathbf{k}_{ij})} e^{i\varphi_{ij}}, \\ \mathbf{D}(\mathbf{r}_{mn}) = \Re \sum_{ij} i \frac{\mathbf{k}_{ij}}{|\mathbf{k}_{ij}|} e^{i\mathbf{k}_{ij} \cdot \mathbf{r}_{mn}} \sqrt{\Gamma(\mathbf{k}_{ij})} e^{i\varphi_{ij}}, \quad (60)$$

where $\Gamma$ is the prescribed spectrum and $\varphi_{ij}$ are random uniform and independent phases on $[0, 2\pi]$. The nonlinear surface is parameterized by the points $(\mathbf{r}_{mn} + \mathbf{D}(\mathbf{r}_{mn}), h(\mathbf{r}_{mn}))$.

<figure><img src="image_9.png" /><figcaption>Figure 8. The 512 × 512 point linear sea surface with $U_{10} = 15 \text{ m s}^{-1}$.</figcaption></figure>

An example is given on Figures 8 and 9 with an Elfouhaily directional spectrum at wind $U_{10} = 15 \text{ m s}^{-1}$. It is a 3 m × 3 m patch of a total 50 m × 50 m sea surface generated with 8192 × 8192 points. The spectrum has been truncated at $k_{\text{max}} = 500 \text{ rad m}^{-1}$, corresponding to a minimal surface wavelength of 1 cm, and the surface is sampled regularly at the Shannon frequency $2k_{\text{max}}$. Since the operation is based on abscissa displacements, the resulting surface is given on a non regular grid. As can be seen on the encircled region of the plot, the crests of the CWM are sharpened while those of the linear surface are smoother.

### 3.4.2. Time Evolution

[39] The time evolution of nonlinear surface is a challenging issue, essentially because of the absence of a simple and well-defined dispersion relation. However, CWM surface at a given time is obtained by the same transformation (1) of a time-dependent linear surface. Therefore, it suffices to let the reference linear surface evolve and to perform the local transformation at the current time. If we, in addition, assume a fully developed time-independent spectrum, the evolution of the linear surface is simply obtained by use of the

<figure><img src="image_10.png" /><figcaption>Figure 9. The 512 × 512 point corresponding to CWM sea surface with $U_{10} = 15 \text{ m s}^{-1}$.</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>Figure 10. Time evolution of a 2-D sea surface. $\Delta t = 0.1$ s.</figcaption></figure>

gravity-waves dispersion relation $\omega^2 = g|\mathbf{k}|$, and amounts to change the original phases $\varphi_{ij}$ by an additional factor $-\omega_{ij}t = \pm\sqrt{g|\mathbf{k}_{ij}|}t$ in the FFT (60), depending on the travel direction of the waves.

[40] Figure 10 exemplifies the time evolution of a 2-D linear surface with one-sided time spectrum and the corresponding nonlinear surface. The undressed spectrum was chosen to be the fully developed Elfouhaily spectrum by a wind $U_{10} = 3$ m s$^{-1}$. The sample surface was taken to be 64 m long, with extreme frequencies $k_{\min} = 1.10^{-3}$ rad m$^{-1}$ and $k_{\max} = 100$ rad m$^{-1}$ and a sampling of 4096 points. The evolution of 4 m patch is represented for both the linear and CWM surface, with a time step $\Delta t = 0.1$ s.

## 4. Comparison With Classical Nonlinear Theories

[41] As mentioned in the Introduction, a certain number of fully nonlinear and numerically efficient solutions of potential flows have been developed in recent years. This makes it, in principle, possible to validate approximate theories. In practice, comparing the latter with various exact numerical solutions raises some difficulties, such as the lack of control of the final spectrum in an evolving nonlinear solution, the sensibility to the initial state or the relevance of sample surfaces comparisons. However, fast numerical schemes now allow the derivation of statistical properties of the surface through the use of extensive Monte-Carlo computations, especially for one-dimensional surfaces [e.g., Chalikov, 2005; Toffoli et al., 2008]. Nevertheless, going through a validation procedure by systematic comparisons of relevant statistical quantities is an important work which goes far beyond the scope of this paper and is left for further investigation.

### 4.1. Stokes Expansion

[42] A perturbative expansion of the implicit function $\tilde{h}$ can be obtained in the case of small displacements. We will here limit the discussion to the 2-D case. Supposed that the profile is obtained by dilation of a single dimensionless template $h_0$:

$$h(x) = ah_0(Kx), \quad D(x) = aD_0(Kx), \qquad (61)$$

where $a$ and $s = Ka$ are height and slope parameters, respectively. Then easy algebra leads to the following expansion, correct at second order in slope

$$\tilde{h} = h - Dh' + DD'h' + \frac{1}{2}D^2h'' \qquad (62)$$

[43] In the case of a single wave $h(x) = -a\cos(Kx)$, this perturbative series can be compared with a Stokes expansion. After rearrangement of the different terms in (62) we obtain:

$$\tilde{h}(x) = a \left( -\frac{1}{2}s + \left[ 1 - \frac{3}{8}s^2 \right] \cos(Kx) - \frac{1}{2}s \cos(2Kx) + \frac{3}{8}s^2 \cos(3Kx) \right), \qquad (63)$$

which coincides with a Stokes expansion at third order in slope. Note, however, that the CWM is more general than a mere superposition of Stokes waves, as frequency and phase coupling between the different modes comes in play through the nonlinear terms of the spatial expansion.

### 4.2. Longuet-Higgins Theory

[44] The classical approach [Hasselmann, 1962; Longuet-Higgins, 1963] to the nonlinear theory of gravity waves is to seek both the elevation $h$ and velocity potential $\Phi$ in a perturbation series,

$$\tilde{h}(\mathbf{r}, t) = h^{(1)}(\mathbf{r}, t) + h^{(2)}(\mathbf{r}, t) + \dots \qquad (64)$$

$$\Phi(\mathbf{r}, t) = \Phi^{(1)}(\mathbf{r}, t) + \Phi^{(2)}(\mathbf{r}, t) + \dots, \qquad (65)$$

where the first terms are given by the linear spectral representation of a Gaussian process,

$$\tilde{h}^{(1)}(\mathbf{r}, t) = \sum_{j=1}^{N} a_j \cos \psi_j, \quad \psi_j = \mathbf{k}_j \cdot \mathbf{r} - \omega_j t - \varphi_j \qquad (66)$$

$$\Phi^{(1)}(\mathbf{r}, t) = \sum_{j=1}^{N} b_j \cos \psi_j, \qquad (67)$$

and the following terms in the expansion involve nth-order multiplicative combinations of these linear spectral components. The perturbative expansions of elevation and velocity potential are identified simultaneously by injecting the successive Fourier expansions in the equations of motion. The leading, quadratic, nonlinear term for elevation was provided by Longuet-Higgins [1963] in the form (the factor 1/2 in the kernels $K_{ij}$ and $K'_{ij}$ is missing in the original paper by Longuet-Higgins [1963], as was later acknowledged by the author himself [see Srokosz and Longuet-Higgins, 1986]:

$$h^{(2)}(\mathbf{r}, t) = \frac{1}{2} \sum_{i,j=1}^{N} a_i a_j \left[ K_{ij} \cos \psi_i \cos \psi_j + K'_{ij} \sin \psi_i \sin \psi_j \right], \qquad (68)$$
where

$$
\begin{align}
K_{ij} &= (|\mathbf{k}_i||\mathbf{k}_j|)^{-\frac{1}{2}} [B_{ij}^{-} + B_{ij}^{+} - \mathbf{k}_i \cdot \mathbf{k}_j] + |\mathbf{k}_i| + |\mathbf{k}_j| \nonumber \\
K'_{ij} &= (|\mathbf{k}_i||\mathbf{k}_j|)^{-\frac{1}{2}} [B_{ij}^{-} - B_{ij}^{+} - |\mathbf{k}_i||\mathbf{k}_j|] \nonumber \\
B_{ij}^{\pm} &= \frac{\Omega_{ij}^{\pm} (\mathbf{k}_i \cdot \mathbf{k}_j \mp |\mathbf{k}_i||\mathbf{k}_j|)}{\Omega_{ij}^{\pm} - |\mathbf{k}_i \pm \mathbf{k}_j|} \tag{69} \\
\Omega_{ij}^{\pm} &= \left(\sqrt{|\mathbf{k}_i|} \pm \sqrt{|\mathbf{k}_j|}\right)^2 \nonumber
\end{align}
$$

To simplify the comparison we will again concentrate on the 2-D case. For long-crested waves we may operate the substitution $\mathbf{k}_i \cdot \mathbf{k}_j \rightarrow \text{sign}(k_i k_j) |k_i k_j|$, leading to simplified expressions of the kernels:

$$
K'_{ij} = -\operatorname{sign}(k_i k_j) \max(|k_i|, |k_j|) \quad (70)
$$

$$
K_{ij} = \min(|k_1|, |k_2|) \qquad (71)
$$

Now, in the perturbative expansion (62) after the CWM we
have:

$$
\begin{align}
h^{(2)}(x, t) &= -D(x, t) \partial_x h^{(1)}(x, t) \nonumber \\
&= -\sum_{i,j=1}^{N} a_i a_j k_j \operatorname{sign}(k_i) \sin \psi_i \sin \psi_j. \tag{72}
\end{align}
$$

Since the former process (68) is centered while the latter
(72) is not, we must rather compare with a recentered right-
hand side:

$$
\frac{1}{2} \sum_{i,j=1}^{N} a_i a_j (\tilde{K}_{ij} \cos(\psi_i) \cos(\psi_j) + \tilde{K}'_{ij} \sin(\psi_i) \sin(\psi_j)) \quad (73) \\
-\frac{1}{2} \sum_{i,j=1}^{N} a_i a_j \bar{K}_{ij} (\cos(\psi_i - \psi_j) - \delta_{ij}) \\
\text{with} \quad \left\{ \begin{array}{l} \tilde{K}'_{ij} = -\operatorname{sign}(k_i k_j) \max(|k_i|, |k_j|) \\ \bar{K}_{ij} = \operatorname{sign}(k_i k_j) \min(|k_i|, |k_j|), \end{array} \right.
$$

where $\delta_{ij}$ is the Kronecker symbol. This last expression is
resembling but not identical to the second-order correction
(68) of Longuet-Higgins. Note, however, that the respective
kernels coincide on the diagonal.

[45] This makes the CWM consistent with Longuet-
Higgins theory, at least for narrow spectra. Passing to the
limit of infinitely many spectral components, Longuet-
Higgins [1963] could also derive general formula for the
first few cumulants of the second-order nonlinear surface.
The mean and RMS of elevation at second-order are found
to be identical to those of linear process and the third
cumulant turns out to be nonvanishing (in the paper by
Longuet-Higgins [1963], the following expression is given
for one-sided spectrum only):

$$
\kappa_3 = 3 \int_{-\infty}^{+\infty} \int_{-\infty}^{+\infty} dk dk' \min(|k|, |k'|) \Gamma(k) \Gamma(k') \approx 3\sigma_0^2 \sigma_1^2 \quad (74)
$$

where as usual $\Gamma$ is the spectrum of the linear process $h^{(1)}$.
The corresponding skewness,

$$
\lambda_3 \approx \frac{3\sigma_1^2}{\sigma_0}, \qquad (75)
$$

has opposite sign with respect to the skewness (14) derived
in the framework of the CWM. However, the absolute
values of these quantities are too small for their sign to be
meaningful. A quick estimation can be performed with a
power law Phillips omnidirectional spectrum, $\Gamma(k) = 0.0025 \times$
$|k|^{-3}$, for $|k| > k_{peak}$, in which case we the skewness predicted
by the two models are found quasi-independent of the peak
wave number, $\lambda_3 \simeq 0.015$ for the Longuet-Higgins theory
and $\lambda_3 \simeq -3.10^{-6}$ for the CWM. Note that some recent
numerical experiments for one-dimensional surfaces after the
so-called ChSh method [Chalikov, 2005] show a unambigu-
ously positive skewness, so that the precision of the CWM
might no be sufficient to capture the latter correctly.

4.3. Weber and Barrick Theory

[46] In their 1977 companion papers Weber and Barrick
[1977] and Barrick and Weber [1977] revisited the nonlin-
ear theory for random seas with continuous spectra. The
adopted methodology is essentially the same as Longuet-
Higgins [1963] but the perturbative expansion is operated
on the continuous Fourier components of the surface. The
time-evolving surface elevation is sought in the form:

$$
\tilde{h}(\mathbf{r}, t) = \int_{\mathbb{R}^2} \mathbf{d}\mathbf{k} \int_{\mathbb{R}} d\omega \, \hat{h}(\mathbf{k}, \omega) e^{i(\mathbf{k} \cdot \mathbf{r} - \omega t)}, \quad (76)
$$

and a perturbative expansion is operated on the Fourier
components.

$$
\hat{h}(\mathbf{k}, \omega) = h_1(\mathbf{k}, \omega) + h_2(\mathbf{k}, \omega) \quad (77)
$$

The first-order term correspond to free waves propagating
with the gravity wave dispersion relation $\omega = \sqrt{g|\mathbf{k}|}$ (the
"linear" term),

$$
h_1(\mathbf{k}, \omega) = h_1^+(\mathbf{k})\delta\left(\omega - \sqrt{g|\mathbf{k}|}\right) + h_1^-(\mathbf{k})\delta\left(\omega + \sqrt{g|\mathbf{k}|}\right) \quad (78)
$$

while the second-order term is found to be:

$$
\begin{equation}
\begin{split}
h_2(\mathbf{k}, \omega) &= \int \mathbf{dk}_1 \mathbf{dk}_2 d\omega_1 d\omega_2 A(\mathbf{k}_1, \mathbf{k}_2) \delta(\mathbf{k}_1 + \mathbf{k}_2 - \mathbf{k}) \\
&\qquad \delta(\omega_1 + \omega_2 - \omega) h_1(\mathbf{k}_1, \omega_1) h_1(\mathbf{k}_2, \omega_2)
\end{split}
\tag{79}
\end{equation}
$$

Here the kernel $A$ is given by:

$$
\frac{1}{2} \left( |\mathbf{k}_1| + |\mathbf{k}_2| + \sqrt{|\mathbf{k}_1||\mathbf{k}_2|} \left(1 - \hat{\mathbf{k}}_1 \cdot \hat{\mathbf{k}}_2\right) \frac{|\mathbf{k}_1 + \mathbf{k}_2| + \Omega_{12}^+}{|\mathbf{k}_1 + \mathbf{k}_2| - \Omega_{12}^+} \right), \quad (80)
$$

where $\Omega_{12}^{+}$ is given by (70), and $A = 0$ whenever $\mathbf{k}_2 = -\mathbf{k}_1$
and $\omega_2 = -\omega_1$. Even through this is not obvious at first
sight, this kernel is consistent with Longuet-Higgins
perturbative theory since:

$$
A(\mathbf{k}_i, \mathbf{k}_j) = \frac{1}{2} (K_{ij} - K'_{ij}) \quad (81)
$$

To make a comparison with the CWM, we will consider the
surface frozen at a given time, say $t = 0$, in which case the
spatial process $h(\mathbf{r}) = h(\mathbf{r}, 0)$ at first- and second-order can be written:

$$ \tilde{h}(\mathbf{r}) = \int_{\mathbb{R}^2} \mathbf{dk} [\hat{h}_1(\mathbf{k}) + \hat{h}_2(\mathbf{k})] e^{i(\mathbf{k} \cdot \mathbf{r})}, \quad (82) $$

with

$$ \hat{h}_1(\mathbf{k}) = [\hat{h}_1^+(\mathbf{k}) + \hat{h}_1^-(\mathbf{k})] \quad (83) $$

and

$$ h_2(\mathbf{k}) = \sum_{s_1, s_2 = \pm 1} \int \mathbf{dk}_1 \mathbf{dk}_2 A(\mathbf{k}_1, \mathbf{k}_2) \delta(\mathbf{k}_1 + \mathbf{k}_2 - \mathbf{k}) h_1^{s_1}(\mathbf{k}_1) h_1^{s_2}(\mathbf{k}_2) \quad (84) $$

Denoting as usual $\Gamma$ and $\tilde{\Gamma}$ the first- and higher-order wave number spectra,

$$ \begin{aligned} \langle h_1(\mathbf{k}_1)h_1(\mathbf{k}_2) \rangle &= \delta(\mathbf{k}_1 - \mathbf{k}_2)\Gamma(\mathbf{k}_1), \\ \langle \tilde{h}(\mathbf{k}_1)\tilde{h}(\mathbf{k}_2) \rangle &= \delta(\mathbf{k}_1 - \mathbf{k}_2)\tilde{\Gamma}(\mathbf{k}_1), \end{aligned} \quad (85) $$

we can easily establish the following relationship:

$$ \tilde{\Gamma}(\mathbf{k}) = \Gamma(\mathbf{k}) + \int \mathbf{dk}_1 \Pi(\mathbf{k}_1, \mathbf{k} - \mathbf{k}_2) \Gamma(\mathbf{k}_1) \Gamma(\mathbf{k} - \mathbf{k}_1), \quad (86) $$

with

$$ \Pi(\mathbf{k}_1, \mathbf{k}_2) = 2|A(\mathbf{k}_1, \mathbf{k}_2)|^2 \quad (87) $$

In the 2-D case, the kernel reduces to:

$$ \Pi(k_1, k_2) = \frac{1}{2}(|k_1| + |k_2|)^2 \quad (88) $$

and thus

$$ \tilde{\Gamma}(k) = \Gamma(k) + \frac{1}{2} \int_{-\infty}^{+\infty} dk' k^2 \Gamma(k') \Gamma(k - k'), \quad (89) $$

which is the first integrand appearing in the low-frequency expansion after the CWM (29). As discussed later by Creamer et al. [1989], retaining this sole term leads to a divergence of the second-order correction at higher wave numbers. This is explained by the fact that the second-order spectrum (fourth order in surface amplitude) is not complete, since it misses the contribution of the $h_1 \times h_3$ term.

## 4.4. Creamer Theory

[47] In order to generate nonlinear sea surfaces Creamer et al. [1989] uses a canonical transformation of physical variables (surface elevation and potential) in order to improve the accuracy of the Hamiltonian expansion. This transformation has the same domain of validity of the CWM in a sense that it can be used for surface gravity waves and reproduces the effects of the lowest-order nonlinearities for the first-order development of the transformation. The 3-D

formulation remains, however, quite involved and its numerical implementation require further approximations [Soriano et al., 2006]. In the 2-D Creamer model, the nonlinear process $\tilde{h}$ is given by:

$$ \tilde{h}(x) = h(x) + \delta h(x), \quad (90) $$

where the corrective term $\delta h$ is expressed by its Fourier transform

$$ \hat{\delta h}(k) = \int_{-\infty}^{+\infty} dx \, e^{-ikx} \left( \frac{e^{ikD} - 1}{|k|} - i \operatorname{sign}(k) D \right) \quad (91) $$

and $D$ is the Hilbert transform of $h$. This expression is unpractical for further analytical investigation. However, at low frequencies ($kD \ll 1$) the exponential may be expanded,

$$ \hat{\delta h}(k) \approx \int_{-\infty}^{+\infty} dx \, e^{-ikx} \left( -\frac{|k|}{2} D^2 - \frac{i}{6} k^2 \operatorname{sign}(k) D^3 \right), \quad (92) $$

leading to the lowest-order approximation for the dressed spectrum [Creamer et al., 1989, equation 6.11]:

$$ \tilde{\Gamma}(k) = \Gamma(k) + \frac{1}{2} \int_{\mathbb{R}} dk' [k^2 \Gamma(k') \Gamma(k-k') - 2k^2 \Gamma(k) \Gamma(k')] \quad (93) $$

This expression is similar to the first integral in the low-frequency expansion (29). Figure 11 displays a comparison of Creamer et al.'s [1989], Weber and Barrick's [1977] and CWM low-frequency expansion, for an omnidirectional $k^{-3}$ spectrum with exponential cutoff at peak frequency $kp = 0.7 \text{ rad m}^{-1}$ (corresponding to a wind of 3 m s$^{-1}$) and upper limit $k_u = 120 \text{ rad m}^{-1}$. The undressed (linear) spectrum is shown together with the corrections brought by the dressed spectrum. Creamer et al.'s [1989] and CWM expansions are extremely close at low frequency but CWM eventually diverges at higher frequency ($k > 100 \text{ rad m}^{-1}$). Weber and Barrick [1977] diverge very early ($k > 3 \text{ rad m}^{-1}$) and is slightly higher than CWM and Creamer et al.'s [1989] corrections.

# 5. Comparison With Experimental Data

[48] The reference data basis for the sea wave slope distribution is the optically derived measurement of Cox and Munk [1954], which has been used to calibrate many models of the literature. Since the CWM in its current state is restricted to gravity waves only, it cannot describe the scales smaller than, say 5 cm, and the related slopes, making the comparison with Cox and Munk data irrelevant. Instead, we will resort to a recent airborne campaign [Vandemark et al., 2004], which has provided laser measurements of the omnidirectional slope statistics of long gravity waves. This amounts to filter out in the slope statistics the contribution of wavelengths smaller than about 2 m and renders the comparison with the CWM possible. The main outcome of this study was an elevated kurtosis for the omnidirectional slope, a result that can be put on the account of either the strong directionality of the wavefield or its non-Gaussian character. We will investigate the respective contributions of
<figure><img src="image_12.png" /><figcaption>Figure 11. Comparison of Creamer et al.'s [1989], Weber and Barrick's [1977], and CWM low-frequency expansions for the corrective term to the undressed spectrum.</figcaption></figure>

these two effects in the framework of the CWM. Denote $P_{2-omni}(S)$ the omnidirectional slope distribution, that is the distribution of absolute magnitude of slope $S = |\nabla h|$. For an isotropic Gaussian distribution with variance $\sigma_2^2$, this is a Rayleigh distribution with parameter $\sigma_2$,

$$P_{2-omni}(S) = \frac{S}{\sigma_2^2} \exp\left(-\frac{S^2}{2\sigma_2^2}\right), \quad (94)$$

whose kurtosis is $\lambda'_4 = 3.245$. For a directional Gaussian slope distribution with upwind/crosswind mean square slope (mss) ratio $\rho^2 = \sigma_{200}^2/\sigma_{020}^2$ and total mss $\sigma_2^2$, the nth moments $M_n = \langle S^n P_{omni}(S) \rangle$ of the omnidirectional slope distributions are found to be:

$$M_n = \frac{\sigma_2^n}{2\pi\rho} \int_0^{2\pi} d\theta \left[ \cos^2\theta + \frac{\sin^2\theta}{\rho^2} \right]^{-(1+n/2)} R_n \quad (95)$$

where $R_n$ is the nth moment of the normalized Rayleigh distribution ($S \exp(-S^2/2)$). The variation of kurtosis with the directionality parameter $\rho$ can be estimated numerically. A maximum value $\lambda'_4 = 4.166$ is reached at $\rho = 3$, while the minimum kurtosis is obtained at $\rho = 1$ for the Rayleigh distribution ($\lambda_4 = 3.245$). It follows that the elevated values of kurtosis reported by Vandemark et al. [2004], which ranges from 4.5 to 6, cannot be explained by mere directional effects of the slope distribution.

[49] The kurtosis has been computed as a function of wind speed for both linear and CWM surfaces generated with a directional Elfouhaily (undressed) spectrum. The fourth moment of the theoretical CWM slope distribution is in principle infinite, but the corresponding integral can be shown to have a slow, logarithmic divergence. Therefore, the slope distribution has been truncated to a maximum value of 1.7, corresponding to a steep wave of about 60 degree. For

small and moderate winds ($U_{10} \le 12 \text{ m s}^{-1}$), the resulting fourth moment is quite insensitive to the chosen threshold. Furthermore, we have checked that the lack of normalization of the slope distribution after truncation has a negligible impact on the computation of the first cumulants. At higher winds, the slope kurtosis is found to increase slightly with the slope threshold. However, we do not expect the CWM to remain meaningful for steep waves. The simulated excess kurtosis is shown on Figure 12 and compared with recorded data. To reproduce the filtering of small waves slopes realized in the paper by Vandemark et al. [2004], the Elfouhaily spectrum has been truncated to a maximum wave number of $k_u = 6 \text{ rad m}^{-1}$. The corresponding surfaces are referred to as “long gravity waves.” The comparison is given with the untruncated gravity waves Elfouhaily spectrum ($k_u = 200 \text{ rad m}^{-1}$). The horizontal line at $\gamma = 0.245$ is the excess kurtosis of the Rayleigh distribution, obtained for Gaussian isotropic slope distribution. The line at $\gamma = 0.7$ is the estimation of Cox and Munk [1954], which is insensitive to wind and identical for slick and clean surfaces. The inclusion of nonlinearities through the CWM drastically increases the excess kurtosis and brings it to values intermediate between Vandemark et al.'s [2004] data and Cox and Munk's [1954] data, while the linear model remains closer to the Rayleigh distribution.

## 6. Conclusion

[50] As reported, CWM provides an analytically tractable, numerically efficient solution to approach the geometrical description of nonlinear surface waves. CWM is also robust to the inclusion of high frequencies. CWM explicitly builds on a phase perturbation method to modify the surface coordinates, and statistical properties can be derived. We establish the complete first- and second-order statistical properties of surface elevations and slopes for long-crested
<figure><img src="image_13.png" /><figcaption>**Figure 12.** Excess kurtosis of the omnidirectional slopes as a function of wind speed.</figcaption></figure>

as well as fully two-dimensional surfaces. As compared to standard approximation, the CWM is shown to be a reasonably accurate model for weak nonlinear gravity-wave interactions. It is based on the local deformation of a reference Gaussian process and the first few cumulants up to fourth order can be expressed in terms of the underlying Gaussian statistics. Relations between dressed and undressed spectra have been established and found to favorably extend the classical low-frequency formulations of Weber and Barrick [1977] and Creamer et al. [1989].

[51] As already pointed out [Elfouhaily et al., 1999], it can be crucial to determined the required input undressed spectrum for which the simulated moments remain consistent with a measured spectrum. The CWM can then be used to define an inversion scheme to consistently evaluate the first-order cumulants (elevation skewness, the elevation and slope cross skewness) to evaluate the predicted long wave geometrical contribution to altimeter sea state bias [Elfouhaily et al., 2000; Vandemark et al., 2005].

[52] Moreover, the nonlinear surface wave geometry with shallow troughs and enhanced crests, implies an excess of both zero and steep slope occurrences. As numerically derived, CWM predictions unambiguously confirm that bound harmonics associated to the simplified surface coordinate changes will indeed lead to non negligible surface slope kurtosis. Compared to measurements, CWM is found to help to bridge the differences between a linear Gaussian model and reported large slope kurtosis.

[53] For short gravity waves, the CWM can also be used to heuristically introduce the skewness of individual slopes. These effects can indeed be subsequently incorporated in the model through a generalization of the horizontal displacement $D(r, t)$ to steepen slightly the forward face of individual

waves, especially when the local steepness exceeds a threshold value.

[54] Acknowledgments. Frédéric Nouguier is funded by the Délégation Générale pour l'Armement (DGA) under the supervision of Y. Hurtaud. Many thanks go to Philippe Forget for his careful reading of the paper.

## References

Aberg, S. (2007), Wave intensities and slopes in Lagrangian seas, *Adv. Appl. Prob.*, **39**(4), 1020–1035.

Barrick, D. E., and B. L. Weber (1977), On the nonlinear theory for gravity waves on the ocean's surface. Part II: Interpretation and applications, *J. Phys. Oceanogr.*, **7**, 11–21.

Chalikov, D. (2005), Statistical properties of nonlinear one-dimensional wave fields, *Nonlinear Process. Geophys.*, **12**(5), 671–689.

Chalikov, D., and D. Sheinin (2005), Modeling extreme waves based on equations of potential flow with a free surface, *J. Comput. Phys.*, **210**, 247–273.

Chapron, B., V. Kerbaol, D. Vandemark, and T. Elfouhaily (2000), Importance of peakedness in sea surface slope measurements and applications, *J. Geophys. Res.*, **105**, 17,195–17,202.

Cox, C., and W. Munk (1954), Statistics from the sea surface derived from the sun glitter, *J. Mar. Res.*, **13**, 198–227.

Creamer, D., F. Henyey, R. Schult, and J. Wright (1989), Improved linear representation of ocean surface waves, *J. Fluid Mech.*, **205**, 135–161.

Dias, F., and T. J. Bridges (2006), The numerical computation of freely propagating time-dependent irrotational water waves, *Fluid Dyn. Res.*, **38**, 803–830.

Elfouhaily, T. (2000), Truncated Hamiltonian versus surface perturbation in nonlinear waves theories, *Waves Random Media*, **10**, 103–116.

Elfouhaily, T., B. Chapron, K. Katsaros, and D. Vandemark (1997), A unified directional spectrum for long and short wind-driven waves, *J. Geophys. Res.*, **102**, 15,781–15,796.

Elfouhaily, T., D. Thompson, D. Vandemark, and B. Chapron (1999), Weakly nonlinear theory and sea state bias estimation, *J. Geophys. Res.*, **104**, 7641–7647.

Elfouhaily, T., D. Thompson, B. Chapron, and D. Vandemark (2000), Improved electromagnetic bias theory, *J. Geophys. Res.*, **105**, 1299–1310.

Fournier, A., and W. T. Reeves (1986), A simple model of ocean waves, *Comput. Graphics*, **20**(4), 75–84.
Fructus, D., D. Clamond, J. Grue, and Ø. Kristiansen (2005), An efficient model for three-dimensional surface wave simulations: Part I: Free space problems, *J. Comput. Phys.*, **205**, 665–685.

Gerstner, F. (1809), Theorie der wellen samt einer daraus abgeleiteten
theorie der deichprofile, *Ann. Phys.*, **2**, 412–445.

Gjosund, S. (2003), A Lagrangian model for irregular waves and wave
kinematics, *J. Offshore Mech. Arctic Eng.*, **125**, 94–102.

Hasselmann, K. (1962), On the nonlinear energy transfer in a gravity-wave
spectrum. Part 1. General theory, *J. Fluid Mech.*, **12**, 481–500.

Joelson, M., and M. Néel (2008), On alpha stable distribution of wind
driven water surface wave slope, *Chaos*, **18**, 033117.

Longuet-Higgins, M. (1963), The effects of non-linearities on statistical distributions in the theory of sea waves, *J. Fluid Mech.*, **17**, 459–480.

Papoulis, A. (1965), *Probability, Random Variables and Stochastic Processes*, pp. 153–162, McGraw-Hill, Boston, Mass.

Pierson, W. (1961), Models of random seas based on the Lagrangian equa-
tions of motion, *Tech. Rep. Nonr-285(03)*, Dep. of Meteorol. and Ocea-
nogr., Coll. of Eng. Res. Div., N. Y. Univ., New York.

Pierson, W. (1962), Perturbation analysis of the Navier-Stokes equations in Lagrangian form with selected linear solutions, *J. Geophys. Res.*, **67**, 3151–3160.

Ruban, V. (2005), Water waves over a time-dependent bottom: Exact de-
scription for 2D potential flows, *Phys. Lett. A*, **340**(1–4), 194–200.

Soriano, G., M. Joelson, and M. Saillard (2006), Doppler spectra from a two-dimensional ocean surface at l-band, *IEEE Trans. Geosci. Remote Sens.*, **44**(9), 2430–2437.

Srokosz, M., and M. Longuet-Higgins (1986), On the skewness of sea-
surface elevation, *J. Fluid Mech.*, **164**, 487–497.

Stokes, G. (1847), On the theory of oscillatory waves, *Trans. Cambridge
Philos. Soc.*, **8**(441), 197–229.

Stokes, G. (1880), Supplement to a paper on the theory of oscillatory
waves, *Math. Phys. Pap.*, **1**, 225–228.

Tayfun, M. (1980), Narrow-band nonlinear sea waves, *J. Geophys. Res.*, **85**, 1548–1552.

Tayfun, M. A., and F. Fedele (2007), Wave-height distributions and non-linear effects, *Ocean Eng.*, **34**, 1631–1649.

Tick, L. (1959), A non-linear random model of gravity waves, *J. Math. Mech.*, **8**(5), 643–651.

Toffoli, A., E. Bitner-Gregersen, M. Onorato, and A. Babanin (2008), Wave
crest and trough distributions in a broad-banded directional wave field,
*Ocean Eng.*, **35**, 1784–1792.

Vandemark, D., B. Chapron, J. Sun, G. Crescenti, and H. Graber (2004), Ocean wave slope observations using radar backscatter and laser altimeters, *J. Phys. Oceanogr.*, **34**, 2825–2842.

Vandemark, D., B. Chapron, T. Elfouhaily, and J. Campbell (2005), Impact of high-frequency waves on the ocean altimeter range bias, *J. Geophys. Res.*, **110**, C11006, doi:10.1029/2005JC002979.

Watson, K., and B. West (1975), A transport-equation description of non-linear ocean surface wave interaction, *J. Fluid Mech.*, **10**, 815–826.

Weber, B., and D. Barrick (1977), On the nonlinear theory for gravity
waves on the ocean's surface. Part I: Derivations, *J. Phys. Oceanogr.*,
*7*, 3–10.

West, B. J., K. A. Brueckner, R. S. Janda, D. M. Milder, and R. L. Milton (1987), A new numerical method for surface hydrodynamics, *J. Geophys. Res.*, **92**, 11,803–11,824.

Zakharov, V. (1968), Stability of periodic waves of finite amplitude on the surface of a deep fluid, *J. Appl. Mech. Tech. Phys. Engl. Transl.*, **2**, 190–194.

Zakharov, V., A. Dyachenko, and O. Vasilyev (2002), New method for numerical simulation of a nonstationary potential flow of incompressible fluid with a free surface, *Eur. J. Mech. B Fluids*, **21**(3), 283–291.

<footnote>B. Chapron, Laboratoire d'Océanographie Spatiale, IFREMER, Zone 1
Industrielle Pointe du Diable, B.P. 70, F-29280 Plouzané, France.</footnote>

<footnote>C.-A. Guérin, LSEET, UMR 6017, Université du Sud-Toulon-Var,
CNRS, B.P. 20132, F-83957 La Garde CEDEX, France.</footnote>

<footnote>F. Nouguier, Institut Fresnel, UMR 6133, Faculté de Saint-Jérôme,
Université Paul Cézanne, CNRS, F-13397 Marseille CEDEX 20, France.
(frederic.nouguier@fresnel.fr)</footnote>