## Abstract

The properties of directional distributions in ocean wave spectra are studied, with an emphasis on sea states with bimodal directional distributions in the high-frequency tails of single-peaked wave systems. A peak-splitting tendency has been a challenge in the interpretation of results from some data-adaptive estimation methods. After a survey of the theory, mathematical and numerical explanations are presented regarding domains of uni- and bimodality for symmetric Burg and Shannon maximum entropy methods. The study finds that both the Burg and Shannon maximum entropy methods have a tendency to split peaks, and that the domains of uni- and bimodality for these two methods depend on the Fourier coefficients input into the algorithms. Comparisons of data-adaptive methods based on data collected near the Ekofisk oil field in the North Sea and from nonlinear wave simulations are presented. The maximum likelihood (ML) method, the iterative maximum likelihood (IML) method, and the Burg and Shannon maximum entropy methods are applied. A large fraction of the directional wave spectra from Ekofisk shows bimodal features for distributions above the spectral peak for all of the abovementioned methods. In particular, strong similarity in bimodal features between the iterative maximum likelihood and the Burg maximum entropy methods are found. In general, the bimodality is consistent with previous observations, and it seems to be associated with wave and spectral development owing to nonlinear wave–wave interactions rather than being associated with the peak-splitting tendency in the estimates from any of the algorithms. The bimodal directional distributions were sometimes persistent and sometimes formed or decayed within the order of hours.

## 1. Introduction

Under second-order stationary and homogeneous conditions, an ocean wave field is stochastically characterized by the three-dimensional spectrum . The complexity involved in obtaining estimates of *χ* has turned the research toward derived quantities like the frequency spectrum and frequency-dependent directional (angular) distributions . However, from here to claim that is *the* directional wave spectrum is physically wrong and often leads to confusion because the operation loses information. Well-known studies such as Mitsuyasu et al. (1975), Hasselmann et al. (1980), and Donelan et al. (1985) observed and provide parameterizations for unimodal directional distributions with directional spread narrowest at the peak frequency and broader toward both lower and higher frequencies. Banner and Young (1994) reported numerical evidence that directional distributions may in fact be bimodal above the spectral peak of a single-peaked (in frequency) wave system. Bimodal distributions have been measured with spatial arrays (Young et al. 1995), buoys (Ewans 1998; Ewans and Van der Vlugt 1999; Wang and Hwang 2000), and HF radars (Kahma et al. 2005, p. 313) in fetch-limited sea states, and recently with a stereo video system (Leckler et al. 2015) for young wind waves. Similar observations with various other wave-measurement systems can be found in Wang and Hwang (2001), Hwang et al. (2000), Hwang and Wang (2001), Wang and Hwang (2003), and Long and Resio (2007).

Data from heave–pitch–roll measurement systems [heave–slope (HS)] will be central in the following analysis. Such systems provide data with few restrictive assumptions about the wave field, apart from homogeneity and stationarity. The standard output from the system is collocated time series of surface elevation and surface slopes in two orthogonal directions. Data from a spatially compact array may be interpolated into an HS setting and analyzed as triplet data.

It is well known, since the pioneering work of Longuet-Higgins (1976), that nonlinear energy transfer tends to move wave energy along directions ±35° away from a narrow peak in the wavenumber vector spectrum. The result has been verified using the full solution to nonlinear interaction equations (Banner and Young 1994) and nonlinear Schrödinger equations based on narrowband Gaussian initial spectra (Dysthe et al. 2003). In a recent study, Toffoli et al. (2010) have investigated the temporal evolution of directional wave spectra based on the potential Euler equations, applying a higher-order spectral method (Dommermuth and Yue 1987) that is free of bandwidth constraints and without any external forcing, such as wind input or breaking dissipation. For a fairly narrowbanded spectrum, it is observed that the energy of short waves redistributes in agreement with Longuet-Higgins (1976). However, for broadbanded waves the energy redistribution direction is shifted, to some extent, away from directions ±35° depending on the shapes of the initial directional distributions. Moreover, Morland (1996) has suggested that energy transfer from the wind to the waves as a result of inviscid critical layer mechanisms may give rise to double-peaked distributions for high frequencies.

