# Ocean Wave Directional Distribution from GPS Buoy Observations off the West Coast of Ireland: Assessment of a Wavelet-Based Method

Daniel Peláez-Zapata, Vikram Pakrashi, Frédéric Dias

► To cite this version:

Daniel Peláez-Zapata, Vikram Pakrashi, Frédéric Dias. Ocean Wave Directional Distribution from GPS Buoy Observations off the West Coast of Ireland: Assessment of a Wavelet-Based Method. Journal of Atmospheric and Oceanic Technology, 2024, 41 (8), pp.749 - 765. 10.1175/jtech-d-23-0058.1. hal-05163980

HAL Id: hal-05163980

https://hal.science/hal-05163980v1

Submitted on 15 Jul 2025

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
Ocean Wave Directional Distribution from GPS Buoy Observations off the
West Coast of Ireland: Assessment of a Wavelet-Based Method

DANIEL PELÁEZ-ZAPATA<sup>a,b</sup> VIKRAM PAKRASHI,<sup>c</sup> AND FRÉDÉRIC DIAS<sup>a,b</sup>

(Manuscript received 23 April 2023, in final form 3 May 2024, accepted 6 June 2024)

**ABSTRACT:** Knowledge of the directional distribution of a wave field is crucial for a better understanding of complex air-sea interactions. However, the dynamic and unpredictable nature of ocean waves, combined with the limitations of existing measurement technologies and analysis techniques, makes it difficult to obtain precise directional information, leading to a poor understanding of this important quantity. This study investigates the potential use of a wavelet-based method applied to GPS buoy observations as an alternative approach to the conventional methods for estimating the directional distribution of ocean waves. The results indicate that the wavelet-based estimations are consistently good when compared to the framework of widely used parameterizations for the directional distribution. The wavelet-based method presents advantages in comparison with the conventional methods, including being purely data-driven and not requiring any assumptions about the shape of the distribution. In addition, it was found that the wave directional distribution is narrower at the spectral peak and broadens asymmetrically at higher and lower scales, particularly sharply for frequencies below the peak. The directional spreading appears to be independent of the wave age across the entire range of frequencies, implying that the angular width of the directional spectrum is primarily controlled by nonlinear wave-wave interactions rather than by wind forcing. These results support the use of the wavelet-based method as a practical alternative for the estimation of the wave directional distribution. In addition, this study highlights the need for continued innovation in the field of ocean wave measuring technologies and analysis techniques to improve our understanding of air-sea interactions.

SIGNIFICANCE STATEMENT: This study presents a wavelet-based technique for obtaining the directional distribution of ocean waves applied to GPS buoy. This method serves as an alternative to conventional methods and is relatively easy to implement, making it a practical option for researchers and engineers. The study was conducted in a highly energetic environment characterized by high wind speeds and large waves, providing a valuable dataset for understanding the dynamics of marine environment in extreme conditions. This research has implications for improving our understanding of directional characteristics of ocean waves, which is crucial for navigation, offshore engineering, weather forecasting, and coastal hazard mitigation. This study also highlights the challenges associated with understanding wave directionality and emphasizes a need for further observations.

**KEYWORDS:** Wind waves; Air-sea interaction; Buoy observations; In situ oceanic observations; Spectral analysis/models/distribution; Wavelets

**1. Introduction**

The directional wave spectrum is fundamental to numerical wave modeling. It has a wide range of applications in science and engineering, such as air-sea interactions, wave climate, sea-state forecasting, microseism prediction, coastal erosion, and wave energy harvesting (e.g., Lenain and Melville 2017; Hauser et al. 2005; Romero et al. 2012). Wave buoys have traditionally been the most extended instrument to measure ocean wave parameters. In recent years, there has been a surge in the use of real-time, GPS wave buoys (Lancaster et al. 2021; Raghukumar et al. 2019), which can measure the three-

dimensional wave motion at relatively low cost and with simple deployments. Nevertheless, extracting directional data from the wave field is not straightforward, particularly when only information from a single point is available (Kuik et al. 1988; Hauser et al. 2005). This is the reason why a number of methods based on different approaches have been proposed and followed over the last 4 decades with no consensus on what might be the most reasonable approach. Here, we present a novel analysis technique and compare its performance with conventional methods using a dataset obtained from a GPS wave buoy.

The directional wave spectrum $E(f, \theta)$ is often expressed as the product of two functions: the frequency spectrum and a directional distribution function, which are dependent on both frequency and direction, i.e., $E(f, \theta) = S(f)D(f, \theta)$, where f is the wave frequency and θ is the wave direction. The primary challenge lies in determining the directional distribution function, which is required to satisfy the following constraint:

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

<footnote><sup>a</sup> Centre Borelli, École Normale Supérieure Paris-Saclay, Gif-sur-Yvette, France</footnote>

<footnote><sup>b</sup> School of Mathematics and Statistics, University College Dublin, Dublin, Ireland</footnote>

<footnote><sup>c</sup> Dynamical Systems and Risk Laboratory, School of Mechanical and Materials Engineering, University College Dublin, Dublin, Ireland</footnote>

<footnote>Denotes content that is immediately available upon publication as open access.</footnote>

<footnote>Corresponding author: Daniel Peláez-Zapata, daniel.pelaez-zapata@ucdconnect.ie</footnote>
$$ \int_{-\pi}^{\pi} D(f, \theta) d\theta = 1. \qquad (1) $$

The most direct way to obtain the directional distribution of ocean waves is by means of Fourier series expansion, i.e.,

$$ D(f, \theta) = \frac{1}{2\pi} + \frac{1}{\pi} \sum_{n=1}^{\infty} (a_n \cos n\theta + b_n \sin n\theta), \quad (2) $$

where

$$ \begin{aligned} a_n(f) &= \int_{-\pi}^{\pi} D(f, \theta) \cos n\theta \, d\theta \quad \text{and} \\ b_n(f) &= \int_{-\pi}^{\pi} D(f, \theta) \sin n\theta \, d\theta \end{aligned} \qquad (3) $$

are the Fourier coefficients or circular moments of the directional distribution $D(f, \theta)$. These coefficients can be obtained from buoy records through the conventional cross-spectral analysis (Kuik et al. 1988). However, the series is generally truncated to the number of terms that are possible to be determined with the available data. For instance, GPS and heave-pitch-roll buoys can provide three perpendicular quantities; hence, only the first five Fourier coefficients can be determined.

This approach tends to produce poor and inconsistent directional distributions (Benoit et al. 1997). Different alternatives to increase the number of Fourier coefficients have been proposed. Some of them include expanding the capabilities of a single buoy to be able to measure the sea surface curvature (Mitsuyasu et al. 1975) or processing the available data to approximate the actual curvature (Gryazin and Gleb 2022). However, in both cases, only nine Fourier coefficients can be determined, which seems to be insufficient to describe the directional distribution, particularly in complex seas. Some attempts have been made to fit the directional distribution $D(f, \theta)$ to unimodal and bimodal parametric models. However, they are not recommended because their estimates are of poor quality (Benoit et al. 1997). More recently, Plant and Donelan (2020) proposed that the directional wave spectrum can be estimated from pitch-roll buoys through an iterative deconvolution method. However, this method requires prior knowledge of a parametric model for the directional distribution, which represents a considerable disadvantage.

The conventional methods that have been routinely used to estimate the directional wave spectrum are based on finding the best fit to the directional data. They can be classified into three families. The first one is the maximum-likelihood-based method (Capon 1969), of which the most representative version is the iterative maximum likelihood (IMLM). This method improves the estimation using an adaptive technique (Oltman-Shay and Guza 1984). It is known to be fast and relatively easy to implement and seems to be an acceptable compromise between accuracy and computational time (Benoit 1992). However, according to Donelan et al. (2015), it generally produces too broad directional spreadings. The second family of methods is the maximum-entropy-based methods, which have been proven to be suitable for the estimation of

the directional distribution function, by using different definitions of entropy, such as Burg (Lygre and Krogstad 1986) or Shannon (Hashimoto et al. 1994). The most outstanding member of this family is the so-called extended maximum entropy principle (EMEP) which seems to give the most reliable estimate but is computationally demanding (Benoit 1993). Furthermore, Hauser et al. (2005) mentioned that Burg-based maximum entropy method (MEM) estimates produce a double-peaked distribution when the Fourier coefficients of a unimodal distribution are used as input. This was further confirmed by Simanesew et al. (2018), who reported that both the Burg and Shannon MEMs have a tendency to split peaks of the directional distribution function. Finally, Hashimoto and Kobune (1988) introduced the Bayesian directional method (BDM) which predicts the directional distribution based on Bayesian inference principles and considers the constraints of the spectral cross-correlation coefficients. This method has been recognized as one of the most promising for estimating the directional wave spectrum from an array of sensors. However, it is unstable for single-point measurements where only few information is available. Further details can be found in Benoit (1992), Benoit (1993), and Benoit et al. (1997), who performed a comprehensive comparison of these methods using realistic numerical simulations and laboratory experiments.

Given the uncertainties in estimating methods, the width of the directional distribution (commonly known as spreading) is poorly understood and difficult to model accurately. It is generally accepted that this parameter depends on the wave frequency. Simanesew et al. (2018) found evidence that frequency-dependent spreading may naturally develop from an initial state of no frequency dependence. The vast majority of observations show narrow distributions near the spectral peak and a broadening as the frequency moves away toward both lower and higher scales (e.g., Mitsuyasu et al. 1975; Hasselmann et al. 1980; Donelan et al. 1996; Ewans 1998; Forristall and Ewans 1998; among others). However, the exact form of this frequency dependence remains a subject of investigation (Latheef et al. 2017). The mechanism that explains this behavior is unclear. The most accepted explanation is attributed to the nonlinear wave-wave interactions (Krogstad 2001; Romero and Lubana 2022). However, Lin et al. (2021) suggest that low signal-to-noise ratios may play an important role in the observed increase in directional spreading, particularly above the spectral peak.

