## 1. Introduction

The characterization of surface gravity waves is important for air–sea interaction processes, such as the exchange of energy, momentum, mass, and heat between the ocean and the atmosphere, as well as for risk assessment for offshore structures and related engineering applications. This study is concerned with the characterization of spatial statistics of the sea surface elevation using airborne lidar observations. Several field investigations have reported on the statistics of various parameters describing surface wind waves such as the surface displacement, wave crest and trough displacement, wave heights, and joint statistics of wave periods and amplitudes. These measured statistics have been shown to compare favorably with theoretical distributions of linear and nonlinear waves. However, the vast majority of wave observations in the field come from temporal series at a point or small spatial wave-gauge arrays. In this study we present an analysis of airborne spatiotemporal measurements of the sea surface elevation that provide directional and slope information that cannot be obtained directly using single point measurements. The measurements were collected from the NSF/NCAR C-130 aircraft during the Gulf of Tehuantepec Experiment (GOTEX) in fetch-limited conditions during the month of February 2004 (Romero and Melville 2010a).

A detailed analysis of several statistical parameters, including the probability distribution function (pdf) of wave heights, wave slopes, the joint pdf of wave amplitudes and wavelengths, as well as the length of contours bounding threshold elevations is presented. The measured statistical distributions are used to test analytical linear and nonlinear models, some of which have not been tested previously against field measurements. Another important parameter, which is not directly available from single point measurements, is the length of crests exceeding threshold wave criteria, such as wave height or slope. In this study, the measured crest length distributions are compared against second-order nonlinear stochastic simulations using computed directional wavenumber spectra (Romero and Melville 2010b). Preliminary results of this work have been reported by Melville et al. (2005) and Romero and Melville (2010c).

This paper is organized as follows. In section 2, we provide background information, including relevant definitions and formulas. In section 3, the data and analysis are described. The results are presented and summarized in sections 4 and 5, respectively.

## 2. Background

### a. Surface elevation

*η*can be represented in space and time (

**x**,

*t*) using a Fourier series:where

*a*is the

_{n}*n*th Fourier amplitude,

**k**

*is the wavenumber vector,*

_{n}*ω*is the frequency, and

_{n}*ϵ*is the phase. This study is concerned with the spatial description of

_{n}*η*(

**x**) at a given instant in time (e.g.,

*t*=

*t*), with a corresponding wavenumber spectrum

_{o}*F*(

**k**) ≡

*F*(

*k*

_{1},

*k*

_{2}) that, according to Parseval’s theorem, is related to the variance of the sea surface displacement 〈

*η*

^{2}〉 bywith the brackets representing an average over space and

*m*

_{00}corresponding to the zeroth spectral moment. The first and second moments are defined according towhere the first moments

*m*

_{10}and

*m*

_{01}normalized by the variance provide an estimate of the spectrally weighted mean wavenumber vector defined asand the second moments

*m*

_{20}and

*m*

_{02}are by definition related to the root-mean-square slope (RMSS) of the sea surface elevation (Cox and Munk 1954):Additional spectral moments about the mean wavenumber vector are related to the statistics of the wave envelope (Longuet-Higgins 1957) and are defined as follows:andThe above spectral moments are introduced for convenience and are referenced throughout the text. Below we provide a brief background on the description of the sea surface elevation.

*η*is Gaussian, which in nondimensional form is given bywhere

*η*′ =

*η*/

*η*

_{rms}with

*η*

_{rms}corresponding to the root-mean-square surface elevation. As is well known, the Gaussian distribution has both zero skewness 〈

*η*

^{3}〉/〈

*η*

^{2}〉

^{3/2}and zero excess kurtosis 〈

*η*

^{4}〉/〈

*η*

^{2}〉

^{2}− 3. Longuet-Higgins (1963) included the effect of nonlinearities on the statistical description of the sea surface elevation in the form of a Gram–Charlier distribution and showed that the skewness of the sea surface elevation is always positive. The Gram–Charlier distribution, which allows for modifications of the Gaussian distribution with parametric dependence on the skewness and kurtosis, has the obvious drawback of occasionally giving negative probabilities. Subsequently, Tayfun (1980) derived a second-order distribution for a narrowbanded energy spectrum of unidirectional waves. The Tayfun distribution was rederived by Socquet-Juglard et al. (2005) in a simpler form, given bywherewith

Following the theoretical framework by Longuet-Higgins (1963), Forristall (2000) compared measured wave crest displacement distributions against second-order stochastic simulations using measured frequency spectra. His results show that second-order simulations give good agreement with observed wave crest distributions in the field, including hurricane wind conditions.

### b. Wave height

*a*for a narrow energy spectrum, which is given bywhere

*H*, assuming

*H*= 2

*a*, Forristall (1978) defines Eq. (14) with parametric dependence on the total variance of the energy spectrum

*m*

_{0}as given byEquations (14) and (15) are identical when

*α*= 2.126 and

*β*= 8.42. The reduction of observed wave heights, when compared to the Rayleigh distribution, was further examined by Longuet-Higgins (1980), who showed that it could be attributed to finite-bandwidth effects. Subsequently, several authors derived theoretical distributions that are consistent with the empirical Weibull distribution (16) (Tayfun 1981; Naess 1985; Tayfun 1990, hereafter T90; Boccotti 2000, hereafter B00; among others). For a recent review, refer to Tayfun and Fedele (2007). Here we consider two asymptotic models of wave height distributions for large waves, see T90 and B00, both of which depend on the autocorrelation function of the sea surface

**Ψ**(

*t*) and are given byandrespectively. The normalized temporal covariance or autocorrelation function is related to the frequency spectrum