Despite all these indications from the theory and observations of bimodal distributions in the short-wave region, uncertainties concerning peak-splitting tendency of some of the estimation methods make definite conclusions about the shape difficult. Estimation of ocean wave spectra mostly applies data-adaptive methods, and the resulting interpretations depend on analysis techniques and their underlying assumptions and principles. Comparisons of several directional analysis methods based on directional wave basin experiments obtained from three different measuring systems were reported by Benoit and Teisson (1994). They observed that the maximum likelihood (ML) method (see, e.g., Isobe et al. 1984; Krogstad 1988) provides broader directional distributions compared to the target, and that it performs poorly for sea states with narrow-peaked or bimodal distributions. On the other hand, the iterative maximum likelihood (IML) method (Oltman-Shay and Guza 1984) provides well-resolved directional distributions, but it is less reliable for estimates obtained from wave probe array when compared to Shannon maximum entropy (SME) method (Hashimoto and Kobune 1986) or Bayesian directional methods (Benoit and Teisson 1994). However, the internal iteration method used in the SME method to determine the Lagrangian multipliers experiences convergence problems if estimated Fourier coefficients are near the edge of the infeasible set (Krogstad 2012). An alternative estimation technique based on the Burg maximum entropy (BME) principle was introduced by Lygre and Krogstad (1986). When applied to certain theoretical distributions, the differences between the BME and SME algorithms range from the ability to match both broad and narrow peaks to a strong tendency for peak splitting (Krogstad 2012). With Fourier coefficients of certain theoretical unimodal distributions as input [e.g., the - distributions], the BME distribution shows two peaks that often have led to discussions about its physical significance. Our results show that this peak-splitting tendency is not limited to the BME method but that it also occurs, though somewhat weaker, in other estimated distributions.

Nevertheless, several recent studies indicate that the tendency to split spectral peaks may not be a problem for the BME method when applied to real wave data. Comparisons of HF radar spectra with the BME spectra estimated from Fourier coefficients determined from the HF radar spectra and with directional Waverider BME estimates are reported in Kahma et al. (2005, p. 313). The shapes of the spectra estimated from these three cases were consistent, including the bimodal features. During the study of frequency–wavenumber spectra of young wind waves, Leckler et al. (2015) compared a BME-processed spectrum from HS time series computed from 3D surface data with the original spectrum. The observed bimodal spectra in both BME and the original estimates are in reasonable agreement, confirming that the bimodality is real and consistent with Ewans (1998) observations. Moreover, since the spectra are measured by a stereo video system, any doubts concerning spurious peaks as a result of the analysis technique may be ruled out. Recently, Simanesew et al. (2016) reported estimates of ocean wave spectra from four different directional analysis methods. In many of the spectra reported, they found a striking agreement between the BME estimates and both the IML and SME estimates, suggesting that the bimodality in the BME distributions indeed is real.

Young (2010) examined both wavenumber and frequency spectra of finite-depth wind-generated waves using the wavelet directional method (Donelan et al. 1996). He observed bimodality in the spectrum, consistent with previous observations. Although the corresponding () spectrum was unimodal, he suggested that the possible bimodality in the spectrum was probably masked in the high-frequency region of the spectrum.

Our investigation, based on Ekofisk laser array data, and comprehensive with respect to sea states, shows that the BME estimates are generally consistent with the IML and SME estimates. In particular, for the bimodal sea states, we find strong similarity in bimodal features between the BME and IML distributions. These two methods simultaneously display the bimodal features in the whole range of frequencies, while the SME directional distributions show these bimodal features for relatively high frequencies.

## 2. Theory review

### a. Wave spectra

The most common mathematical model of a wavy ocean surface is a zero-mean weakly stationary (in time *t*) and homogeneous [in space ] stochastic field . This ensures the existence of a covariance function, , and, by the Wiener–Khinchin theorem, the Fourier identity,

Here, is the wavenumber vector and *ω* is the angular frequency. The function is the 3D wavenumber–frequency spectrum. An alternative derivation of *χ*, based on the stochastic integral representation for , is also sometimes used (Kahma et al. 2005, section 2.2). In general, *χ* will be a nonnegative generalized function.

The leading-order solution to the equations and boundary conditions for free surface waves leads to the dispersion relation connecting and possible *ω* values. We shall be concerned with gravity waves on finite depth for which , where

with . Here, *g* is the acceleration of gravity and *h* is the depth. We refer to this as linear wave theory (LWT) and the unique positive solution of Eq. (2) for *k* we denote as .

The general three-dimensional spectrum may also be written as

where is the frequency spectrum and is the wavenumber vector distribution at frequency *ω*, normalized so that for all .

For LWT, the has nonzero contributions only on the circle and may be written as

