Role of Nonlinear Four-Wave Interactions Source Term on the Spectral Shape

: The goal of this paper is to investigate the importance of the four-wave nonlinear interactions (SNL4) on the shape of the power spectrum of ocean waves. To this end, the following results are discussed: a number of authors have conducted modern experimental measurements of ocean waves over the past decades and found that the measured power spectrum has (a) a very high central peak (characterized by the parameter γ , developed in the 1970s in the JONSWAP program) and (b) enhanced high-frequency channels which lead to the phenomenon of “bimodality”, also a well-known phenomenon. We discuss how a numerical hindcast of the Draupner storm (1995) with the standard code WAVEWATCH-III with full Boltzmann interactions also reﬂects these previously experimentally determined spectral shapes. Our results suggest that the use of the full Boltzmann interactions (as opposed to the discrete interaction approximation often employed for forecasting / hindcasting) is important for obtaining this characteristic physical spectral shape of the power spectrum.


Introduction
Recently, third-generation spectral wave models have been applied to investigate extreme waves in wind [1], waves under explosive cyclones [2] and on wave hindcasts performed with the WAVEWATCH-III model (WW3). Based on experimental wave tank results, several investigators [3,4] corroborate that the effects of four-wave nonlinear interactions (SNL4) are stronger in wave fields with narrower spectra. The goal of this paper is to make a relative comparison of hindcasts of the Draupner storm using WAVEWATCH-III with two nonlinear wave interaction models: (1) the discrete interaction approximation (DIA) model [5] and (2) the Webb-Resio-Tracy (WRT) model [6][7][8], which computes the full nonlinear Boltzmann interactions. Furthermore, because we seek to use a standard forecast/hindcast computer code, we have implemented WAVEWATCH-III [9,10], which contains both the DIA and WRT algorithms.
We recognize a well-known problem in real-time forecasting, brought about by the use of the two methods for computing the SNL4. DIA was developed to approximate the nonlinear transfer of wave energy between pairs of four wave numbers in spectral space, having the same number of degrees of freedom as the discrete spectrum. Consequently, DIA is relatively fast and useful for real-time forecasts. However, WRT computes the full Boltzmann integral and is therefore a much more detailed algorithm that consumes considerable computer time, and is therefore not amenable to real-time forecasting. Our goal here in hindcasting of the Draupner storm is to compare the results of the two hindcasts, not only as an exercise in hindcasting, but also as an exercise in wave physics. What information does the WRT spectrum have that the DIA spectrum does not? It has been a long dream, indeed a hope, that extreme waves could be accounted for in wind/wave models. To this end, the instability of a uniform narrow-banded wave train to sideband perturbations, or the so-called Benjamin-Feir or modulational instability (BFI, [11]), has been studied in recent years by many authors. A number of authors have associated the modulation instability with the occurrence of rogue waves in the ocean: this is because of the spontaneous formation of breather states in nonlinear random wave trains (see [12,13] for a full list of references). Breather states, born of the BFI, are nonlinear packets (nonlinear beats) that can rise to quite high amplitudes (three or four times the background sea state) during their evolution. The BFI parameter can be used to characterize the possible presence of rogue waves in a sea state [14]. At present, the modulational instability mechanism is one of the theories for understanding the generation of many types of freak waves.
However, we do not discuss in the present study what physical mechanisms might have generated the huge wave known as the Draupner wave. The main reason is that the third-generation wind/wave models we use are not able to address the extreme wave field, because they are phase-averaged models. The fact that phase locking cannot be addressed therefore excludes the presence of coherent structures (Stokes waves, breathers) in the wave field. Phase averaging destroys coherency and therefore, using present generation wind/wave models to study rogue waves is not a tenable line of research. Perhaps future wind/wave models will address the issue of nonlinear phase locking and the problem can then be revisited at that time. Some attempts have been made to correct this situation [15,16]; however, explicit representation of rogue waves in the spectrum of wind/waves remains a challenge.
Nevertheless, we have decided to reflect on what possible nonlinear properties of ocean waves that modern wind/wave models can produce, primarily from the point of view of the BFI. Because of this focus on nonlinear interactions, in the present study, we use the full Boltzmann interactions for SNL4, but we also address the possibility that the modulation instability might play a role in spectral properties for two-dimensional waves. We have focused on the idea that the shape of the power spectrum naturally comes from highly nonlinear processes. We ask further: Can the Benjamin-Feir instability have an influence on the shape of the power spectrum? While the results that we present are qualitative, they are intriguing and supported by previous work [17], suggesting that several features of power spectra determined from measured data reflect properties of the modulational instability. In two dimensions, one normally uses directional spreading, the cosine to the nth power, to modify the PM and JONSWAP spectra. This study suggests that the actual physical spectral shapes that occur in nature differ from these characteristic forms due to the BFI. It may well be that a generic spectral shape in three dimensions can be derived in future studies, beginning from considerations given herein.
The three-dimensional shape of the power spectrum of nonlinear sea states has not yet been treated in detail. Thanks to several careful investigations [18][19][20][21], there are excellent experimental results showing that the spectrum of surface waves is not a JONSWAP spectrum with cosine to the nth power spreading. Indeed, the spectrum seen experimentally by these and other authors is characterized by a high central peak with two enhanced channels at low frequency and another two channels at high frequency. The physical interpretation of spectral shapes of this type has remained open. This subject was partially addressed in [17], who obtained a theoretical shape of the power spectrum for a standard Pierson-Moskowitz spectrum, which is consistent with the measurements just discussed. We develop this idea further herein.
In the present paper, we shed light on this intriguing problem by also studying the Draupner storm [22]. We provide a hindcast of this storm using WW3 with both the discrete interaction approximation (DIA) and the Webb-Resio-Tracy (WRT) full nonlinear Boltzmann interactions, and find a spectral shape similar to those experimentally obtained and discussed in the previous paragraph. We find, however, that the WRT spectra have shapes more qualitatively consistent with the experimental results, as opposed to those obtained from DIA. This paper is structured as follows. Section 2 discusses the basis of wave modeling as used herein. Section 3 gives details of the hindcast of the Draupner storm, Section 3.1 gives an overview of the actual hindcast, Section 3.2 gives a parameter study of the Draupner storm and Section 3.3 gives an overview of the shape of the Draupner spectrum as computed from the hindcasts. Section 4 provides our perspective on the two-dimensional spectrum of ocean waves based upon the influence of the Benjamin-Feir instability. Section 5 provides conclusions, the most important of which is that one should use the WRT nonlinear interactions, rather than the DIA method, to ensure that the nonlinear physical shape of the computed power spectrum is correct, to within the phase-averaged approximation made in the model.