*φ*(

*ω*) byIn the above development,

**Ψ*** ≡ |

**Ψ**(

*t**)|, with

*t** corresponding to the first minimum from the center at

*t*= 0, and

*r*(

*t*) evaluated at half of the spectrally weighted mean wave period

*T*, and

_{m}**Ψ**, defined by

It has been shown that, for a narrowband spectrum, second-order nonlinearities produce negligible effects on the statistical distribution of heights of large waves (Tayfun 1983). This is because both the troughs and the crests are raised by comparable amounts in a narrowband system. However, recent higher-order numerical simulations by Socquet-Juglard et al. (2005) have shown that, in a narrowband spectrum with narrow directional spreading, corresponding to long-crested waves, the wave height probabilities can significantly deviate from the Raleigh distribution, giving rise to extreme wave heights due to modulational instabilities. This has also been shown in laboratory experiments of long-crested (one dimensional) waves by Mori et al. (2007). However, for field conditions, the wind sea spectrum has finite directional spreading that inhibits the modulational instability typical of long-crested waves (Socquet-Juglard et al. 2005). A similar conclusion was drawn in a recent theoretical study by Tayfun and Fedele (2007).

#### Local wave parameters

*p*(

*a*′,

*λ*′), which was shown to compare well with laboratory measurements of one-dimensional waves. However, to our knowledge their distribution has not been tested against spatial data in the field. The distribution by Xu et al. (2004) is given bywhere the normalized wave amplitude

*λ*′ = 1/

*k*′ with

*ν*is the bandwidth of the wavenumber spectrum given bywherecorresponds to the spectral moments of the one-dimensional wavenumber spectrum, defined as the omnidirectional spectrum,with

*θ*being the azimuthal angle between −

*π*and

*π*. In this study, the directional wavenumber spectrum is defined between

^{1}that arises from the two-dimensional Fourier transform of “stationary” spatial images of the sea surface elevation. It is well known that most of the wind wave energy at scales near the spectral peak is mostly contained within ±90° from the peak (see, e.g., Fig. 4 of Hwang and Wang 2001). Thus it is expected that the inherent ambiguity of the measured spectra does not significantly affect the results of this study.

## 3. Data

We present an analysis of airborne lidar observations collected from the NSF/NCAR C-130 aircraft during GOTEX in February 2004. The measurements were collected in offshore wind conditions with wind speeds between 10 and 20 m s^{−1} and fetches between 20 and 500 km. By definition, fetch-limited wind wave conditions require limited fetch, homogeneous and steady winds blowing off an infinite and straight coastline and, most importantly, a steady wave field. Such idealized conditions are only approximately found in nature. In particular, the Gulf of Tehuantepec is well known for having strong offshore winds that last from 3 to 5 days, which could in principle generate steady wave fields. However, since the winds are forced by pressure differences across a narrow mountain gap, the resulting wind field has an inhomogeneous jetlike structure, with maximum speeds and a mean direction due south close to shore and gradually decreasing and veering due west as the wind jet fans out over the Pacific Ocean (Steenburgh et al. 1998). However, the wave data collected during GOTEX have been shown to agree with the classical fetch relations (Romero and Melville 2010a). This implies that the wave measurements were collected in approximately fetch-limited conditions. The evolution of the directional wavenumber spectra was characterized in detail and compared against numerical simulations using a state-of-the-art wind wave model by Romero and Melville (2010a,b). Below, we provide a brief summary of the data.

### a. Sea surface topography

The Airborne Topographic Mapper (ATM) is a conical scanning lidar, developed by the National Aeronautics and Space Administration and EG&G, which has been shown to successfully measure the directional properties of ocean surface waves (Hwang et al. 2000a,b; Romero and Melville 2010a). During GOTEX, the ATM conical scanning angle was set to 15° with a pulse repetition and sampling frequency *f _{s}* = 5 kHz and a scanning frequency

*f*

_{sc}= 20 Hz. Because the ATM was operated primarily at a nominal altitude of 400 m above the mean sea level, the radius of the scanning pattern on the sea surface,

*R*, is approximately 100 m. Assuming no pulse return dropouts, the maximum horizontal separation between consecutive measurements in the along-flight direction is given by the ratio of the horizontal aircraft speed

*V*to the scanning frequency

_{a}*f*

_{sc}. For the typical speed of

*V*= 100 m s

_{a}^{−1}, the horizontal resolution in the along-flight direction is 5 m. The cross-track resolution is approximately 2.5 m, given by the ratio of the perimeter of the circular scan at the surface (

*P*= 2π

*R*≈ 628 m) and the number of pulses along the scan,

*N*=

*f*/

_{s}*f*

_{sc}= 250. According to Krabill and Martin (1987), the calibrated absolute error per pulse in the elevation measured by the ATM is 8 cm, which includes a 3 cm (rms) range error, 5 cm for positioning through differential GPS, and 5 cm for attitude-induced errors.

The ATM data along the circular scanning pattern on the surface of the ocean were separated into forward and rear scans, providing two partially overlapping swaths of the sea surface elevation separated in time by about 1 s. Each spatiotemporal data subset was binned and interpolated to a regular horizontal grid with a resolution of 5 m. All of the data analyzed in this study correspond to data swaths approximately aligned with the dominant waves, including both upwind and downwind legs.

### b. Spectral analysis

After the data were gridded and interpolated on a regular grid, the directional wavenumber spectrum was computed from the two-dimensional FFT over typical data segments 200 m wide by 4000 m long. Before Fourier transforming, each data segment was detrended and padded with zeros so that the spectral resolutions were approximately *dk _{y}* = 2π/320 = 0.0196 rad m