(Glad and Krogstad 1992). Here, defines the distribution of located on the intersection of the dispersion manifold and the plane *ω* equal to a constant, and *θ* is the direction of . Apart from large-area remote sensing systems, direct field measurements of are out of reach. However, estimates of certain moments are possible from the standard HS systems. The moments are integral expressions defined by

The HS systems provide the following properties: (Glad and Krogstad 1992). Nevertheless, the leading-order moments are far from determining the actual shape of . An example of slices of wavenumber distribution, , is shown in Fig. 1 for a set of positive frequencies. The figure shows numerical data from a higher-order modified nonlinear Schrödinger equation. The LWT dispersion circle and higher-order contributions are observed near the most energetic parts of wave harmonics. Above the spectral peak, the first-order harmonic shows a deviation from the dispersion relation and it is situated within the dispersion circle. This deviation is due to nonlinear evolution, as explained by Krogstad and Trulsen (2010) and Taklo et al. (2015).

Working with the full 3D spectrum is awkward, and it is necessary to introduce a set of simplified functions by integrating over and *ω*. The simplest is of course the well-known frequency spectrum introduced above,

The magnitude is symmetric with respect to the origin, , and has an integral equal to which explains the factor 2 in Eq. (6) and in the following.

We let be the integral over positive frequencies keeping fixed,

where . Term will in general have contributions from several distributions for each , and there is no way to distinguish between linear free- and bound-wave contributions without further information about the functions. Inversion from back to *χ* is thus only possible assuming the dispersion relation,

The transformation from to the commonly used directional spectrum, , is carried out by a straightforward substitution of the dispersion relation and a change in variables ,

### b. Directional distributions

The directional distribution may be viewed as a probability distribution function defined over the direction *θ*. Neglecting the *ω* dependence for simplicity, is a real, nonnegative, and periodic function with Fourier series,

where for , , and . We refer to Kahma et al. (2005) regarding properties of *D*. Considering an HS system, the cross-spectrum from , , and provides the moments from which we determine the leading Fourier coefficients,

(Glad and Krogstad 1992). Here, .

The main directional parameters—that is, the mean wave direction and the directional spread (circular standard deviation of *D*; Mardia 1972) —are usually expressed in terms of the Fourier coefficients,

where . The angular skewness and kurtosis of the distribution may also be used to describe the shape of the directional distributions (Mardia 1972; Fisher 1995; Kuik et al. 1988; Berens 2009). However, contrary to the linear case, the circular skewness and kurtosis do not bring in new information and they will not be considered further in this study.

We shall restrict ourselves to the leading complex Fourier coefficients, , which may be estimated directly by the HS systems. To ensure the existence of a nonnegative distribution with an integral equal to 1, it is necessary and sufficient that the Fourier coefficients form a positive semidefinite sequence for which the following matrix inequalities are valid:

(Kahma et al. 2005). These inequalities are easily seen to impose the following feasibility conditions on and :

for the first two Fourier coefficients. The second inequality may alternatively be written as

where and . See Kahma et al. (2005, section 2.5) for a complete algorithm.

#### 1) Some properties of ME distributions

We shall in the following consider two families of maximum entropy (ME) directional distributions obtained from estimates of feasible Fourier coefficients. The solutions of both families are formulated as a strictly concave functional for *D* with linear constraints. The BME was first used by Burg (1975) in connection with signal processing and later by Lygre and Krogstad (1986) for the estimation of directional distributions in ocean wave spectra. The solution, when it exists, is unique. The idea is to maximize the entropy

under linear constraints,

which leads to a directional distribution of the form

where

(Kay 1988). By introducing the representation from Eq. (16) we obtain the following simplified expressions: , , and .

The SME method (Hashimoto and Kobune 1986; Hashimoto 1997) is based on the entropy

which, apart from the periodicity, is identical to the entropy known from information theory. Maximizing the entropy in Eq. (21) under the constraints in Eq. (18), one obtains an expression for *D* containing the Lagrangian multipliers ,

where . It is obvious from Eq. (22) that is a strictly positive function for all finite choices of . A numerical algorithm, used to determine the multipliers satisfying the constraints in Eq. (18), is described in Krogstad (2012).

#### 2) Unimodal and bimodal regions for symmetric ME distributions

In the following we shall limit ourselves further to the important class of *symmetric* distributions with respect to some direction. To ensure that all Fourier coefficients are real numbers, we change the orientation of *θ* such that the symmetry is relative to .