In an effort to address the limitations of conventional methods, Donelan et al. (1996) introduced the wavelet directional method (WDM). They argued that the conventional methods require assumptions on the stationarity of the wave field that becomes questionable if the sea state is changing rapidly as in the case of a storm passage. The WDM leverages the wavelet time-frequency decomposition obtaining the wave direction from the phase difference in the wavelet transforms associated with the surface elevation at three or more spatially distributed wave gauges. The resulting wave direction is a function of time and frequency; thus, the directional distribution is calculated by estimating the density along the time dimension associated with each frequency. In general, the WDM has been proven to be superior in comparison with the
<figure><img src="image_3.png" /><figcaption>FIG. 1. Location of the Spotter buoy deployed off the west coast of Ireland.</figcaption></figure>

traditional methods previously mentioned (Donelan et al. 1996; Krogstad et al. 2006; Donelan et al. 2015; Plant and Donelan 2020). However, it is essential to note that WDM is applicable exclusively to spatial arrays of wave elevation records and not to wave buoys.

Furthermore, Krogstad et al. (2006) proposed to extend WDM to triplet data (such as heave-pitch-roll buoys) by replacing the Fourier coefficients in the conventional cross-spectral analysis by the corresponding wavelet coefficients. This analysis produces local estimates of the circular moments and wave direction as a function of frequency and time; hence, it is possible to obtain the directional distribution function. They compared their wavelet-based analysis with the original WDM using an array of laser measurements, from which they derived an equivalent triplet dataset of heave/slope through linear spatial interpolation. However, this comparison was never carried out using buoys. Hence, there is a strong motivation to extend and apply this methodology to single-point triplet data, in particular to GPS buoys. The widespread adoption of these buoys for standard wave measurements in recent years highlights the importance of improving and refining directional estimation methods to enhance the accuracy and reliability of wave data.

In this article, we present a complementary extension of the WDM, originally proposed by Donelan et al. (1996), for the estimation of the directional wave spectrum using wave buoy data. In particular, we used a GPS wave buoy deployed off the west coast of Ireland. Moreover, our approach, unlike Krogstad et al. (2006) extension, uses wavelet-based correlations of surface elevation and eastward and northward components of the horizontal velocity, thereby providing a more intuitive estimation of the directional distribution function. Additionally, our method is presented in a manner that enhances clarity and understanding, enabling easier adoption and implementation within the physical oceanography community. The paper is structured as follows: Section 2

introduces the methodology and describes the field observations. Section 3 presents a comparison of the wavelet-based directional distribution function with the conventional methods. In addition, an analysis on the shape of the directional distribution, wave direction, and spreading parameters is presented. Finally, section 4 summarizes the main findings and presents the concluding remarks.

## 2. Methods and dataset description

### a. Field observations

In this work, data from a GPS wave buoy, commercially known as Spotter and manufactured by SOFAR, were used (Raghukumar et al. 2019; Lancaster et al. 2021). GPS wave buoys typically measure the Doppler-shifted radial velocity signal between the buoy and the satellite as the buoy moves with the passing waves. The satellite orbit position and speed are known at all times from the ephemerides; therefore, the difference between the velocity derived from the Doppler shift and the satellite radial velocity is caused by the buoy motion. As these measurements typically yield only radial velocities, a satellite positioned near the nadir can provide an estimate of the buoy vertical velocity, while a satellite positioned closer to the horizon can offer insight into the buoy horizontal motion. In addition, a minimum of four satellites are required simultaneously in order to compute the velocity of the buoy in longitude, latitude, and vertical (Krogstad et al. 1999). Spotter buoys measure continuously at a sampling rate of 2.5 Hz. The buoy was deployed off the coast of Inis Mór, Aran Islands, in the west coast of Ireland (53.1486°N, 9.8872°W), approximately 1.6 km from the shoreline, in a water depth of around 50 m (Fig. 1).

The Spotter buoy continuously recorded data during 22 days, from 13 August to 5 September 2020. The analysis is performed over consecutive 30-min blocks (resulting in 1056 blocks being analyzed). We employed a standard
criterion to spike removal and data quality control: If the surface elevation exceeds four times the standard deviation of the 30-min block, that value is labeled as invalid. Linear interpolation was employed to fill gaps in the time series. Blocks with gaps exceeding 5 s or containing more than 10% of invalid data were disregarded and removed from the analysis (12 out of 1056 blocks were excluded). Wind speed data were obtained from the hourly ERA5 reanalysis at a nearby point. In addition, we obtained directional wave spectra from the high-resolution hindcast ResourceCODE. This database consists of 28-yr hourly wave bulk parameters and directional wave spectra generated using the third-generation wave model WAVEWATCH III (WW3) for the North Atlantic European waters (Accensi et al. 2021, 2022; Raillard et al. 2023).

During the measurement period (Fig. 2), the average sea-state conditions were characterized by significant wave heights around 1.6 m, peak wave periods ranging from 4 to 12 s, and peak wave directions coming mostly from the west. During the passage of Storm Ellen, between 19 and 21 August, the significant wave height hit a maximum of 4.5 m. In addition, a rapid increase in wind speed was recorded, going from 9 to 24 m s⁻¹ in 1 h. This site is well known for storm waves with significant wave heights reaching up to 12 m (Fedele et al. 2019). These waves are even capable of moving boulders weighing tens or hundreds of tons (Cox et al. 2020; Kennedy et al. 2021). Additionally, the typical wave field is characterized by complex mixed seas, with several swell systems coexisting along with locally generated wind seas, making the analysis of the wave directionality even more challenging. The complex nature of the study site and the rapidly evolving sea-state conditions provide a valuable opportunity to compare the effectiveness of wavelet analysis with conventional spectral analysis.

b. Wavelet-based cross correlation

The continuous wavelet transform (CWT), which maps a 1D time series onto a 2D time-frequency domain, is well suited for describing nonstationary processes such as ocean surface waves (Massel 2001; Liberzon et al. 2019; Ewans 2015; Donelan et al. 1996). This transformation can be seen in the time domain as the convolution of a time series x(t) with a set of special functions, known as wavelet function or “mother” wavelet. These functions must satisfy the following conditions: finite energy, localization in both time and frequency, and zero mean (Farge 1992; Torrence and Compo 1998). The complex wavelet coefficients can be written as

$$W(\tau, \sigma) = \int x(t) \frac{1}{\sqrt{\sigma}} \psi^* \left( \frac{t-\tau}{\sigma} \right) dt, \quad (4)$$

where σ and τ are the scale (dilation) factor and translation (shift) factor, respectively. The star denotes complex conjugate. One of the most commonly employed wavelet functions is the Morlet wavelet, which has been found useful in the analysis of ocean waves (Donelan et al. 1996; Krogstad et al. 2006). It is defined as a sinusoidal wave localized by a Gaussian window:

$$\psi\left(\frac{t-\tau}{\sigma}\right) = \pi^{-1/4} e^{i\omega_0(t-\tau)/\sigma} e^{-(t-\tau)^2/2\sigma^2}, \quad (5)$$

<figure><img src="image_4.png" /><figcaption>FIG. 2. Ocean wave conditions during the measurement period: (a) significant wave height, (b) peak period, (c) peak wave direction, (d) peak directional spreading, and (e) 10-m wind speed from ERA5.</figcaption></figure>

where ω₀ is the nondimensional frequency parameter. We used ω₀ = 6 as a convenient value to satisfy the admissibility condition, ensuring that the integral of the mother wavelet is close to zero (Farge 1992; Torrence and Compo 1998). Furthermore, this value has been widely adopted in the literature as it provides a good balance between time and frequency localization and is well suited for wavelet-based feature extraction (Grinsted et al. 2004). In the context of ocean waves, selecting ω₀ = 6 ensures that the wavelet scale closely approximates the corresponding Fourier wavelength.

The cross-wavelet transform W<sub>xy</sub> of two time series, x(t) and y(t), respectively, can be defined as the conventional approach for obtaining the Fourier-based cross spectrum, i.e., as “the product between the transform of the first signal and the complex conjugate of the transform of the second signal,” i.e.,

$$W_{xy}(\tau, \sigma) = W_x(\tau, \sigma)W_y^*(\tau, \sigma). \quad (6)$$

This quantity reveals valuable information about the common power and relative phase of both signals in the time-frequency space (Grinsted et al. 2004). For practical purposes, we can write the shift factor τ as the time dimension t. Similarly, we can write the scale factor σ in terms of the wave frequency f as σ = 1/(fω<sub>λ</sub>), where ω<sub>λ</sub> ≈ 1.033 (with ω₀ = 6) is a factor related to the Fourier wavelength of the Morlet wavelet (Torrence and Compo 1998). This allows us to write the wavelet transform quantities in terms of the local time dimension and the wave frequency.

Selecting a frequency range that covers the relevant scales of variability in the data is crucial. The maximum resolvable frequency is determined by the Nyquist–Shannon sampling theorem, which specifies that it cannot exceed half of the sampling rate, while the minimum resolvable frequency is
<figure><img src="image_5.png" /><figcaption>FIG. 3. Representative case of wavelet spectrogram and local wave direction computed for the buoy records corresponding to 1000 UTC 20 Aug 2020. (a) Sea surface elevation. (b) Wavelet power where the color scale represents the spectral power density and the black arrows correspond to the local wave direction in terms of the wave frequency (vertical axis) and local time (horizontal axis). To enhance clarity, only the arrows representing energy that exceeds 10% of the maximum energy are shown. (c) Time-integrated wavelet power density (black line) compared to the Fourier power density (gray line).</figcaption></figure>