^{−1}and

*dk*= 2π/4000 = 0.0016 rad m

_{x}^{−1}in the cross-track and along-track directions, respectively.

*h*are the wavenumber and water depth, respectively. Thus, in deep water the relative shift

*δk*/

_{x}*k*=

*c*/

*V*is greater for the longer waves. The “true” wavenumber component in the along-flight direction

_{a}*k*was calculated iteratively from Eq. (26) as follows: 1)

_{x}*ω*is obtained from

*k*and then used to calculate

_{y}*δk*; 2) an initial value of

_{x}*k*is obtained according towhere the ± sign depends on the propagation direction of the waves relative to the aircraft, being positive for waves traveling in the flight direction and vice versa; and 3) steps 1 and 2 are repeated twice using the newly calculated values of

_{x}*k*to calculate

_{x}*ω*and

*δk*. This Doppler correction was applied to the measured directional spectra and the local wavenumber vectors calculated in the spatial domain using the Hilbert transform as described in appendix A.

_{x}### c. Data summary

In this statistical study all of the results are presented in nondimensional form. Here a brief summary of the data without normalization is provided. Figure 1 shows the 10-m wind speed *U*_{10}, significant wave height *H _{s}*, dominant wavelength

*λ*, and the alignment factor between the dominant wave direction

_{p}*θ*and the wind direction

_{p}*θ*, defined as cos(

_{w}*θ*−

_{w}*θ*), versus the average offshore distance along the dominant wave direction, or effective fetch

_{p}**X**

*(Romero and Melville 2010a). The data from research flights (RFs) 5, 7, 9, and 10 are shown with circles, diamonds, stars, and squares, respectively. The 10-m wind speed was corrected from the measurements collected between 30 and 50 m above the mean sea level, including corrections due to atmospheric stability as described by Romero and Melville (2010a). For all four flights the wind speed shows a decreasing trend with increasing fetch, decaying on average from about 18 to 15 m s*

_{e}^{−1}within roughly 400 km. The significant wave height on average increases from 1.5 to 3.5 m, with most of the growth occurring within the first 100 km off the coast. The dominant wavelength increases from about 50 to 125 m. As discussed in Romero and Melville, the observed winds and the dominant waves at short fetches were typically oriented due south while rotating toward the west with increasing fetch. As shown in Fig. 1d, both wind and waves are mostly aligned, with the exception of a few data points from RF 5 that were collected near the western edge of the wind jet.

### d. Data averaging

*g*is gravity,

**X**

*the effective fetch, and*

_{e}*χ*such that the total number of dominant waves per record

*N*was approximately the same in all cases analyzed. Following Socquet-Juglard et al. (2005),

_{w}*N*can be estimated bywhere

_{w}*L*and

_{x}*L*are the record lengths in the two orthogonal directions,

_{y}*x*corresponding to the along-flight direction and

*m*

_{20},

*m*

_{02}, and

*m*

_{00}are the second and zeroth moments of the directional wavenumber spectrum, respectively. Equation (29) was calculated from the measured wavenumber spectra with

*k*<

_{s}*k*

_{1}< 0.35 rad m

^{−1}and −0.35 rad m

^{−1}<

*k*

_{2}< 0.35 rad m

^{−1}corresponding to the resolved range of wavenumbers for the windsea part of the spectrum. The wavenumber component

*k*

_{1}is aligned with respect to the peak of the directional spectrum and

*k*is a low wavenumber cutoff separating the windsea from swell, as described by Romero and Melville (2010a).

_{s}## 4. Results

### a. Surface elevation and wave heights

The ATM measurements were bin averaged according to the nondimensional fetch *χ* and the total number of waves in a record as described in section 3, yielding several datasets, each covering an average area of 200 m by 70 km with approximately 2600 dominant waves, and four robust datasets with an average area of 200 m by 330 km and approximately 12 000 dominant waves at intermediate and large nondimensional fetches. Figures 2a,b show the steepness of the dominant waves 〈*η*^{2}〉^{1/2}*k _{p}*, also referred to as the significant slope (Huang et al. 1981), and the RMSS

*k*corresponding to the spectral peak wavenumber, the steepness of the dominant waves decreases with increasing

_{p}*χ*, whereas the RMSS remains nearly constant. This is because the RMSS depends on the second moment of the wavenumber spectrum, being sensitive to the high wavenumber part of the spectrum, whose contribution to RMSS increases with increasing fetch as the peak of the spectrum shifts to lower wavenumbers. The significant slope is a measure of the nonlinearity, which is related to the skewness 〈

*η*

^{3}〉/〈

*η*

^{2}〉

^{3/2}(Fig. 2c) showing a decreasing trend with respect to

*χ*, in good agreement with the second-order narrow band distribution by Tayfun (1980). The skewness is a measure of the vertical asymmetry for which positive values describe higher crests and shallower troughs. This effect is evident in the measured pdfs of the sea surface elevation, plotted in Figs. 3a–d, sorted with increasing nondimensional fetch. The data are well described by the narrowband second-order distribution (Tayfun 1980). However, the agreement degrades for very low probabilities, particularly in Fig. 3c. Another feature shown in Fig. 3c is a slight bias in the empirical distributions with negative displacements consistently exceeding Tayfun’s distribution.

*η*

^{4}〉/〈

*η*

^{2}〉