##### (i) Burg ME distributions

Referring to Eq. (19), the BME distribution becomes the *Poisson distribution* (or *wrapped Cauchy distribution*) for . The opposite limit is the double-peaked delta distribution obtained when is set to 0. All nonsingular symmetric BME distributions are strictly positive. Thus, the stationary points of *D* are equivalent to the stationary points of the denominator in Eq. (19). By setting , one finds . The last pair of solutions requires that

Equation (23) is identical to Eq. (2.129) in Kahma et al. (2005). From Eq. (20) we have and . With only the first two solutions valid and a dominant peak of *D* at , the minimum of *D* has to be located at 180°, and the distribution is unimodal (see Fig. 2a). The last pair of solutions is symmetric with respect to and therefore both are maxima or both are minima. If both are maxima, then the stationary points 0 and 180° become local minima, and the new solution has two maxima and hence is bimodal (Fig. 2b). Finally, if both solutions in the last pair are minima, then the distribution has local peaks at 0 and 180° and the distribution is again bimodal (Fig. 2c).

Figure 3 gives an overview of and for symmetric Burg ME distributions, and the shape of typical distributions in the different regions. This figure, which allows or rather to become negative, is an expanded version of Fig. 2.5 in Kahma et al. (2005). To keep the resemblance with that graph, the ordinate axis shows for and for . The infeasible part is bounded by the coordinate axes and the black curve. The rest of the domain is divided into three parts: one for unimodal and two for bimodal distributions. Typical distributions are inserted for each region. However, the analysis has limited value when it comes to discriminating between uni- and bimodal distributions in general, as demonstrated in Fig. 4.

The green curve is the relationship for the unimodal Poisson distribution, whereas the magenta curve is the - distribution (Longuet-Higgins et al. 1963) for which the Burg ME distribution always has two peaks. The cyan curve is the relationship for the unimodal von Mises distribution, which is a special form of the SME functions. It leaves the Burg ME unimodal region slightly below .

With the von Mises Fourier coefficients as input, the Burg ME distribution shows two peaks for distributions with . A popular example is shown in Fig. 5 for the Poisson and - distributions. The Burg ME algorithm maintains the shape of the unimodal Poisson distribution, while it produces two false peaks for the unimodal - distribution.

##### (ii) Shannon ME distributions

The Shannon ME distributions are symmetric if and only if and are equal to zero, from which it immediately follows that the Fourier coefficients and are equal to 0. Since is only a normalization, one may consider

and define

We observe that *D* is bimodal for with peaks at and when is negative, and at 0 and *π* when is positive. The widths of the peaks decrease as moves away from 0. The signs of and may vary but in general . Assuming , if , then *D* will be a von Mises distribution centered around *π* when and around 0 when . Rewriting Eq. (24) as

one may compute the stationary points of from

which implies that

Hence, and *π* are always stationary points. Moreover, additional stationary points exist for

that is, . So far this mimics the BME case, with replacing the Fourier coefficients. It is likely that there is a mapping between and ; assuming that this is the case, the unimodal domain in turns out to be somewhat different from the previous. Note that here, instead of a distribution with Fourier coefficients, we operate with Lagrangian multipliers; see Fig. 6.

The von Mises distribution is completely within the Shannon ME unimodal region. It is interesting to observe, however, that both the Poisson and the - distributions are not completely within the unimodal region. The Poisson distribution crosses the Shannon ME upper unimodal boundary at , while the - distribution is expected to show two peaks below . This means that with the Poisson Fourier coefficients as input, the Shannon ME distribution shows two peaks for input distributions with : one main peak at the mean direction and one small sidelobe at about 180° offset from the mean direction. With the - Fourier coefficients as input, the Shannon ME distribution shows two but less distinct peaks on either sides of the mean direction for distributions with .

### c. Estimation of directional distributions

There exists a multitude of methods for estimating directional distributions, which in turn produce pages of intercomparison plots (Donelan et al. 2015). These algorithms take time series input and follow it up with multivariate spectral analysis. The spectral terms often relate to Fourier coefficients or other parameters from . However, early attempts using truncated Fourier series were not very successful (Longuet-Higgins et al. 1963). It was not before the 1970s that direction of arrival algorithms caused a significant step forward (Capon 1969). From the mid-1980s the maximum likelihood and maximum entropy formulas have turned out to be a good general choice for HS data. There are also several ways of bringing in Bayesian and inverse problem formalisms. A classic work is Long and Hasselmann (1979).

