{'fraction_non_alphanumeric': 0.03879858037944179, 'fraction_numerical': 0.022170617359681025, 'mean_word_length': 4.456251494142959, 'pattern_counts': {'":': 0, '<': 10, '<?xml version=': 0, '>': 4, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Systematic Study of Rogue Wave Probability Distributions in a Fourth-Order Nonlinear Schr\"odinger Equation

Nonlinear instability and refraction by ocean currents are both important mechanisms that go beyond the Rayleigh approximation and may be responsible for the formation of freak waves. In this paper, we quantitatively study nonlinear effects on the evolution of surface gravity waves on the ocean, to explore systematically the effects of various input parameters on the probability of freak wave formation. The fourth-order current-modified nonlinear Schr\"odinger equation (CNLS4) is employed to describe the wave evolution. By solving CNLS4 numerically, we are able to obtain quantitative predictions for the wave height distribution as a function of key environmental conditions such as average steepness, angular spread, and frequency spread of the local sea state. Additionally, we explore the spatial dependence of the wave height distribution, associated with the buildup of nonlinear development.


Introduction
Oceanic rouge waves, or freak waves, are surface gravity waves whose wave heights are extremely large compared to the typical wave in a given sea state. Freak waves have been reported throughout maritime history, as they are a hazard to mariners, cargo ships, and even to large cruise liners. However, serious scientific investigation commenced more recently, due to the availability of wave state measurement methods including satellite remote sensing, and the introduction of stochastic process theory to ocean wave forecasting [Bates, 1952;Kinsman, 1965].
The scientific community has studied the topic experimentally (in water tanks) [Onorato et al., 2004] and [Onorato et al., 2006], observationally (with satellite imaging, for example) [Dankert et al., 2003] and [Schulz-Stellenfleth and Lehner , 2004], and theoretically [Kharif and Pelinovsky, 2003]. A review of recent progress in freak wave forecasting can be found in an article by Dysthe et al. [2008]. Given the chaoticity of ocean wave dynamics, it is not possible to predict individual freak wave events. Instead, we study wave behavior probabilistically. The simplest theory assumes a random superposition of a large number of monochromatic waves with different frequencies and propagating directions, as in the Longuet-Higgins random seas model [Longuet-Higgins, 1957]. By the central limit theorem, this model leads to a Rayleigh probability distribution of wave heights, where 2H is the wave height, σ 2 is the variance of the surface elevation, and √ 2πσ is the mean wave height. For linear waves, the wave height 2H (crest to trough) equals twice the crest height due to a symmetry of the wave equation. This symmetry is broken for nonlinear waves, and throughout the paper, we focus on the wave height 2H, where the crest height is H to leading order.
Conventionally, freak waves are defined as having wave height 2H ≥ 8.8σ, or alternatively 2H ≥ 2.2Hs, where Hs is the significant wave height [Dean, 1990]. The significant wave height Hs was traditionally defined as the average of the tallest one-third of waves in a wave train, which comes to Hs ≈ 4.004σ when the height distribution is Rayleigh. Since the Rayleigh distribution (1) serves as an excellent zerothorder approximation for realistic wave height distributions in the ocean, Hs = 4σ may be used today as another definition of the significant wave height [Forristall , 1978]. Similarly, extreme freak waves are defined as having wave height 2H ≥ 12.0σ or 2H ≥ 3Hs. According to the simple Rayleigh distribution (1), the occurrence probabilities of a freak wave and an extreme freak wave are 6.3 × 10 −5 and 1.5 × 10 −8 , respectively. Compared with the observational data [Dankert et al., 2003], the stochastic Rayleigh model greatly underestimates the probability of freak waves. Recent advances in technology have allowed multiple wave tank experiments and field observations to be conducted, confirming the need for a more realistic theory to explain the results [Forristall , 2000;Onorato et al., 2004].
Nonlinear instability and refraction by ocean currents are both important mechanisms that go beyond the Rayleigh approximation and may be responsible for the formation of freak waves.
The nonlinear instability, also known as the Benjamin-Feir instability in the one-dimensional case [Benjamin and Feir , 1967], describes how a monochromatic wave can be unstable under a class of small perturbations. In the twodimensional case, however, the spread in the wave vector happens not only in magnitude but also in wave direction. The nonlinear evolution involves both a nonlinear selffocusing effect at short times, which leads to an increase in the occurrence rate of extreme events, and also a dispersion in wave vector on larger time scales, which leads to a decrease in the probability of extreme events as the initial spectrum becomes less unidirectional [Onorato et al., 2002].
We usually use the steepness to describe the nonlinearity of a water wave. Steepness is defined as ε = 2Hk0/2, where 2H is the mean wave height and k0 is the peak wave number. The nonlinear evolution of a surface-gravity water wave is believed to be well described by the nonlinear Schrödinger equation when the steepness is small, 2Hk0 ≪ 1, and the wave bandwidth is narrow, |∆ k|/k0 ∼ ε ≪ 1. This is supported by Trulsen and Dysthe's analysis of the famous New Year's wave event, which occurred at the Draupner platform in the North Sea in 1995 [Trulsen and Dysthe, 1997].
Clearly the exact evolution is always nonlinear to some extent, but the key is to introduce nonlinearities at the right moment and in an insightful and computable way. Realistic fully nonlinear computations wave by wave over large areas are very challenging, but initial attempts have been made to simulate the ocean surface using the full Euler equation both on large scales [Tanaka, 2001] and over smaller areas [Gibson and Taylor , 2005;Gibson and Swan, 2007]. However, in this paper, we rely on series expansions, resulting in a sequence of approximations. Pioneering work using a third-order equation in ε (the original Nonlinear Schrödinger Equation (NLS)) was performed by Zakharov [1968], and then extended to fourth order by Dysthe [1979] (for water waves propagating in an ocean of infinite depth). In nature, the small steepness constraint is usually satisfied, since the steepness under storm conditions is normally less than 0.09 [Dysthe et al., 2008]. However, the narrow bandwidth constraint is often not complied with in naturally occurring conditions. To overcome this problem, a modified nonlinear Schrödinger equation (MNLS) has been obtained by Trulsen and Dysthe [1996] for broader bandwidth, only requiring |∆ k|/k0 ∼ √ ε, ε ≪ 1.
The effect of currents interacting with an incoming sea state has been studied by Heller et al. [2008] using ray dynamics and by Janssen and Herbers [2009] using Monte Carlo simulations. These studies suggest that the interaction between incoming wave and current may serve as a triggering mechanism for the formation of freak waves during nonlinear wave evolution. In the work by Heller et al. [2008], a freak index is defined in terms of the mean wave speed, mean current speed, and the angular spread of the incoming wave, and a quantitative relationship is predicted between this freak index and the occurrence probability of freak waves. However, in this work ray dynamics is used in place of the real wave equation for ocean waves, and nonlinearity is not included in the model. The calculations done by Janssen and Herbers [2009] do include nonlinearity, and serve as a proof of principle for the importance of the interaction between nonlinearity and refraction by currents. Still, results are obtained only for specific values of the input parameters. Thus, a systematic quantitative analysis of the freak wave probability as a function of the strength of nonlinear and refractive effects has not previously been undertaken.
Here, we use the current-modified nonlinear Schrödinger equation (CNLS) developed by Stocker and Peregrine [1999] to describe the evolution of the wave envelope, with the aim of integrating both nonlinear effects and wave-current interaction in one model. This approach allows for a full quantitative analysis of both these mechanisms of freak wave formation. In the present paper we focus on numerically exploring the dependence of the wave height distribution on nonlinear effects only, by setting the currents to zero. We will see that interesting behavior is obtained already in the current-free regime. The joint dependence on steepness and current strength for nonzero random currents will be presented in a future publication. This paper is organized as follows: In section 2 we review the CNLS equations and present the basic computational setup. In section 3 we discuss the expected form of the wave height distribution in the presence of deviations from Rayleigh statistics due to nonlinearity and currents.
In section 4 we systematically study the dependence of the wave height distribution on the parameters describing the incoming sea state, allowing the probability of freak wave formation to be calculated as a function of these parameters. Finally, in section 5 we summarize our results and discuss the outlook for the future.

Model
The current-modified nonlinear Schrödinger equation equation (CNLS) extends the NLS to include a random surface current varying around a mean value u0. The timeindependent current field can be expressed as: where the random current velocity fluctuations δ u are assumed to be O(ε 2 ), and slowly varying on the scale of a wavelength. For surface gravity waves in deep water (i.e., when the water depth is much larger than the wavelength), the dispersion relation is given by: where ω is the angular frequency of the incoming wave whose wavelength is given by λ = 2π/k, and the wave velocity is ∇ k ω = g/4kk + u. For convenience, we take the peak wave vector of the incoming wave to be in the positive x direction. The wave vector of each wave component is then expressed as: where k0 > 0.
To simplify the equations and simulations, we work in the frame of reference moving with the peak wave velocity g/4k0x + u0, and study the evolution of the wave envelope instead of the wave function itself. The complex wave envelope a( r, t) is defined by where ζ( r, t) is the surface elevation. For a narrow-banded spectrum with small angular spread (i.e., |∆ k|/k0 ∼ √ ε ≪ 1), the wave envelope varies on a scale long compared to the wavelength. In that regime, the probability distribution of the wave envelope a( r, t) closely approximates the crest height distribution. Finally, it is convenient to employ dimensionless variables A = k0a, δ U = k0/gδ u, and (X, Y, Z, T ) = (k0x − √ gk0t 2 , k0y, k0z, gk0t) . (6) In these variables, the third-order CNLS equation is [Stocker and Peregrine, 1999]: where subscripts denote partial derivatives. For larger steepness, i.e., larger A, the third-order equation does not sufficiently account for nonlinear effects, and we need the fourthorder CNLS (CNLS 4 ) equation, also given by Stocker and Peregrine [1999]: On the left hand side of (8) are the third-order terms (for comparison with (7)), while the right hand side contains all terms appearing at fourth order. Here ∇H = (∂/∂X, ∂/∂Y ) is the gradient in the horizontal plane, andΦ, Φc, and B represent the mean flow, surface current, and oscillatory parts, X -3 respectively, of the velocity potential φ: where ϕ = k0x − √ gk0t = X −T /2 is the phase. The surface elevation, which is the quantity of interest for our purposes, is similarly expanded as where the expansion coefficients may be obtained from the velocity potential as We have performed calculations using both third-order and fourth-order equations, but all results shown below are obtained using CNLS 4 . A split-operator Fourier transform method developed by Weidman and Herbst [1986] was employed to solve the PDE numerically. The evolution equation is separated into two parts: a free evolution part and a part containing nonlinear and current terms. In each time step, the wave envelope is transformed into momentum space where the free evolution operator is applied, and then transformed back into position space and acted on by the nonlinear and current operator. In the case of CNLS 4 , the free evolution part is: and the nonlinear and current part is: Without loss of generality, we choose the mean wavelength of the incoming wave to be 156 m, corresponding to a mean wave speed of 7.8 m/s. Considering both the computation time and statistical accuracy, the calculation The 10 km by 20 km area is discretized using a 512 by 1024 grid. The incoming wave is prepared at the left side of the grid, developing along the x axis from left to right, and is eventually absorbed in the right boundary.
was conducted on a 512 by 1024 grid, corresponding to a 10 km by 20 km ocean area or 64 by 128 wavelengths. The initial incoming wave is a random linear superposition of N monochromatic plane waves with different frequencies and propagating directions: Each wave vector ki in (12) can be expressed as where the random wave number variation k ′ follows a Gaussian distribution of half height width ∆k, and the angle θ from the mean propagation direction is uniformly dis- (To leading order, the results depend only on the standard deviation, ∆θ, and not on the higher moments of the angular distribution.) Each coefficient ci is drawn from a normal distribution, and the final sum is normalized to have the desired mean wave height. We have confirmed that our results are Nindependent for sufficiently large N ; specifically, N = 2 · 10 5 is used for obtaining the data shown in this paper.
The basic setup of the calculation is shown in Figure 1. The incoming wave is prepared as a random superposition of plane waves at the left side of the grid, developing along the x axis from left to the right, and is eventually absorbed at the right side boundary. The split-operator Fourier transform method used in the calculation automatically imposes periodic boundary conditions on the system. Absorption of the outgoing wave on the right is required to prevent it from re-entering the grid on the left side. However, to avoid reflection or other boundary effects, we allow the outgoing wave on the right to decay gradually using a tanh() multiplicative factor, and similarly another tanh() multiplicative factor allows the incoming wave to build up gradually in the initialization area on the left side of the grid, as indicated by the dotted lines in Figure 1. The boundary conditions in the transverse (y) direction remain periodic. By varying the system size in both x and y, we have confirmed that for the parameters used below, boundary conditions have a negligible effect on the wave height statistics.

Form of the Wave Height Distribution
In the Longuet-Higgins random seas model [Longuet-Higgins, 1957], the sea state is given by a random superposition of many plane waves with differing directions and frequencies, and by the central limit theorem, the surface elevation function ζ( r) is distributed as a complex Gaussian random variable with standard deviation σ. Furthermore, for a narrow-banded spectrum with small angular spread (|δ k| ≪ | k|), the crest height H is equal to the wave function amplitude |ζ|, and the probability distribution of wave height 2H is given by the Rayleigh distribution of (1). The corresponding cumulative probability of encountering a wave of height 2H or larger is We now consider the effect of nonlinearity and currents, where the wave envelope a( r, t) defined by (5) is evolving in accordance with the equations presented in the previous section. Since we still assume a narrow spectrum and a small angular spread, the envelope evolves slowly in space and time, on the scale of the mean wavelength and mean wave period, respectively. In analogy with the situation for linear ray dynamics in the presence of random currents [Ying et al., 2011], in the neighborhood of any spacetime point ( r, t), we have a wave intensity proportional to H 2 = |a( r, t)| 2 = |ζ( r, t)| 2 . Thus, in contrast with the original Longuet-Higgins model, the locally-averaged wave intensity is not uniform over all space and time but instead exhibits "hot spots" and "cold spots" in space-time associated with energy focusing and defocusing. Hence the wave height distribution is given locally in space and time by a Rayleigh distribution around the local mean height (corresponding to a locally random superposition of plane waves), while the local mean height itself varies slowly on the scale of the mean wavelength and mean period.
At each point in space-time, the central limit theorem and thus the Rayleigh distribution still apply, and we have where I( r, t) is the local energy density, normalized by the average energy density of the incoming sea state (so that I = 1), and σ 2 is the variance of the surface elevation in the incoming sea state, before the nonlinear evolution begins to act. Now averaging over space and over time, we obtain a total cumulative wave height distribution In (16), the full cumulative distribution of wave heights for a given sea state has been expressed as a convolution of two factors: (i) the local density distribution g(I), which at least for small deviations from Rayleigh is conjectured to be χ 2 distributed in analogy to the ray dynamics approximation, and (ii) the universal Longuet-Higgins distribution of wave heights for a given local density. Similar decompositions of chaotic wave function statistics into non-universal and universal components have found broad applicability in quantum chaos, including for example in the theory of scars [Kaplan, 1999;Smith and Kaplan, 2009;Bäcker and Schubert, 2002]. In the context of rogue waves, a similar approach was adopted by Regev et al. to study wave statistics in a onedimensional inhomogeneous sea, where the inhomogeneity arises from the interaction of an initially homogeneous sea with a (deterministic) long swell [Regev et al., 2008]. Taking the local mean intensity to be χ 2 distributed with N degrees of freedom, and convolving the χ 2 distribution of the mean intensity with the Rayleigh distribution around the mean intensity, we can obtain as in the linear case [Höhmann et al., 2010] a K-distribution for the total distribution of wave heights: where Kn(y) is a modified Bessel function. Defining the dimensionless variable x = 2H/Hs ≈ 2H/(4σ), so that a rogue wave is given by x = 2.2 and an extreme rogue wave by x = 3.0, we find the probability of a wave height exceeding x significant wave heights: to be compared with the random seas prediction in the same dimensionless units. Notably, the N parameter describing deviations from Rayleigh statistics may be directly related to the excess kurtosis of the sea state, γ2 = η 4 /(η 2 ) 2 − 3, which measures deviations from Gaussian statistics for the surface elevation η [Onorato et al., 2002;Janssen and Herbers, 2009]. For a narrow-banded spectrum, we have γ2 = (3/2)(2H) 4 /(2H) 2 2 − 3, and from (18)

Wave Height Probability Distributions
In our simulations, we choose the mean wave number to be k0 = 40 km −1 (a typical value in a normal sea state), corresponding to a group velocity v0 = 7.8 m/s, wavelength λ0 = 156 m, and period T0 = 10 seconds. Note that there is no loss of generality in this choice, as it merely sets the fundamental spatial and time scales for the wave dynamics. Each run simulates wave evolution for 4 · 10 5 seconds, corresponding to 4 · 10 4 periods. In the setup pictured in Figure 1, the wave front takes approximately t = 2000 seconds to travel from the initialization area to the absorbing area taking the straight line path, and the system fully equilibrates at around t = 7000 seconds. The time interval 7000 s < t < 400000 s is used to obtain wave height statistics in the wave evolution region defined by 4 km < x < 16 km, which indicated by the shaded area in Figure 1. All wave heights are normalized by the mean wave height 2H in a given simulation run, obtained using the time interval 7000 s < t < 20000 s. In the end, we calculate the probability distribution of wave heights in each run, and by repeating this process for different input parameters we obtain the dependence of the wave height distribution on the sea parameters.
At present, we set the current to zero, and focus on the nonlinear effects. In this case, we find that the wave height distribution for given sea parameters is reliably obtained from a single run given a sufficiently long run time, as discussed above. The two important properties of the incoming wave are the steepness and the wave vector spread. The , and compared with the random seas model of (20) (lowest solid curve). In each case, the solid curve is a best fit to the K-distribution of (19). Here the we fix the angular spread ∆θ = 2.6 • and wave number spread ∆k/k0 = 0.1 of the incoming sea.

X -5
steepness ε = k02H/2 controls the relative importance of the nonlinear term 1 2 B|B| 2 in (8) compared with the linear terms; the steepness is adjusted by varying the mean height 2H of the incoming sea, since the mean wave number k0 is fixed throughout.
Two parameters are needed to describe the wave vector spread: the variation in magnitude (wave number) is characterized by the ratio ∆k/k0 (or equivalently by the frequency spread ∆ω/ω0 ≈ (1/2)∆k/k0), and the angular spread is given by ∆θ. The objective is to find the response of the freak wave occurrence probability to these three input parameters.
First, we keep the wave vector spread constant, and vary the steepness in the range 0.016 < ε < 0.07. The steepness as defined above is normally less than 0.09 under storm conditions [Dysthe et al., 2008]. Since we are focused on the normal sea state, ignoring for example tsunamis (large tidal waves caused by underwater earthquake), we are not interested in the behavior at anomalously large values of the steepness. The lower limit is determined by two considerations. Physically, when the steepness becomes very small, the distribution of wave heights approaches the Rayleigh limit, where the wave height distribution is already well understood, and the number of extreme waves is very small. Numerically, at small values of the steepness, the finite simulation time makes it difficult to quantify deviations from the baseline Rayleigh distribution, due to the paucity of events in the tail of the wave height distribution.
Typical results are represented by dashed or dotted curves in Figure 2, where we fix ∆k/k0 = 0.1 and ∆θ = 2.6 • (the values of ∆θ required to see very strong effects from nonlinear focusing are typically smaller than those needed to observe significant deviations from Rayleigh by linear scattering [Ying et al., 2011]). The cumulative probability distribution of the wave height 2H, in units of the significant wave height Hs, is shown for three nonzero values of the wave steepness ε. From bottom to top, the three thick dashed or dotted lines show results for steepness ε = 0.019, 0.032, and 0.042. As expected, the Rayleigh probability distribution of (20) is recovered in the limit ε → 0, and ever stronger enhancement in the tail is observed as the steepness of the incoming sea increases. The occurrence probability of extreme rogue waves, 2H/Hs = 3.0, is enhanced by one to three orders of magnitude for the parameters shown.
The shape of the probability distribution and the typical enhancements in the probability compared with the Rayleigh baseline are consistent with the results of previous numerical simulations. For example, Onorato et al. [2006] have obtained an enhancement of up to one order of magnitude in the probability of freak wave formation (2H/HS > 2.2) and up to 2.5 orders of magnitude in the probability of extreme freak wave formation (2H/Hs > 3), for ε = 0.08 (k0Hs/2 = 0.125), ∆k/k0 = 0.18 (∆ω/ω0 = 0.09), and ∆θ = 0.
In Figure 2, each data set is fit to the K-distribution of (19), arising from the local Rayleigh approximation. We see that the fits, indicated by solid lines, perform adequately for probabilities down to 10 −6 , where statistical noise begins to dominate. The best-fit values are N = 53.2, 12.4, and 8.2, with smaller steepness corresponding to larger N . The Rayleigh distribution (N → ∞) is also shown here, corresponding to ε = 0.
In particular, we clearly observe the crossover between the Gaussian behavior (20) at small to moderate heights and asymptotic exponential behavior at large heights. However, systematic deviations do exist, which are especially visible at larger values of ε, corresponding to smaller values of the N parameter. These systematic deviations are in large part due to the fact that the true wave height distribution for any given set of input parameters exhibits spatial dependence, evolving from the original Rayleigh distribution imposed by incoming boundary conditions to the broader Kdistribution, and then gradually back to a Rayleigh distribution as the wave energy is transferred to longer wavelengths and the steepness decreases [Janssen and Herbers, 2009]. An example of this spatial dependence appears below in Figure 3. Thus, a more accurate model for the total wave height distribution consists of a sum of several K-distributions, or equivalently the tail of the full distribution may be modeled by a K-distribution multiplied by a prefactor C < 1. Nevertheless, as seen in Figure 2, (19) correctly describes wave height probabilities at the ±20% level of accuracy, allows for an extremely simple one-parameter characterization of the wave height distribution, and facilitates easy comparison between the effects of linear and nonlinear focusing.

Effect of Several Sea Parameters on the Wave Height Distribution
Based on Figure 2, we can clearly qualitatively conclude that larger steepness will lead to a higher freak wave occurrence probability. Now by fitting the distribution to (19), in Figure 4, we examine quantitatively how increasing the  Figure 4. The best-fit N value in (19) describing the wave height probability distribution is shown as a function of the steepness ε for several values of the initial angular spread ∆θ. The initial wave number spread is fixed at ∆k/k0 = 0.1, as in Figure 2. The line shows the best-fit scaling N ∼ ε −2.9 for ∆θ = 2.6 • . steepness will affect the N value. In the figure, each point corresponds to one particular sea state, and points of a given shape correspond to the same initial angular spread but different values of the initial steepness. The horizontal axis indicates the steepness, and the vertical axis shows the N value. From bottom to top, the four curves correspond to angular spread ∆θ = 1 • , 2.6 • , 3.6 • , and 5.2 • . A power law is observed on the log-log scale for each angular spread. Taking ∆θ = 2.6 • (indicated by squares in Figure 4) as an example, we obtain the scaling where c ≈ −3. The ∆θ = 3.6 • and ∆θ = 5.2 • curves lie very close to each other, indeed, at larger values of the steepness (not shown), saturation occurs, indicating that the occurrence rate of freak waves does not fall below a certain level when the steepness is sufficiently high, independent of the angle ∆θ. However, we also note that the CNLS equations are igure 6. The horizontal axis indicates the wave number spread ∆k/k0, and the vertical axis shows the best-fit N value, which determines the wave height probability distribution. Again, the steepness is fixed at ε = 0.032. The best-fit line shows the scaling of N with ∆k/k0, in this example, N ∼ (∆k/k0) 1.15 for ∆θ = 2.6 • .
based on a small steepness approximation, and thus the behavior obtained at very high values of the steepness may not be reliable. In the more realistic low-to-moderate steepness regime shown in the figure, N grows with increasing initial angular spread, i.e., the occurrence rate of freak waves falls. Now that we have an understanding of how N changes with steepness, the next step is to fix the steepness and vary the initial angular and wave number spread. At the end, we will be able to bring these results together, to obtain the dependence of the freak wave occurrence rate on all three key environmental parameters. Figure 5 shows the dependence of N on the initial angular spread for several different values of the frequency spread, while the steepness is fixed at 0.032. In this figure, points of a given shape represent sea states with the same wave number spread. The initial angular spread is varied on the horizontal axis.
From bottom to top, the four data sets represent wave number spreads ∆k/k0 = 0.2, 0.15, 0.1, and 0.08, and the dependence on initial angular spread is well represented by N ∼ (∆θ) α , where α = 0.96, 1.04, 1.07, and 1.03, for the four data sets, respectively. Thus, N scales as the first power of the angular spread at this value of the steepness, independent of the wave number spread.
Similarly, we can choose several typical values of the initial angular spread, and obtain the dependence of N on the  (19) describing the wave height probability distribution is shown as a function of wave vector spread (keeping ∆k/k0 = 7 · ∆θ/π while varying ∆k/k0 for several values of the steepness ε). The line shows the best-fit scaling with power 2.32 for ε = 0.032. initial wave number spread, as shown in Figure 6. Again, the steepness is fixed at 0.032, and the other physical conditions are the same as in the previous analysis. All sea states having the same initial angular spread are represented by the same symbol. From bottom to top, the angular spread is ∆θ = 1 • , 2.6 • , 3.6 • , and 5.2 • . All the initial sea states in this figure have steepness ε = 0.032. The power law is consistent for different initial angular spread, showing that N scales approximately as the first power of the wave number (or frequency) dispersion.
We notice deviations from the power-law behavior near the left edge of Figures 5 and 6, corresponding to very small angular spread ∆θ ≈ 1 • and very small wave number spread ∆k/k0 ≈ 0.08 (∆ω/ω0 ≈ 0.04), respectively. At even smaller values of ∆θ and ∆k (not shown), the CNLS equations we rely on, which are based on a perturbative expansion in powers of the surface elevation, begin to break down. This breakdown is evidenced, for example by large discrepancies appearing between the wave height distributions obtained using the fourth-order CNLS (8) and thirdorder CNLS (7), which make the CNLS expansion untrustworthy. Furthermore, we recall that the functional form of the K-distribution (18) is obtained assuming modest fluctuations of the local intensity around the mean intensity, i.e., this functional form makes sense only for N > 1, where the wave height distribution may be considered to be a modified Rayleigh distribution. We note in particular that the Benjamin-Feir modulational instability is known to be present in the NLS equations for (∆k/k0)/ε < 32/π [Alber , 1978], and for sufficiently small values of (∆k/k0)/ε, this instability will dominate the wave evolution so that extreme waves become commonplace rather than rare events. This regime is outside the range of validity of the present analysis. Nevertheless, as we will see below, the moderate parameter values for which our analysis is applicable, which are also the parameters most likely to occur in nature, can given rise to enhancement factors as large as 10 3 in the likelihood of occurrence of extreme freak waves (Table 2).

A Unified Model
In the previous subsection, we examined separately the dependence of the wave height distribution on steepness, angular spread, and wave number spread, keeping the other two variables fixed. We would now like to understand the combined dependence on these three parameters describing the initial sea state.
Since both angular and wave number spread control the area in wave vector space in which wave modes are permitted, we examine a scenario in which the angular spread  Figure 9. The initial wave number spread ∆k/k0 is varied for several values of the steepness ε, keeping the initial angular spread ∆θ = 5.2 • fixed. From bottom to top, the four data sets correspond to ε = 0.07, 0.042, 0.032, and 0.026. and wave number spread change at a constant ratio for a given value of the steepness. We arbitrary choose the ratio as ∆k/k0 = 7 · ∆θ/π. The results are shown in Figure 7, again plotted on a log-log scale, with the frequency spread ∆k/k0 on the horizontal axis. Each symbol represents all sea states with a given value of the steepness. The line shows the best-fit power law N ∼ (∆k/k0) 2.32 for ε = 0.032. From bottom to top, the other three curves correspond to steepness ε = 0.051, 0.042, and 0.026, corresponding to a scaling with power 2.3, 2.0, and 2.5, respectively.
First, we focus on the case of steepness ε = 0.032, where we have a clear understanding of how N depends on the angular and frequency spread. We notice that the power law exponent 2.32 appearing in Figure 7 when ∆k/k0 and ∆θ are varied together is close to the sum of the exponent obtained when ∆k/k0 and ∆θ are varied separately. Thus, we propose a model to describe the dependence of N on angular and frequency spread: where B, α, and β may depend on the steepness. Now, if we fix the ratio between the angular and wave number spread, d∆θ π = A · ∆k k 0 , (23) yields Comparing (23) and (24) with the data shown in Figure 5 and Figure 7, and setting A = 1/7, we find the results are consistent and give B = 5279, α = 1, and β = 1. Thus, we obtain the following predictive model (when ε = 0.032): The next step is to find out how B, α, and β in (23) may evolve when the steepness changes. For several typical values of the steepness ε, we vary the initial angular spread ∆θ, while keeping the frequency spread fixed at a typical value, ∆k/k0 = 0.1. The results are shown in Figure 8. Here each symbol represents sea states with one value of the steepness: 1.1 · 10 2 5.2 · 10 4 3 5 37 7.3 · 10 3 1.2 10 16 1.3 · 10 3 0.6 20 6.8 2.2 · 10 2 0.3 50 2.9 27 0.12 100 1.8 7.8 0.06 a The enhancement in the probability of rogue wave formation (wave height 2H = 2.2 Hs) as well as the enhancement of the probability of extreme rogue wave formation (wave height 2H = 3.0 Hs) and the excess kurtosis γ 2 (defined in (21)) are calculated for several values of the N parameter. Here E(x) = P total (x)/P Rayleigh (x). from bottom to top, the data sets correspond to steepness ε = 0.051, 0.042, 0.032, and 0.026. Fitting each data set to a power law relating N and ∆θ, we find that the power (α in (23)) is 0.93, 1, 1.15, and 1.19 for these four values of ε, respectively. Thus, α is negatively correlated to the steepness, i.e., the dependence on angular spread becomes slightly weaker as the steepness increases (consistent with the saturation effect noted earlier). Nevertheless, the effect is quite small for moderate steepness values, so to a good approximation we may treat α as a constant near unity.
Similarly, in Figure 9 we fix the initial angular spread fixed at a typical value ∆θ = 5.2 • , and examine the N dependence on the wave number spread ∆k/k0 for different values of the steepness ε. Each wave height probability distribution is fitted to obtain the N parameter, and a power law relating N and ∆k/k0 (the β exponent in (23)) is obtained for each value of ε. From bottom to top, the data sets correspond to steepness ε = 0.052, 0.042, 0.032, and 0.026, and yield best fit exponents β = 1.07, 0.98, 1.03, and 1.39, respectively. Since the variation in β is non-systematic and weak, we conclude that β, like α, may be approximated as an ε-independent constant near unity.
Thus, we conclude that in (23), both α and β are independent of steepness ε, and only B is a function of steepness: In Table 1, we present the value of B corresponding to values of the steepness ε. From Figure 4 and (22), we conclude that that N ∼ ε −3 , If we use this scaling to fit the data in Table 1, we obtain B = 0.203ε −3 , and thus the final model is described by: As noted previously, the scaling breaks down for very small ∆k/k0 or ∆θ, or very large ε, for example, when the twodimensional Benjamin-Feir index [Mori et al., 2011] is large. We now examine the implications for the occurrence probability of a freak wave. According to the conventional definition of a freak wave, the wave height for a freak wave is 3.51 2H, where 2H is the mean wave height. By integrating the wave height distribution over heights larger than 3.51 2H, we obtain the total freak wave occurrence probability. For example, for ε = 0.032, ∆θ = 5.2 • , ∆k/k0 = 0.1, we have N = 21.4. Comparing the result with a Gaussian theory, the freak wave occurrence probability is 6 times higher. Table 2 below shows the enhancement rate of the probability of freak wave occurrence in the nonlinear theory as compared with the simple Gaussian theory. We note that even at N values between 50 and 100, corresponding to the upper range of values in Figures 4 through 9, the occurrence of extreme rogue waves is enhanced by an order of magnitude. Exponentially larger enhancement is predicted for parameters associated with smaller values of N .

Spatial Dependence and the Kurtosis
The distribution data in Figures 2 to 9 is taken over the entire area of wave development thus does not give information about the spatial build up of rogue waves. To study the spatial distribution of the wave heights, we can divide the wave development area into slices of thickness around 0.3 km (approximately two wavelengths) in the wave propagation direction, and calculate the fourth moment of the wave height distribution, defined as (2H) 4 , in units where 2H = 1, for each slice. For zero steepness, the waves are linear, and the wave height distribution is Rayleigh, giving a value of 3.242 for the fourth moment of the wave heights. Figure 3 shows the fourth moment of the wave height as a function of position for different values of the steepness, where the angular spread and wave number variation are fixed at ∆θ = 1 • and ∆k/k0 = 0.1.
We can see from Figure 3 that the fourth moment of the wave height grows initially over a distance scale of several wavelengths, and then returns to the value given by a Rayleigh distribution. The result agrees with the Monte Carlo study conducted by Janssen and Herbers [2009].
In the following, we take a closer look at the wave height distribution in each slice. In this case we choose an example corresponding to the solid red curve in Figure 3, with ∆θ = 1 • , ∆k/k0 = 0.1, and steepness ε = 0.054. To reduce statistical noise, we divide the wave development area into larger slices of thickness around 1.3 km (approximately 8 wavelengths). The cumulative probability distribution of the wave height 2H, in units of the significant wave height Hs, is shown in Figure 10 for three different areas. The uppermost dot-dashed curve shows the result for the slice  Figure 10. The wave height distribution, in units of the significant wave height, is calculated for three different spatial regions (dashed or dot-dashed lines), and compared with the Rayleigh random seas model of (20) (lowest solid curve). In each case, the solid curve is a best fit to a single K-distribution of (19). Here the we fix the angular spread ∆θ = 1 • , steepness ε = 0.054, and wave number spread ∆k/k0 = 0.1 of the incoming sea.  Figure 11. The excess kurtosis γ2 = (3/2)(2H) 4 /(2H) 2 2 − 3 is plotted versus distance along the direction of wave development. The upper data set (squares) is obtained directly from simulation results, and the lower data set(circles) is calculated using (21) from a K-distribution fitted to the wave height distribution in each slice. The parameters are the same as in Figure 10. centered at the 7 km position; from Figure 3 we see that this is the peak region for wave height fluctuations. The short-dashed line shows the result for the slice centered at 16 km, which is near the end of the wave development area, where wave statistics are gradually reverting to the Rayleigh limit. Finally, the long-dashed line in the middle is the wave height distribution for the whole wave development area, as indicated in Figure 1.
In Figure 10, each data set is fit to the K-distribution of (19), arising from the local Rayleigh approximation. As noted in section 4.1, when wave height data are collected over a large spatial field that includes some areas of very strong deviations from Rayleigh statistics and other areas where such deviations have not yet had an opportunity to develop, the full distribution may not be well approximated by a single K-distribution. That is exactly what we observe in Figure 10: the wave height distribution for a single slice is better approximated by the K-distribution than that for the whole area.
We can also observe that the wave height distribution in a slice in the ending area is better approximated by a single Kdistribution than the distribution in the peak area. This is due to the fact that the distribution is evolving much faster in the peak area. Expanding to second order in slice thickness ∆x for large N , the deviation from a K-distribution due to finite slice thickness scales as (∆x/N ) 2 · (dN/dx) 2 , where dN/dx is the rate at which N evolves with position in the middle of that slice. This is consistent with the observation in Figure 3: the size of each slice is the same, but the slope dN/dx is larger and N is smaller in the peak area. Therefore, the deviation of the wave height distribution from a K-distribution due to finite slice thickness is larger in the peak area than in the ending area.
In order to further investigate the deviation of the exact wave height distribution from a K-distribution, we compare for each slice the excess kurtosis γ2 calculated directly from the simulation with the excess kurtosis obtained by fitting the wave height distribution to a K-distribution and then using (21). Of course, the answers would agree exactly if the distribution within each slice were perfectly described by a single K-distribution. The results are presented in Figure 11, which includes data from 8 slices, ranging from the peak area (7 km) to the ending area (16 km). The initial nonlinear focusing area between 4 km and 7 km is omitted because the N value is varying too dramatically within a slice, making a fit of the wave height distribution to a single K-distribution meaningless. In the figure, we note that from right to left the difference between the exact kurtosis calculation and the analytic formula of (21) increases systematically, but decreases in the leftmost slice, centered on the peak. If we compare with the solid curve in Figure 3, we see that |dN/dx| increases systematically from right to left and then decreases in the peak area. That is consistent with the results in Figure 11, given the scaling ∼ (∆x/N ) 2 (dN/dx) 2 of the error due to finite slice thickness ∆x.

Conclusions and Outlook
In this paper, we study the dependence of the wave height distribution on the steepness, initial angular spread, and frequency of the incoming sea. We use the CNLS 4 nonlinear wave equation to simulate the ocean wave development, and obtain quantitative predictions for the wave height distribution as a function of steepness, and angular and frequency spread. We then fit the wave height probability distribution to a K-distribution, governed by a single parameter N that quantifies deviations from the Rayleigh distribution predicted by a simple random seas model. Furthermore, we show that N may be related to the excess kurtosis. We find simple power laws for the dependence of the N parameter on the environment parameters as long as the sea conditions are not extreme, and as a result, we propose a simple model to predict N as a function of the initial sea state. By obtaining N for each sea state, we can quantitatively predict the enhancement of the freak wave occurrence probability for that sea state.
However, this is not a completed investigation. After fully understanding the quantitative consequences of the nonlinear effect, the ultimate goal of the project is to combine nonlinearity and deflection by random currents in a single model, which will allow the probability of rogue wave formation to be predicted for a wide range of realistic sea conditions. Previous linear results show that the effect of current is well characterized by a single parameter: the freak index [Heller et al., 2008;Ying et al., 2011]. Preliminary results of combining nonlinearity and deflection by random currents show that the freak index is still a dominant parameter.
When we introduce a nonzero current field, additional length scales come in to play, including the characteristic eddy size (the scale on with the random current varies) and the typical distance required for the first focal point to appear. It is still unclear how these parameters will interact with the wavelength of the incoming wave and the typical distance scale for the buildup of the nonlinear effect. An in-depth investigation is required to understand the underlying mechanism through which the formation of hot and cold spots is aided by nonlinear focusing.
Also, there is a clear need to compare the model simulations with observations and experiments. Although comprehensive global data are not available at this point, it may be possible to compare the results with local observations where data are more readily available, e.g. in the North Sea.