^{2}− 3, which is a measure of non-Gaussianity, as well as the incidence of extreme wave heights (Mori and Janssen 2006). Figure 2d shows the excess kurtosis against the nondimensional fetch. The kurtosis shows substantial variability, decreasing with increasing nondimensional fetch at a rate comparable to that predicted by the theoretical narrowband second-order statistics (Tayfun 1980). However, the measured kurtosis is always larger than the narrowband approximation. In an attempt to incorporate higher order nonlinearities into the prediction of the extreme wave heights, Mori and Janssen (2006) derived an expression that relates the kurtosis of the sea surface displacement to the Benjamin–Feir index (BFI) based on third-order nonlinear theory of unidirectional waves. The BFI is a measure of the slope normalized by the bandwidth, which Janssen (2005) defines aswhere the bandwidth parameter Δ is given byandFigures 2e,f show the measured BFI parameter and the normalized maximum wave height observed (

*H*

_{max}/

*H*) plotted against the nondimensional fetch. Both the BFI parameter and excess kurtosis decrease with increasing nondimensional fetch. However, contrary to the predictions by Mori and Janssen (2006), the BFI and maximum wave height observed are not correlated. We find waves that meet the typical “rogue wave” criterion (

_{s}*H*> 2

*H*) (Dysthe et al. 2008) only at longer fetches for which BFI is small.

_{s}*H*is defined as the zero-crossing distance from trough to crest, or crest to trough. In this study, the wave height is calculated from the two-dimensional data using the zero-up crossing definition along parallel lines in the along-track direction, which was approximately aligned with the dominant waves. Figure 4 shows the measured distribution of wave heights normalized by the significant wave height

*H*= 4〈

_{s}*η*

^{2}〉

^{1/2}, ordered with increasing nondimensional fetch in Figs. 4a–d, with

*H*estimated from the entire spectrum. With the exception of the shortest fetch, the data are in good agreement with the distributions by T90, B00, and the empirical distribution by Forristall (1978), with the Rayleigh distribution serving as an upper bound on the distribution of large wave heights. Of the three non-Rayleigh models considered, B00 gives the lowest probabilities of finding large waves. Here, the analytical models by T90 and B00 [Eqs. (17) and (18)] have been evaluated using the spatial two-dimensional autocorrelation function

_{s}*r*with

_{m}*x** corresponds to the first minimum of

*θ*,

_{p}### b. Local wave properties

Several studies in the literature have reported on empirical estimates of the joint probability distribution function of wave amplitudes (or heights) and periods from temporal field measurement collected at a point, for example, Longuet-Higgins (1975); Cavanié et al. (1976); Shum and Melville (1984), among others. A more recent study by Xu et al. (2004) derived the joint distribution of wave amplitudes and wavelengths, which compared favorably with laboratory observations. To our knowledge, field measurements of the joint distribution of wave amplitudes and wavelengths have not been made before. Here we present the statistics of local wave properties over a horizontal area, including wave elevation, wavelength, and wave slope calculated from the bandpassed data using the Hilbert transform (for details, see appendix A). The measurements are found to agree well with analytical linear models available in the literature.

Figure 5 shows the measured joint pdfs of dimensionless wavelengths *λ*′ = 1/*k*′ and wave amplitudes *a*′, where *a*′ ≈ 0.5, *λ*′ ≈ 0.5) waves.

The wave slope characterizes the nonlinearity of the wave field. Here, the wave slope vector *s*_{1}, *s*_{2}) = (*ak*_{1}, *ak*_{2}) calculated from the Hilbert transform of the bandpassed ATM data (for details, see appendix A). The directional pdfs of measured wave slopes *p*(*s*_{1}, *s*_{2}) are shown in Fig. 6. The empirical distributions show a narrowing trend with increasing nondimensional fetch, comparing favorably with the analytical expression (A14) derived in appendix A by following Longuet-Higgins (1957). The agreement between the empirical and theoretical distributions is best at short fetches (Figs. 6a,b). At long fetches the observed distributions are narrower than the theory with reduced probabilities of finding steep wave slope components orthogonal to the dominant waves. A similar bias was found between the empirical and theoretical distributions of wavenumber vectors shown in appendix A.

*s*from bandpassed ATM data are shown in Fig. 7. The data show a significant reduction in the probability of finding steep dominant waves with increasing nondimensional fetch, in good agreement with the theory. The empirical distributions are compared to three other analytical expressions, that depend on different moments of the spectrum: 1)

*p*

_{LH}(

*s*) corresponding to the directional slope distribution

*p*(

*s*

_{1},

*s*

_{2}), derived in Eq. (A14), following the work of Longuet-Higgins, and integrated numerically in azimuth between

*p*(

_{X}*s*) estimated numerically from the distribution of wavelengths and amplitudes (Xu et al. 2004) in Eq. (22); and 3)

*p*

_{R}(

*s*) corresponding to the Rayleigh distribution for wave amplitudes (Cartwright and Longuet-Higgins 1956; Longuet-Higgins 1957) scaled with the mean wavenumber as given byAll three distributions, plotted in Fig. 7, show very similar trends, with

*p*

_{LH}giving the best agreement with the measured distributions and

*p*

_{R}doing a better job than

*p*. Both

_{X}*p*

_{R}and

*p*serve as upper bounds on the distribution of steep wave slopes.

_{X }### c. Lengths of contours and crests

The two-dimensional spatiotemporal data permit direct estimates of spatial information such as the lengths of crests exceeding threshold wave heights and slopes. Here the statistics of the length of contours surrounding large crests and the distributions of crest length per unit area are presented.

The crest lengths from the two-dimensional wave data were determined as follows. 1) The crests are identified by finding the areas that meet a threshold criterion, *η* > *η _{o}* or

*s*>

*s*, with an additional constraint that the local wave phase of the bandpassed ATM data are within ±10°, where Θ = 0° corresponds to a local maximum of

_{o}*η*, thus avoiding negative elevations and zero crossings while restricting the analysis to the crests of the dominant waves. 2) The data are binarized where the thresholded areas are assigned to be equal to one and zero elsewhere. 3) Each binary wave front is contoured, the contour length is calculated, and each contour is fitted to an ellipse. 4) The crest length is approximated using the major axis of the ellipse, and the local wavenumber vector is determined from the mean wavenumber vector within the thresholded area. Figure 8 shows an example of the bandpassed filtered data with contours for different thresholds and sample ellipses fitted to the contours where

*η*> 2.2

*η*

_{rms}.

*p*(

*η*,

*α*,

*θ*) =

*p*(

*η*)

*p*(

*α*,

*θ*) is the joint distribution of the surface elevation

*η*and the surface gradient

**∇**

*η*, with

*α*and

*θ*corresponding to the magnitude and direction of

**∇**

*η*, respectively. After integration, Eq. (36) becomeswhere

*γ*= (

*m*

_{02}/

*m*

_{20})

^{1/2}is a measure of the long crestedness of the wave field with

*γ*= 0 and

*γ*= 1 corresponding to long- and short-crested waves, respectively, and

*f*(

*γ*) is an elliptic integral between 1 and 1.2207 for long- and short-crested waves, respectively, with a weak singularity near the origin (

*γ*= 0). Here, the distribution of

*C*is the perimeter of the

_{n}*n*th contour around the area where

*η*>

*η*, and

_{o}*A*

_{tot}the total surface area. The empirical estimates of

*s*can be approximated bywhich neglects the modulations of the local wavenumber and includes a factor of ½ to exclude the contributions from the wave troughs. Equation (40) was obtained by following the derivation of Eq. (37) from the joint pdf of wave amplitudes and gradient of the wave envelope

_{o}*p*(

*a*,

*a*,

_{x}*a*) =

_{y}*p*(

*a*)

*p*(

*a*,

_{x}*a*), with the amplitudes scaled with mean wavenumber

_{y}*p*(

*s*,

*a*,

_{x}*a*) =

_{y}*p*(

*a*,

*a*,

_{x}*a*)|

_{y}*da*/