determined by the duration of the time series being analyzed. It is convenient to write the frequency scales as fractional powers of 2. For many applications in ocean waves, the minimum and maximum frequencies can be set to 2⁻⁵ and 1 Hz, respectively. The frequency resolution is determined by the number of voices per octave, which is the number of subdivisions between two powers of 2. In this study, we used 16 voices per octave to achieve a suitable frequency resolution while maintaining computational efficiency.

### c. Estimation of the local wave direction

The concept of directional wave spectrum is commonly used to describe, in a statistical sense, the randomness of the sea surface in the presence of ocean waves. It is widely accepted to represent the sea surface elevation as a superposition of cosine waves with different wavelengths, directions, and random phases (Hauser et al. 2005). In specific scenarios, such as waves in very shallow waters, waves in the presence of strong currents, or wave interaction with structures, nonlinear effects could be potentially important. However, for the vast majority of applications, the use of linear theory offers a pragmatic and feasible approach to estimate wave characteristics with reasonable accuracy.

The linear wave theory predicts that a floating particle, such as a wave buoy, will move in a closed, circular, or elliptical orbit in response to the wave-induced velocity field. Specifically, the orbital velocities in the wave troughs always align with the up-wave direction, while the velocities in the crests align with the down-wave direction of wave propagation (see

Holthuijsen 2010, p. 121). In simpler terms, the horizontal velocity and surface elevation signals are in phase. Exploiting this relationship allows us to calculate the local wave direction using the velocity signal, which is covariant with the surface elevation, i.e.,

$$ \theta(f, t) = \tan^{-1} \left[ \frac{\Re\{W_{v\eta}(f, t)\}}{\Re\{W_{u\eta}(f, t)\}} \right], \quad (7) $$

where $\theta(f, t)$ represents the wave direction at a specific frequency $f$ and time $t$. The term $\Re\{\cdot\}$ denotes the real part of the cross-wavelet transform between the sea surface elevation $\eta$ and the wave-induced horizontal velocity components, denoted as $W_{v\eta}(f, t)$ and $W_{u\eta}(f, t)$, respectively.

It is worth mentioning that the same result is achieved by replacing the Fourier-based coefficients of the conventional cross-spectral analysis of triplet data from single-point measurements (Kuik et al. 1988) by those equivalent wavelet-based coefficients. This substitution is possible considering the statistical properties of the Morlet wavelet as it was previously demonstrated by Krogstad et al. (2006). It would produce a version of the directional moments that depends on both frequency and time. The directional distribution function $D(f, \theta)$ is then obtained by estimating the density of the wave direction along the time dimension per each frequency scale. For the purpose of this study, a kernel density estimation (KDE) has been implemented (see appendix).

Figure 3 shows a representative case of the wavelet power density spectrogram and the local wave direction for a 30-min
block corresponding to 1000 UTC 20 August 2020. The recorded sea surface elevation is shown in the upper panel (Fig. 3a). The crest-to-trough wave height was around 4.5 m, and it varied according to the passage of wave groups, which can also be observed in the wavelet power density spectrogram (Fig. 3b) as the brightest colors. Most of the energy is concentrated in a frequency band centered around $f = 2^{-3}$ Hz. The black arrows show that the waves in this frequency band are traveling primarily toward the northeast. Wind-driven shorter waves, with frequencies on the order of $f = 2^{-2}$ Hz, turn their propagation direction approximately toward the north-northwest, following the wind direction, which is consistent with the ERA5 reanalysis that reports an average wind speed of 13 m s⁻¹ with a mean direction of 98° (the angles reported throughout the text follow a Cartesian convention; i.e., they are measured counterclockwise and indicate the direction of wind or wave propagation). Additionally, the wavelet spectrogram reveals some energy in the infragavity wave frequency band (the order of $f = 2^{-4}–2^{-5}$ Hz), localized around 1015–1030 UTC. It is not completely clear whether this energy corresponds to ongoing/reflected infragavity waves or it is only contamination of the buoy measurements due to the influence of the mooring, which may affect the lowest frequencies of the wave spectrum (Draycott et al. 2022). Finally, as shown in Fig. 3c, the time-integrated wavelet spectrum is smoother in comparison with the Fourier-based spectrum. This aligns with the trade-off between time and frequency resolution that is inherent of the wavelet transform. As wavelets provide better localization in the time domain, it necessarily leads to a decrease in frequency resolution. In particular, the frequency resolution decreases as the frequency increases. This consequently produces smoother spectra. Nevertheless, the integral of both curves remains the same.

### d. Implementation of the wavelet-based method

The step-by-step implementation of the wavelet-based method to obtain the directional distribution function from buoy observations can be summarized as follows:

1) Collect the time series of the sea surface elevation $\eta(t)$ and the horizontal velocities: the eastward component $u(t)$ and the northward component $v(t)$. The duration of the time series should be sufficient to capture the relevant temporal variability, typically ranging between 10 and 30 min.

2) Estimate the complex wavelet coefficients of each time series, specifically $W_\eta(f, t)$, $W_u(f, t)$, and $W_v(f, t)$. For a practical guide on the implementation of the CWT, see Torrence and Compo (1998). This step involves the choice of the mother wavelet, central frequency, and frequency range. These parameters will determine the spectral resolution and ultimately the level of smoothing in the directional distribution function.

3) Estimate the wavelet-based cross correlation of the velocity components and the surface elevation, $W_{u\eta}(f, t)$ and $W_{v\eta}(f, t)$, using Eq. (6). These values will provide an idea of the amount of waves propagating eastward and northward, respectively.

4) Compute the local wave direction as a function of frequency and time $\theta(f, t)$ using Eq. (7).

5) Calculate the probability density function using a histogram or KDE for wave directions per each frequency scale and construct the directional distribution function. The choice of angular resolution is made at this step. Typical angular resolutions range between $\Delta\theta = 5°–15°$.

## 3. Results and discussion

### a. Directional distribution function

For comparison, the directional distribution function $D(f, \theta)$ was computed using the most widely used methods designed for buoy data, namely, direct Fourier transform method (DFTM), IMLM, EMEP, and Bayesian direct method (BDM), all of those implemented in the DIrectional WAve SPectrum analysis (DIWASP) toolbox provided by Johnson (2002). The number of iterations for the iterative methods (IMLM, EMEP, and BDM) was set to 10 to achieve a balance between result quality and computational efficiency. The Fourier-based cross spectra were estimated over 30-min blocks (4500 data points) using Welch's approach. These blocks were divided into smaller segments, each containing 256 data points. This results in 17 nonoverlapping segments, leading to 34 degrees of freedom. Consequently, the frequency bandwidth was determined to be $\Delta f = 0.0098$ Hz, providing an indication of the frequency resolution.

The wavelet-based directional distribution was computed by estimating the density of the wave direction along the time dimension per each frequency scale, as described in section 2. The angular resolution was set to $\Delta\theta = 5°$ for all cases. The frequency scales are nonlinearly spaced increasing in powers of 2 (octaves), from $2^{-5}$ to 1 Hz, with 16 voices per octave. This adds up to 81 frequency components. Unlike Fourier-based cross spectra, the wavelet frequency bandwidth is not a constant value but rather a function of the wavelet scale. Within the chosen frequency range, the wavelet frequency bandwidth, understood as the difference between consecutive frequency points, varies from $\Delta f = 0.0014$ Hz for the larger scales to $\Delta f = 0.0424$ Hz for the shorter scales, with an average value of $\Delta f \approx 0.0121$ Hz. Although this variation in bandwidth with scale makes it difficult to directly compare with Fourier-based approaches, the average wavelet bandwidth closely matches the Fourier bandwidth. This enables a meaningful comparison, in terms of frequency resolution and level of smoothing, between the wavelet-based approach and the conventional methods.

Figure 4 illustrates a typical example of the polar representation of the directional distribution function $D(f, \theta)$ obtained by the wavelet-based method (Fig. 4a), the conventional Fourier-based methods (Figs. 4b–e), and the output of a numerical simulation (Fig. 4f). All methods present a similar structure, i.e., a narrow wave system with a frequency around $f = 0.1$ Hz traveling northeast and a secondary, much broader, wind-driven wave system, which travels mainly toward the northwest and becomes broader at shorter scales.
<figure><img src="image_6.png" /><figcaption>FIG. 4. Polar representation of the directional spreading function $D(f, \theta)$ corresponding to 1000 UTC 20 Aug 2020 (same run as in Fig. 3). (a) The results of the wavelet-based method, (b)–(e) the conventional Fourier-based methods, and (f) the corresponding directional distribution function obtained from the ResourceCODE (WW3) hindcast. Black arrows point toward the wind direction.</figcaption></figure>

The wavelet-based method presents a well-defined and smoother directional distribution in comparison with the conventional methods. It shows a smooth transition from the wave direction associated with the peak to the waves at shorter scales traveling toward the wind direction, as well as a consistent tendency to increase the spreading as the frequency increases. This characteristic was also noticed by Donelan et al. (1996) when comparing WDM and MLM using wave staff data. Clearly, the DFTM, which is based on the Fourier expansion truncated up to the second term in Eq. (2), has a tendency to overspread the energy in direction. As discussed by Benoit et al. (1997), this is most likely a consequence of the insufficient number of terms used in the Fourier expansion. Similarly, the IMLM, EMEP, and BDM present comparable results. However, they present important random fluctuations and seem to contain some spurious spikes at high frequencies, showing waves that appear to propagate in the opposite direction of wind. Although BDM agrees well with the other methods in this particular case, we will not include it in our analysis because, as pointed out by Benoit (1992), its use is not advisable for buoy data, where only information in a single point is available. The numerical results of WW3 exhibit poor results regarding the directional spreading, showing a very narrow distribution in comparison with the observations.