Theoretical Background of Wave Model
The evolution in time and space of the two-dimensional wave spectrum F( f , θ, x, y, t) as function of frequency f, direction θ, position x, y and time t is described by the energy balance equation of ocean waves (deep water case) as: ∂F ∂t where C g is the wave group velocity vector at which information travels across the spectral domain. On the right-hand side of the equation is the source function that consists of all the physical processes that allow the growth and the decay of the waves. These processes are the wind input from the atmosphere to the water surface (S in ), the dissipation of wave energy due to white capping (S dis ) and the nonlinear four-wave interactions implicit in the waves themselves (S nl , SNL4). Details about the wind and dissipation source terms (S in and S dis ) can be found in [23,24] The S nl source term (1), describing the nonlinear energy transfer among the waves due to resonant interactions, is given by the Boltzmann integral [25]: where n i is the action density at wavenumber k i , n i = n(k i ), for n(k) = P(k)/ω where P(k) is the wavenumber spectrum. G(k 1 , k 2 , k 3 , k 4 ) is the coupling coefficient [6] and δ κ and δ σ are delta functions representing the resonance conditions for wavenumber and frequency. Of all source terms, the one for the S nl it is the only source term for which a closed theoretical solution based on first principles exists, although it is it is expressed in terms of a complicated six-fold integral in a three-dimensional wave number manifold, making it very time-consuming for operational purposes.
In most of the operational wave forecasts, the peak of the spectrum is eroded when the DIA (discrete interaction approximation of Hasselmann et al., 1985) method is used, rather than the full Boltzmann integral. However, in addition to DIA, there are several other methods for computing the S nl , such as we mentioned earlier: the WRT, the generalized multiple DIA [26], the two-scale approximation (TSA) [27], the SRIAM (Simplified Research Institute of Applied Mechanics) [28] and the generalized kinetic equation of Gramstad and Babanin [16].
In [29], we conducted several numerical experiments to help better understand the physics in the kinetic equation. We make a number of observations with regard to the behavior of the dynamics of the spectrum that surprisingly behave similarly to certain aspects of the modulation instability. In the present work, we essentially discuss the WW3 results obtained by using DIA and the WRT methods for the case of the Draupner storm. The WRT method was implemented in WW3 using the portable subroutines developed by Van Vledder [30,31].