*ds*|, where

*p*(

*s*) is given in Eq. (35), and

*p*(

*a*,

_{x}*a*) corresponds to Eq. (2.8.15) in Longuet-Higgins (1957, p. 372).

_{y}*n*th crest and

*η*>

*η*, which is estimated as the major axis of the ellipse fitted over the area exceeding the threshold criterion.

_{o}*L*can be analytically approximated byand nondimensionalized according toSimilarly the dimensionless crest-length distribution using wave slope threshold criteria is defined and approximated bywith

_{η}*η*and

_{o}*s*normalized by

_{o}*η*

_{rms}and

## 5. Second-order stochastic simulations

In this section, we investigate the performance of second-order nonlinear Monte Carlo simulations to predict the statistics of crest lengths based on computed directional wavenumber spectra corresponding to research flight 5. This flight was selected because of its strong winds (see Fig. 1) and also because the modeled directional spectra gave the best agreement with the observations (Romero and Melville 2010b). The observed evolution of the significant wave height and dominant wavelength versus the fetch are shown in Figs. 1b,c, respectively. For reference, Figs. 5 and 8 in Romero and Melville (2010a) show examples of the observed directional and one-dimensional wavenumber spectra measured at various fetches.

The numerical wind wave model included forcing by the wind (Snyder et al. 1981), wave dissipation (primarily due to wave breaking), and nonlinear four-wave resonant interactions with the so-called exact method described in van Vledder (2006), based on the work by Tracy and Resio (1982) and Resio and Perrie (1991). The computed directional spectra were validated against the ATM observations, giving good agreement particularly for the lower moments of the spectra, but the directional spectra were in general narrower than the observations by about 10°. Figure 4 in Romero and Melville (2010b) shows a direct comparison between measured and simulated directional wavenumber spectra. For more details on the numerical simulations, see Romero and Melville (2010b).

The approach used for the stochastic simulations follows Forristall (2000). In this study, several linear two-dimensional spatial wave fields, of 5 m by 5 m resolution, approximately covering the same area as the ATM observations, were generated using the computed directional spectra and a uniformly random distribution of wave phases. Then, the second-order nonlinear correction was calculated from the kernels derived by Longuet-Higgins (1963), including a missing factor of ½ (Sharma and Dean 1979). The statistical analysis of the crest length per unit area

*k*<

_{p}*k*

_{1}< 1.75

*k*and |

_{p}*k*

_{2}| < 1.75

*k*) linear spectrum, shown with open squares in Fig. 11. Both narrowband and broadband linear simulations give nearly identical distributions of

_{p}Means and standard errors of the moments of the spectrum where *i* = 0, 1, 2) and *m _{ij}* (

*i*= 0, 1, 2;

*j*= 0, 1, 2), correspond to the moments of the omnidirectional and directional wavenumber spectra given in Eqs. (24) and (3), respectively. The moments about the mean wavenumber

*μ*(

_{ij}*i*= 0, 1, 2;

*j*= 0, 1, 2) are defined in Eqs. (6)–(8). All spectra were bandpassed in the range 0.5

*k*<

_{p}*k*

_{1}< 1.75

*k*and −1.75

_{p}*k*<

_{p}*k*

_{2}< 1.75

*k*, with

_{p}*k*