Slicing the directional distribution across directions at several values of the nondimensional frequency $f/f_p$ offers a clearer perspective of the main differences associated with

each method. As illustrated by Fig. 5, all distributions show similar patterns across the entire range of frequencies, except those estimated by DFTM, which exhibit excessively broad spreadings. Near the spectral peak, the estimates of EMEP and wavelets are in good agreement. However, for $f/f_p > 2.5$, there are some discrepancies in the mean wave direction. Additionally, both EMEP and IMLM tend to develop an inconsistent double-peak distribution (separated by approximately 180°), with one peak in the direction of the wind and another peak in the opposite direction. This behavior was also noticed by Donelan et al. (1996) and seems to be most likely an artifact of the statistical fitting due to the limited number of data.

Furthermore, we have noticed that EMEP occasionally fails to accurately determine the directional distribution at certain values of $f/f_p$, producing excessively too narrow and asymmetric shapes and sometimes even spreading the energy uniformly across all directions. We have not been able to identify a clear relationship between these results and the underlying wave conditions. However, it might be linked to the poor tolerance of the method to the occurrence of spurious data. In that sense, the wavelet method, being a purely database method, appears to be more resilient to the influence of outliers and spurious data. This makes it a potentially useful tool for analyzing signals with possibly noisy or unreliable data as in the case of moored wave buoy data.

A reasonable strategy to compare the main differences between the analyzed directional distributions is to examine the directional parameters such as wave direction and directional
<figure><img src="image_7.png" /><figcaption>FIG. 5. Directional distribution $D(f, \theta)$ at various nondimensional frequencies $f/f_p$ corresponding to 1000 UTC 20 Aug 2020 (same run as in Fig. 3). The vertical axis spans from 0 to 0.02. The light blue arrows mark the average wind direction reported by ERA5 reanalysis.</figcaption></figure>

spreading. These parameters provide a summary of the directional distribution characteristics and can be used to identify trends or patterns in the data that may not be immediately apparent from the raw directional distributions. The main results are discussed in the following subsections.

## b. Wave direction parameters

Wave direction can be considered as a measure of the central tendency of the directional distribution. Physically, it can be interpreted as the direction $\theta$ in which waves, of a particular frequency $f$, propagate. Two quantities can be defined considering the first and second pairs of Fourier coefficients:

$$ \theta_1(f) = \tan^{-1} \left\{ \frac{b_1(f)}{a_1(f)} \right\}, \quad \text{and} \quad (8) $$

$$ \theta_2(f) = \frac{1}{2} \tan^{-1} \left\{ \frac{b_2(f)}{a_2(f)} \right\}. \quad (9) $$

Following Ewans (1998) and Kuik et al. (1988), the angles $\theta_1(f)$ and $\theta_2(f)$ are referred to as the “mean wave direction” and the “dominant wave direction,” respectively. Note that when using the arc-tangent function to determine wave direction parameters, it is important to select the appropriate quadrant to avoid 180° ambiguity and ensure accurate results.

By definition, the Fourier coefficients $a_1(f)$, $a_2(f)$, $b_1(f)$, and $b_2(f)$ represent the circular moments of the directional

distribution $D(f, \theta)$—as defined by Eq. (3). Hence, any estimation of $D(f, \theta)$ allows obtaining a set of Fourier coefficients, as well as other parameters including wave direction and directional spreading. The Fourier coefficients can also be determined directly from wave buoy observations through cross-spectral analysis between the surface elevation and the eastward and northward components of the horizontal velocity or wave displacements (Kuik et al. 1988). These parameters are referred to as “observed” as they are derived directly from cross-spectral analysis of the buoy data. The circular moments obtained from the estimated directional distribution may not always align with those estimated from the cross-spectral analysis. These are considered as a baseline as they are closer to direct observations and are not influenced by the limitations of specific methods. Consequently, any difference could be attributed to potential errors in the specific method. In this section, we specifically evaluate the differences between the wave direction parameters estimated from the directional distribution for each method with those obtained from the observed Fourier coefficients.

Quantifying relative errors for circular data is not straightforward. Multiple metrics, depending on the specific context and applications, can be defined. We employed a measure of the distance between two angles defined by Mardia and Jupp (2000):

$$ d_{\theta} = \frac{1}{N} \sum_{i=1}^{N} [1 - \cos(\theta_{i}^{*} - \theta_{i})], \quad (10) $$
<table><caption>TABLE 1. Mean angular deviation for the estimation of the mean wave directions $\theta_1$ and the dominant wave direction $\theta_2$ as a function of the nondimensional frequency $ff_p$. The observed wave direction parameters, i.e., those based on the cross-spectral circular moments, are taken as reference and are compared to those derived from the directional distributions for each specific method.</caption><thead><tr><th></th><th colspan="11">ff<sub>p</sub></th></tr><tr><th></th><th>0.5</th><th>0.7</th><th>0.9</th><th>1.0</th><th>1.2</th><th>1.5</th><th>2.0</th><th>2.5</th><th>3.0</th><th>3.5</th><th>4.0</th></tr></thead><tbody><tr><th scope="row">θ<sub>1</sub></th><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><th scope="row">Wavelet</th><td>0.24</td><td>0.09</td><td>0.01</td><td>0.01</td><td>0.01</td><td>0.03</td><td>0.05</td><td>0.05</td><td>0.07</td><td>0.08</td><td>0.10</td></tr><tr><th scope="row">DFTM</th><td>0.37</td><td>0.09</td><td>0.01</td><td>0.01</td><td>0.01</td><td>0.02</td><td>0.05</td><td>0.05</td><td>0.06</td><td>0.08</td><td>0.11</td></tr><tr><th scope="row">IMLM</th><td>0.38</td><td>0.10</td><td>0.01</td><td>0.01</td><td>0.01</td><td>0.03</td><td>0.06</td><td>0.06</td><td>0.06</td><td>0.08</td><td>0.11</td></tr><tr><th scope="row">EMEP</th><td>0.37</td><td>0.09</td><td>0.01</td><td>0.01</td><td>0.01</td><td>0.04</td><td>0.06</td><td>0.06</td><td>0.07</td><td>0.09</td><td>0.12</td></tr><tr><th scope="row">θ<sub>2</sub></th><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><th scope="row">Wavelet</th><td>0.21</td><td>0.18</td><td>0.02</td><td>0.01</td><td>0.03</td><td>0.04</td><td>0.10</td><td>0.17</td><td>0.26</td><td>0.27</td><td>0.31</td></tr><tr><th scope="row">DFTM</th><td>0.35</td><td>0.25</td><td>0.02</td><td>0.01</td><td>0.03</td><td>0.04</td><td>0.13</td><td>0.14</td><td>0.24</td><td>0.28</td><td>0.30</td></tr><tr><th scope="row">IMLM</th><td>0.45</td><td>0.23</td><td>0.02</td><td>0.02</td><td>0.03</td><td>0.05</td><td>0.15</td><td>0.16</td><td>0.25</td><td>0.30</td><td>0.33</td></tr><tr><th scope="row">EMEP</th><td>0.28</td><td>0.19</td><td>0.02</td><td>0.02</td><td>0.05</td><td>0.07</td><td>0.15</td><td>0.19</td><td>0.29</td><td>0.37</td><td>0.43</td></tr></tbody></table>

where $N = 1056$ is the number of observations, $\theta$ is considered as the observed direction, and $\theta^*$ is the one estimated from the directional distribution for each specific method. This metric takes values between $d_\theta = 0$, when the two angles are equal, and $d_\theta = 2$, if one angle points in the opposite direction to the other. Similarly, $d_\theta = 1$ for two perpendicular angles.

Table 1 shows the mean angular differences, based on Eq. (10), in the observed and estimated wave direction parameters $\theta_1$ and $\theta_2$ at different nondimensional frequencies $ff_p$. The errors associated with $\theta_1$ are typically smaller in comparison with $\theta_2$. The larger differences are generally presented at both the lowest frequencies ($ff_p \approx 0.5$) and the highest frequencies ($ff_p > 3$). Near the spectral peak ($ff_p = 1$), the lowest differences are observed. A value of $d_\theta = 0.01$ would be equivalent to an angular difference between $6^\circ$ and $9^\circ$, which is the order of the angular resolution of the directional distribution ($\Delta\theta = 5^\circ$). The wavelet-based estimates are slightly closer to the observed directions in comparison with the other methods. However, the mean angular differences for all of the analyzed methods appear to be reasonably low, indicating that none of the methods would likely have difficulties estimating the wave direction parameters.

### c. Directional spreading parameters

In this section, the parameters that define the width of the directional distribution produced by the investigated methods are analyzed. We introduce the spreading factors based on parametric models and contrast our results with other studies. Finally, a discussion about the dependence of the spreading factors on the wave development stage is provided.

#### 1) PARAMETRIC MODELS

There are several analytical models that have been developed to represent the directional distribution function $D(f, \theta)$. These models are typically based on statistical analysis of wave data collected from buoys, wave gauges, or other wave measurement instruments. One of the first efforts was undertaken by Mitsuyasu et al. (1975), who collected directional wave information using a cloverleaf buoy at several locations in the Pacific Ocean off the coast of Japan. On the basis of

their data, they proposed an idealized shape of the directional distribution of the form:

$$D(f, \theta) = A(s) \cos^2s \left(\frac{\theta - \bar{\theta}}{2}\right), \quad (11)$$