To obtain estimates of directional distributions from in situ measurements, we shall apply the following four directional analysis methods: ML, IML, and the two families of ME distributions. Both families are constructed based on HS-estimated Fourier coefficients (Glad and Krogstad 1992; Ochi 2005). Unlike the ML/IML algorithms, the HS systems additionally provide an estimate for the RMS wavenumber , and hence the ratio between and the LWT wavenumber, often called the check ratio:

The ML method was introduced by Capon (1969) in signal processing and was later modified in the context of ocean wave spectra (see, e.g., Davis and Regier 1977; Isobe et al. 1984; Krogstad 1988; Krogstad et al. 1988). This method has a known deficiency in that it tends to fail to reproduce the cross-spectra when is used to recompute the original spectra. As a consequence, it results in a smeared directional distribution. An appealing improvement may be achieved using the IML method introduced by Pawka (1983), which was later linked to standard inverse problem algorithms (Krogstad 1988; Ochi 2005). The method establishes an iterative improvement to the conventional ML estimate in order to deconvolve the smearing behavior of the ML method. We refer to Simanesew et al. (2016) and the references therein for further discussion about the ML and IML methods.

## 3. Ekofisk laser array measurements

The Ekofisk laser array is a research system mounted at the Ekofisk oil field in the North Sea since February 2003. It consists of four down-looking Optech lasers mounted on a bridge connecting the Kilo and Bravo platforms situated about 1 n mi (1 n mi = 1.852 km) northwest of the main Ekofisk complex. The lasers are placed at the four corners of a 2.6 m × 2.6 m square structure located approximately 20 m above the mean surface of a 70-m-deep sea. Raw data are collected continuously at 5 Hz with 1-mm accuracy. The data stream is stored in 20-min files with 1-min overlap.

The system was designed by the Norwegian Meteorological Institute [Meteorologisk Institutt (MET Norway)] in cooperation with the University of Miami, and the data collection is carried out under the operational responsibility of ConocoPhillips Inc. Extensive testing of the array has been performed, some of which is published (Krogstad et al. 2006; Krogstad and Trulsen 2010) and some documented in internal reports (e.g., Machado and Krogstad 2004) that are available on request from ConocoPhillips Inc.

Because of operational human health requirements on laser light power, the system has been forced to run with insufficient power, resulting in frequent dropouts in unfavorable conditions with a weak return signal. The raw data have been through a careful data check and restoration. Moreover, each record contains a summary of the data restoration and is downsampled to 1.7 Hz. The present analysis uses a fraction of the data from several tens of thousands of records, focusing on a set of 20-min records with and a unidirectivity index (UI) greater than 0.98. The size of the UI is used here to select records with one dominant wave field. The index UI is defined as

Figure 7 shows scatterplots of and for the abovementioned records, where curves of constant steepness are also plotted using the definition . Here denotes the wave period corresponding to the frequency at the spectral peak. It appears that the overall steepness lies below 0.1, which is a fairly extreme value.

Even if the full array consists of four lasers, the results are virtually unchanged when using only three lasers because of the array’s compact size relative to the size of the waves. This is actually convenient, since we are going to simulate an HS system by interpolating elevation and slope in the center of the triangle formed by the three remaining sensors and their elevation recordings. Both the reduced spatial array consisting of three lasers and the interpolated HS system have their usual limitations. The former suffers from spatial aliasing, since the shortest leg in the array is 2.6m, resulting in a limiting wavelength . The corresponding frequency , which is well below the temporal Nyquist limit

## 4. Data analysis

The ML and IML analyses are performed on spatial-array data, whereas the Burg and Shannon ME methods are applied to interpolated HS triplet data. As far as the Ekofisk compact array is concerned, it does not matter whether the ML/IML analysis is performed on spatial-array data or on the corresponding HS triplet data because the results are virtually identical. Both the ML and IML methods provide estimates of directional distributions from which we may determine the four leading Fourier coefficients. However, the HS system provides these Fourier coefficients directly and independent of any additional assumptions beyond the existence of a 3D spectrum. Apart from applying the same data series, the ML/IML and HS algorithms have no direct connections. Figure 8 shows plots of ML/IML directional spread against the HS triplet directional spread at four different frequencies for 1000 records. It is observed that the directional spread at the peak varies within the range 20°–65°. The IML spread is in fairly good agreement with the HS spread, while the ML spread shows a significant positive bias relative to HS.