_{1}corresponding to the direction of propagation of the dominant waves and

*k*being the wavenumber at the spectral peak. The data were aggregated from multiple research flights (5, 7, 9, and 10) with the nonoverlapping bin intervals shown in Fig. 2a along the dimensionless fetch

_{p}*χ*.

The direct comparison between the measured statistics of the crest lengths and the nonlinear stochastic simulations is shown in Fig. 12 with ^{−1}, or about 18-m wavelengths. Since the length of a given crest depends on the slope of the surface, it is expected that *w* = [0.04 0.20 0.04; 0.2 1 0.2; 0.04 0.2 0.04] oriented relative to the flight track, which reduced the noise levels at high wavenumbers while preserving most of the variance (80% on average). Figure 12 shows that the second-order nonlinear simulation of

Means and standard errors of the moments of the measured and simulated spectra corresponding to RF 5, where *m _{ij}* (

*i*= 0, 2;

*j*= 0, 2) are the moments of the directional wavenumber spectrum defined in Eq. (3). The moments of the simulated data include the linear and second-order nonlinear wave fields computed from the full bandwidth and the bandpassed linear wave field with wavenumbers in the ranges 0.5

*k*<

_{p}*k*

_{1}< 1.75

*k*and −1.75

_{p}*k*<

_{p}*k*

_{2}< 1.75

*k*, with

_{p}*k*

_{1}corresponding to the dominant wave direction and

*k*is the wavenumber at the peak of the spectrum. The measured ATM spectra include a low-pass filter applied in the spatial domain to reduce the noise as described in section 5, as well as the full spectrum with the high wavenumber noise included.

_{p}## 6. Discussion and conclusions