where $A(s)$ is a normalization factor to satisfy Eq. (1), $s$ is a parameter that defines the width of the distribution, and $\bar{\theta}$ is the mean wave direction at each frequency $f$. The width of the distribution $s$ is well known to be a function of the nondimensional frequency $ff_p$. Additionally, Mitsuyasu et al. (1975) provided evidence that there is a certain dependence of the parameter $s$ on the dimensionless fetch and, thus, with the growth of the wind waves. Consequently, they suggested the following parameterization:

$$s = \begin{cases} s_p(f/f_p)^5, & f/f_p < 1 \\ s_p(f/f_p)^{-2.5}, & f/f_p \ge 1 \end{cases}$$

$$s_p = 11.5(U_{10}/c_p)^{-2.5}. \quad (12)$$

The ratio $U_{10}/c_p$ (where $U_{10}$ is the wind speed at 10 m above sea level and $c_p$ is the wave phase speed at the spectral peak) is known as the inverse wave age. It is a measure of the state of development of the wave field. As the ratio $U_{10}/c_p$ increases, waves are considered to be younger. Typically, long-period swells present values under 0.8, mature wind seas fall between 0.8 and 1, and young wind seas are characterized by higher values of $U_{10}$ relative to $c_p$. According to Mitsuyasu et al. (1975), developing young waves should present broader spreading than well-developed mature waves or swells.

During the so-called Joint North Sea Wave Project (JONS-WAP) experiment in the North Sea, Hasselmann et al. (1980) analyzed directional wave observations collected by a heave-pitch-roll buoy. The observed directional distributions by Hasselmann et al. (1980) fitted well to the $\cos^2s$ shape. Nevertheless, in contrast to Mitsuyasu et al. (1975), they observed no dependence on $U_{10}/c_p$ for frequencies below the spectral peak, producing the following parameterization:
<figure><img src="image_8.png" /><figcaption>FIG. 6. Spreading parameters (left) $\beta$ and (right) $s$, from the sech$^2\beta$ and $\cos^{2s}$ distributions, respectively. The circles represent the median of each parameter. The error bars indicate one standard deviation. In the right panel, the lower and upper bounds of both Mitsuyasu and Hasselmann models correspond to $U_{10}/c_p = 1$ and $U_{10}/c_p = 1.4$, respectively.</figcaption></figure>

$$ s = \begin{cases} 6.97(f/f_p)^{4.06}, & f/f_p < 1.05 \\ 9.77(f/f_p)^{\mu}, & f/f_p \ge 1.05, \end{cases} \\ \mu = -2.33 - 1.45(U_{10}/c_p - 1.17). \quad (13) $$

On the other hand, Donelan et al. (1985) investigated the directional characteristics of the wave field using data recorded with an array of 14 wave staffs in Lake Ontario. They concluded that the directional distribution given by Eq. (11) did not adequately describe their data. Considering the theoretical propagation of second-order Stokes wave groups, they proposed that the directional distribution can be represented by

$$ D(f, \theta) = \frac{\beta}{2} \operatorname{sech}^2 \beta(\theta - \bar{\theta}). \quad (14) $$

The parameter $\beta$ controls the width of the distribution and is a function of the nondimensional frequency $f/f_p$. In contrast to the previous studies, the directional distributions obtained by Donelan et al. (1985) did not exhibit any evidence of dependence on wave age. They only examined data up to a range of $f/f_p = 1.6$ and therefore concluded that above this value, the spreading remained constant. Banner (1990) refuted this conclusion using high-frequency stereophotography data and suggested that the spreading should broaden as the frequency increases. The final parameterization of the so-called Donelan-Banner distribution is

$$ \beta = \begin{cases} 2.61(f/f_p)^{1.3}, & 0.56 < f/f_p < 0.95 \\ 2.28(f/f_p)^{-1.3}, & 0.95 < f/f_p < 1.60 \\ 10^p, & f/f_p > 1.60 \end{cases} \\ p = -0.4 + 0.8393 \cdot e^{-0.567 \ln[(f/f_p)^2]} \quad (15) $$

Roughly speaking, both $\beta$ and $s$ are inversely related to wave spreading. The higher their values, the narrower the directional distribution. This predicts that the spreading tends to be narrower in the vicinity of the spectral peak and it broadens as the frequency moves away from the peak toward both lower and higher frequencies. Figure 6 shows the directional spreading parameters $\beta$ and $s$ for the directional distributions obtained by the investigated methods. The parameter $\beta$ was determined by using the nonlinear least squares method to find the best fit of the assessed directional distributions to the model described in Eq. (14). The parameter $s$ was calculated according to Eq. (5.9) in Ewans (1998). The data were grouped into several values of $f/f_p$, and the median and the standard deviation are presented.

The directional distribution obtained by the wavelet-based method was found to adequately fit both $\beta$ and $s$ spreading parameters, particularly for values greater than $f/f_p \approx 1.5$. These results are in good agreement with the observations and conclusions of Banner (1990). The EMEP method performed better close to the spectral peak, but had a higher
<figure><img src="image_9.png" /><figcaption>FIG. 7. (a) Wavelet-based estimates of circular rms spreading $\sigma_1$ (°) and (b) directional spreading factor $\sigma_2$. The circles represent the median for the entire dataset at each nondimensional frequency. The error bars indicate one standard deviation. In this case, the lower and upper bounds of the Mitsuyasu and Hasselmann parameterizations correspond to $U_{10}/c_p = 1.4$ and $U_{10}/c_p = 0.7$, respectively.</figcaption></figure>

standard deviation around the median across the entire range of frequencies, especially for the *s* parameter. The IMLM method consistently produced broader distributions in comparison with EMEP and wavelet. However, it presents less overall variability, particularly considering the $\cos^{2s}$ model. For frequencies below $f_p$, none of the methods follows the $\sec^2\beta$ parameterization. They were found to produce much broader distributions than the ones predicted by the Donelan et al. (1985) and Banner (1990) observations. Nevertheless, looking at the *s* parameter, all methods appear to adequately represent the directional spreading in the region $f/f_p < 1$, in accordance with the observations of Hasselmann et al. (1980). Finally, it is difficult to intercompare these two parameters because their values are specific to the models they represent. However, they can be translated into a universal spreading measure as shown in the next section.

## 2) CIRCULAR RMS SPREADING AND DIRECTIONAL SPREADING FACTOR

Multiple parameters can be defined to measure angular spreading. For example, Mardia and Jupp (2000) described various quantities for the variance, standard deviation, and angular dispersion, while Ewans (1998) proposed an angular width based on a symmetric double Gaussian directional distribution. In this study, we focused on two spreading parameters. First, the circular rms spreading $\sigma_1(f)$, which has been routinely used by operational hindcasts and observation programs (Ewans 1998; Forristall and Ewans 1998; Kuik et al. 1988). This parameter represents the average width of the directional distribution for each frequency. The second parameter is the directional spreading factor $\sigma_2(f)$. This parameter is defined as the square root of the ratio of the wave variance in the direction of the dominant wave direction to the total wave variance (Forristall and Ewans 1998; Latheef et al. 2017). The directional spreading factor varies from $\sqrt{1/2}$ for maximal

spreading, to 1 for unidirectional waves. Note that unlike $\sigma_1$, which increases with spreading, $\sigma_2$ decreases with spreading. We adopted the definitions provided by Ewans (1998) and Forristall and Ewans (1998) to enable direct comparisons with their respective studies. Similarly to the wave direction parameters, the spreading parameters can be obtained in terms of the first and second pairs of Fourier coefficients:

$$ \sigma_1(f) = \sqrt{2\left[1 - \sqrt{a_1^2(f) + b_1^2(f)}\right]}, \quad \text{and} \quad (16) $$

$$ \sigma_2(f) = \sqrt{\frac{1}{2}\left[1 + \sqrt{a_2^2(f) + b_2^2(f)}\right]} \quad (17) $$

Figure 7 compares the wavelet-based estimates of spreading parameters with the parametric models discussed earlier. This comparison aims to evaluate the consistency and agreement between the wavelet-based estimates and the parametric models in representing the spreading characteristics of ocean waves. The circular rms spreading $\sigma_1$ exhibits the lowest values in the vicinity of the spectral peak and increases for both lower and higher frequencies (Fig. 7a). The standard deviation typically varies between 6° and 10° and also seems to be the lowest at this frequency scale. The wavelet-based estimates of circular rms spreading at $f/f_p = 1$ are approximately 35°, whereas parametric models predict values between 15° and 33°, suggesting that our results have a tendency to produce broader spreadings near the peak frequency. This behavior is evidenced in Fig. 6 by all methods investigated, indicating that this might be either a limitation of the buoy measurements or a particular feature of this geographic region. The latter could be a solid explanation according to the findings of Forristall and Ewans (1998), who examined worldwide measurements of directional spreading and found that at sites exposed to extratropical storms, there is a tendency for
the waves to be more widely spread, particularly at higher lat-
itudes. At $f/f_p = 2$, we observe, on average, a spreading of
$\sigma_1 = 45^\circ$, in accordance with most of the investigated parame-
terizations and recent satellite observations (Le Merle et al.
2021). Conversely, the wavelet-based estimates reveal a slight
underestimation of the circular rms spreading for $f/f_p > 2.5$ in
comparison with the analyzed parameterizations. However,
the spreading value of $\sigma_1 = 55^\circ$ seems to be within the range
of $50^\circ-60^\circ$ reported by Ewans (1998) and Forristall and Ewans
(1998) at $f/f_p = 4$.