Figure 9 shows color plots of the directional spectrum and the directional distribution of a typical unimodal sea state based on various analysis methods. For the directional spectra, shown in the upper row, there is only a minor difference among the methods, all showing clear unimodality, which is also reflected in the directional distributions shown in the lower row. Multimodality is a persistent feature in the Ekofisk data, where a large fraction of the directional spectra have evidence of bimodality for distributions above the spectral peak. Figure 10 shows occurrences of bimodal distributions for a set of records with . We observe bimodality for high frequencies in more than 80% of the spectra analyzed from nearly 1000 records. The angular separation of the two peaks mostly starts in the frequency range: . The spectrum in Fig. 11, which in this case is bimodal, is compared by applying the various analysis methods. The ML analysis of the data gives a rather smeared directional distribution, whereas the bimodal features are clear with ME and the bimodality is less pronounced with ML. On the other hand, the estimates from IML and both families of ME distributions show fairly good correspondence. For the IML and BME directional distributions, the bimodality shows up immediately above the spectral peak, while for the SME distributions, it shows up for moderately high frequencies. In numerous records, the angular separation of the two peaks for SME distributions occurs at relatively higher frequencies compared to both IML and BME distributions. The Burg ME algorithm provides relatively narrow directional distributions with pronounced peaks and a wide bottom as expected.

In the following section, we analyze selected data from storms that occurred during the period of 18–23 January 2007. During this period the local wind condition was unstable with the wind speed fluctuating between 1 and 25 m s^{−1}, and the significant wave height reaching up to 9 m. Figure 12 shows the development of the bimodal directional distribution from initially a unimodal sea state, observed in the evening of 22 January 2007. Estimates of directional spectra from BME are shown in the upper row with estimates of the directional distributions from various methods in the lower three rows. The IML estimates are in the upper-middle row, the BME estimates are in the lower-middle row, and the SME estimates are in the last row. Our results are qualitatively similar to Wang and Hwang (2001). At 1800 UTC, the wave energy propagates along the mean wave direction; however, as the wave field develops, the short waves are aligned along two directions oblique to the mean wave direction. Such a phenomenon in the ocean wave spectra was already observed by Wang and Hwang (2001) and attributed to nonlinear wave–wave interactions (Banner and Young 1994; Ewans 1998; Longuet-Higgins 1976; Dysthe et al. 2003; Toffoli et al. 2010). From the estimates of directional spectra in the upper row it is evident that, in this case, the bimodality does not seem to be associated with either the presence of two distinct wave systems or cases of new wave development alongside old ones. A 10% increase in steepness is observed between the spectra at 1800 and 2200 UTC.

In many of the spectra from the Ekofisk oil field, the ML algorithm often provides broader peaks, while the BME provides a narrow distribution compensated by a wide bottom, in agreement with previous observations. The SME distributions, while being nonsmeared, look more like the ML distributions around the spectral peak. However, when there is bimodality, the SME distributions shows more detailed bimodal features for moderate to high frequencies. On the other hand, IML distributions are quite similar to the BME distributions and lie somewhere in between the Shannon and Burg ME distributions.

To further confirm that the IML and BME methods provide virtually identical directional features, we consider nearly fully developed bimodal distributions from a number of neighboring records. The scatterplots in Fig. 13 show locations of peaks for bimodal directional distributions from 31 records. The peaks are computed as the local maxima of smoothed distributions. When there are several peaks, which is often the case for ML and IML distributions, we simply choose the two dominant ones. The overall shape of the distributions is in agreement, with the main difference being the existence of bimodality in the low-frequency range for the ME estimates and not for the ML and IML estimates. We have no explanation for this discrepancy; however, numerical evidence suggests that the distributions widen below the spectral peak as a result of nonlinear interactions (Simanesew et al. 2016). Near the spectral peak, the two bimodal arms coincide with the two dashed lines obtained by linear approximation of a parametric curve with respect to *ω* and *θ*, namely, the quartet resonance of the “figure of 8” for gravity waves derived by Phillips (1960). The relevance of these dashed lines for the directional redistribution of energy is captured by the two-dimensional cubic Schrödinger equation as discussed by Dysthe et al. (2003) in the space and is also discussed by Longuet-Higgins (1976). Here *θ* is the angle between two of the interacting waves.