The analysis of the measured surface waves showed that the significant slope and elevation skewness are weakly dependent on the dimensionless fetch, with the skewness nearly in agreement with that estimated from the Tayfun distribution (Socquet-Juglard et al. 2005) calculated from the measured significant slope. The kurtosis also shares a similar trend, generally decreasing with increasing nondimensional fetch. However, the measured kurtosis is always larger than that predicted by the narrowband Tayfun distribution. Contrary to the predictions by Mori and Janssen (2006), the maximum wave height observed *H*_{max}/*H _{s}* is not correlated with the Benjamin–Feir index, with

*H*

_{max}/

*H*being nearly independent of fetch. The pdfs of the sea surface elevation calculated from the data show good agreement with the second-order narrowband distribution by Tayfun (1980), except for very large displacements and particularly for negative displacements, which may be attributed to higher order nonlinear effects and the lack of data at such low probabilities. In general, the exceedance probability of wave height calculated in the direction of the dominant waves using standard zero-crossing techniques shows a good agreement with the distributions by T90 and B00, with the Rayleigh distribution (Cartwright and Longuet-Higgins 1956) overestimating the observed wave height probabilities.

_{s}We presented empirical field measurements of the joint distribution of wave amplitudes and wavelengths from bandpassed data, giving a good agreement with the analytical solution derived by Xu et al. (2004), which to our knowledge had not been previously tested against measurements in the field. Following Longuet-Higgins (1957), a linear theoretical distribution of wave slope vectors *p*(*s*_{1}, *s*_{2}) was derived and shown to compare favorably with the measurements. However, the measured distributions of wave slope at large nondimensional fetches gave slightly reduced probabilities of finding steep waves in the direction orthogonal to the dominant waves. The measured distributions of wave slope moduli were shown to compare favorably with theoretical distributions. In particular, the derived distribution of wave slope vectors numerically integrated in azimuth gave excellent agreement with the observations, while the Rayleigh distribution served as an upper bound on the data.

The analytical solution describing the statistics of the length of contours per unit area derived by Longuet-Higgins (1957) was tested for the first time by field measurements. The results are in good agreement with the theory. A similar expression was derived for the lengths of contours bounding areas exceeding wave steepness thresholds. The data and the linear analytical approximation are also in good agreement. Based on an empirical relationship between the length of contour and the crest length, analytical approximations describing the distributions of the crests length with elevation and wave slope thresholds were obtained and shown to be in good agreement with the ATM observations.

The distribution of crest length per unit area *η _{o}*, were simulated using stochastic second-order numerical simulations from computed directional wavenumber spectra with a third-generation wind–wave model using exact computations of the nonlinear four-wave resonant interactions as described by Romero and Melville (2010b). The second-order nonlinear simulations showed good agreement with the observations, in general increasing the crest-length distribution of large waves by about a factor of 2 or more. The simulations also showed that the crest-length distribution is less affected by the bandwidth than by second-order nonlinearities.

## Acknowledgments

We acknowledge the collaboration of Carl A. Friehe and Djamal Khelif at the University of California, Irvine, in planning and conducting the GOTEX experiments and in help with the analysis of the atmospheric boundary layer data. We are grateful to Allen Schanot, Henry Boynton, Lowell Genzlinger, Ed Ringleman, and the support staff at the NCAR Research Aviation Facility. We thank Bill Krabill, Bob Swift, Jim Yungel, John Sonntag, and Robbie Russell at NASA/EG&G for access to the ATM, its deployment, and initial data processing. LR is thankful to Miguel Onorato for useful comments and discussions during preliminary stages of this work. This research was supported by grants to WKM from the National Science Foundation (OCE), the Office of Naval Research (Physical Oceanography).

## APPENDIX A

### Distributions of Local Wave Parameters

*Z*(

**x**) over a horizontal infinite area is defined bywhere

*η*(

**x**) is the sea surface elevation and

*ξ*is the Hilbert transform of

*η*, which is defined bywhere

*P*is the Cauchy principal value. The two-dimensional Hilbert transform has been previously described by Stark (1971), Zhu et al. (1990), and Bülow and Sommer (2001), among others. It can be efficiently obtained from discrete data using the Fourier transform

**k**is the wavenumber vector in Fourier space. Equation (A3) has a 180° ambiguity that, by analogy to the one-dimensional case, can be removed by defining a unit vector

**e**such that the positive wavenumbers correspond to

**e**·

**k**> 0 andwhich neglects the contribution from the removed mean elevation

*a*and phase

*ϵ*are defined byandThe local wavenumber vector

**k**= (

*k*

_{1},

*k*

_{2}) is given by the gradient of the local phase function according towhere the spatial derivatives of

*ϵ*(

**x**) are modulo 2

*π*. Because this use of the Hilbert transform is limited to narrowband processes, the data analyzed were bandpassed in the range 0.5

*k*<

_{p}*k*

_{1}< 1.75

*k*and −1.75

_{p}*k*<

_{p}*k*

_{2}< 1.75

*k*, where

_{p}*k*is the wavenumber at the peak of the wind wave spectrum. The limits where chosen in an attempt to preserve the most variance while excluding the harmonics of the dominant waves. Additional analysis showed that small changes in the limits of the filter did not affect the results.

_{p}*δ*=

*μ*

_{20}

*μ*

_{02}and

*μ*

_{11},

*μ*

_{20}, and

*μ*

_{20}are the first and second spectral moments about the mean wavenumber vector; andis the Jacobian of the transformation. Figure A1 shows

*ak*

_{1},

*ak*

_{2}) = (

*s*

_{1},

*s*

_{2}) and amplitudes, defined aswhereis the Jacobian of the transformation. The integration of Eq. (A13) with respect to

*a*from zero to infinity gives the pdf of wave slope vectorswith

*δ*=

*μ*

_{20}

*μ*

_{02}. The integration was carried out using Mathematica, and it was assumed that the cross-spectral moment

*μ*

_{11}= 0.

## APPENDIX B

### Crest and Contour Lengths

The relationship between the estimates of the crest length *λ _{c}* and the length of contour

*C*for several surface elevation thresholds is investigated from their joint distribution, which is shown in Fig. B1. The data show that the population is dominated by both small crests and small contour lengths, with data distributed about a slope of 0.37, which is shown with a solid black line.

## REFERENCES

Boccotti, P., 2000:

*Wave Mechanics for Ocean Engineering*. Elsevier Science, 496 pp.Bülow, T., , and G. Sommer, 2001: Hypercomplex signals-a novel extension of the analytic signal to the multidimensional case.

,*IEEE Trans. Signal Process.***49**, 2844–2852.Cartwright, D. E., , and M. S. Longuet-Higgins, 1956: The statistical distribution of the maxima of a random function.

,*Proc. Roy. Soc. London***237A**, 212–232.Cavanié, A., , M. Arthan, , and R. Ezraty, 1976: A statistical relationship between individual heights and periods of storm waves.

*Proc. Conf. on Behavior Offshore Structures,*Trondheim, Norway, Norwegian Institute of Technology, 354–360.Cox, C., , and W. Munk, 1954: Statistics of the sea surface derived from sun glitter.

,*J. Mar. Res.***13**, 198–227.Dysthe, K. B., , H. E. Krogstad, , and P. Muller, 2008: Oceanic rogue waves.

,*Annu. Rev. Fluid Mech.***40**, 287–310.Forristall, G., 1978: On the statistical distribution of wave heights in a storm.

,*J. Geophys. Res.***83**, 2353–2358.Forristall, G., 2000: Wave crest distributions: Observations and second-order theory.

,*J. Phys. Oceanogr.***30**, 1931–1943.Huang, N., , S. R. Long, , and L. Bliven, 1981: On the importance of the significant slope in empirical wind-wave studies.

,*J. Phys. Oceanogr.***11**, 569–573.Hwang, P. A., , and D. W. Wang, 2001: Directional distributions and mean square slopes in the equilibrium and saturation ranges of the wave spectrum.

,*J. Phys. Oceanogr.***31**, 1346–1360.Hwang, P. A., , D. W. Wang, , E. J. Walsh, , W. B. Krabill, , and R. N. Swift, 2000a: Airborne measurements of the wavenumber spectra of ocean surface waves. Part I: Spectral slope and dimensionless spectral coefficient.

,*J. Phys. Oceanogr.***30**, 2753–2767.Hwang, P. A., , D. W. Wang, , E. J. Walsh, , W. B. Krabill, , and R. N. Swift, 2000b: Airborne measurements of the wavenumber spectra of ocean surface waves. Part II: Directional distribution.

,*J. Phys. Oceanogr.***30**, 2768–2787.Janssen, P. A. E. M., 2005: Nonlinear four-wave interaction and freak waves.

*Rogue Waves: Proc. ‘Aha Huliko’a Hawaiian Winter Workshop,*Honolulu, HI, University of Hawaii at Manoa, 85–90Krabill, W. B., , and C. F. Martin, 1987: Aircraft positioning using global positioning system carrier phase data.

,*Navigation***34**, 1–21.Longuet-Higgins, M. S., 1952: On the statistical distributions of sea waves.

,*J. Mar. Res.***11**, 245–265.Longuet-Higgins, M. S., 1957: The statistical analysis of a random, moving surface.

,*Proc. Roy. Soc. London***249A**, 321–387.Longuet-Higgins, M. S., 1963: The effect of non-linearities on statistical distributions in the theory of sea waves.

,*J. Fluid Mech.***17**, 459–480.Longuet-Higgins, M. S., 1975: On the joint distribution of wave periods and amplitudes of sea waves.

,*J. Geophys. Res.***80**, 2688–2694.Longuet-Higgins, M. S., 1980: On the distribution of the heights of sea waves: Some effects of nonlinearity and finite band width.

,*J. Geophys. Res.***85**(C3), 1519–1523.Longuet-Higgins, M. S., 1983: On the joint distribution of wave periods and amplitudes in a random wave field.

,*Proc. Roy. Soc. London***389A**, 241–258.Melville, W. K., , L. Romero, , and J. M. Kleiss, 2005: Extreme waves in the Gulf of Tehuantepec.

*Rogue Waves: Proc. ‘Aha Huliko‘a Hawaiian Winter Workshop,*Honolulu, HI, University of Hawaii at Manoa, 23–28.Mori, N., , and P. A. E. M. Janssen, 2006: On kurtosis and occurrence probability of freak waves.

,*J. Phys. Oceanogr.***36**, 1471–1483.Mori, N., , M. Onorato, , P. A. E. M. Janssen, , A. R. Osborne, , and M. Serio, 2007: On the extreme statistics of long-crested deep water waves: Theory and experiments.

,*J. Geophys. Res.***112**, C09011, doi:10.1029/2006JC004024.Naess, A., 1985: On the distribution of crest to trough wave heights.

,*Ocean Eng.***12**, 221–234.Resio, D., , and W. Perrie, 1991: A numerical study of nonlinear energy fluxes due to wave-wave interactions. Part 1: Methodology and basic results.

,*J. Fluid Mech.***223**, 609–629.Romero, L., , and W. K. Melville, 2010a: Airborne observations of fetch-limited waves in the Gulf of Tehuantepec.

,*J. Phys. Oceanogr.***40**, 441–465.Romero, L., , and W. K. Melville, 2010b: Numerical modeling of fetch-limited waves in the Gulf of Tehuantepec.

,*J. Phys. Oceanogr.***40**, 466–486.Romero, L., , and W. K. Melville, 2010c: Observations and modeling of linear and nonlinear spatio-temporal surface wave statistics.

*Proc. 29th Int. Conf. on Offshore Mechanics and Arctic Engineering,*Shanghai, China, ASME, 1–13.Sharma, J. N., , and R. G. Dean, 1979: Development and evaluation of procedure for simulating a random directional second order sea surface and associated wave forces. University of Delaware Ocean Engineering Rep. 20, 112 pp.

Shum, K. T., , and W. K. Melville, 1984: Estimates of the joint statistics of amplitudes and periods of ocean waves using an integral transform.

,*J. Geophys. Res.***89**(C4), 6467–6476.Snyder, R. L., , F. Dobson, , J. Elliott, , and R. B. Long, 1981: Array measurements of atmospheric pressure fluctuations above surface gravity waves.

,*J. Fluid Mech.***102**, 1–59.Socquet-Juglard, H., , K. Dysthe, , K. Trulsen, , H. E. Krodgstad, , and J. D. Liu, 2005: Probability distributions of surface gravity waves during spectral changes.

,*J. Fluid Mech.***542**, 195–216.Stark, H., 1971: An extension of the Hilbert transform product theorem.

,*Proc. IEEE***59**, 1359–1360.Steenburgh, W. J., , D. M. Schultz, , and B. A. Colle, 1998: The structure and evolution of gap outflow over the Gulf of Tehuantepec, Mexico.

,*Mon. Wea. Rev.***126**, 2673–2691.Tayfun, M. A., 1980: Narrow-band nonlinear sea waves.

,*J. Geophys. Res.***85**(C3), 1543–1552.Tayfun, M. A., 1981: Distribution of crest-to-trough wave heights.

,*J. Waterw. Port Coastal Ocean Div.***107**, 149–158.Tayfun, M. A., 1983: Nonlinear effects on the distribution of crest-to-trough wave heights.

,*Ocean Div.***10**, 97–106.Tayfun, M. A., 1990: Distribution of large wave heights.

,*J. Waterw. Port Coast. Ocean Div.***116**, 686–707.Tayfun, M. A., , and F. Fedele, 2007: Wave height distributions and nonlinear effects.

,*Ocean Eng.***34**, 1631–1649.Tracy, B. A., , and D. T. Resio, 1982: Theory and calculation of the nonlinear energy transfer between sea waves in deep water. Army Engineer Waterways Experiment Station Vicksburg MS Hydraulics Lab Tech. Rep. 11, 54 pp.

van Vledder, G. P., 2006: The WRT method for the computation of non-linear four-wave interactions in discrete spectral wave models.

,*Coastal Eng.***53**, 223–242.Walsh, E. J., , D. W. Hancock, , D. E. Hines, , R. N. Swift, , and J. F. Scott, 1985: Directional wave spectra measured with the surface contour radar.

,*J. Phys. Oceanogr.***15**, 566–592.Xu, D., , X. Li, , L. Zhang, , N. Xu, , and H. Lu, 2004: On the distributions of wave periods, wavelengths, and amplitudes in a random wave field.

,*J. Geophys. Res.***109**, C05016, doi:10.1029/2003JC002073.Zhu, Y. M., , F. Peyrin, , and R. Goutte, 1990: The use of a two-dimensional Hilbert transform for Wigner analysis of 2-dimensional real signals.

,*Signal Process.***19**, 205–220.

^{1}

This ambiguity is often resolved in real data by additional environmental information, including the wind direction and the source of the swell.