Hindcast of the Draupner Storm
The hindcasts of the Draupner storm of January 1995 were conducted by applying the WW3 model ( The shallow water model was here used for the North Sea nested grid. For the computation of the nonlinear interactions, the discrete interaction approximation (DIA) and the full Boltzmann integral (WRT) methods were activated in separate hindcast exercises. The WRT method is based on the original work of the six-dimensional Boltzmann integral formulation of [32,33], and additional considerations by [6][7][8].
For the wind input and the wave energy dissipation, the ST3 WW3 module was adopted [34]. The JONSWAP bottom friction formulation and refraction were activated. For the propagation of the spectral wave energy through the geographical domain, a third-order propagation scheme was chosen using Tolman's [35] averaging technique. Reflection, ice and currents were not considered for the present study.
The hindcast was performed for 49 days starting on 15 November 1994 at 00 UTC and ending on 3 of January 1995 at 00 UTC by activating the DIA method (hereafter WW3-DIA) for the creation of the initial conditions for the hindcasts (WW3-WRT) for reproducing the sea state conditions of the Draupner storm. The WW3-WRT hindcasts began on 31 December 1994 at 00 UTC and ended on 3 January 1995 at 00 UTC.
With the spectral wave model configuration described above, the WW3-WRT hindcast took almost 5 days of wall clock time, in contrast with the WW3-DIA calculation that took only 27.8 h. Computations were performed in the Fionn Supercomputer of the Irish Centre for High-End Computing (ICHEC).
The WW3 model was driven by the 1 hourly the NCEP (National Centers for Environmental Prediction) Climate Forecast System Reanalysis (CFSR) of NOAA (National Oceanic and Atmospheric Administration) [36].
According to the NOAA's reanalysis, the wind speed on 1 January 1995 at 15 UTC reached a maximum of 22 m/s blowing from the NW along the main axis of the North Sea, where the fetch is largest (Figure 1, left panel), whereas the maximum significant wave height reached almost 11 m (middle panel of Figure 1) in the Draupner platform, represented with a circle in accordance with the value reported in [22,37]. In addition, from the WW3-WRT hindcast, peak periods higher than 15 s where observed at the moment of the big wave distributed in the large part of the basin.
The sea state conditions computed here indicated no crossing seas at the moment of the Draupner wave event (and before this date).
Frequency wave spectra from the observed data and from the two different hindcasts were estimated and compared to each other ( Figure 2). Good agreement was obtained; however, a slight shift can be seen in the peaks. Qualitatively, the WRT-based simulation appears to provide a better match to the observations, especially in the context of the energy level at the peak. The reason for this shift can be associated with the spectral resolution employed in the settings of the performed hindcasts, which for the present exercise had a configuration which consists of 36 directions, 40 frequencies, an increment factor of 1.070, frequency range from 0.0350 up to 0.4898 Hz and a directional increment of 10 • .  Comparison of the observed power spectrum (red), using Haver's data from the Draupner platform, and the spectra from the WW3-discrete interaction approximation (DIA) (blue dotted) and WW3-WRT (black dotted) hindcasts at the Draupner location. The data spectrum and the spectrum from WW3-WRT are in good agreement. The WW3-DIA spectrum is substantially lower than the data and WW3-WRT spectra.
Our WW3 results show that using the WRT for computing the SNL4 values of the main averaged parameters are higher than those obtained using the DIA approximation ( Figure 3). As can be seen,   Comparison of the observed power spectrum (red), using Haver's data from the Draupner platform, and the spectra from the WW3-discrete interaction approximation (DIA) (blue dotted) and WW3-WRT (black dotted) hindcasts at the Draupner location. The data spectrum and the spectrum from WW3-WRT are in good agreement. The WW3-DIA spectrum is substantially lower than the data and WW3-WRT spectra.
Our WW3 results show that using the WRT for computing the SNL4 values of the main averaged parameters are higher than those obtained using the DIA approximation ( Figure 3). As can be seen, at the peak of the Draupner storm, the H s computed with WW3-DIA is H s = 9.94 m, whereas with WW3-WRT, the H s = 10.8 m [22] reported a H s = 11.9 m on 1 January 1995 at 15 UTC and 20 min; [38] reported H s = 10.5 m at the time of the freak wave, about 1 m lower than the measured value.
Our results indicate some improvement of the H s since the difference between the observed value and that computed with WW3-WRT is 1.1 m. Regarding the directional spreading, at the time of the Draupner wave (1 January 1995 at 15 UTC), we obtained 18.4° on the case of the WW3-WRT Figure 2. Comparison of the observed power spectrum (red), using Haver's data from the Draupner platform, and the spectra from the WW3-discrete interaction approximation (DIA) (blue dotted) and WW3-WRT (black dotted) hindcasts at the Draupner location. The data spectrum and the spectrum from WW3-WRT are in good agreement. The WW3-DIA spectrum is substantially lower than the data and WW3-WRT spectra.
Our WW3 results show that using the WRT for computing the SNL4 values of the main averaged parameters are higher than those obtained using the DIA approximation ( Figure 3). As can be seen, at the peak of the Draupner storm, the H s computed with WW3-DIA is H s = 9.94 m, whereas with WW3-WRT, the H s = 10.8 m [22] reported a H s = 11.9 m on 1 January 1995 at 15 UTC and 20 min; [38] reported H s = 10.5 m at the time of the freak wave, about 1 m lower than the measured value.
Our results indicate some improvement of the H s since the difference between the observed value and that computed with WW3-WRT is 1.1 m. Regarding the directional spreading, at the time of the Draupner wave (1 January 1995 at 15 UTC), we obtained 18.4 • on the case of the WW3-WRT simulation and 22.1 • with WW3-DIA, which shows a difference of about 4 • . In [39], the difference obtained in the directional spreading was not detectable when they compared the hindcast by the SRIAM nonlinear source terms and a hindcast using WW3-DIA.  [39], the difference obtained in the directional spreading was not detectable when they compared the hindcast by the SRIAM nonlinear source terms and a hindcast using WW3-DIA. Simultaneously, a noticeable decrease of the peak period (WW3-WRT) was found (up to 8.1699 s) on 31 December 1994 at 06 UTC, which was not observed in the WW3-DIA hindcast ( Figure 3). The T p decreases as H s increases and this fact can be explained with a pre-existing swell (with a long period), which begins to be dominated by "sea" as H s increases.
More recently, several studies reported the bimodal directional distribution of extreme events based on measurements of a stereo video system, Doppler radar and HF (High Frequency) radar [40][41][42][43]. In our simulations of the Draupner wave, we see the bimodality distribution of the spectrum ( Figure 4). We give the cross-section obtained perpendicular to the frequency axis (at f / f p = 2 ). From the latter, it can be seen that the DIA result over time (not all times shown here) gives a rather muddled shape of the spectrum as time increases. However, the WRT evolves systematically by robustly filling in the channels as time increases (Figure 4, compare middle and right panels). Note the bimodal shape Simultaneously, a noticeable decrease of the peak period (WW3-WRT) was found (up to 8.1699 s) on 31 December 1994 at 06 UTC, which was not observed in the WW3-DIA hindcast ( Figure 3). The T p decreases as H s increases and this fact can be explained with a pre-existing swell (with a long period), which begins to be dominated by "sea" as H s increases.
More recently, several studies reported the bimodal directional distribution of extreme events based on measurements of a stereo video system, Doppler radar and HF (High Frequency) radar [40][41][42][43]. In our simulations of the Draupner wave, we see the bimodality distribution of the spectrum (Figure 4).  [39], the difference obtained in the directional spreading was not detectable when they compared the hindcast by the SRIAM nonlinear source terms and a hindcast using WW3-DIA. Simultaneously, a noticeable decrease of the peak period (WW3-WRT) was found (up to 8.1699 s) on 31 December 1994 at 06 UTC, which was not observed in the WW3-DIA hindcast (Figure 3). The T p decreases as H s increases and this fact can be explained with a pre-existing swell (with a long period), which begins to be dominated by "sea" as H s increases.
More recently, several studies reported the bimodal directional distribution of extreme events based on measurements of a stereo video system, Doppler radar and HF (High Frequency) radar [40][41][42][43]. In our simulations of the Draupner wave, we see the bimodality distribution of the spectrum (Figure 4). We give the cross-section obtained perpendicular to the frequency axis (at f / f p = 2 ). From the latter, it can be seen that the DIA result over time (not all times shown here) gives a rather muddled shape of the spectrum as time increases. However, the WRT evolves systematically by robustly filling in the channels as time increases (Figure 4, compare middle and right panels). Note the bimodal shape We give the cross-section obtained perpendicular to the frequency axis (at f / f p = 2). From the latter, it can be seen that the DIA result over time (not all times shown here) gives a rather muddled shape of the spectrum as time increases. However, the WRT evolves systematically by robustly filling in the channels as time increases (Figure 4, compare middle and right panels). Note the bimodal shape for the case with the WRT algorithm. Similar results for a Pierson-Moskowitz spectrum were obtained in [29].

A Parameter Study of the Draupner Storm
Moreover, we estimated integral parameters that characterize the shape of the spectrum, based on [44], and compared their magnitudes for both hindcasts (WW3-DIA and WW3-WRT). From Table 1, it can be seen that spectral bandwidth in the frequency axis represented with the parameter of Goda (Q p ) is higher (4.1993) in the case of WW3-WRT. In addition, the spectral steepness was higher (0.0579), as well as the BFI (0.6097) and the Kurtosis (0.0222) for WW3-WRT. The spreading, as mentioned above, is lower for the case of the WW3-WRT (18.39 • ). This suggests that the full nonlinear interactions are important for the estimation of the integral extreme wave parameters. The deviations from Gaussianity are measured in terms of the kurtosis C 4 of the surface elevation probability density function [44]. A fit was found for the maximum of the kurtosis (based on results from numerical simulations of the Nonlinear Schrödinger equation), which follows from the narrow-band limit of the Zakharov equation (see [44]) where BFI is the Benjamin-Feir index, which is given by: where ε = k p √ m o is the integral steepness parameter, k p is the peak wavenumber and m o is the spectral zero moment. The BFI is the ratio between nonlinearity and spectral bandwidth. Its relation to the kurtosis has been found in the limit of large times and narrowband spectra neglecting directional spreading [15]. Including the contribution from the shape of the waves, the total kurtosis becomes: Q p above is the Goda peakedness parameter, The Mean Directional Spreading σ θ in Table 1 is: Additionally, we compute the nonlinear Benjamin-Feir instability index for wind waves (I BF ) (Equation (8)) [12,13] that estimates the number of rogue wave unstable modes in the nonlinear wave spectrum (it is the integer part of I BF ), based on the application of the theory of the nonlinear Schrödinger NLS equation to random wave trains. In this case, for the WRT hindcast, we obtained a higher I BF (77.6) than for the DIA simulation (71.4), considering a sea state duration of 3 h.
While the BFI of Janssen is a measure of how peaked the wave spectrum is, the IBF of Osborne is an estimate of the number of breathers trains in the spectrum, where breathers can be seen as pairwise phase-locked Stokes waves for sufficiently steep and narrow-banded wave motion, such that I BF > 1 [13]. From a parametric point of view, it is easy use Equation (8) to estimate I BF , but the reader is forewarned that this value has dubious value in modern wave theories that are not phase resolving.

The Shape of the Draupner Spectrum from the Hindcasts
Here, we emphasize an important feature of the WW3-WRT power spectrum: the peak is about 40 percent higher than for WW3-DIA (with a ratio of their peaks near 7/5). The inference is that the higher peak in WW3-WRT characterizes the importance of the full Boltzmann nonlinear interaction computation as compared to the "less nonlinear" discrete interaction approximation. If we compare the bidimensional wave spectra obtained from the WW3-WRT and WW3-DIA ( Figure 5), we can note further that the two channels in the WW3-WRT spectrum have considerably greater definition with respect to the channels in the WW3-DIA spectrum. Our results show a lower peak and reduced channels for the power spectrum of WW3-DIA ( Figure 5) with respect to WW3-WRT. This suggests that the WRT nonlinear interactions produce more realistic power spectra for large waves where high and peaked spectra occur. The enhanced power spectral features (high central peak and channels) suggest that the WRT algorithm produces a more physically meaningful directional spectral shape. Of course, the downside is that the WRT method requires substantial amounts of computer time, much more than the DIA approach.
The additional, wide-angle channels at low frequency, as observed in ocean measurements (see discussion in Section 4), are also possible in the Draupner storm, but in the present case, the exponential drop off at low frequency is so abrupt that the low-frequency channels are suppressed (see discussion in [17] for the physics of this process).
We analyzed term by term the source magnitudes for the WW3-WRT and WW3-DIA hindcasts and compared them in the whole spectral domain, see Figure 6. As could be expected, magnitudes are stronger for the case of WW3-WRT (top panels), which determines the shape of the power spectrum shown in Figure 5 (right panel), exhibiting a high and central peak with two channels. The Draupner spectrum obtained from the hindcast WW3-DIA (left panel in Figure 5) has smaller magnitudes for the source terms. Moreover, the DIA spectrum is broader in direction than the WRT spectrum. We have strong evidence that the S nl are the responsible process for the formation of the robust channels in the high-frequency part of the power spectrum (Figures 6 and 7). In particular, the shape that the S nl spectrum has at this time (1 January 1995 at 15:00 UTC), appears somewhat similar to the shape that the power spectrum has in the high frequencies with two strong, well-defined channels ( Figure 5). The full nonlinear interactions play an important role in the nonlinear wave dynamics and in the definition of the shape of the power spectrum.

On the Shape of a Two-Dimensional Spectrum Based on the Benjamin-Feir Instability
We now discuss the shape of the two-dimensional ocean wave power spectrum from experimental and theoretical points of view, with the measurements and analysis of data obtained by a number of authors [18][19][20][21]. One of the important features of these works is that the wave spectrum does not have the traditional shape of the typical JONSWAP shape with cosine squared spreading. To address this issue, [19] analyzed data, based on the ideas of [18], and developed a 2D spectral fitting procedure that was then used by [21] to represent the average spectrum measured in Currituck Sound (see also [13]. We have graphed this average spectrum in Figure 8. Note that the directional spectrum is characterized by a high central peak and channels that move off to high and low frequencies, frankly at angles and orientations not at first obvious from the point of view of the physics of nonlinear waves (see Figure 9). The goal of this section is to give a brief physical interpretation of this spectral shape in terms of the modulational instability.
We begin with the study of the 2D nonlinear Schrödinger equation (NLS) as conducted by [45]. This paper is a classical one for understanding both deterministic and stochastic behavior of We have strong evidence that the S nl are the responsible process for the formation of the robust channels in the high-frequency part of the power spectrum (Figures 6 and 7). In particular, the shape that the S nl spectrum has at this time (1 January 1995 at 15:00 UTC), appears somewhat similar to the shape that the power spectrum has in the high frequencies with two strong, well-defined channels ( Figure 5). The full nonlinear interactions play an important role in the nonlinear wave dynamics and in the definition of the shape of the power spectrum.

On the Shape of a Two-Dimensional Spectrum Based on the Benjamin-Feir Instability
We now discuss the shape of the two-dimensional ocean wave power spectrum from experimental and theoretical points of view, with the measurements and analysis of data obtained by a number of authors [18][19][20][21]. One of the important features of these works is that the wave spectrum does not have the traditional shape of the typical JONSWAP shape with cosine squared spreading. To address this issue, [19] analyzed data, based on the ideas of [18], and developed a 2D spectral fitting procedure that was then used by [21] to represent the average spectrum measured in Currituck Sound (see also [13]. We have graphed this average spectrum in Figure 8. Note that the directional spectrum is characterized by a high central peak and channels that move off to high and low frequencies, frankly at angles and orientations not at first obvious from the point of view of the physics of nonlinear waves (see Figure 9). The goal of this section is to give a brief physical interpretation of this spectral shape in terms of the modulational instability.
We begin with the study of the 2D nonlinear Schrödinger equation (NLS) as conducted by [45].

On the Shape of a Two-Dimensional Spectrum Based on the Benjamin-Feir Instability
We now discuss the shape of the two-dimensional ocean wave power spectrum from experimental and theoretical points of view, with the measurements and analysis of data obtained by a number of authors [18][19][20][21]. One of the important features of these works is that the wave spectrum does not have the traditional shape of the typical JONSWAP shape with cosine squared spreading. To address this issue, [19] analyzed data, based on the ideas of [18], and developed a 2D spectral fitting procedure that was then used by [21] to represent the average spectrum measured in Currituck Sound (see also [13]. We have graphed this average spectrum in Figure 8. Note that the directional spectrum is characterized by a high central peak and channels that move off to high and low frequencies, frankly at angles and orientations not at first obvious from the point of view of the physics of nonlinear waves (see Figure 9). The goal of this section is to give a brief physical interpretation of this spectral shape in terms of the modulational instability.
This expression says that the modulational frequency Ω = Ω(K x , Ky ) is a function of the x and y modulational wavenumbers, K x , Ky . Please note that the usual wavenumbers are given by k x = k o + K x , k y = K y , where k o is the carrier or peak frequency of the spectrum. Figure 9 shows the boundaries of the modulationally unstable region, which is shown by the color cyan. Inside the cyan region, the waves are unstable and so the frequency becomes imaginary: this is the region where a nonlinear wave component can grow exponentially, leading to rogue wave packets. Outside the cyan region, the waves are stable, with real frequency, and so we can at most form sine wave or Stokes wave components in the spectrum. Figure 8. Average directional spectrum estimated from Currituck Sound data as found by [21]. The parameterization is due to [19]. Graphed is the spectrum P( f, θ ) times frequency to the forth power: P( f ,θ ) f 4 , thus emphasizing directional properties at high frequency.
This contribution of Longuet-Higgins is quite important for describing waves governed by the modulational instability, and the cyan regions in Figures 9 and 10 tell us where the waves are modulationally unstable in the spectrum. In the lower panel of Figure 9, we show the Longuet-Higgins instability diagram after it has been transformed from wavenumber ( k x , ky ) to frequencydirection ( f , θ ) coordinates. Figure 9 also shows that the center of the power spectrum lies where the two straight lines cross in the upper panel and where the two curved lines cross in the lower panel. Also note the single contours in the power spectrum (ellipses), here shown for reference, in the upper and lower panels. Figure 8. Average directional spectrum estimated from Currituck Sound data as found by [21]. The parameterization is due to [19]. Graphed is the spectrum P( f , θ) times frequency to the forth power: P( f , θ) f 4 , thus emphasizing directional properties at high frequency.
We begin with the study of the 2D nonlinear Schrödinger equation (NLS) as conducted by [45]. This paper is a classical one for understanding both deterministic and stochastic behavior of nonlinear waves in two dimensions. Longuet-Higgins derived the modulational dispersion relation for the 2D NLS equation, which is given by: This expression says that the modulational frequency Ω = Ω(K x , K y ) is a function of the x and y modulational wavenumbers, K x , K y . Please note that the usual wavenumbers are given by k x = k o + K x , k y = K y , where k o is the carrier or peak frequency of the spectrum. Figure 9 shows the boundaries of the modulationally unstable region, which is shown by the color cyan. Inside the cyan region, the waves are unstable and so the frequency becomes imaginary: this is the region where a nonlinear wave component can grow exponentially, leading to rogue wave packets. Outside the cyan region, the waves are stable, with real frequency, and so we can at most form sine wave or Stokes wave components in the spectrum. As a final step, we take the lower panel in Figure 9 and place it above a graph of the power spectrum from Figure 8. The composite picture is given in Figure 10. It is then easy to see how the modulational channels in the upper panel correspond to the enhanced channels in the Currituck Sound spectrum, both at low and high frequency. These channels are seen to reflect the modulational instability of the 2D NLS equation. It is now physically and geometrically clear how "bimodality" occurs in the 2D spectrum of ocean waves. The enhanced regions in the power spectrum cause the channels to occur because of the presence of many small breathers in the unstable, cyan regions, oscillating up and down in space and time, whose average amplitude is larger than a single sine or Stokes mode. Of course, one expects large breathers near the peak of the spectrum, therefore giving rise to the high central peak.
We have suggested that the channels of the Longuet-Higgins instability diagram, after transforming from wavenumber coordinates ( k x , ky ) to frequency-direction coordinates ( f , θ ), correspond to the requisite channels seen in the power spectra of measured data. However, a deeper understanding is required to fully understand the actual shapes of the channels. Much of the machinery for carrying out these computations can be found in [46], albeit in brisk form. This accounts for the averaging over the modulational channels for each component in the spectrum, smoothing out the gentle shape of the channels. A future paper gives details. It is worth noting that any correct theoretical formulation would give spectra higher at the central peak than WRT because phase locking would be taken into account, thus physically giving rise to the presence of large breathers in the spectrum, normally centered about the peak [12,13]. This contribution of Longuet-Higgins is quite important for describing waves governed by the modulational instability, and the cyan regions in Figures 9 and 10 tell us where the waves are modulationally unstable in the spectrum. In the lower panel of Figure 9, we show the Longuet-Higgins instability diagram after it has been transformed from wavenumber (k x , k y ) to frequency-direction ( f , θ) coordinates. Figure 9 also shows that the center of the power spectrum lies where the two straight lines cross in the upper panel and where the two curved lines cross in the lower panel. Also note the single contours in the power spectrum (ellipses), here shown for reference, in the upper and lower panels.
As a final step, we take the lower panel in Figure 9 and place it above a graph of the power spectrum from Figure 8. The composite picture is given in Figure 10. It is then easy to see how the modulational channels in the upper panel correspond to the enhanced channels in the Currituck Sound spectrum, both at low and high frequency. These channels are seen to reflect the modulational instability of the 2D NLS equation. It is now physically and geometrically clear how "bimodality" occurs in the 2D spectrum of ocean waves. The enhanced regions in the power spectrum cause the channels to occur because of the presence of many small breathers in the unstable, cyan regions, oscillating up and down in space and time, whose average amplitude is larger than a single sine or Stokes mode. Of course, one expects large breathers near the peak of the spectrum, therefore giving rise to the high central peak.  [45]. The proximity of the channels in the modulational instability diagram with the channels on the average Currituck Sound spectrum are connected by arrows. This result suggests a one-to-one relationship between the shape of the instability diagram (upper) and the spectrum (lower). Then, one presumes that the "bimodality" in the tail of the spectrum (see Figure 4) is caused by the modulational instability.

Conclusions
We have employed WAVEWATCH-III to hindcast the Draupner storm of 1995, using the Webb-Resio-Tracy (WRT) method to compute the exact Boltzmann integrations, and also the discrete interaction approximation (DIA).
We found a spectral shape that is not generically a spread JONSWAP spectrum, but instead a spectrum that often has a high peak and enhanced channels at high and low frequencies. These results are similar to those found experimentally by [18][19][20][21]. We discussed the average shape of the power spectra of the experimental data from Currituck Sound [21], seen in Figure 8, where the peaked spectrum and low-and high-frequency channels seem to be characteristic of the nonlinear physics. Our point of view is that these channels are related to the physics of the modulational instability [13,17,46], and our results are supported by theoretical and experimental analyses. We have discussed how the spectrum should naturally have a high peak, together with "channels" at both high and low frequency. We note that Figure 5 gives the WW3-DIA spectrum in the lower panel, while the upper panel is the spectrum for WW3-WRT. The DIA hindcast does show some peakedness and highfrequency channels. However, the WRT hindcast shows a much higher central peak and more robust, enhanced high-frequency channels. Neither hindcast gives an indication of the lower frequency channels, an effect that does not always show up in data (see [17] for a discussion).
We suggest that a search for a fast algorithm for the WRT method should be pursued. Future work might include increasing the resolution for the spectrum for the frequency and direction bins in order to get additional features of the spectrum. It is hoped that the two-scale algorithm of Resio and colleagues [47][48][49] will provide users with the needed speed to do fully nonlinear forecasting. Additional work should make attempts to introduce wave coherency into wind/wave models. An important step has been the development of the generalized kinetic equation of [16].  [45]. The proximity of the channels in the modulational instability diagram with the channels on the average Currituck Sound spectrum are connected by arrows. This result suggests a one-to-one relationship between the shape of the instability diagram (upper) and the spectrum (lower). Then, one presumes that the "bimodality" in the tail of the spectrum (see Figure 4) is caused by the modulational instability.
We have suggested that the channels of the Longuet-Higgins instability diagram, after transforming from wavenumber coordinates (k x , k y ) to frequency-direction coordinates ( f , θ), correspond to the requisite channels seen in the power spectra of measured data. However, a deeper understanding is required to fully understand the actual shapes of the channels. Much of the machinery for carrying out these computations can be found in [46], albeit in brisk form. This accounts for the averaging over the modulational channels for each component in the spectrum, smoothing out the gentle shape of the channels. A future paper gives details. It is worth noting that any correct theoretical formulation would give spectra higher at the central peak than WRT because phase locking would be taken into account, thus physically giving rise to the presence of large breathers in the spectrum, normally centered about the peak [12,13].

Conclusions
We have employed WAVEWATCH-III to hindcast the Draupner storm of 1995, using the Webb-Resio-Tracy (WRT) method to compute the exact Boltzmann integrations, and also the discrete interaction approximation (DIA).
We found a spectral shape that is not generically a spread JONSWAP spectrum, but instead a spectrum that often has a high peak and enhanced channels at high and low frequencies. These results are similar to those found experimentally by [18][19][20][21]. We discussed the average shape of the power spectra of the experimental data from Currituck Sound [21], seen in Figure 8, where the peaked spectrum and low-and high-frequency channels seem to be characteristic of the nonlinear physics. Our point of view is that these channels are related to the physics of the modulational instability [13,17,46], and our results are supported by theoretical and experimental analyses. We have discussed how the spectrum should naturally have a high peak, together with "channels" at both high and low frequency. We note that Figure 5 gives the WW3-DIA spectrum in the lower panel, while the upper panel is the spectrum for WW3-WRT. The DIA hindcast does show some peakedness and high-frequency channels. However, the WRT hindcast shows a much higher central peak and more robust, enhanced high-frequency channels. Neither hindcast gives an indication of the lower frequency channels, an effect that does not always show up in data (see [17] for a discussion).
We suggest that a search for a fast algorithm for the WRT method should be pursued. Future work might include increasing the resolution for the spectrum for the frequency and direction bins in order to get additional features of the spectrum. It is hoped that the two-scale algorithm of Resio and colleagues [47][48][49] will provide users with the needed speed to do fully nonlinear forecasting. Additional work should make attempts to introduce wave coherency into wind/wave models. An important step has been the development of the generalized kinetic equation of [16].