Figure 14 shows scatterplots of angular separations between the two peaks for frequencies above the spectral peak, in the range . For ML, IML, and BME distributions, peak separation begins slightly above the spectral peak, whereas for the SME distributions, this happens at about twice the spectral peak. We see that the separation of the two peaks increases with frequency reaching up to 140° at about 4 times the spectral peak. Previous observations show a clear variation of the angular separation as a function of dimensionless frequency, based on data (Ewans 1998) and exact nonlinear simulation (Young et al. 1995). As shown in the figure, Ewans’s symmetric double Gaussian parameterization (Ewans 1998) clearly predicts the variation with dimensionless frequency; however, it appears to be a lower limit in this case. In Fig. 15, the ML, IML, and SME angular separations are compared with results from the BME. The ML is biased low and the SME is biased high, while the IML shows remarkable agreement with the BME angular separation.

Examples of the relationship between and are shown in Fig. 16 for symmetric Burg/Shannon ME distributions and for the Ekofisk data based on estimates of Fourier coefficients from HS systems. The decomposition of the scatter diagram into three subsets— similar to Fig. 3, which consists of one unimodal and two bimodal regions—is valid only for symmetric BME and SME distributions. One cannot prove in general that the distribution is uni- or bimodal by looking at the locations of in the diagram. They are indicators only in the sense that one may anticipate the shape of the distribution if certain Fourier coefficients are applied to either the BME or SME method. We have used data from four situations where and the main wave direction are reasonably steady, and data with below 3 m have also been removed.

Moreover, a stringent criterion has been imposed on the data based on the directions of and . One may rotate the coordinate system so that is real and is still complex. We shall redefine the Fourier coefficients as such that and compute , where and . The parameter CP has the following implications for symmetric distributions where and are either parallel or antiparallel: and are nearly parallel, and and are nearly antiparallel. The *x* axis now shows , and the *y* axis shows sign ()—positive for nearly parallel and negative for nearly antiparallel. For high values of (corresponding to small spread), the scatters are seen to cluster within the unimodal region and lie somewhere above the boxcar distribution (green curve). As the values of decrease, the cluster of points drifts into the bimodal type 2 region, densely populating within the boxcar distribution. It is also clear that part of the scatters for positive are outside the BME or SME distribution’s unimodal region, suggesting false peak splitting (see Fig. 17 for a better view).

In summary, the data analysis and the numerical experiments in section 2 suggest that there is peak-splitting tendency in both the Burg and Shannon maximum entropy algorithms for certain Fourier coefficients. However, many of the spectra, based on data from the Ekofisk array, show evidence of bimodality for distributions above the spectral peak. Since the bimodality is consistently revealed in all four estimation methods, it is less likely to be associated with spurious peaks generated by the analysis methods.

In addition, estimates of the directional spread obtained from IML and HS algorithms are virtually identical as shown in Fig. 8. For the IML and BME distributions, the bimodal features are also nearly identical. These two methods are dissimilar, yet they produce quite similar distributions. This similarity grants a certain confidence in the results. Thus, it may be recommendable to include both methods in the estimation procedure when additional information is needed.

## 5. Numerical simulations

We have employed the Dysthe equation, also known as the modified nonlinear Schrödinger (MNLS) equation, for the temporal evolution of a directional wave field (Dysthe 1979; Trulsen and Dysthe 1996; Trulsen et al. 2000). We aim to study the evolution of the wave spectrum, similar to the work of Dysthe et al. (2003), with a focus on the development of bimodal distributions. The MNLS equation accounts for nonlinear interactions that can lead to frequency-dependent directional distributions (Simanesew et al. 2016).

We employ the numerical method of Lo and Mei (1985, 1987) with periodic boundary conditions in both the *x* and *y* directions. Evolution in the temporal direction is achieved with a splitting scheme in which the linear part is integrated exactly in Fourier space and the nonlinear part is integrated by finite differences. A spatial grid with points is used for simulating the wave field. The computational domain in space is set to 50 characteristic wavelengths in both directions. For the directional analysis, the solutions are extracted at selected times.

All the simulations presented here are initialized by a two-dimensional Gaussian spectrum, as given in Eq. (3) of Dysthe et al. (2003), with given values of spectral bandwidths in the *x* and *y* directions denoted by and , respectively. The evolution is shown in Fig. 18 for and for . At , energy is concentrated around the spectral peak and propagates along the direction . As *t* increases, nonlinear effects subsequently move wave energy into oblique directions, thus developing a bimodal wave field above the spectral peak. Moreover, the energy redistribution seems to be significantly faster with a broader initial spectrum, .