The circular rms spreading exhibits a particularly sharp in-
crease for frequencies below the spectral peak. Krogstad (2001)
attributed this behavior to the contribution of the second-order
bound nonlinear effects. Due to the proximity of coastlines, re-
flected and refracted energy could also contribute to the overall
spreading in the low-frequency band. The wavelet-based rms
spreading grows from $\sigma_1 = 35^\circ$ at $f/f_p = 0.9$ to $\sigma_1 = 70^\circ$ at
$f/f_p = 0.5$. This is clearly consistent with the predictions of the
Hasselmann model and to some extent of the Mitsuyasu model.
However, it becomes evident that there is an important deviation
between our results and the prediction of the Donelan–Banner
parameterization, which indicates a spreading of $\sigma_1 = 38^\circ$ near
$f/f_p = 0.5$. Our results are consistent with the values of $\sigma_1$ re-
ported by similar buoy-based observations (Mitsuyasu et al.
1975; Hasselmann et al. 1980; Rogers and Wang 2007; Ewans
1998; Forristall and Ewans 1998; Ewans 2001). However, looking
at observations made with other techniques, such as high-
frequency stereophotography (Banner 1990), spatial array of
wave staffs (Donelan et al. 1985), satellite-mounted real-aperture
radars (Le Merle et al. 2021), or airborne-based laser arrays
(Romero and Melville 2010), we observe that the reported val-
ues are of the same order of magnitude as those predicted by
the Donelan–Banner model, approximately half the magni-
tude of our estimates. Moreover, the estimated rms spreading
at this range of frequencies is similar among the investigated
methods. This would indicate that the discrepancies observed
at low frequencies ($f/f_p < 1$) are likely due to the measure-
ment technique rather than the methods themselves. Further-
more, Pettersson et al. (2003) compared directional spreading
from three different sensors and found significant disagree-
ments among them. They suggested that the mooring might
cause some effect in the directional spreading at lower fre-
quencies. Nevertheless, determining which measurement tech-
nique is more suitable to describe the directional spreading is
not immediately obvious. Although spatial measurements may
be more accurate to represent the wave field, they might be lim-
ited in the long waves. Likewise, the size of the image or the
distance between wave staffs will determine the maximum re-
solvable wavelength and thus the directional spreading accuracy
at this frequency range.

The directional spreading factor $\sigma_2$ reaches a maximum at the spectral peak and decreases toward both sides (Fig. 7b). This is consistent with the observations of Ewans (1998) and Forristall and Ewans (1998). However, they observed a minimum at $f/f_p = 2.8$ and then a rebound toward higher frequencies. This behavior is weakly reproduced by the Mitsuyasu and Hasselmann models at different $f/f_p$ depending on $U_{10}/c_p$ but is not shown by the Donelan–Banner model. This

rebound is not clearly observed for the wavelet-based esti-
mates of the directional spreading factor.

3) DEPENDENCE ON WAVE AGE

The relationship between the directional spreading and the
state of development of ocean waves has been a matter of de-
bate for several years. According to Mitsuyasu et al. (1975)
and Babanin and Soloviev (1998), wave age defines the width
of the directional distribution across the entire range of fre-
quencies. The Mitsuyasu model predicts values of the circular
rms spreading parameter ranging from $\sigma_1 = 15^\circ$ to $\sigma_1 = 33^\circ$
for $U_{10}/c_p = 0.7$ and $U_{10}/c_p = 1.4$, respectively. Hasselmann
et al. (1980) found that such dependency was only important
for $f/f_p > 1$. Conversely, Donelan et al. (1985) and Ewans
(1998) found no evidence of any wave age dependency.
Figure 8a shows the median values of $\sigma_1$ and $\sigma_2$ grouped
in three different ranges of the inverse wave age parameter
$U_{10}/c_p$ that correspond to swell (0–0.8), mature (0.8–1.0), and
young wind seas (1.0–2.0). Our results support the conclusions
of Hasselmann et al. (1980) in the sense that the directional
spreading clearly does not depend on wave age for $f/f_p < 1$.
The range of frequencies above the spectral peak shows a
weak yet inconclusive variation of the spreading for different
values of $U_{10}/c_p$. This finding would imply—as suggested by
Donelan et al. (1985), Ewans (1998), Romero et al. (2012),
Le Merle et al. (2021), and others—that the directional distri-
bution, and thus the shape of the directional wave spectrum,
is primarily determined by nonlinear wave–wave interactions
rather than by wind forcing. However, it is worth noting that
the influence of changing wind conditions on wave spreading
has not been addressed in this study. Variations in wind direc-
tion, gustiness, and accelerating wind conditions can have a
significant impact on wave spreading. This aspect remains to
be explored in future research.

4) DEPENDENCE ON SIGNIFICANT WAVE HEIGHT

Figure 9 suggests that near the spectral peak, the higher the
waves the narrower their directional distribution. The median
circular rms spreading $\sigma_1$ was about 10° larger for the smallest
wave category than for the largest. The directional spreading
factor $\sigma_2$ shows a maximum value around the spectral peak
and closely mirrors the inverted behavior of $\sigma_1$. The spreading
as a function of frequency for $H_s < 1$ m flattens at about
$f/f_p = 2$. For $f/f_p > 2.5$, the spreading in the smallest wave cat-
egory is lower than that observed for larger waves. This seems
to confirm the findings of Forristall and Ewans (1998), Ewans
(2001), and Latheef et al. (2017), who showed that with in-
creasing the significant wave height, the sea state becomes
less spread. However, these results are highly inconclusive as
the estimated standard deviations are around 5°–13°, signifi-
cantly high in comparison with the difference between the
median values reported for the different $H_s$ categories. The
flattening of spreading for low wave heights occurring when
$f/f_p > 2.5$ could be reasonably attributed to a low signal-to-
noise ratio. However, it is not completely clear how the
spreading would respond to a noisy signal. Guimarães et al.
<figure><img src="image_10.png" /><figcaption>FIG. 8. As in Fig. 7, but for various intervals of inverse wave age $U_{10}/c_p$.</figcaption></figure>

(2018) suggest that wave measurements from GPS data have more noise than those from inertial measurement units (IMUs), especially for low wave heights. However, this noise was observed to impact only low frequencies. By using Monte Carlo simulations in the context of ADCP measurements, Herbers and Lentz (2010) found that directional spreading values increase with decreasing wave height (i.e., decreasing signal-to-noise ratio). This has been also suggested by laboratory experiments (Lin et al. 2021) and wave buoy observations (Herbers et al. 2012). The overall relationship between directional spreading, significant wave height, and the role of signal-to-noise ratio remains uncertain.

## 4. Conclusions

In summary, we investigated the potential use of a wavelet-based method as an alternative approach to the conventional methods for estimating the directional distribution of ocean waves. We employed observations of a GPS buoy deployed off the west coast of Ireland. The wavelet-based estimates of the wave directional distribution are consistently good when compared to the conventional methods, namely, DFTM, ILM, and EMEP. The comparison reveals that DFTM is not suitable to describe the directional distribution, particularly for complex

multimodal wave systems; ILM produces broader directional spreading in comparison with the other methods; and EMEP performance is good overall; however, it appears to be highly unstable, in particular with the presence of spurious data. The wavelet-based method produces smoother directional distributions in comparison with the Fourier-based methods. This is consistent with the fundamental difference in frequency resolution between wavelets and Fourier analysis. The results of this study are in line with previous research and support the use of the wavelet-based method as a practical alternative for the estimation of the wave directional distribution.

The study also examined the overall behavior of the wavelet-based directional spreading parameters using the available dataset. A strong dependency on the nondimensional frequency was found. It was found that the spreading is narrower at the spectral peak and broadens asymmetrically at higher and lower scales, particularly sharply for frequencies below the peak. The analysis revealed that the spreading is weakly influenced by wave age across the entire range of frequencies, suggesting that nonlinear wave-wave interactions may play a primary role in controlling the angular width of the directional spectrum, rather than wind forcing.

The implications of these findings are significant for improving our methods of analysis of wave observations as well as our

<figure><img src="image_11.png" /><figcaption>FIG. 9. As in Fig. 7, but for various intervals of significant wave height $H_s$.</figcaption></figure>
understanding of the directional characteristics of the wave field, which play a crucial role in weather forecasting, oceanography, coastal engineering, and climate modeling. The wavelet-based method shows promise as a purely data-driven approach that does not require any assumptions about the shape of the distribution and can be used to complement conventional methods. In addition, it is relatively easy to implement in comparison with other methods, making it a practical option for researchers and engineers. The GPS buoy offers a practical and cost-effective means to obtain directional information of ocean waves, particularly in remote or hard-to-access locations. Moreover, improved estimation of the directional wave spectrum from these buoys can serve as valuable input for the calibration and validation of numerical wave models, enhancing the accuracy of wave forecasts and contributing to improved predictive capabilities.

In terms of future research directions, the findings of this study provide a basis for expanding the scope of the investigation to different kinds of buoys, such as heave-pitch-roll buoys. In addition, this method could also be potentially extended to acoustic Doppler current profiler data. Further research is also needed to perform extensive comparison between the different methods and techniques, in particular under controlled scenarios such as wave tanks to investigate their sensitivity to parameters such as the smoothing level.

*Acknowledgments.* This work was funded by the European Research Council (ERC) under the EU Horizon 2020 research and innovation program (Grant Agreement 833125-HIGHWAVE). The authors acknowledge Arnaud Disant for his assistance with the logistics of the buoy deployment. We extend our gratitude to the anonymous reviewers for their helpful comments that led to the improvement of the quality of this manuscript. DSPZ would like to thank the maintainers and developers of open-source scientific Python ecosystem for making data processing and visualization more efficient. DSPZ is grateful to the Navier (ENS Paris-Saclay) and Stokes (University College Dublin) teams for the fruitful discussions during the analysis of the data.

