Four-wave interactions are shown to play an important role in the evolution of the spectrum of surface gravity waves. This fact follows from direct simulations of an ensemble of ocean waves using the Zakharov equation. The theory of homogeneous four-wave interactions, extended to include effects of nonresonant transfer, compares favorably with the ensemble-averaged results of the Monte Carlo simulations. In particular, there is good agreement regarding spectral shape. Also, the kurtosis of the surface elevation probability distribution is determined well by theory even for waves with a narrow spectrum and large steepness. These extreme conditions are favorable for the occurrence of freak waves.
At present there is a considerable interest in understanding the occurrence of freak waves. The notion of freak waves was first introduced by Draper (1965), and this term is applied for single waves that are extremely unlikely as judged by the Rayleigh distribution of wave heights (Dean 1990). In practice this means that when one studies wave records of a finite length (say of 10–20 min), a wave is considered to be a freak wave if the wave height H (defined as the distance from crest to trough) exceeds the significant wave height HS by a factor 2.2. It is difficult to collect hard evidence on such extreme wave phenomena because they occur so rarely. Nevertheless, observational evidence from time series collected over the past decade does suggest that for large surface elevations the probability distribution for the surface elevation deviates substantially from the one that follows from linear theory with random phase, namely the Gaussian distribution (e.g., Wolfram and Linfoot 2000).
There are a number of reasons why freak wave phenomena may occur. Often, extreme wave events can be explained by the presence of ocean currents or bottom topography that may cause wave energy to focus in a small area because of refraction, reflection, and wave trapping. These mechanisms are well understood and may be explained by linear wave theory (e.g., Lavrenov 1998).
Trulsen and Dysthe (1997) argue, however, that it is not well understood why exceptionally large waves may occur in the open ocean away from nonuniform currents or bathymetry. As an example they discuss the case of an extreme wave event that happened on 1 January 1995 in the Norwegian sector of the North Sea. Their basic premise is that these waves can be produced by nonlinear self modulation of a slowly varying wave train. An example of nonlinear modulation or focusing is the instability of a uniform narrowband wave train to side-band perturbations. This instability—known as the side-band, modulational, or Benjamin–Feir instability (Benjamin and Feir 1967)—will result in focusing of wave energy in space and/or time as is illustrated by the experiments of Lake et al. (1977).
To a first approximation, the evolution in time of the envelope of a narrowband wave train is described by the nonlinear Schrödinger equation. This equation, which occurs in many branches of physics, was first discussed in the general context of nonlinear dispersive waves by Benney and Newell (1967) (see also Ostrowskii 1967). For water waves it was first derived by Zakharov (1968) using a spectral method and by Hasimoto and Ono (1972) and Davey (1972) using multiple-scale methods. The nonlinear Schrödinger equation in one-space dimension may be solved by means of the inverse scattering transform. For vanishing boundary conditions, Zakharov and Shabat (1972) found that for large times the solution consists of a combination of envelope solitons and radiation modes, in analogy with the solution of the Korteweg–de Vries equation. However, for two-dimensional propagation, Zakharov and Rubenchik (1974) discovered that envelope solitons are unstable to transverse perturbations, and Cohen et al. (1976) found that a random wave field would break up envelope solitons. This meant that solitons could not be used as building blocks of the nonlinear evolution of gravity waves.
For periodic boundary conditions, the solution of the nonlinear Schrödinger equation is more complex. Linear stability analysis of a uniform wave train shows that close side bands grow exponentially in time in good qualitative agreement with the experimental results of Benjamin and Feir (1967) and Lake et al. (1977). For large times there is a considerable energy transfer from the carrier wave to the side bands. In one-space dimension (1D) if there is only one unstable side band, Fermi–Pasta–Ulam recurrence occurs (Yuen and Ferguson 1978) in qualitative agreement with the experiments of Lake et al. (1977). In the presence of many unstable side bands, the evolution of a narrowband wave train becomes much more complex. No recurrence is then found (Caponi et al. 1982), and these authors have termed this confined chaos in a nonlinear wave system because most of the energy resides in the unstable modes. Also, in two-space dimensions (2D) the phenomenon of recurrence is the exception rather than the rule. In addition, in 2D the instability region is unbounded in the perturbation wave vector space, resulting in energy leakage to high wave-number modes; hence there is no confined chaos in 2D (Martin and Yuen 1980). This suggests that the 2D nonlinear Schrödinger equation is inadequate to describe the evolution of weakly nonlinear waves. This was pointed out already by Longuet-Higgins (1978) who performed a stability analysis on the exact equations and found that the instability region is finite in extent. More-realistic evolution equations such as the fourth-order evolution equation of Dysthe (1979) or the Zakharov equation (1968) are needed to give an appropriate description of nonlinear gravity waves in two-space dimensions.
Nevertheless, studies of the properties of the nonlinear Schrödinger equation have been vital in understanding the conditions under which freak waves may occur. This was discussed in detail by Osborne et al. (2000). For periodic boundary conditions, the one-dimensional nonlinear Schrödinger equation may be solved by the inverse scattering method as well. The role of the solitons is then replaced by unstable modes. In the linear regime, these modes just describe the evolution in time according to the Benjamin–Feir instability, whereas, by means of the inverse scattering transform, the fate of the unstable mode may be followed right into the nonlinear regime. Using the inverse scattering transform, the solution of the 1D nonlinear Schrödinger equation may be written as a “linear” superposition of stable modes, unstable modes, and their mutual nonlinear interactions. Here, the stable modes form a Gaussian background wave field from which the unstable modes occasionally rise up and subsequently disappear again, repeating the process quasi-periodically in time. Making use of the inverse scattering transform, these authors readily construct a few examples of giant waves from the one-dimensional nonlinear Schrödinger equation. The question now is what happens in the case of two-dimensional propagation. The notion of solitons is no longer useful, because solitons are unstable in two dimensions. Osborne et al. (2000) show that unstable modes do indeed still exist and that in the nonlinear regime they can take the form of large-amplitude freak waves. Furthermore, the notion of unstable modes seems to be a generic property of deep-water wave trains as the authors find nonlinear unstable modes in both the one- and two-dimensional versions of Dysthe's fourth-order evolution equation. To summarize this discussion, it seems that freak waves are likely to occur as long as the wave train is subject to nonlinear focusing. In addition, we only need to study the case of one-dimensional propagation, because it captures the essentials of the generation of freak waves.
Therefore, in the context of the deterministic approach to wave evolution, there seems to be a reasonable theoretical understanding of why in the open ocean freak waves occur. In ocean wave forecasting practice one follows, however, a stochastic approach; that is, one attempts to predict the ensemble average of a spectrum of random waves because knowledge on the phases is not available. The main problem then is to what extent one can make statements regarding the occurrence of freak waves in a random wave field. Of course, in the context of wave forecasting only statements of a probablistic nature can be made. Because freak waves imply considerable deviations from the normal Gaussian probability distribution function (pdf) of the surface elevation, the main question therefore is whether one can determine in a reliable manner the pdf of the surface elevation. Because the wave spectrum plays a central role in the stochastic approach, the question therefore is whether for given wave spectrum the probability of extreme events may be determined.
Present-day wave forecasting systems are based on the energy balance equation (Komen et al. 1994), including a parameterized version of Hasselmann's four-wave nonlinear transfer (Hasselmann 1962). Resonant four-wave interactions for a random, homogeneous sea play an important role in the evolution of the spectrum of wind waves, because on the one hand they determine the high-frequency part of the spectrum, giving rise to an ω−4 tail (Zakharov and Filonenko 1968) while on the other hand the peak of the spectrum is shifted toward lower frequencies. The homogeneous nonlinear interactions give rise to deviations from the Gaussian pdf for the surface elevation, because the third-order nonlinearity generates fourth cumulants of the pdf while the finite fourth cumulant results in spectral change. An important issue is, however, whether the standard homogeneous theory can properly describe the generation of freak waves, simply because it does not seem to incorporate the Benjamin–Feir instability mechanism (Alber 1978; Alber and Saffman 1978; Crawford et al. 1980; Janssen 1983b). This follows from simple scaling considerations applied to the Hasselmann evolution equation for four-wave interactions. Because the rate of change of the action density N is proportional to N3, the nonlinear transfer occurs on the time scale TNL = O(1/ε4ω0). Here, ε is a typical wave steepness, which is assumed to be small, and ω0 is a typical angular frequency of the wave field. In contrast, the Benjamin–Feir instability occurs on the much faster timescale of O(1/ε2ω0).
The Benjamin–Feir instability is an example of a nonresonant four-wave interaction in which the carrier wave is phase-locked with the side bands. This process cannot be described by a theory that assumes that the Fourier amplitudes are not correlated (i.e., a homogeneous wave field) and in which only resonant four-wave interactions are considered. For an inhomogeneous, Gaussian narrowband wave train, Alber and Saffman (1978) and Alber (1978) derived an evolution equation for the Wigner distribution of the sea state. Inhomogeneities gave rise to a much faster energy transfer, comparable with the typical timescale of the modulational instability. In fact, these authors discovered the random version of the Benjamin–Feir instability: a random narrowband wave train is unstable to side-band perturbations provided the width of the spectrum is sufficiently narrow. Therefore, one would expect the Alber and Saffman approach to be an ideal starting point for treating freak waves in a random wave context. However, it is emphasized that this approach has it limitations because deviations from normality have not yet been taken into account. In this paper it will be shown, using numerical simulations of an ensemble of ocean waves, that non-Gaussian effects are important while inhomogeneities play only a minor role in the evolution of the ensemble-averaged wave spectrum.
On the other hand, nonresonant interactions appear to be relevant. Hasselmann's treatment of four-wave interactions is extended by including the effects of nonresonant interactions. As a consequence, the resonance function is, for short times, broader than the usual δ function and depends on the angular frequency resonance conditions and on time. The standard nonlinear transfer is based on the assumption that the action density spectrum is a slowly varying function of time. It is then argued that the resonance function may be replaced by its large time limit, giving the usual delta function. However, the time span required for the resonance function to evolve toward a delta function is so large that considerable changes in the action density function may have occurred in the meantime. This will be shown for the special case of one-dimensional propagation of surface gravity waves. In those circumstances, the standard approach to nonlinear wave–wave interactions would not give rise to nonlinear transfer, whereas considerable changes of the wave spectrum occur in the new approach. In fact, there is close agreement between results on the ensemble-averaged spectrum and the kurtosis of the pdf of the surface elevation, as obtained from numerical simulations of an ensemble of ocean waves. Because time series from the numerical simulations indicate the occurence of freak waves when the waves are sufficiently steep [see also Trulsen and Dysthe (1997) or Osborne et al. (2000)], the implication is that an approach to nonlinear transfer that includes nonresonant interactions seems to capture freak wave events. However, it is strongly emphasized that such an approach can only give statements of a probablistic nature on the occurrence of extreme wave events.
The structure of this paper is as follows. In section 2 I review developments regarding the evolution of a random wave field, but I discuss only the ideas needed for understanding results in the remainder of this paper. In particular, I extend the standard theory of four-wave interactions by including effects of nonresonant interactions and derive an explicit expression for the kurtosis in terms of the action density spectrum. I also discuss Alber and Saffman's key result, that according to lowest-order inhomogeneous theory there is only Benjamin–Feir instability when the wave spectrum is sufficiently narrow. In section 3 I present results from Monte Carlo simulations of the nonlinear Schrödinger equation following similar work by Onorato et al. (2000). Only one-dimensional wave propagation is discussed. Apart from reasons of economy (runs are typically done with 500-member ensembles), the main reason for this choice is that for one dimension the nonlinear transfer according to the standard homogeneous theory of four-wave interactions vanishes identically. The ensemble-averaged evolution of the wave spectrum clearly shows that there is an irreversible energy transfer resulting in a broadening of the spectrum while the pdf of the surface elevation has considerable deviations from the Gaussian distribution. These deviations from normality may be described, as expected from four-wave interactions, by means of the fourth cumulant. In case of nonlinear focusing, the correction to the pdf is such that there is an enhanced probability of extreme events, while in the case of nonlinear defocusing (this was achieved by changing the sign of the nonlinear term) the opposite occurs, namely the probability of extreme events is reduced. This is in agreement with results by Tanaka (1992) who found an increase in groupiness in the case of nonlinear focusing while in the opposite case of a stable wave train groupiness decreases. Defocusing of surface gravity wave trains occurs in shallow waters when the parameter k0D (with k0 being a typical wave number and D the depth) is less than 1.36 (Mori and Yasuda 2002a). In principle, the approach could be extended to the case of shallow water to study what happens with the probability of extreme events when k0D < 1.36. However, this would introduce an extra complication. Nonlinear defocusing will therefore only be discussed in a qualitative sense, by changing the sign of the nonlinear term.
Both the spectral broadening and the fourth cumulant (or kurtosis) are found to depend on a single parameter characterizing the narrowband wave train, namely the ratio of mean square slope to the normalized width of the (frequency) spectrum. It is suggested to call this ratio the Benjamin–Feir index (BFI). If the BFI is larger than 1, then according to Alber and Saffman (1978) the random wave field is modulationally unstable. This result would suggest that if the BFI is less than 1 no changes in the spectrum occur, whereas in the opposite case the unstable side bands would give rise to a broadening of the wave spectrum. Hence, BFI = 1 is a bifurcation point. The numerical simulations provide no convincing evidence of a bifurcation at BFI = 1. Rather, there is already a considerable broadening of the wave spectrum around BFI = 1, while the dependence of the broadening on the BFI appears to be smooth rather then abrupt (cf. Tanaka 1992).
I continue in section 3 by presenting results from Monte Carlo simulations of the Zakharov equation (Zakharov 1968). Results are similar in spirit to those obtained with the nonlinear Schrödinger equation, except that the modulational instability seems to occur for larger BFI. For the nonlinear Schrödinger equation, the spectral change owing to nonlinear transfer is symmetrical with respect to the spectral maximum, but this is not the case for Zakharov equation. In the latter case, the nonlinear transfer coefficients and the angular frequency are asymmetrical with respect to the spectral peak and as a consequence there is a downshift of the peak of the spectrum. It is emphasized that this downshift occurs in the absence of dissipation, whereas quantities such as action, wave momentum, and total wave energy are conserved.
In section 4, an interpretation of the numerical results of section 3 is given. First, it is shown that inhomogeneities only play a minor role in the evolution of the wave spectrum, and deviations from normality are more relevant. Second, results from the numerical solution of the extended version of Hasselmann's wave–wave interaction approach are presented and compared with the results from Monte Carlo simulations. A good agreement is obtained. Apart from the fact that I have given a direct validation of Hasselmann's four-wave theory, it also shows that even in extreme conditions such as occur during the generation of freak waves, reliable estimates of deviations from normality can be made.
In section 5 a summary of conclusions is given. Much to my surprise, effects of inhomogeneity only play a minor role in understanding the ensemble-averaged evolution of surface gravity waves. Homogeneous four-wave interactions, albeit extended by allowing for a time-dependent resonance function, seem to capture most essential features of the averaged nonlinear wave evolution. It now seems possible to estimate the enhanced occurrence of extreme waves and freak waves on the open ocean, because the kurtosis may be estimated directly from the wave spectrum.
2. Review and extension of the theory of a random wave field
The starting point is the Zakharov equation, which is a deterministic evolution equation for surface gravity waves in deep water. It is obtained from the Hamiltonian for water waves, first found by Zakharov (1968). Consider the potential flow of an ideal fluid of infinite depth. Coordinates are chosen in such a way that the undisturbed surface of the fluid coincides with the x–y plane. The z axis is pointed upward, and the acceleration of gravity g is pointed in the negative z direction. Let η be the shape of the surface of the fluid and let ϕ be the potential of the flow. Hence, the velocity of the flow follows from u = −∇ϕ.
By choosing as canonical variables
Zakharov (1968) showed that the total energy E of the fluid may be used as a Hamiltonian. Here,
The x integrals extend over the total basin considered. If an infinite basin is considered, the resulting total energy is infinite, unless the wave motion is localized within a finite region. This problem may be avoided by introducing the energy per unit area by dividing (2) by the total surface L × L, where L is the length of the basin, and taking the limit of L → ∞ afterward. As a consequence, integrals over wavenumber k are replaced by summations while δ functions are replaced by Kronecker δs. [For a more complete discussion, see Komen et al. (1994).] I will adopt this approach implicitly in the remainder of this paper.
The boundary conditions at the surface, namely the kinematic boundary condition and Bernoulli's equation, are then equivalent to Hamilton's equations,
where δE/δψ is the functional derivative of E with respect to ψ, and so on. Inside the fluid, the potential ϕ satisfies Laplace's equation,
with boundary conditions
If one is able to solve the potential problem, then ϕ may be expressed in terms of the canonical variables η and ψ. Then the energy E may be evaluated in terms of the canonical variables, and the evolution in time of η and ψ follows at once from Hamilton's equations (3). This was done by Zakharov (1968), who obtained the deterministic evolution equations for deep-water waves by solving the potential problem (4)–(6) in an iterative fashion for small steepness ε. In addition, the Fourier transforms of η and ϕ were introduced, and results could be expressed in a concise way by use of the action variable A(k, t). For example, in terms of A, the surface elevation η becomes
Here, k is the wavenumber vector, k is its absolute value, and ω = gk denotes the dispersion relation of deep-water gravity waves. Substitution of the series solution for ϕ into the Hamiltonian (2) gives an expansion of the total energy E of the fluid in terms of wave steepness,
Retaining only the second-order term of E corresponds to the linear theory of surface gravity waves, the third-order term corresponds to three-wave interactions, and the fourth-order term corresponds to four wave interactions. Because resonant three-wave interactions are absent for deep-water gravity waves, a meaningful description of the wave field is only obtained by going to fourth order in ε. In fact, Krasitskii (1990) has shown that in the absence of resonant three-wave interactions there is a nonsingular, canonical transformation from the action variable A to the new variable a that allows elimination of the third-order contribution to the wave energy. In loose terms, the new variable a describes the free-wave part of the wave field. Apart from a constant factor, the energy of the free waves becomes
where a1 = a(k1), and so on, δ is the Dirac delta function, and the interaction coefficient T is given by Krasitskii (1990). The interaction coefficient enjoys a number of symmetry conditions, of which the most important one is T1,2,3,4 = T3,4,1,2, because this condition implies that E is conserved. Hamilton's equations now become the single equation
and, evaluating the functional derivative of E with respect to a*, the evolution equation for a becomes
known as the Zakharov equation. Apart from the free-wave energy (9), the Zakharov equation admits conservation of action and of wave momentum, respectively, as
a. Comments on the Zakharov equation
The properties of the Zakharov equation have been studied in great detail by, for example, Crawford et al. (1981) [for an overview see Yuen and Lake (1982)]. Thus the nonlinear dispersion relation, first obtained by Stokes (1947), follows from (11), and also the instability of a weakly nonlinear, uniform wave train (the so-called Benjamin–Feir instability) is described well by the Zakharov equation; the results on growth rates, for example, are qualitatively in good agreement with the results of Longuet-Higgins (1978). However, these results were obtained with a form of the interaction coefficient T that did not result in a Hamiltonian form of (11). Krasitskii (1990) found the correct canonical transformation to eliminate the cubic interactions, which resulted in a T that satisfied the appropriate symmetry conditions for (11) to be Hamiltonian. Krasitskii and Kalmykov (1993) studied the differences between the Hamiltonian and the non-Hamiltonian forms of the Zakharov equation, but only for large amplitude were differences in the solution found.
In this paper, I initially use a narrowband approximation to the Zakharov equation, because the main impact of the Benjamin–Feir instability is found near the spectral peak. This approximate evolution equation is obtained by means of a Taylor expansion of angular frequency ω and the interaction coefficient T around the carrier wavenumber k0. The nonlinear Schrödinger equation is then obtained by using only the lowest-order approximation to T given by k30, and angular frequency ω is expanded to second order in the modulation wavenumber p = k − k0. The main advantage of the use of the nonlinear Schrödinger equation is that many properties of this equation are known and that it can be solved numerically in an efficient way. The draw-back, however, is that it overestimates the growth rates of the Benjamin–Feir instability and that the nonlinear energy transfer is symmetrical with respect to the carrier wavenumber. For this reason, I study solutions of the complete Zakharov equation as well, using the Krasitskii (1990) expression for the interaction coefficient T. Similarly, one could study higher-order evolution equations such as the one by Dysthe (1979), but I found that spectra may become so broad that the narrowband approximation becomes invalid.
Another reason for studying the nonlinear Schrödinger equation is that it allows one to introduce an important parameter that will be used to stratify the numerical and theoretical results. From the physical point of view, we are basically studying a problem that concerns the balance between dispersion of the waves and its nonlinearity. For the full Zakharov equation, it will be difficult to introduce a unique measure of, for example, nonlinearity because the nonlinear transfer coefficient T is a complicated function of wavenumber. However, in the narrowband approximation, giving the nonlinear Schrödinger equation, this is more straightforward to do. Balancing the nonlinear term and the dispersive term in the narrowband version of (11) therefore gives the dimensionless number
where ω0 is the angular frequency at k0. Because my interest is in the dynamics of a continuous spectrum of waves, the slope parameter ε and the relative width σ′ω of the frequency spectrum relate to spectral properties, hence ε = (k20〈η2〉)1/2, with 〈η2〉 being the average surface elevation variance, and σ′ω = σω/ω0. For positive sign of the dimensionless parameter (13), there is focusing (modulational instability) while in the opposite case there is defocusing of the weakly nonlinear wave train. Based on this, I introduce the Benjamin–Feir index, which, apart from a constant, is the square root of the dimensionless number (13). Using the dispersion relation for deep-water gravity waves and the expression for the nonlinear interaction coefficient, T0 = k30, the BFI becomes
The BFI turns out to be very useful in ordering the theoretical and numerical results presented in the following sections. For simple initial wave spectra (defined in terms of the modulation wavenumber p) that only depend on the variance and on the spectral width, it can be shown that for the nonlinear Schrödinger equation the large-time solution is completely characterized by the BFI. For the Zakharov equation this is not the case, but the BFI is still expected to be a useful parameter for narrowband wave trains. The BFI plays a key role in the inhomogeneous theory of wave–wave interactions (Alber 1978), and a similar parameter has been introduced and discussed by Onorato et al. (2001) in the context of freak waves in random sea states.
b. Stochastic approach
The Zakharov equation [(11)] predicts amplitude and phase of the waves. For practical applications such as wave prediction, the detailed information regarding the phase of the waves is not available. Therefore, at best one can hope to predict average quantities such as the second moment,
where the angle brackets denote an ensemble average. Here, we briefly sketch the derivation of the evolution equation for the second moment from the Zakharov equation, assuming a zero mean value, 〈a1〉 = 0. It is known, however, that because of nonlinearity, the evolution of the second moment is determined by the fourth moment, and so on, resulting in an infinite hierarchy of equations (Davidson 1972). To obtain a meaningful truncation of this hierachy, it is customary to assume that the sea surface is close to a Gaussian state. This means that the amplitudes a1 are simultaneously Gaussian, an assumption that is a reasonable one for small wave steepness ε. In that event, higher-order moments can be expressed in lower-order moments. For a zero-mean stochastic variable a, the fourth moment can be written as (Hasselmann 1962; Crawford et al. 1980)
Here, it is assumed that, when compared with Bj,k, the moment 〈ajak〉 is small (note that for homogeneous fields this moment vanishes). In addition, D is the so-called fourth cumulant, which vanishes for a Gaussian sea state. Resonant nonlinear interactions, however, will tend to create correlations in such a way that a finite fourth cumulant results. However, for small steepness, D is expected to be small, so that an approximate closure of the infinite hierarchy of equations may be achieved.
Let us now sketch the derivation of the evolution equation for the second moment 〈aia*j〉 from the Zakharov equation [(11)]. To that end, I multiply (11) for ai by a*j, add the complex conjugate with i and j interchanged, and take the ensemble average:
where c.c. denotes complex conjugate, and i j denotes the operation of interchanging indices i and j in the previous term. Because of nonlinearity, the equation for the second moment involves the fourth moment. Similarly, the equation for the fourth moment involves the sixth moment. It becomes
So far, no approximations have been made. In the next section, I discuss the implications of the assumptions of a homogeneous weakly nonlinear wave field. Homogeneity of the wave field, however, does not allow a description of the Benjamin–Feir instability, and therefore in the following section I discuss the consequences for spectral evolution when the wave field is allowed to be inhomogeneous.
c. Evolution of a homogeneous random wave field
A wave field is considered to be homogeneous if the two-point correlation function 〈η(x1)η(x2)〉 depends only on the distance x1 − x2. Using the expression for the surface elevation, (7), it is then straightforward to verify that a wave field is homogeneous provided that the second moment Bi,j satisfies
where Ni is the spectral action density, which is equivalent to a number density because ωiNi is the spectral energy density, while kiNi is the spectral momentum density (apart from a factor ρw, the water density).
For weakly nonlinear waves, the fourth cumulant D is small when compared with the product of second-order cumulants [this may be verified afterward; it follows immediately from (18)]. Now, invoking the random-phase approximation [i.e. (16)] with D = 0 on (17), combined with the assumption of a homogeneous wave field, results in constancy of the second moment Bi,j. Hence, the need to go to higher order; that is, the fourth moment has to be determined through (18).
Application of the random phase approximation to the sixth moment (which implies that the sixth cumulant is ignored) and solving (18) for the fourth cumulant D, subject to the initial condition D(t = 0) = 0, gives
where Δω is shorthand for ωi + ωj − ωk − ωl, and I have made extensive use of the symmetry properties of the nonlinear transfer coefficient T, in particular the Hamiltonian symmetry. In addition, I used the property that, according to (17), the action density N only evolves on the slow timescale. The function G is defined as
The function G develops for large time t into the usual generalized functions P/Δω (where P means the Cauchy principal value), and δ(Δω), since,
The limit in (24) is a limit in the sense of generalized functions and is, in a strict sense, only meaningful inside integrals over wavenumber when multiplied by a smooth function.
where now Δω = ω1 + ω2 − ω3 − ω4. This evolution equation is usually called the Boltzmann equation.
Two limits of the resonance function Ri(Δω, t) are of interest to mention. For small times one has
and for large times one has
Hence, according to (25), for short times the evolution of the action density N is caused by both resonant and nonresonant four-wave interactions, and for large times, when the resonance function evolves toward a δ function, only resonant interactions contribute to spectral change.
In the standard treatment of resonant wave–wave interactions (e.g., Hasselmann 1962; Davidson 1972) it is argued that the resonance function Ri(Δω, t) may be replaced by its time-asymptotic value [(27)], because the action density spectrum is a slowly varying function of time. However, the time required for the resonance function to evolve toward a delta function may be so large that in the meantime considerable changes in the action density may have occurred. For this reason I will keep the full expression for the resonance function.
An important consequence of this choice concerns the estimation of a typical time scale TNL for the nonlinear wave–wave interactions in a homogeneous wave field. With ε being a typical wave steepness and ω0 being a typical angular frequency of the wave field, one finds from the Boltzmann equation [(25)] that for short times TNL = O(1/ε2ω0), and for large times TNL = O(1/ε4ω0). Hence, although the standard nonlinear transfer, which uses (27) as resonance function, does not capture the physics of the modulational instability (which operates on the fast timescale 1/ε2ω0), the full resonance function does not suffer from this defect.
It is also important to note that according to the standard theory there is only nonlinear transfer for two-dimensional wave propagation. In the one-dimensional case there is no nonlinear transfer in a homogeneous wave field. The reason for this is that only those waves interact nonlinearly that satisfy the resonance conditions k1 + k2 = k3 + k4 and ω1 + ω2 = ω3 + ω4. In one dimension these resonance conditions can only be met for the combinations k1 = k3, k2 = k4 or k1 = k4, k2 = k3. Then, the rate of change of the action density, as given by (25) and (27), vanishes identically because of the symmetry properties of the term involving the action densities. This contrasts with the Benjamin–Feir instability, which has its largest growth rates for waves in one dimension. On the other hand, using the complete expression for the resonance function, there is always an irreversible nonlinear transfer even in the case of one-dimensional propagation.
The Boltzmann equation, (25), admits just as the deterministic Zakharov equation, conservation of total action, wave momentum, and the ensemble average of the Hamiltonian [(9)] is conserved as well. [The last conservation law follows from (25) by consistently utilizing the assumption of a slowly varying action density.] It is emphasized that the Hamiltonian consists of two parts, the energy according to linear wave theory, and a nonlinear interaction term. Therefore, unlike the standard theory of four-wave interactions, the linear expression for the wave energy is not conserved. The exception occurs for large times when the resonance function Ri has evolved toward a δ function, and then just as in the standard theory the linear wave energy is conserved. This follows also from the numerical simulations presented in section 3, which show that the ensemble average of the Hamiltonian is conserved but, in particular for short times, not the linear wave energy. Furthermore, it should be mentioned that the Boltzmann equation [(25)] has the time reversal symmetry of the original Zakharov equation, because the resonance function changes sign when time t changes sign. Also, as Ri vanishes for t = 0, the time derivative of the action density spectrum is continuous around t = 0 and does not show a cusp (Komen et al. 1994). Nevertheless, despite the fact that there is time reversal, (25) has the irreversibility property: the memory of the initial conditions gets lost in the course of time owing to phase mixing.
The standard nonlinear transfer in a homogeneous wave field has been studied extensively in the past four decades. The Joint North Sea Wave Project (JONSWAP) study (Hasselmann et al. 1973) has shown the prominent role played by four-wave interactions in shaping the wave spectrum and in shifting the peak of the spectrum toward lower frequencies. Modern wave forecasting systems therefore use a parameterization of the nonlinear transfer (Komen et al. 1994).
The main interest in this paper is in the statistical aspects of random, weakly nonlinear waves in the context of the Zakharov equation. In particular I am interested in the relation between the deviations from the Gaussian distribution and four-wave interactions. Because of the symmetries of the Zakharov equation, the first moment of interest is then the fourth moment and the related kurtosis. The third moment and its related skewness vanish: information on the odd moments can only be obtained by making explicit use of Krasitskii's (1990) canonical transformation. Now, the fourth moment 〈η4〉 may be obtained in a straightforward manner from (16) and the expression for the fourth cumulant [(20)] as
Denoting the second moment 〈η2〉 by m0, deviations from normality are then most conveniently established by calculating the kurtosis
because for a Gaussian pdf C4 vanishes. The result for C4 is
where Rr is defined by (22), and the integral should be interpretated as a principal value integral. For large times, unlike the evolution of the action density, the kurtosis does not involve a Dirac δ function but rather depends on P/Δω. Therefore, the kurtosis is determined by the resonant and nonresonant interactions. It is instructive to apply (29) to the case of a narrowband wave spectrum in one dimension. Hence, performing the usual Taylor expansions around the carrier wavenumber k0 to lowest significant order, one finds for large times
where p = k − k0 is the wavenumber with respect to the carrier. It is seen that the sign of the kurtosis is determined by the ratio T0/ω″0, which is the same parameter that determines whether a wave train is stable to side-band perturbations. Note that numerically the integral is found to be negative, at least for bell-shaped spectra. Hence, from (30) it is immediately plausible that for an unstable wave system, which has negative T0/ω″0, the kurtosis will be positive and thus will result in an increased probability of extreme events. On the other hand, for a stable wave system there will a reduction in the probability of extreme events.
A further simplification of the expression for the kurtosis may be achieved if it is assumed that the wavenumber spectrum F(p) = ω0N(p)/g only depends on two parameters, namely, the variance m0 and the spectral width σk. Introduce the scaled wavenumber x = p/σk and the correspondingly scaled spectrum m0H(x)dx = F(p)dp. Then, using the deep-water dispersion relation and T0 = k30, (30) becomes
where ε is the significant steepness k0m1/20 and σ′ω is the relative width in angular frequency space σω/ω0 = 0.5σk/k0. The parameter J is given by the expression
d. Evolution of an inhomogeneous random wave field
The Benjamin–Feir instability is the result of a nonlinear interaction of waves that are phase locked, as the carrier wave is phase locked with the side bands, and therefore this process cannot be described by a theory that assumes that the Fourier amplitudes are not correlated, as expressed by the assumption of homogeneity of the wave field [(19)]. Therefore, this suggests that local nonlinear events such as freak waves could be beyond the scope of the standard description of ocean waves.
The investigation of the effect of inhomogeneities on the nonlinear energy transfer started with the work of Alber (1978) and Alber and Saffman (1978), whereas Crawford et al. (1980) combined the effects of inhomogeneity and nonnormality on the evolution of weakly nonlinear water waves. A review of this may be found in Yuen and Lake (1982). I will only discuss the lowest-order effects of inhomogeneity, disregarding any effects resulting from deviations from normality, and I only discuss one-dimensional wave propagation.
Hence, I do not impose the condition of a homogeneous wave field [(19)]. Now invoking the Gaussian approximation on the fourth moment [(16) with D = 0] and substituting the result in the evolution equation for the second moment, (17), gives
Here, we used the property that the second moment B is hermitian, Bi,j = B*j,i, and we made use of the symmetry properties of T.
In principle, (32) could be used to study the (in)stability of a homogeneous wave spectra, but to my knowledge this has not been done so far. Instead of this, Alber (1978) and Alber and Saffman (1978) studied the stability of a narrowband homogeneous wave spectrum. Following Crawford et al. (1980) and Yuen and Lake (1982), a considerable simplification of the evolution equation for Bi,j may be achieved by expanding angular frequency ω and interaction coefficient T around the carrier wavenumber k0. At the same time, one introduces the sum and difference wavenumbers
and one introduces the relative wavenumber p = n − k0. The correlation function B is from now on regarded as a function m and n. Realizing that in the narrowband approximation n is close to k0 while m is small, one obtains from (32) the following approximate evolution equation for B:
Here, a prime denotes differentation with respect to the carrier wavenumber k0, and T0 = k30. A key role in the work of Alber and Saffman is played by the envelope spectral function W(p, x, t), which is in fact a Wigner distribution (Wigner 1932). It is related to the Fourier transform of B(n, m, t) with respect to m,
and a homogeneous sea state simply has a Wigner distribution that is independent of the spatial coordinate x, in agreement with the definition of homogeneous sea given in (19). In terms of the Wigner distribution, (34) becomes a transport equation in x, p, and t, which bears a similarity with the Vlasov equation from plasma physics. This transport equation is obtained by means of a Taylor expansion of the difference term in the right-hand side of (34) with respect to l, giving an infinite sum. The result is
where ρ(x, t) = 2〈η2〉 is the mean-square envelope variance, given by
and the dots on the right-hand side of (36) represent the remaining terms of the Taylor series expansion. Note that all terms of the series are required to recover properly the random version of the Benjamin–Feir instability.
Alber and Saffman (1978) and Alber (1978) studied the stability of a homogeneous spectrum and found that it is unstable to long-wavelength perturbations if the width of the spectrum is sufficiently small. In other words, in case of instability, inhomogeneities would be generated by what is termed the random version of the Benjamin–Feir instability, therefore violating the assumption of homogeneity made in the standard theory of wave–wave interactions.
To see whether a homogeneous spectrum W0(p) is stable, one proceeds in the usual fashion by perturbing W0(p) slightly according to
Linearizing the evolution equation for W around the equilibrium W0 and considering normal mode perturbations, one obtains a dispersion relation between the angular frequency ω and the wavenumber k of the perturbation. Instability is found for Im(ω) > 0. Alber (1978) considered as special case the Gaussian spectrum
where 〈a20〉 is a constant envelope variance and σk is the width of the spectrum in wavenumber space. Stability of the random wave train was found when the relative width of the spectrum, σk/k0, exceeds a measure of mean-square slope. In terms of the relative width σω/ω0 of the frequency spectrum, which is just one-half of the relative width for the wavenumber spectrum, one finds stability if
in the opposite case there is instability of the random wave train. Note that in terms of the BFI the stability condition (40) simply becomes BFI < 1.
As a consequence, one should expect to find in nature spectra with a width larger than the right-hand side of (40), because for smaller width the random version of the Benjamin–Feir instability would occur, resulting in a rapid broadening of the spectral shape. For a random narrowband wave train this broadening is an irreversible process because of phase mixing (Janssen 1983b). The broadening of the spectrum is associated with the generation of inhomogeneities in the wave field. To appreciate this point, it is mentioned that the evolution equation [(36)] satisfies a number of conservation laws. Using the already introduced envelope surface elevation variance ρ(x), the first few conservation laws are given by
assuming periodic boundary conditions in x space and the vanishing of W for large p. The first equation expresses conservation of wave variance, the second one implies conservation of wave momentum, and the last one is the most interesting one in the present discussion because it relates the rate of change of spectral width to the inhomogeneity of the wave field. If the wave field is homogeneous, then ρ(x) is independent of x and the second integral in (41c) is then, because of the first conservation law, independent of time. Therefore, for a homogeneous wave field there is, as expected, no change in spectral width with time; only inhomogeneities will give rise to spectral change according to lowest-order inhomogeneous theory of wave–wave interactions.
Note that the first two conservation laws of (41) may also be obtained immediately from the ensemble average of (12), and the last conservation law follows from the expression of the free-wave energy given in (9) by performing ensemble averaging and by invoking the narrowband approximation. Let us give some of the details of this last derivation. Thus, in the first term the angular frequency is expanded around the carrier wave number k0 up to second order; in the second term the interaction coefficient is replaced by its value at k0. For one-dimensional propagation, we therefore get
Now, the first two terms are already conserved because of conservation of action and momentum, so I will omit them. Performing ensemble averaging while invoking the assumption of a Gaussian state that is, (16) with D = 0, and renaming of the integration variable gives
To summarize the present discussion, note that the central role of the BFI is immediately evident in the context of the lowest-order inhomogeneous theory of wave–wave interactions. According to the stability criterion (40) there is change of stability for BFI = 1. In other words, BFI is a bifurcation parameter: on the short timescale, spectra will be stable and therefore do not change if BFI < 1 while in the opposite case inhomogeneities will be generated, giving rise to a broadening of the spectrum. However, this prediction follows from an approximate theory that neglects deviations from normality. In general, considerable deviations from normality are to be expected, in particular in the case of Benjamin–Feir instability. It is therefore of interest to explore the consequences of nonnormality. This will be done in the next section by means of a numerical simulation of an ensemble of surface gravity waves.
3. Numerical simulation of an ensemble of waves
It is important to determine the range of validity of both the homogeneous and inhomogeneous theories of four-wave interactions. Both theories assume that the wave steepness is sufficiently small, the homogeneous theory ignores effects of inhomogeneity, and the inhomogeneous theory assumes that deviations from normality are small. To address these questions, I simulate the evolution of an ensemble of waves by running a deterministic model with random initial conditions. Only wave propagation in one dimension will be considered from now on.
For given wave number spectrum F(k), which is related to the action density spectrum through F = ωN/g, initial conditions for the amplitude and phase of the waves are drawn from a Gaussian probability distribution of the surface elevation. The phase of the wave components is then random between 0 and 2π while the amplitude should also be drawn from a probability distribution (Komen et al. 1994). Regarding each wave component as independent, narrowband wave trains, a Rayleigh distribution seems to be appropriate for the amplitude (Tucker et al. 1984). However, because the surface elevation is only determined by a finite number of waves, extreme states are not well represented. As a consequence, the kurtosis of the pdf is underestimated. For example, for linear waves it was checked that, even with 51 wave components and a wavenumber resolution of 0.2σk, the kurtosis was underestimated by more than 5%. The size of the ensemble was varied between 500 and 5000. On the other hand, drawing random phases but choosing the amplitudes of the waves in a deterministic fashion, as is common practice, gave an underestimation of kurtosis by only 0.1%. Because the main interest is in the proper representation of extreme events, and because computer resources are limited, it was therefore decided to take only the initial phase as random variable; hence,
where θ(k) is a random phase = 2πxr, xr is a random number between 0 and 1, and Δk is the resolution in wavenumber space.
Each member of the ensemble is integrated for a long enough time to reach equilibrium conditions, typically on the order of 60 dominant wave periods. At every time step of interest, the ensemble average of quantities such as the correlation function B, the pdf of the surface elevation, and integral parameters such as wave height, spectral width, and kurtosis is taken. Typically, the size of the ensemble Nens is 500 members. This choice was made to ensure that quantities such as the wave spectrum were sufficiently smooth and that the statistical scatter in the spectra, which is inversely proportional to Nens, is small enough to give statistically significant results. I now apply this Monte Carlo approach to the nonlinear Schrödinger equation and to the Zakharov equation.
a. Nonlinear transfer according to the nonlinear Schrödinger equation
As a starting point I choose the Zakharov equation [(11)] with transfer coefficients and dispersion relation appropriate for the nonlinear Schrödinger equation. The action variable is written as a sum of δ functions
I have solved this set of differential equations with a Runge–Kutta method with variable time step. Relative and absolute error of the solution have been chosen in such a way that conserved quantities such as action, wave momentum, and wave energy are conserved to at least five significant digits.
In the case of the nonlinear Schrödinger equation, I expand the angular frequency around the carrier wavenumber k0 up to second order. Again using the difference wave number p = k − k0, I find
and I eliminate the contribution of the first term by a change of variable of the form a = a′ exp(−iω0t) and dropping of the prime while the contribution of the second term is removed by transforming to a frame moving with the group velocity at k0. Furthermore, the interaction coefficient T is replaced by its value at k0. As a result we obtain
Before results on the evolution of the spectral properties of the system (45)–(46) are presented, I mention that the nonlinear Schrödinger equation admits a straightforward scaling relation. To see this, let us remove the dependence of the initial condition on the variance 〈η2〉 and the width σk by introducing dimensionless variables
where ε is the wave steepness defined below (13) and c0 is the phase speed corresponding to the carrier wavenumber, c0 = ω0/k0. Writing the nonlinear Schrödinger equation in terms of these new variables, it is immediately evident that for large times its solution can only depend on a single parameter, namely, k0ε/σk, which apart from a constant is just the BFI as defined in (13).
Initial results obtained from the ensemble average of Monte Carlo forecasting did not show the simple scaling behavior with respect to the BFI, until it was realized that results should only be compared for the same dimensionless time t′, which depends in a sensitive manner on the spectral width σk. I therefore integrated the system of equations in (45) until a fixed dimensionless time t′ = 15. A spectral width σk = 0.2k0 was chosen and without loss of generality the carrier wavenumber k0 = 1 was taken. The integration interval then corresponds to about 60 wave peak periods. Furthermore, the resolution in wavenumber space was taken as Δk = σk/3 while the total number of wave components was 41, therefore covering a wide range in wavenumber space. As already noted, this choice gave for linear waves a reasonable simulation of the pdf of the surface elevation.
Note that the specification of a random initial phase has important consequences for the evolution of a narrowband wave train. This is immediately evident when one compares in Fig. 1 time series for the surface elevation from a run with a fixed initial phase θ(k) = 0 with results from a run with a random choice of the initial phase. With a deterministic choice of initial phase, the nonlinear Schrödinger equation generates in an almost periodic fashion extreme events (Fig. 1a), satisfying the criteria for freak waves; with a random choice of initial phase (Fig. 1b), this is much less evident. Comparing the timeseries from the two cases in detail, it is clear that for fixed phase small waves and large waves occur more frequently than in the random-phase case. This impression is confirmed by the pdf of the surface elevation shown in Fig. 2. For reference we have also shown the Gaussian probability distribution. In both cases there are considerable deviations from normality, but in particular for deterministic phase the deviations are large. Similar deviations from the normal distribution were found by Janssen and Komen (1982). Their approach was entirely analytical, and they started from the assumption that for large time the solution of the nonlinear Schrödinger equation would evolve toward a series of envelope solitons, described by an elliptic function. Although they only considered the pdf of the envelope (which under normal conditions is given by the Rayleigh distribution), one may obtain the pdf of the surface elevation as well. The resulting analytical pdf has characteristics similar to the pdf for the case of deterministic phase.
The Monte Carlo approach was adopted because it is not evident that for the system under discussion the ergodic hypothesis applies. This hypothesis implies replacement of the ensemble average by a time average. However, if one happens to choose initial phases in a way that is favorable for the generation of envelope solitons, then there is a high probability that the solution stays close to the envelope soliton branch and will hardly ever visit other parts of phase space. To guarantee a representative picture, we therefore decided to perform Nens runs, where for each run amplitude and phase are drawn in an independent manner. In the remainder, only ensemble-averaged results will be discussed. These are presented graphically in the figures. To limit the number of plots, each figure shows results from the numerical simulations [usually denoted by Monte Carlo forecasting of waves (MCFW)], and corresponding results from the stochastic approach (labeled “theory”). The theoretical results are discussed in section 4.
In Fig. 3 I show the evolution of the spectral width σk with dimensionless time t′ for several values of BFI. Here, σk is defined using integrals of the wave spectrum F over wavenumber p:
Note that for simulations with the nonlinear Schrödinger equation this turned out to give a remarkably stable estimate of the width of the spectral peak, because the spectra vanish sufficiently rapidly for large p. According to this simulation, there is a considerable broadening of the spectrum, which occurs on a fairly short timescale of about 10 peak wave periods. In this case of one-dimensional propagation, the standard theory of nonlinear transfer would give no spectral change. Note that σk shows in the early stages of time evolution an overshoot followed by a rapid transition toward an equilibrium value. The number of oscillations around this equilibrium value depends on the precise details of the discretization scheme. In particular, more, larger-amplitude oscillations are found for coarser spectral resolution. The overshoot is in agreement with results of Janssen (1983b), who studied the evolution of a single unstable mode in the context of inhomogeneous theory of wave–wave interactions. For sufficiently narrow spectra, overshoot in the amplitude of the unstable mode was found followed by a damped oscillation toward an equilibrium value. The damping timescale was found to depend on the width of the spectrum, vanishing for small width. In physical terms, the damping is caused by phase mixing (or destructive interference), and its effect depends on wavenumber resolution. Note that other parameters, such as the kurtosis, show a similar time behavior.
As an example of spectral evolution, Fig. 4 shows for BFI = 1.40 the initial and final time–wavenumber spectrum. To give an idea about the representativeness of the results, 95% confidence limits, based on Nens − 1 degrees of freedom, are shown as well. The broadening of the spectrum as caused by the nonlinear interactions is statistically significant. Although the spectral change should be symmetrical with respect to the carrier wavenumber, that is, p = 0, it is clear that there are asymmetries present in the ensemble average of the numerical results. However, these deviations are within the statistical uncertainty. To make sure of this, I redid the case for BFI = 1.40 but now with an ensemble size of 2000. As expected, statistical uncertainty was reduced by a factor of 2 while asymmetries were reduced as well.
To examine whether the Monte Carlo results show evidence of a bifurcation at BFI = 1, I plot in Fig. 5 the relative increase in spectral width, defined as [σk(t′∞) − σk(0)]/σk(0), as a function of the BFI evaluated with the initial value for spectral width. The results suggest that there is only evolution of the spectrum for sufficiently large BFI, but in contrast to inhomogeneous theory of wave–wave interactions, BFI = 1 does not appear to be a bifurcation point, given that considerable changes in the wave spectrum already start to occur for BFI = 1/2. Although from inhomogeneous theory one would expect a sudden transition from no spectral change to spectral change, Fig. 5 seems to suggest that the transition is gradual. I attribute this discrepancy to the assumption in inhomogeneous theory that deviations from normality are small, because these may give rise to irreversible changes of the spectrum as well. This will be discussed more extensively in the next section. It is illuminating to plot the information on spectral width in a slightly different manner, namely, by relating the final time value of BFI with its initial value. This is done in Fig. 6, and it clearly shows that for large times BFI hardly exceeds the value of 1. This seems to agree with the conjecture given in Section 2b that according to inhomogeneous theory [(40)] one should not expect spectra to have a BFI much larger than 1. According to the Monte Carlo results (Fig. 3) the timescale of change for large BFI is, on average, only a few wave periods.
Nonlinear effects give rise to considerable changes in the probability distribution of the surface elevation from the Gaussian distribution (see also Onorato et al. 2000; Mori and Yasuda 2002a), although the deviations are of course much less than in the cases discussed in Fig. 2. This is shown in Fig. 7 for BFI = 1.40. To emphasize the occurrence of extreme events, I have plotted the logarithm of the pdf as function of the surface elevation normalized with the wave variance. The Gaussian distribution then corresponds to a parabola. The simulated pdf, in the range of 2–4, shows an almost linear behavior that suggests an exponential decay of the pdf. In Fig. 8 I summarize the results on the deviations from normality by plotting the final time value of the kurtosis C4 = 〈η4〉/3m20 − 1 as a function of the final time BFI. Here, the fourth moment 〈η4〉 was determined from the pdf of the surface elevation, which was obtained by sampling the second half of the time series for the surface elevation at an arbitrarily chosen location. Alternatively, the fourth moment may be obtained from (28), giving very similar answers. For small nonlinearity one would expect a vanishing kurtosis, but the simulation still underestimates, as already mentioned, the kurtosis by 2%. The kurtosis depends almost quadratically on the Benjamin–Feir index up to a value close to 1. This quadratic dependence will be explained in the next section, when an interpretation of results is provided. Near BFI = 1, on the other hand, the kurtosis behaves in a more singular fashion, because, in agreement with the discussion of Fig. 6, the Benjamin–Feir index cannot pass the barrier near 1.
The nonlinear Schrödinger equation [(45)] for deep-water waves is an example in which nonlinearity leads to focusing of wave energy and therefore counteracts the dispersion by the linear term, which is proportional to ω″0. The results from the numerical simulation do indeed suggest that, when nonlinearity is sufficiently strong, focusing of energy occurs, giving considerable enhancements to the probability of extreme events, at least when compared with the normal distribution. In the opposite case, when the nonlinear term has opposite sign, defocusing of wave energy occurs and one would expect a reduction in the number of extreme events. To show this, I performed simulations with the nonlinear Schrödinger equation [(45)] but now with negative nonlinear transfer coefficient (T0 = −k30). Results of this case are shown in Figs. 5, 6, and 8; the logarithm of the pdf of the surface elevation is shown in Fig. 9. These plots show that in the case of defocusing the broadening of the spectrum is less dramatic. Furthermore, the final time Benjamin–Feir index does not have a limiting value of about 1. On the other hand, the kurtosis for this case is negative, resulting, as can be seen from Fig. 9, in a large reduction of the probability of extreme events. The dependence of the kurtosis on BFI is different from the case of focusing, because there are clear signs of saturation beyond BFI = 1, whereas only in the range BFI < 0.5 is there a quadratic dependence of kurtosis C4 on BFI.
b. Nonlinear transfer according to the Zakharov equation
The nonlinear Schrödinger equation gives the lowest-order effects of finite bandwidth on the evolution of a weakly nonlinear wave train. Dysthe (1979) investigated the consequences of next order in bandwidth and he found a surprisingly large impact on the results for the growth rates of the modulational instability. Similarly, Crawford et al. (1981) studied the stability of a uniform wave train using the complete Zakharov equation, which retains all of the high-order dispersion effects. In general, growth rates are reduced when compared with results from the nonlinear Schrödinger equation; therefore, according to the Zakharov and Dysthe equations, a uniform wave train is less unstable. In fact, growth rates and thresholds for instability were in better agreement with experimental results of Benjamin and Feir (1967) and Lake et al. (1977) (see also Janssen 1983a). The Zakharov and Dysthe equations have, in addition, the interesting property that the nonlinear transfer coefficient and the angular frequencies are not symmetrical with respect to the carrier wavenumber. It will be seen that this has important consequences for the spectral shape.
The Dysthe equation follows from the Zakharov equation by expanding angular frequency to third order in the modulation wave number p while the interaction coefficient T is expanded up to first order in p. For narrowband wave trains it gives an accurate description of the sea state. However, wave spectra may become so broad that the narrowband approximation becomes invalid, and therefore I have chosen to study numerical results from the Zakharov equation.
The Zakharov equation I solved was given by (44), where the nonlinear transfer coefficient was from Krasitskii (1990), and the exact dispersion relation for deep-water gravity waves was taken. The initial condition was provided by (46). The discretization details were identical to those of the numerical simulations with the nonlinear Schrödinger equation. Because the Zakharov equation contains all higher-order terms in the modulation wave number p, it is not possible to prove that the large-time solution of the initial value problem is determined completely by the BFI, but in good approximation the BFI can still be used for this purpose as long as the spectra are narrow-banded.
In Fig. 10 is plotted the ensemble-averaged wavenumber spectrum for BFI = 1.4, and it shows a clear downshift of the peak of the spectrum while also considerable amounts of energy have been pumped into the high-wavenumber part of the spectrum. The wavenumber downshift is caused by the asymmetries in the nonlinear transfer coefficient and to the same extent by the asymmetries in the angular frequency with respect to the carrier wavenumber. This was checked by running (44) with constant nonlinear transfer coefficient, and similar-looking ensemble mean spectra, but with one-half of the wavenumber downshift, were obtained. There is also a noticable broadening of the spectrum. However, because of the increased spectral levels at high wavenumbers, use of the second moment of the wave number spectrum, as was done for the nonlinear Schrödinger equation [(48)], to measure the width of the spectral peak is not appropriate. Instead I use here the width as obtained from fitting the peak of the spectrum with a Gaussian shape function.
The relation between the final time BFI versus the initial value of BFI is shown in Fig. 11 and is compared with the corresponding one from the nonlinear Schrödinger equation. Also, Fig. 12 shows the normalized kurtosis versus the final time BFI. Results from the Zakharov equation are in qualitative agreement with the ones from the nonlinear Schrödinger equation. However, because growth rates are smaller, the broadening of the spectrum is less, the final time BFI is higher by about 10%, and the normalized kurtosis is smaller. A unique feature of the Zakharov equation is the downshift of the peak of the wavenumber spectrum. This is shown in Fig. 13 where is plotted the final time value of the peak wavenumber, normalized with its initial value versus the initial BFI. For large values of BFI reductions in peak wavenumber of more than 10% are found from the results of the numerical simulations, but the dependence of the downshift in peak wavenumber on BFI is not smooth. This is caused by the fact that the ensemble-averaged spectra do not always have a well-defined spectral peak.
4. Interpretation of numerical results
In the previous section I have discussed results from the Monte Carlo simulation of the nonlinear Schrödinger equation and the Zakharov equation. These results show that on average there is a rapid broadening of the wave spectrum while nonlinearity gives rise to considerable deviations from Gaussian statistics. The question now is whether the average of the Monte Carlo results may be obtained in the framework of a simple theoretical description. In section 2 I discussed two attempts to achieve this. The first one is the standard theory of wave–wave interactions, extended with the effects of nonresonant four-wave interactions. This approach assumes a homogeneous wave field but allows for deviations from the Gaussian sea state. The second theory is the inhomogeneous theory of wave–wave interactions, which assumes that effects of inhomogeneity in the wave field are dominant while deviations from normality only play a minor role. This approach seems to be an ideal starting point for treating inhomogeneous and nonstationary phenomena such as freak waves because it describes the random version of the Benjamin–Feir instability. Let us therefore first discuss the validity of inhomogeneous theory using results from the Monte Carlo forecasting of ocean waves obtained from the nonlinear Schrödinger equation.
According to inhomogeneous theory the broadening of the wave spectrum is caused by the inhomogeneity of the wave field. This is clearly expressed by the conservation law [(41c)] and is explained in the discussion that follows. From this conservation law, one may therefore obtain a measure of inhomogeneity of the wave field, namely,
Here, the integrals over space are weighted by the size of the domain, and using the definitions for ρ [(37)] and the Wigner distribution [(35)] one may express the inhomogeneity measure I in terms of the correlation function B(n, m, t) as
The correlation functions B(p, m) may be readily obtained from the numerical results for the complex amplitude a(k, t) after ensemble averaging. For a homogeneous wave field B(p, m) = N(p)δ(m); hence I2 = I21, or I = 1.
The initial conditions used in the numerical simulation of waves have been chosen in such a way that the sea state corresponds to a Gaussian one. As a consequence, because the complex amplitudes a(k, t) are not correlated, this implies that initially the sea state is homogeneous as well (Komen et al. 1994). However, the wave ensemble consists of a finite number of members, and this means that the initial probability distribution is not a perfect Gaussian (e.g., the kurtosis is slightly underestimated) but it also means that the initial conditions are slightly inhomogeneous. According to the inhomogeneous theory, the perturbations should grow exponentially in time, resulting in, for example, a broadening of the wave spectrum.
For BFI = 1.4, the evolution in time of the inhomogeneity I is shown in Fig. 14. Initially, inhomogeneity is small but grows rapidly in the course of time, which is then followed by an oscillation around the level 1.004. This level of inhomogeneity and the variation with time is, however, extremely small (note that for the case of Fig. 1a, I is on the order of 3), and it cannot explain the large changes in the wavenumber spectrum seen in the numerical simulations.
In addition, according to the inhomogeneous theory, the conservation law of (41c) should be satisfied. In section 3b it was explained that this conservation law follows from the conservation of Hamiltonian, assuming that deviations from normality may be ignored. It is of interest, of course, to test whether it is justified to ignore effects of the fourth cumulant. To that end, Fig. 15 compares the evolution in time of the Hamiltonian as obtained from the numerical simulation (this will be called the “exact” Hamiltonian from now on) with the Hamiltonian according to inhomogeneous theory. Although the exact Hamiltonian is a constant (at least up to five significant digits), it is clear that the approximate Hamiltonian is not conserved when evaluated using the numerical results. In fact, there are large deviations as the approximate Hamiltonian becomes negative, whereas the exact Hamiltonian is positive definite. The disagreement between the approximate and the exact Hamiltonian is caused by the neglect of the higher-order cumulants. This is immediately clear from Fig. 15 in which the nonlinear contribution to the Hamiltonian according to lowest-order inhomogeneous theory [called approximate; NL (c4 = 0)] is compared with the corresponding nonlinear contribution that includes higher-order cumulants (called exact; NL + c4). The approximate nonlinear contribution hardly varies with time, which is in agreement with the results from Fig. 14 that effects of inhomogeneity are small. The exact nonlinear contribution shows a significant variation with time. The differences between exact and approximate are considerable, and therefore it is not justified to ignore effects of deviations from normality in a simple theoretical description of the evolution of the sea state. As a matter of fact, the deviations from normality are the main reason for the spectral broadening as the time-varying nonlinear contribution to the Hamiltonian, including effects of the fourth-order cumulant, just compensates for the changes with time of the linear part of the wave energy. Clearly, according to the Monte Carlo simulations the linear wave energy is not conserved.
The previous discussion therefore suggests that effects of the deviations from normality are important, which may explain why there is no bifurcation at BFI = 1 as is suggested by inhomogeneous theory (Fig. 5). Alternatively, there could be another possible reason for the discrepancy between the inhomogeneous theory and the numerical results. The evidence for a bifurcation at BFI = 1 follows from a linear stability analysis; therefore the perturbation of the homogeneous spectrum is assumed to be small. In the numerical simulations the perturbation is not prescribed a priori and its amplitude is determined by the full nonlinear evolution equation. The initial size of the perturbation is, however, determined in an indirect way by the ensemble size, because the larger the ensemble is, the closer one approaches to conditions for a homogeneous sea state. This is evident in Fig. 14 in which the inhomogeneity I is smaller for larger ensemble size. I therefore performed Monte Carlo simulations with different ensemble size, namely, for Nens = 50 and Nens = 2000, but almost identical results for, for example, the relative increase in spectral width as a function of BFI were found. Although this does not present conclusive evidence, it does suggest that the initial size of the perturbations is small enough and that the main reason for the absence of a bifurcation at BFI = 1 is the neglect of deviations from normality.
In summary, it has been shown that, in the inhomogeneous theory of four-wave interactions, effects of the generation of the fourth cumulant cannot be ignored. At the same time it has been shown that the numerical ensemble of ocean waves may be regarded to good approximation as a homogeneous ensemble. Hence, the standard theory of four-wave interactions (extended by including nonresonant interactions), which assumes a homogeneous wave field, may be a good candidate to explain the results of the numerical simulations in section 3.
Therefore, the Boltzmann equation [(25)] was used to evolve the action density N(k) for the same cases as presented in section 3. The differential equation was solved with a Runge–Kutta method with variable time step, and the continuous problem was discretized in the same way as was done in the case of the solution of the Zakharov equation. Run times using the homogeneous theory are typically two orders of magnitudes faster than when following the ensemble approach.
In contrast to inhomogeneous theory, the standard theory gives a much better approximation to the exact Hamiltonian as shown in Fig. 15. There are, as should be, small differences because the standard theory is an approximation as well, given that both effects of inhomogeneity and the sixth cumulant have been neglected. Further results from the discretized version of the homogeneous theory are compared with the ones from the simulations with the nonlinear Schrödinger equation in Figs. 3–8.
From Fig. 3, which shows the evolution in time of the spectral width for several values of the BFI, it is seen that for large times there is good agreement between homogeneous theory and the ensemble-averaged results from the Monte Carlo simulations. For short times it is, however, evident that (25) does not show the overshoot found in the numerical simulations. A likely reason for the absence of overshoot found in the theoretical calculations is the assumption that the action density varies slowly as compared with the timescale implied by the resonance function Ri(Δω, t). Both the numerical simulations and homogeneous theory show for short times a rapid broadening of the wave spectrum, which for large times is followed by a transition towards a steady state. The evolution toward a steady state can be understood as follows: First, it should be noted that, according to section 2c, for one-dimensional propagation there is no nonlinear transfer because of resonant nonlinear interactions. Now, initially the resonance function Ri(Δω, t) will be wide so that nonresonant wave–wave interactions are allowed to modify the action density spectrum. After about 5–10 wave periods, however, the resonance function becomes progressively narrower until it becomes approximately a δ function; hence only resonant waves are selected. In that event there is no change of the action density spectrum possible anymore so that for large times a steady state is achieved.
An example of the comparison between theoretical and simulated spectra is given in Fig. 4. There is a fair agreement between the two. However, it should be mentioned that typically the simulated spectrum is slightly more peaked than the theoretical one despite the fact that there is a close agreement in spectral width. This agreement in spectral width between theory and simulation is also very much evident in the Figs. 5 and 6 over the full range of the initial value of BFI. In particular, note that there is close agreement between the upper limit of the final time BFI from theory and the simulation. Hence, homogeneous theory also provides an explanation of why there is a lower bound to spectral width.
We therefore have the curious situation in which both homogeneous and inhomogeneous theory explain why there is a lower bound to spectral width as found in the numerical simulations. However, because it has been shown that inhomogeneities only play a minor role in the numerical simulations it follows that only homogeneous theory provides a proper explanation. In situ observations from the North Sea seem to indicate the presence of a lower bound to spectral width as well (see, e.g., Janssen 1991). Although it is impossible to prove at present that at sea inhomogeneities do not play a role, homogeneous theory even seems to give a plausible explanation of the lower bound found at sea.
As discussed in section 2c, nonlinearity gives rise to deviations from the normal distribution. I determined the normalized kurtosis using (29), which is obtained from the fourth cumulant D. Introducing the normalized height x = η/m0, the pdf of the normalized surface elevation x is then given by (see also Mori and Yasuda 2002b)
where f0 is given by the normal distribution
Equation (54) follows from an expansion of the pdf p in terms of orthogonal functions (d/dx)nf0. Here, n is even because of the symmetry of the Zakharov equation. The expansion coefficients are then obtained by determining the first, second, and fourth moments. For the range of BFI studied in this paper, it was verified that higher moments only gave a small contribution to the shape of the pdf p(x). The pdf according to theory is compared in Fig. 7 with the simulated one, and a good agreement is obtained, even for extreme sea state conditions. Clearly, in the case of nonlinear focusing, the probability of extreme states is, as expected, larger when compared with the normal distribution. In Fig. 8, theoretical and simulated final-time kurtosis is plotted as function of the final-time BFI. A good agreement between the two results is obtained, even close to the limiting value of the final-time BFI. For BFI < 1, both simulated and theoretical kurtosis depend in an almost quadratic fashion on BFI, in agreement with the simple estimates of C4 given in section 2c [(31)].
In the case of nonlinear focusing, a good agreement between the numerical simulations and the homogeneous theory has been obtained, even for extreme values of the Benjamin–Feir index. Here, it should be emphasized that at sea BFI has typical values of 0.5 or less, and only occasionally are values on the order of 1 reached. We have performed simulations to values of BFI of up to 3 and, even for these extreme conditions (having large values of kurtosis, for example), a reasonable agreement is obtained. This is surprising because the homogeneous theory was derived under the assumption of small deviations from normality.
In the case of nonlinear defocusing, the range of validity of homogeneous theory is much more restricted. This is made plainly clear in Figs. 5, 6, 8, and 9, where results from the simulations and homogeneous theory are compared for the case of nonlinear defocusing. To be able to interpret this comparison, note that homogeneous theory does not distinguish between focusing and defocusing, because the nonlinear transfer is independent of the sign of the interaction coefficient T. Only the kurtosis depends on the sign of T. Judging from Figs. 5, 6, and 8, the range of validity of homogeneous theory is restricted to BFI < 0.6. In particular, for large BFI there are large qualitative deviations in the kurtosis of the pdf: from the numerical simulations there are clear signs of saturation in kurtosis, but there is no indication of saturation in the results from homogeneous theory. In addition, in the case of nonlinear defocusing the kurtosis is negative so that for large normalized elevation x the pdf given in (54) may become negative. This is clearly unrealistic, and to correct this undesirable property of homogeneous theory one needs, for large BFI, to take into account the effects of higher-than-fourth-order cumulants as well. It is believed this is the main reason why homogeneous theory has such a restricted validity in this case. In the opposite case of nonlinear focusing, the kurtosis is positive, giving for large x a positive correction to the normal distribution. The pdf of the surface elevation is therefore positive, at least for the cases that have been studied here, and as a consequence homogeneous theory has a much larger range of validity.
Last, I applied homogeneous theory to the Zakharov equation. Results are compared with the numerical simulations in Figs. 10, 11, 12, and 13. There is a fair agreement between simulated and theoretical spectrum (Fig. 10), between simulated and theoretical final-time BFI index (Fig. 11), and in normalized kurtosis (Fig. 12). Less favorable is the agreement between simulated and theoretical peak wavenumber, as shown in Fig. 13. The theoretical results show a smooth dependence of wavenumber downshift on BFI, giving shifts of 20% or more for large values of BFI. However, the simulation shows more scatter while the downshift is at most 15%. The reason for the scatter in the simulated results is probably that the peak of the wavenumber spectrum is not always well defined. In contrast, homogeneous theory gives a smoother spectrum and a well-defined peak of the spectrum.
Present-day wave forecasting systems are based on a description of the ensemble-averaged sea state. In this approach, the wavenumber spectrum plays a central role and its evolution follows from the energy balance equation. The question discussed here is whether it is possible to make statements, necessarily of a statistical nature, on the occurrence of extreme events such as freak waves.
To show that this is possible, the following approach has been adopted. The starting point is a deterministic set of equations, namely, the Zakharov equation or its narrowband limit, the nonlinear Schrödinger equation. There is ample evidence that these equations admit freak wave–type solutions. These freak waves occur when the waves are sufficiently steep, because nonlinear focusing may then overcome the spreading of energy by linear dispersion. For the same reason, the Benjamin–Feir instability occurs. As shown in Fig. 1, the occurrence of freak waves depends in a sensitive manner on the choice of the initial phase of the waves. In addition, on the open ocean, waves propagate from different locations toward a certain point of interest and may therefore be regarded as independent. Hence, for open-ocean applications, the random-phase approximation seems to be appropriate. I therefore simulated the evolution of an ensemble of ocean waves by running a deterministic model with random initial phase. These Monte Carlo simulations are expensive (typically, the size of the ensemble is 500) so that the study is restricted to the case of one-dimensional propagation only.
The ensemble average of the results from the Monte Carlo simulations shows that when the Benjamin–Feir index is sufficiently large (as occurs for the combination of narrow spectra and steep waves) the wave spectrum broadens while at the same time considerable deviations from the Gaussian pdf are found. In the case of the Zakharov equation, the spectral change is even asymmetrical with respect to the peak of a symmetrical spectrum, giving a downshift of the peak wavenumber while as a consequence of the conservation of wave action and wave energy considerable amounts of energy are being pumped into the high-wavenumber part of the spectrum. These spectral changes occur on a short timescale, typically on the order of 10 wave periods. This timescale is comparable to the Benjamin–Feir timescale.
The standard homogeneous theory of resonant nonlinear transfer does not give spectral change in the case of one-dimensional propagation. This theory was therefore extended to allow for nonresonant interactions as well, because the nonlinear focusing related to the Benjamin-Feir instability is an example of a nonresonant four-wave interaction. This nonlinear four-wave transfer is associated with the generation of higher-order cumulants such as the kurtosis. Deviations of the surface elevation pdf from the normal distribution can therefore be expressed in terms of a six-dimensional integral involving the cube of the action spectrum [(29)]. In the case of nonlinear focusing, there is, for a large range of values of the BFI, a good agreement between the ensemble-averaged results from the numerical simulations and homogeneous theory. This is in particular true for the broadening of the spectrum, the spectral shape, and the estimation of the kurtosis. In comparison with the simulations, theory overestimates, however, the peak wavenumber downshift.
Homogeneous theory also explains why for one-dimensional propagation the wave spectrum evolves toward a steady state. The resonance function evolves rapidly toward a δ function; hence for large times only resonant wave–wave interactions are selected that in one dimension do not give rise to spectral change. This is in sharp contrast with the case of two-dimensional propagation. No trend toward a steady state is expected in that case because resonant four-wave interactions do contribute to a change in the spectrum.
The extended version of the homogeneous four-wave theory has two timescales, a fast one on which the nonresonant interactions take place and a long timescale on which the resonant interactions occur. The nonresonant interactions play a role similar to that of transients in the solution of an initial-value problem. They are simply generated because initially there is a mismatch between the choice of the probability distribution of the waves, a Gaussian, and between the initial choice of the wave spectrum, representing a sea state with narrowband, steep nonlinear waves. For example, if one could choose a probability distribution function which is in equilibrium with the nonlinear sea state (theoretically one can, by the way), then nonresonant interactions would not contribute. Only resonant wave–wave interactions will then give rise to nonlinear transfer. In the general case for which there is a finite mismatch between pdf and wave spectrum, the nonresonant contribution will die out very quickly owing to phase mixing but will, nevertheless, as we have seen, result in considerable changes in the wave spectrum. The question therefore is whether there is a need to include effects of nonresonant interactions. This depends on the application. In wave-tank experiments, where one can program a wave maker to produce the initial conditions used in this paper, it seems that effects of nonresonant transfer need to be taken account for. For the open-ocean case this is not clear. The point is that in nature the combination of steep waves and a strictly Gaussian distribution most likely does not occur. Changes in nature are expected to be more gradual so that the mismatch between pdf and wave spectrum is small. Only when a wind starts blowing suddenly, hence for short fetches and duration, are effects of nonresonant interactions expected to be relevant. More research in this direction is, however, required.
The results from Monte Carlo simulations do not provide evidence that there are significant deviations from homogeneity of the ensemble of waves. Deviations from normality are found to be much more important. Nevertheless, it cannot be concluded from the present study that the inhomogeneous approach of Alber and Saffman (1978) is not relevant for real ocean waves. For example, effects of inhomogeneity might be relevant in fetch-limited, rapidly varying circumstances. However, its effects are expected to be small, and therefore only the lowest-order approximation, as given explicitly in (36), needs to be retained.
For real ocean waves, not only four-wave interactions determine the evolution of the wave spectrum. Wind input and dissipation due to white-capping are relevant, and these processes may affect the kurtosis of the sea surface as well. However, it has been shown in this paper that in particular the nonresonant nonlinear transfer acts on a short timescale that is much shorter then the timescales associated with wind input and dissipation source functions used in wave forecasting (Komen et al. 1994). Hence, one would expect that the expression for kurtosis found in this paper [(29)] should be relevant in nature. Nevertheless, nonlinear focusing may result in steep waves. If their steepness exceeds a certain threshold, one would expect a significant amount of wave breaking, thus limiting the height of these waves, and affecting the extreme statistics. A realistic, deterministic model of wave breaking is needed to assess the importance of wave breaking in these circumstances. It may be more effective, however, to try to compare results from the present approach directly with observations of extremes collected over a long period.
The author acknowledges useful and stimulating discussions with Miguel Onorato. His presentation during the 2001 Wise meeting in Canada gave stimulus to have a fresh look at wave–wave interactions and freak waves. Furthermore, the author thanks Theo van Steyn for providing the Runge–Kutta software for accurate integration of a set of ordinary differential equations. Last, the author thanks Anthony Hollingsworth and Martin Miller for critically reading the manuscript.
Corresponding author address: Peter A. E. M. Janssen, European Centre for Medium-Range Weather Forecasts, Shinfield Park, Reading RG2 9AX, United Kingdom. Email: firstname.lastname@example.org