The corresponding directional distributions have been extracted from the spectra above and are shown in Fig. 19. In addition to generating the bimodal directional distributions seen above the spectral peak, nonlinear interactions are also responsible for the widening of the directional distributions below the spectral peak and hence the increase in the directional spread toward the low frequencies. This is observed as a gradual widening of the directional distributions from minimum near the spectral peak to maximum both at lower- and upper-frequency tails. While the simulation directly produces the wavenumber spectrum , conversion to the directional spectrum requires an assumption of the linear dispersion relation. The bandwidth for application of the Schrödinger equation is limited to and, consequently, conclusions about the bimodality cannot be made for frequencies above .

The maximum entropy methods are also applied here to construct the directional spectrum from the three-dimensional MNLS simulation within the time interval from to . The estimates are based on the Fourier coefficients defined in section 2. The resulting spectra are shown in Fig. 20 along with the two principal directions for energy transfer. Both the BME and SME estimates show that the directional spectra extend mainly along these two directions as expected. The BME spectrum shows more detailed bimodal features, consistent with the shape of the spectra shown in Fig. 18.

## 6. Conclusions

In this study we present ocean surface wave directional analysis results from Ekofisk laser array measurements and from nonlinear numerical simulations. We employ four directional analysis methods for the estimation of directional spectra, namely, maximum likelihood (ML), iterative maximum likelihood (IML), Burg maximum entropy (BME), and Shannon maximum entropy (SME).

Previous studies suggest that BME directional estimates have a tendency to split peaks, which, to the best of our knowledge, has not been reported for the SME estimates. Our data analysis and numerical experiments, based on synthetic Fourier coefficients, suggest that there is a peak-splitting tendency in both the BME and SME estimates.

A large fraction of the directional wave spectra from the Ekofisk array shows evidence of bimodality at the higher frequencies for all four estimation methods. The consistency in the results between the various estimation methods suggests that the bimodality is real and that it is associated with the waves’ own development rather than being an artifact generated by the estimation methods. The angular separation of the two peaks starts mostly within the range of frequencies between the peak and 3 times the spectral peak and increases with frequency, reaching up to about 140° at about 4 times the spectral peak.

## Acknowledgments

We would like to acknowledge the constructive comments by the referees. The field data are from the laser array at the Ekofisk oil field in the North Sea and are used by permission of ConocoPhillips Inc. This research has been funded by the Research Council of Norway (RCN) and the University of Oslo through Projects RCN 225933, RCN 214556, and RCN 256466.

## REFERENCES

*Coastal Engineering 1994: Proceedings of the Twenty-Fourth International Conference*, B. L. Edge, Ed., Vol. 1, American Society of Civil Engineers, 42–56, https://doi.org/10.1061/9780784400890.004.

*Statistical Analysis of Circular Data.*Cambridge University Press, 277 pp., https://doi.org/10.1002/bimj.4710380307.

*Advances in Coastal and Ocean Engineering*, P. L.-F. Liu, Ed., Vol. 3, World Scientific, 103–143.

*Proceedings of the Fifth International Offshore Mechanics and Arctic Engineering Symposium*, J. S. Chung et al., Eds., Vol. 1, American Society of Mechanical Engineers, 80

**–**85.

*Proc. Symp. on Description and Modelling of Directional Seas*, Kongens Lyngby, Denmark, Danish Hydraulic Institute and Danish Maritime Institute, A-6, 15 pp.

*Modern Spectral Estimation: Theory and Application.*Pearson Education, 543 pp.

*Marine Technology and Engineering: CENTEC Anniversary Book*, C. G. Soares et al., Eds., Vol. 1, CRC Press, 109–124.

*Ocean Wave Spectra: Proceedings of a Conference*, Prentice-Hall, 111–136.

*Statistics of Directional Data.*Academic Press, 357 pp.

*Ocean Waves: The Stochastic Approach*. Cambridge Ocean Technology Series, Cambridge University Press, 332 pp.

*J. Geophys. Res.*,

**115**, C03006, https://doi.org/10.1029/2009JC005495.

*Coastal Engineering 2000*, B. L. Edge, Ed., American Society of Civil Engineers, 1254–1267, https://doi.org/10.1061/40549(276)97.

## Footnotes

© 2018 American Meteorological Society. For information regarding reuse of this content and general copyright information, consult the AMS Copyright Policy (www.ametsoc.org/PUBSReuseLicenses).