*Data availability statement.* The codes for computing the wavelet-based directional distribution can be accessed through the Python package Extended Wavelet Directional Method (EWDM) available at the GitHub repository (https://github.com/dspelaez/extended-wdm/). The Spotter buoy measurements used in this study were collected during the HIGH-WAVE project and are available upon reasonable request from the authors. The ResourceCODE data can be freely accessed and downloaded from https://resourcecode.ifremer.fr/. The ERA5 wind data reanalysis is available upon registration from https://cds.climate.copernicus.eu/.

# APPENDIX

## Kernel Density Estimation

Determining the most accurate representation of the density of a set of data has been an interesting problem in statistics for many years. An overview of different alternatives

can be found in Silverman (1998). The simpler approach would be to estimate the histogram along the time dimension. However, choosing a bin width that is too small will result in an unsteady or “jittery” directional distribution, and for bin widths that are too large, some information can be lost as the resulting directional resolution would be too low (Fig. A1). Even though histograms are an excellent tool for data description, in this work we explore KDE as an alternative method for density estimation, considering the “smooth” nature of the directional spreading function. KDE is a nonparametric method used to estimate the density of a dataset. It infers the overall density summing the contribution of the individual data points (Sheather 2004). For our particular application, the density of the local wave direction along the time dimension for a given frequency scale can be written as

$$ \widehat{D}_j(\theta) = \frac{1}{N} \sum_{i=1}^{N} K(\theta - \theta_i), \quad (A1) $$

where $j$ represents the frequency index, $N$ is the number of data points along the time dimension, and $K$ is the kernel function. Considering that the wave direction is measured in a circle, a kernel function which is wrapped between $[-\pi, \pi]$ is required. The von Mises distribution is a suitable choice since it has been widely used to model the probability distribution of data in a circle and particularly the distribution of wave directions (e.g., Hauser et al. 2005). The kernel function for a von Mises distribution is given by

$$ K(\theta) = \frac{1}{2\pi I_0(\kappa)} \exp\{\kappa \cos(\theta - \bar{\theta})\}, \quad (A2) $$

where $\bar{\theta}$ is the mean wave direction, $I_0(\kappa)$ is the zero order modified Bessel function, and $\kappa$ is known as the concentration parameter which is an inverse measure of the variance. The kernel function $K$ determines the shape of the individual contribution to the overall density, while $\kappa$ determines their bandwidth, and in turn, it controls the smoothness of the final directional spreading function.

The main problem of histograms persists: If the bandwidth is too small, the variance becomes large, resulting in an undersmooth spreading, whereas if the bandwidth is too large, the bias becomes large and the spreading will be oversmooth. The bandwidth selection in KDE has been a difficult problem for many years as it generally involves the prior knowledge of the curvature of the underlying theoretical distribution (Sheather 2004; Silverman 1998; Oliveira et al. 2012). In this case, in order to define a value of $\kappa$ that represents a good compromise between bias and variance, we take as a reference the value of $\Delta\theta = 5^\circ$, which is close to the uncertainty in the wave direction found by Wang and Freise (1997) when they analyzed the error of an NDBC 3-m pitch-roll discus buoy. In Fig. A1, the KDE estimation of directional spreading with $\kappa = 50$ yields similar results as the histogram estimation with $\Delta\theta = 5^\circ$, but with a smoother distribution. However, it is important to note that identifying the optimal value of $\kappa$ is beyond the scope of
<figure><img src="image_12.png" /><figcaption>FIG. A1. Comparison of the directional spreading obtained (top) by estimating the histogram for different directional resolutions and (bottom) by using KDE for different concentration parameters of the von Mises distribution.</figcaption></figure>

this study and should be explored further by comparing directional spreading observations from multiple sources.

## REFERENCES

Accensi, M., and Coauthors, 2021: ResourceCODE framework: A high-resolution wave parameter dataset for the European shelf and analysis toolbox. *EWTEC 2021—14th European Wave and Tidal Energy Conf.*, Plymouth, United Kingdom, University of Plymouth, 2182-1-2182-9, https://archimer.ifremer.fr/doc/00736/84812/.

—, M. F. Alday Gonzalez, and C. Maisondieu, 2022: RESOURCECODE – resource characterization to reduce the cost of energy through coordinated data enterprise. Database user manual. Tech. Rep. 86306, 46 pp., https://doi.org/10.13155/86306.

Babanin, A. V., and Y. P. Soloviev, 1998: Variability of directional spectra of wind-generated waves, studied by means of wave staff arrays. *Mar. Freshwater Res.*, **49**, 89–101, https://doi.org/10.1071/MF96126.

Banner, M. L., 1990: Equilibrium spectra of wind waves. *J. Phys. Oceanogr.*, **20**, 966–984, https://doi.org/10.1175/1520-0485(1990)020:0966:esoww:2.0:co2.

Benoit, M., 1992: Practical comparative performance survey of methods used for estimating directional wave spectra from heave-pitch-roll data. *Proc. 23rd Coastal Engineering Conf.*, Venice, Italy, American Society of Civil Engineers, 62–75, https://doi.org/10.1061/9780872629332.005.

—, 1993: Extensive comparison of directional wave analysis methods from gauge array data. *Proc. Second Int. Symp. on Ocean Wave Measurement and Analysis (ASCE)*, New Orleans, LA, American Society of Civil Engineers, 740–754, https://www.researchgate.net/publication/236846987_Extensive_comparison_of directional wave_analysis_methods_from_gauge_array_data.

—, P. Frigaard, and H. A. Schäffer, 1997: Analyzing multidirectional wave spectra: A tentative classification of available methods. *Proc. IAHR Seminar on Multidirectional Waves and Their Interaction with Structures*, 27th IAHR Congress, San Francisco, CA, National Research Council of Canada, 131–158, https://cir.nii.ac.jp/crid/1570572699752618112.

Capon, J., 1969: High-resolution frequency-wavenumber spectrum analysis. *Proc. IEEE*, **57**, 1408–1418, https://doi.org/10.1109/proc.1969.7278.

Cox, R., and Coauthors, 2020: Systematic review shows that work done by storm waves can be misinterpreted as tsunami-related because commonly used hydrodynamic equations are flawed. *Front. Mar. Sci.*, **7**, 4, https://doi.org/10.3389/fmars.2020.00004.

Donelan, M., A. Babanin, E. Sanina, and D. Chalikov, 2015: A comparison of methods for estimating directional spectra of surface waves. *J. Geophys. Res. Oceans*, **120**, 5040–5053, https://doi.org/10.1002/2015JC010808.

Donelan, M. A., J. Hamilton, and W. H. Hui, 1985: Directional spectra of wind-generated ocean waves. *Philos. Trans. Roy. Soc.*, **A315**, 509–562, https://doi.org/10.1098/rsta.1985.0054.
—, W. M. Drennan, and A. K. Magnusson, 1996: Nonstationary analysis of the directional properties of propagating waves. *J. Phys. Oceanogr.*, **26**, 1901–1914, https://doi.org/10.1175/1520-0485(1996)026%3C1901:NAOTDP%3E2.0.CO;2.

Draycott, S., A. C. Pillai, R. Gabl, P. K. Stansby, and T. Davey, 2022: An experimental assessment of the effect of current on wave buoy measurements. *Coastal Eng.*, **174**, 104114, https://doi.org/10.1016/j.coastaleng.2022.104114.

Ewans, K. C., 1998: Observations of the directional spectrum of fetch-limited waves. *J. Phys. Oceanogr.*, **28**, 495–512, https://doi.org/10.1175/1520-0485(1998)028<0495:OOTDSO>2.0.CO;2.

—, 2001: Directional spreading in ocean swell. *Ocean Wave Measurement and Analysis*, American Society of Civil Engineers, 517–529, https://doi.org/10.1061/40604(273)54.

—, 2015: A wavelet-based test for swell stationarity. *Appl. Ocean Res.*, **51**, 255–267, https://doi.org/10.1016/j.apor.2015.01.005.

Farge, M., 1992: Wavelet transforms and their applications to turbulence. *Annu. Rev. Fluid Mech.*, **24**, 395–458, https://doi.org/10.1146/annurev.fl.24.010192.002143.

Fedele, F., J. Herterich, A. Tayfun, and F. Dias, 2019: Large nearshore storm waves off the Irish coast. *Sci. Rep.*, **9**, 15406, https://doi.org/10.1038/s41598-019-51706-8.

Forristall, G. Z., and K. C. Ewans, 1998: Worldwide measurements of directional wave spreading. *J. Atmos. Oceanic Technol.*, **15**, 440–469, https://doi.org/10.1175/1520-0426(1998)015<0440:wmodws>2.0.co;2.

Grinsted, A., J. C. Moore, and S. Jevrejeva, 2004: Application of the cross wavelet transform and wavelet coherence to geophysical time series. *Nonlinear Processes Geophys.*, **11**, 561–566, https://doi.org/10.5194/npg-11-561-2004.

Gryazin, D., and K. Gleb, 2022: A new method to determine directional spectrum of sea waves and its application to wave buoys. *J. Ocean Eng. Mar. Energy*, **8**, 269–283, https://doi.org/10.1007/s40722-022-00228-z.

Guimarães, P. V., and Coauthors, 2018: A Surface Kinematics Buoy (SKIB) for wave-current interaction studies. *Ocean Sci.*, **14**, 1449–1460, https://doi.org/10.5194/os-14-1449-2018.

Hashimoto, N., and K. Kobune, 1988: Directional spectrum estimation from a Bayesian approach. *Proc. 21st Int. Conf. on Coastal Engineering*, Costa del Sol, Malaga, Spain, American Society of Civil Engineers, 62–76, https://doi.org/10.1061/9780872626874.005.

—, T. Nagai, and T. Asai, 1994: Extension of Maximum Entropy Principle for directional wave spectrum estimation. *Proc 24th Int Conf. on Coastal Engineering*, Kobe, Japan, American Society of Civil Engineers, 232–246, https://doi.org/10.1061/9780784400890.019.

Hasselmann, D. E., M. Dunckel, and J. A. Ewing, 1980: Directional wave spectra observed during JONSWAP 1973. *J. Phys. Oceanogr.*, **10**, 1264–1280, https://doi.org/10.1175/1520-0485(1980)010<1264:DWSODJ>2.0.CO;2.

Hauser, D., K. Kahma, H. Krogstad, S. Monbaliu, S. Lehner, and L. Wyatt, Eds., 2005: Measuring and analysing the directional spectra of ocean waves. Tech. Note COST Action 714, Working Group 3, 485 pp., https://repository.oceanbestpractices.org/handle/11329/1303.

Herbers, T. H. C., and S. J. Lentz, 2010: Observing directional properties of ocean swell with an acoustic Doppler current profiler (ADCP). *J. Atmos. Oceanic Technol.*, **27**, 210–225, https://doi.org/10.1175/2009JTECHO681.1.

—, P. F. Jessen, T. T. Janssen, D. B. Colbert, and J. H. MacMahan, 2012: Observing ocean surface waves with GPS-tracked buoys. *J. Atmos. Oceanic Technol.*, **29**, 944–959, https://doi.org/10.1175/JTECH-D-11-00128.1.

Holthuijsen, L. H., 2010: *Waves in Oceanic and Coastal Waters*. Cambridge University Press, 387 pp.

Johnson, D., 2002: Directional Wave Spectra Toolbox Version 1.3. MetOcean Solutions Ltd.

Kennedy, A., R. Cox, and F. Dias, 2021: Storm waves may be the source of some “tsunami” coastal boulder deposits. *Geophys. Res. Lett.*, **48**, e2020GL090775, https://doi.org/10.1029/2020GL090775.

Krogstad, H. E., 2001: Second order wave spectra and heave/slope wave measurements. *Ocean Wave Measurement and Analysis*, American Society of Civil Engineers, 288–296, https://doi.org/10.1061/40604(273)30.

—, S. F. Barstow, S. E. Aasen, and I. Rodriguez, 1999: Some recent developments in wave buoy measurement technology. *Coastal Eng.*, **37**, 309–329, https://doi.org/10.1016/S0378-3839(99)00031-9.

—, A. K. Magnusson, and M. Donelan, 2006: Wavelet and local directional analysis of ocean waves. *Int. J. Offshore Polar Eng.*, **16**, 97–103.

Kuik, A. J., G. P. van Vledder, and L. H. Holthuijsen, 1988: A method for the routine analysis of pitch-and-roll buoy wave data. *J. Phys. Oceanogr.*, **18**, 1020–1034, https://doi.org/10.1175/1520-0485(1988)018%3C1020:AMFTRA%3E2.0.CO;2.

Lancaster, O., R. Cossu, S. Boulay, S. Hunter, and T. E. Baldock, 2021: Comparative wave measurements at a wave energy site with a recently developed low-cost wave buoy (Spotter), ADCP and pressure loggers. *J. Atmos. Oceanic Technol.*, **38**, 1019–1033, https://doi.org/10.1175/JTECH-D-20-0168.1.

Latheef, M., C. Swan, and J. Spinneke, 2017: A laboratory study of nonlinear changes in the directionality of extreme seas. *Proc. Roy. Soc., A473*, 20160290, https://doi.org/10.1098/rspa.2016.0290.

Le Merle, E., D. Hauser, C. Peureux, L. Aouf, P. Schippers, C. Dufour, and A. Dalphinet, 2021: Directional and frequency spread of surface ocean waves from SWIM measurements. *J. Geophys. Res. Oceans*, **126**, e2021JC017220, https://doi.org/10.1029/2021JC017220.

Lenain, L., and W. K. Melville, 2017: Measurements of the directional spectrum across the equilibrium saturation ranges of wind-generated surface waves. *J. Phys. Oceanogr.*, **47**, 2123–2138, https://doi.org/10.1175/JPO-D-17-0017.1.

Liberzon, D., A. Vreme, S. Knobler, and I. Bentwich, 2019: Detection of breaking waves in single wave gauge records of surface elevation fluctuations. *J. Atmos. Oceanic Technol.*, **36**, 1863–1879, https://doi.org/10.1175/JTECH-D-19-0011.1.

Lin, Z., T. A. A. Adcock, and M. L. McAllister, 2021: Estimating ocean wave directional spreading using wave following buoys: A comparison of experimental buoy and gauge data. *J. Ocean Eng. Mar. Energy*, **8**, 83–97, https://doi.org/10.1007/s40722-021-00218-7.

Lygre, A., and H. E. Krogstad, 1986: Maximum entropy estimation of the directional distribution in ocean wave spectra. *J. Phys. Oceanogr.*, **16**, 2052–2060, https://doi.org/10.1175/1520-0485(1986)016%3C2052:MEEOTD%3E2.0.CO;2.

Mardia, K. V., and P. E. Jupp, 2000: *Directional Statistics*. Vol. 2, Wiley Online Library, 432 pp.

Massel, S. R., 2001: Wavelet analysis for processing of ocean surface wave records. *Ocean Eng.*, **28**, 957–987, https://doi.org/10.1016/S0029-8018(00)00044-5.
Mitsuyasu, H., F. Tasai, T. Suhara, S. Mizuno, M. Ohkusu, T. Honda, and K. Rikiishi, 1975: Observations of the directional spectrum of ocean waves using a cloverleaf buoy. *J. Phys. Oceanogr.*, **5**, 750–760, https://doi.org/10.1175/1520-0485(1975)005<0750:OOTDSO>2.0.CO;2.

Oliveira, M., R. M. Crujeiras, and A. Rodriguez-Casal, 2012: A plug-in rule for bandwidth selection in circular density estimation. *Comput. Stat. Data Anal.*, **56**, 3898–3908, https://doi.org/10.1016/j.csda.2012.05.021.

Oltman-Shay, J., and R. T. Guza, 1984: A data-adaptive ocean wave directional-spectrum estimator for pitch and roll type measurements. *J. Phys. Oceanogr.*, **14**, 1800–1810, https://doi.org/10.1175/1520-0485(1984)014<1800:ADAOWD>2.0.CO;2.

Pettersson, H., H. C. Graber, D. Hauser, C. Quentin, K. K. Kahma, W. M. Drennan, and M. A. Donelan, 2003: Directional wave measurements from three wave sensors during the FETCH experiment. *J. Geophys. Res.*, **108**, 8061, https://doi.org/10.1029/2001JC001164.

Plant, W. J., and M. A. Donelan, 2020: Directional surface wave spectra from point measurements of height and slope. *J. Atmos. Oceanic Technol.*, **37**, 67–83, https://doi.org/10.1175/JTECH-D-19-0128.1.

Raghukumar, K., G. Chang, F. Spada, C. Jones, T. Janssen, and A. Gans, 2019: Performance characteristics of “Spotter,” a newly developed real-time wave measurement buoy. *J. Atmos. Oceanic Technol.*, **36**, 1127–1141, https://doi.org/10.1175/JTECH-D-18-0151.1.

Raillard, N., S. Chabot, C. Maisondieu, D. Darbynian, G. Payne, and L. Papillon, 2023: RESOURCECODE: A python package for statistical analysis of sea-state hindcast data. *J. Open Source Software*, **8**, 4366, https://doi.org/10.21105/joss.04366.

Rogers, W. E., and D. W. C. Wang, 2007: Directional validation of wave predictions. *J. Atmos. Oceanic Technol.*, **24**, 504–520, https://doi.org/10.1175/JTECH1990.1.

Romero, L., and W. K. Melville, 2010: Airborne observations of fetch-limited waves in the Gulf of Tehuantepec. *J. Phys. Oceanogr.*, **40**, 441–465, https://doi.org/10.1175/2009JPO4127.1.

—, and K. Lubana, 2022: On the bimodality of the wind-wave spectrum: Mean square slopes and azimuthal overlap integral. *J. Phys. Oceanogr.*, **52**, 1549–1562, https://doi.org/10.1175/JPO-D-21-0299.1.

—, W. K. Melville, and J. M. Kleiss, 2012: Spectral energy dissipation due to surface wave breaking. *J. Phys. Oceanogr.*, **42**, 1421–1444, https://doi.org/10.1175/JPO-D-11-072.1.

Sheather, S. J., 2004: Density estimation. *Stat. Sci.*, **19**, 588–597, https://doi.org/10.1214/088342304000000297.

Silverman, B. W., 1998: *Density Estimation for Statistics and Data Analysis*. Routledge, 176 pp.

Simanesew, A. W., H. E. Krogstad, K. Trulsen, and J. C. N. Borge, 2018: Bimodality of directional distributions in ocean wave spectra: A comparison of data-adaptive estimation techniques. *J. Atmos. Oceanic Technol.*, **35**, 365–384, https://doi.org/10.1175/JTECH-D-17-0007.1.

Torrence, C., and G. P. Compo, 1998: A practical guide to wavelet analysis. *Bull. Amer. Meteor. Soc.*, **79**, 61–78, https://doi.org/10.1175/1520-0477(1998)079<0061:APGTWA>2.0.CO;2.

Wang, H. T., and C. B. Freise, 1997: Error analysis of the directional wave spectra obtained by the NDBC 3-m pitch-roll discus buoy. *IEEE J. Oceanic Eng.*, **22**, 639–648, https://doi.org/10.1109/48.650830.