Nonlinear Processes in Geophysics The effect of third-order nonlinearity on statistical properties of random directional waves in finite depth

It is well established that third-order nonlinearity produces a strong deviation from Gaussian statistics in water of infinite depth, provided the wave field is long crested, narrow banded and sufficiently steep. A reduction of thirdorder effects is however expected when the wave energy is distributed on a wide range of directions. In water of arbitrary depth, on the other hand, third-order effects tend to be suppressed by finite depth effects if waves are long crested. Numerical simulations of the truncated potential Euler equations are here used to address the combined effect of directionality and finite depth on the statistical properties of surface gravity waves; only relative water depth k greater than 0.8 are here considered. Results show that random directional wave fields in intermediate water depths, kh=O(1), weakly deviate from Gaussian statistics independently of the degree of directional spreading of the wave energy.


Introduction
The statistical description of the surface elevation and, in particular, the probability of occurrence for extreme waves is an important input for the design and operation of marine and coastal structures.Because waves are on average only weakly nonlinear, solutions of the water wave problem can be written as a power series, where the small parameter of the expansion is the wave steepness (ε); in the case of arbitrary water depth, a general expression for the nonlinear parameter can be found in Whitham (1974).At the lowest Correspondence to: A. Toffoli (toffoli.alessandro@gmail.com)order of approximation, the water wave problem is linear and random wave fields can be approximated as a superposition of sinusoidal waves with random amplitudes and phases (linear wave theory).In this case, according to the central limit theorem, the surface elevations is Normally (Gaussian) distributed.
In the ocean, however, waves tend to behave differently as crests are higher and troughs are shallower than linear theory predicts.In this respect, it is a common practice to approximate the surface elevation by including the second-order bound contribution for each free wave mode, i.e. second-order wave theory (Hasselmann, 1962;Longuet-Higgins, 1963).A number of probability density functions describing second-order waves have been derived by several authors (see, e.g., Tayfun, 1980;Forristall, 2000;Prevosto et al., 2000;Tayfun and Fedele, 2007, among others).
Despite the fact that the second-order approximation agrees with field measurements reasonably well (see, for example, Forristall, 2000;Toffoli et al., 2007), it does not include effects related to the dynamics of free waves.At thirdorder in wave steepness, though, there is a substantial change in the description of water waves.Whereas bound modes are still present, resonant and non-resonant interactions between free waves are also possible.In this respect, a number of numerical and theoretical works (see, e.g., Janssen, 2003;Onorato et al., 2001Onorato et al., , 2005Onorato et al., , 2006;;Mori and Yasuda, 2002;Mori and Janssen, 2006;Mori et al., 2007;Toffoli et al., 2008b) have shown that deep water third-order nonlinearity can lead to focusing of wave energy with the consequent formation of extreme events, provided waves are long crested (i.e.unidirectional), narrow banded and sufficiently steep.In a random wave field, this results in a strong deviation from Gaussian and second-order statistics.However, the effect related to A. Toffoli et al.: Third-order nonlinearity in water of finite depth deep water third-order nonlinearity is substantially reduced when non-collinear wave components (short crested or directional wave fields) are considered (Onorato et al., 2002).In particular, experimental and numerical studies revealed that short crested waves deviate from Gaussian statistics only due to bound wave contributions despite third-order nonlinear evolution (see, e.g., Socquet-Juglard et al., 2005;Gramstad and Trulsen, 2007;Waseda, 2006;Toffoli et al., 2008a).
In water of finite depth, waves induce currents consequently reducing the energy available for nonlinear focusing.For sufficiently small relative water depths (kh<1.36,where k is the wavenumber and h is the water depth), in particular, the third-order effect does not lead to wave focusing at least for long crested waves (Benjamin, 1967;Mori and Yasuda, 2002;Janssen and Onorato, 2007;Whitham, 1974).Nevertheless, unlike in deep water, nonlinear focusing can be triggered by transverse (three dimensional) perturbations (see, for example, Davey and Stewartson, 1974;Francius and Kharif, 2006).For shallow water waves (kh<0.8),Toffoli et al. (2008c) have shown that the occurrence of extreme events in random wave fields may depend on the properties of the directional distribution.At present, however, it is not yet clear whether transverse perturbations may significantly contribute to the statistical description of directional wave fields in the transition region from deep to shallow water conditions (arbitrary water depth), where a number of engineering activities usually takes place.
Here we use numerical simulations of the random sea surface to address the combined effect of finite depth and directionality on the evolution of third-order nonlinear wave fields; both spectral and statistical properties will be discussed.The simulations were performed by integrating numerically the truncated potential Euler equations with a higher order spectral method (see Dommermuth and Yue, 1987;West et al., 1987).A wide range of wave fields with different relative water depth and directional spreading were considered.We mention that, although several methods exist to simulate the truncated potential Euler equations, the selected approach has the advantage to allow a large number of random simulations within a reasonable computational time.Moreover, it does not have any constrains on the spectral bandwidth but its limitation is related to the order of nonlinearity implemented in the numerical method.
The present paper is organized as follows.In Sect.2, we first present a description of the numerical model and the initial conditions for the simulations.In Sect.3, the evolution of the spectral shape is analyzed and discussed.In Sect.4, we present a description of the evolution of the statistical properties of the free surface elevation; in particular, the probability density function and its high order moments (skewness and kurtosis) are discussed.Concluding remarks are presented in the last section.

The model
For the modeling of the dynamics of free surface waves, we adopt the assumption of an irrotational, inviscid and incompressible fluid flow.In this case there exists a velocity potential φ(x,y,z,t) which satisfies the Laplace's equation everywhere in the fluid.We restrict ourselves to the case of domains with constant water depth.At the bottom (z= − h) the boundary condition is such that the vertical velocity is zero (φ z |h=0).At the free surface (z=η(x,y,t)), the kinematic and dynamic boundary conditions are satisfied for the free surface elevation and the velocity potential at the free surface ψ(x,y,t)=φ(x,y,η(x,y,t),t).Using the free surface variables these boundary conditions are as follows (Zakharov, 1968): where the subscripts denote partial derivatives, and W(x,y,t)=φ z | η represents the vertical velocity evaluated at the free surface.
The time evolution of the surface elevation can be calculated directly from Eqs. (1 and 2).Here, we use a higher order spectral method (HOSM), which was independently proposed by Dommermuth and Yue (1987) and West et al. (1987).A comparison of these two approaches (Clamond et al., 2006) has shown that the formulation proposed by Dommermuth and Yue (1987) is less consistent than the one proposed by West et al. (1987).The latter, therefore, was applied for the present study.
HOSM is a pseudo-spectral method, which uses a series expansion in the wave steepness ε of the velocity potential of the form: where each φ (m) is a quantity of order O(ε m ).In the above expansion M is the order of approximation in nonlinearity.A Taylor expansion around z=0 is then performed for each φ (m) term and combined with the above expansion for the potential.After collecting all terms at each order in wave steepness we obtain a system of the form (see, e.g.West et al., 1987): Following West et al. (1987), the vertical velocity at the free surface W(x,y,t) is similarly expanded in series of terms of order O(ε m ): where the terms W (m) are computed from the φ (m) terms: For the case of a rectangular domain in space with dimensions L x and L y in x and y, assuming periodicity in both directions for the wave field, we can use the following expression based on a double Fourier series for each φ (m) term in finite water depth (see, e.g.Whitham, 1974): k,l (t) of the potentials φ (m) can be computed from Eq. ( 4) by using two-dimensional (fast) Fourier transform when the free surface elevation and the free surface velocity potentials are given as input.
Here, we considered a third-order expansion (i.e.M=3) so that three and four waves interaction is included (see Tanaka, 2001a,b).After evaluating the vertical velocity at the free surface at order M, the free surface velocity potential ψ(x,y,t) and the surface elevation η(x,y,t) can be integrated in time from Eqs. (1 and 2).The time integration is performed by means of a four-stage fourth-order Runge-Kutta method with a constant time step.All aliasing errors generated in the nonlinear terms are removed (see West et al., 1987;Tanaka, 2001b, for details).
A concise review of HOSM can be found in Tanaka (2001a), while application of this method to study the occurrence of extreme waves can be found in e.g.Mori and Yasuda (2002), Ducrozet et al. (2007), Toffoli et al. (2008a) and Toffoli et al. (2008b).Note that other numerical methods have also been proposed by Annenkov and Shrira (2001), Clamond andGrue (2001), andZakharov et al. (2002).A comparative analysis between the performance of the HOSM and other numerical approaches can be found in Clamond et al. (2006).

Initial conditions and simulations
In order to prepare the initial wave field, it is necessary to generate a directional, frequency spectrum E(ω, ϑ)=S(ω)D(ϑ), where S(ω) is the frequency spectrum and D(ϑ) is the directional function, and then to transform it into the associated wavenumber spectrum E(k x ,k y ) (details on the transformation from (ω, ϑ) to (k x ,k y ) coordinates can be found in e.g.Holthuijsen, 2007).As it is frequently used for many practical applications, the JONSWAP spectrum (see, e.g., Komen et al., 1994) was used to express S(ω), while a frequency-independent cos N (ϑ) function was then applied to model the directional function.The spectrum in wavenumber coordinates (k x , k y ) can be written as follows: k p =2π/L p and ϑ=arctan(k y /k x ); the parameter σ is equal to 0.07 if k k p and 0.09 if k>k p ; (N ) is the normalizing factor of the directional distribution function.For the present study, for convenience, we described the wave field with a dominant wavelength λ p =156 m (this corresponds to a peak period T p =10 s in water of infinite depth), Phillips parameter α=0.014 and peak enhancement factor γ =3.Such a configuration corresponds to a significant wave height H s =6.36 m and wave steepness k p a=0.13, where k p is the wavenumber associated to the dominant wavelength and a is half the significant wave height.
In order to consider different degrees of the directional spreading, different values of the spreading coefficient N were used, ranging from fairly long crested (large N ) to fairly short crested (small N ) waves.The following values were selected: N=840, 200, 90, 48 and 24 (see Fig. 1, where the values of the angle ϑ are reported in degree for convenience).Note that we only considered cases with different directional spreading, while the wave steepness k p a was kept unchanged, i.e., k p and a were constant.Hereafter, the coefficient N will be used to identify the wave fields.
In order to consider finite water depth effects, the evolution of each directional wave field was then simulated at different water depths.The following relative depths were used: k p h=20, 2.5, 1.8, 1.36, 1 and 0.8.Note that, for each depth, the peak period can be calculated from the dominant wavelength by using the linear dispersion relation.It is important to mention that the adopted model does not describe effects related to bottom topography (a uniform water depth is in fact assumed) and wave breaking.Because bottom topography and breaking dissipation may have a significant effect on the temporal evolution of shallow water waves (see, e.g., Herbers et al., 2007;Toffoli et al., 2008c), cases with k p h<0.8 were not considered in this study.
From the wavenumber spectrum, E(k k,l ), an initial twodimensional surface η(x, y, t=0) was computed using the inverse Fourier transform with the random amplitude and phase approximation.The random phases were assumed to be uniformly distributed over the interval (0, 2π), while  the random amplitudes were Rayleigh distributed (the initial wave field is therefore Gaussian).The velocity potential ψ(x,y,t=0) was obtained from the input surface using linear theory (see Eq. 7).The wave field was contained in a square domain of 1920 m with spatial mesh of 256×256 nodes.With this resolution, the spatial domain contains about 12 dominant waves; each wave was thus discretized with about 21 grid points.Theoretical and numerical studies (Janssen, 2003;Socquet-Juglard et al., 2005) have shown that deviations from the Gaussian statistics due to third-order nonlinearity occur on a short timescale, typically on the order of 30-40 wave periods.In the present study, the total duration of the simulation was set equal to 60 T p (here T p is the peak period in water of infinite depth).A small time step, t=T p /200=0.05 s, was used to minimize the energy leakage; we mention that the selected time step was much smaller than the period of the shortest waves considered in this study.The accuracy of the computation was checked by monitoring the variation of the total energy (see, e.g., Tanaka, 2001a).Despite the fact that the energy content showed a decreasing trend throughout the temporal evolution, its variation is negligible as the relative error in total energy does not exceed 0.4% over the simulation time (cf.Toffoli et al., 2008a).
The model output consists in the surface elevation, η(x,y,t), which was archived every six wave periods.For each simulated condition, we performed 100 repetitions with the same input spectrum but different random amplitudes and random phases in order to have enough samples to achieve statistically significant results.The stability of the statistical moments is discussed in e.g.Toffoli et al. (2008a).

Spectral evolution
Due to nonlinear wave-wave interaction, a part of the spectral energy is transferred between wave components.Therefore, the input spectrum does not retain its initial shape as the wave field evolves in time.More specifically, a fraction of the energy is moved from high to low wavenumbers, generating the downshift of the spectral peak (Hasselmann, 1962).In arbitrary water depth, however, the wave-induced currents have a notable impact on the energy transfer around the spectral peak.This was observed by Janssen and Onorato ( 2007), who performed Monte Carlo simulations based on the Zakharov equations of long crested wave fields in water of arbitrary depth.Their results showed that the spectral downshift did not occur for k p h=1.36, where no significant spectral changes were observed at all, while an upshift took place for shallower depths.
In order to verify whether the simulations of the truncated potential Euler equations recover the results in Janssen and Onorato (2007), the wavenumber spectrum was extracted from the output surfaces as an inverse Fourier transform; an ensemble average over the 100 random repetitions was performed too.In Fig. 2, we present the temporal evolution of the energy spectrum for wave fields with N=90 (the spectra are here plotted as the integral over the k y coordinate); similar results were obtained with other directional spreadings.The spectrum was scaled so that the spectral standard deviation is 1; the wavenumbers were normalized with the initial peak wavenumber denoted as k p(initial) .For deep water conditions, k p h=20, the downshift of the spectral peak was clearly evident, as the position of the peak was estimated at k p(final) /k p(initial) =0.93.As the relative water depth was decreased, though, the downshift became less pronounced and eventually disappeared for k p h=1.36 (spectral peak at k p(final) /k p(initial) =1.01) in agreement with Janssen and Onorato (2007).Part of the spectral energy was moved towards lower wavenumbers, nevertheless.For shallower depths, k p h<1.36, the spectral peak slightly migrated towards higher wavenumbers (k p(final) /k p(initial) =1.03).This weakly pronounced upshift is also consistent with results in Janssen and Onorato (2007).
A small fraction of the spectral energy is also transferred towards high wavenumbers and redistributed over the directional domain (Longuet-Higgins, 1976).As a consequence, the spectrum becomes broader (especially in the upper tail) as the wave field evolves.In order to provide an overall description of the evolution of the directional spreading, the directional properties of the energy spectrum were summarized into a mean directional spread factor.The latter was calculated from the wave spectrum of the output surfaces as the wavenumber-average of the second-order moment of the directional distribution expressed in (k, ϑ) coordinates (see, e.g., Hwang et al., 2000, for details): .
Hereafter, the wavenumber-average of Eq. ( 9) is referred to as σ 2a and converted from radians to degrees for convenience.
In Fig. 3, as an example, we present the temporal evolution of the mean directional spread factor at different relative water depths for wave fields with initial spreading coefficient N =90.It is evident that the directional spectrum changed its initial form, broadening the directional spreading.This result is consistent with numerical simulations of the modified nonlinear Schrödinger equation presented in Dysthe et al. (2003).However, we observed that the variation of the directional spreading was weakly affected by finite depth effects.As k p h was reduced, in fact, the broadening of the directional spreading occurred slightly less rapidly than in conditions of deep water depth.

Skewness and kurtosis
In the following section, we discuss the combined effect of directionality and finite depth on the third and the fourth order moments of the probability density function of the free surface elevation, which are known as the skewness (λ 3 ) and kurtosis (λ 4 ), respectively.Whereas the first is a measure of the vertical asymmetry of the wave profile, the second provides information on the probability of occurrence for extreme waves (see, e.g., Mori and Janssen, 2006).Note that λ 3 =0 and λ 4 =3 in Gaussian random systems (linear waves).However, some definitions of kurtosis assign value of zero to the Normal distribution.Such a convention was used in e.g.Janssen and Onorato (2007), where the kurtosis was calculated as λ 4 /3−1.
In Fig. 4, the temporal evolution of the skewness is presented for different degrees of directional spreading and relative water depths.As the linear (input) wave field evolves, the wave profile becomes more vertically asymmetric with the sharpening of the wave crests and the flattening of the wave troughs as a result of the bound wave contribution.Thus, the skewness increased, departing from the value expected in Gaussian random systems.Such a deviation occurred almost immediately as it was already visible after six peak periods.Because the skewness is weakly dependent upon the dynamics of free waves (see, e.g., Mori and Yasuda, 2002;Mori et al., 2007;Toffoli et al., 2008b), it remained rather constant for the rest of the temporal evolution.
As expected, substantial variation in the skewness were observed with the decrease of k p h, as finite depth effects become more pronounced.In Fig. 5, we illustrate the variation of the maximum value of the skewness over the duration of the temporal evolution of the wave field.On the whole, the directional properties of the wave field only had little influence on the vertical asymmetry, especially for k p h 1.36.It is interesting to note, however, that the coexistence of directional components increases the effect of the bound wave contribution for k p h<1.36 with a consequent increase of the skewness in fairly short crested wave fields.Similar results were also obtained from second-order simulations of the surface elevation (Forristall, 2000;Toffoli et al., 2006).
On the other hand, the kurtosis is influenced by the dynamics of the free wave components, especially for fairly long crested conditions.In deep water, this results in strong deviations from Gaussian statistics, which occur within a fairly short time scale of a few tens of peak periods (Janssen, 2003).In Fig. 6, the temporal evolution of the kurtosis is presented for different degrees of directionality and relative water depths.It is evident that the coexistence of different directional components notably reduces the effect of thirdorder nonlinearity and hence the deviation from Gaussian statistics in deep water conditions.This result is consistent with previous experimental and numerical studies (see, e.g., Onorato et al., 2002;Socquet-Juglard et al., 2005 et al., 2008a;Waseda, 2006).As expected, furthermore, the effect of third-order nonlinearity is also suppressed by finite depth effects.However, whereas substantial deviation of the kurtosis from the Gaussian statistics were still observed for long crested wave fields with k p h>1.36, third-order nonlinearity had no longer a significant effect for k p h 1.36, as the kurtosis weakly deviates from the value of 3 throughout the considered temporal evolution (the variation of the maximum kurtosis over the duration of the simulation as a function of k p h is illustrated in Fig. 7).This finding is consistent with previous investigation of unidirectional wave fields in Mori and Yasuda (2002) and Janssen and Onorato (2007).Our simulations, moreover, showed that third-order nonlinear interaction between non-collinear components in water of finite depth did not produce notable effects on the kurtosis.Within the mentioned water depth regime, these results are to some extent consistent with field measurements (see, for example, Toffoli et al., 2007).However, the kurtosis of the most short crested wave field, i.e.N=24, was substantially higher than in fairly long crested conditions for k p h=0.8 (see Fig. 7).For this relative depth, furthermore, the kurtosis of directional sea states weakly grew throughout the temporal evolution of the wave fields.
It is interesting to mention that Janssen and Onorato (2007) found negative values of kurtosis (or kurtosis lower than 3 with the convention used in the present study) for k p h<1.36; similar results were also obtained in shallow water regime using the Korteweg-de Vries equation by Pelinovsky and Sergeeva (2006).Although our numerical model recovered similar results for fairly long crested conditions (for example, see the case k p h=0.8 in Fig. 6), the simulations of directional sea states showed that the kurtosis was higher than the value expected for Gaussian random systems.

Probability density function of free surface elevation
The skewness and kurtosis provide an indication of the departure from Gaussian statistics.For the sake of completeness, it is interesting to observe the deviation of the tails of the probability density function of the free surface elevation from the Normal distribution.In Figs. 8, 9 and 10, in this respect, we present some examples of the probability density function of the simulated surface elevation; only the most long and short crested wave fields are shown (N=840 and N=24, respectively).In the case of deep water (k p h=20), the upper tail of the probability density function strongly deviated from the Normal distribution for long crested wave fields.We mention, however, that the lower tail of the distribution (i.e.wave troughs) fitted the Normal distribution reasonably well.On the other hand, the coexistence of transverse components reduced this deviation, which became no more significant than in second-order theory (cf.Socquet-Juglard et al., 2005;Toffoli et al., 2008a).Whereas the statistical properties of short crested wave fields remained rather similar as the relative water depth was reduced, the transition between deep to shallow water waves had a substantial effect on the statistics of long crested wave fields.As a result, the deviation of the upper tail of the probability density function from the Normal distribution was gradually reduced.For k p h=1.36, in particular, the upper tails were similar in both long and short crested wave fields.Note that the wave troughs remained Gaussian distributed for this relative depth.For lower water depth conditions, e.g.k p h=0.8, the probability of occurrence for extreme waves notably increased when directional components were taken into account (see Fig. 10).It is also important to mentions that the lower tail of the probability density function showed a substantial deviation from the Normal distribution for this relative depth, indicating shallower troughs than in Gaussian sea states.

Conclusions
For random wave fields in water of infinite depth, thirdorder nonlinearity is responsible for strong deviations from Gaussian statistics, provided waves are long crested, narrow banded and sufficiently steep.However, third-order effects tend to be suppressed by directional spreading and finite depth effects.In order to study the combined influence of wave directionality and water depth on third-order nonlinearity, a large number of simulations of the random sea surface was performed with the truncated potential Euler equations, considering a wide range of directional spreading and relative water depth.The higher order spectral method proposed by West et al. (1987) was used to this end; a third-order expansion was implemented so that the three and four waves interaction was taken into account.Under deep water conditions, our simulations recovered the strong deviation from Gaussian statistics in long crested wave fields as well as the suppression of third-order effects in directional wave fields.As the water depth is decreased, however, the deviation from Gaussian statistics of long crested, deep water waves was gradually reduced.At k p h=1.36, in particular, third-order effects were observed to be negligible on wave statistics.It is interesting to mention, nevertheless, that third-order nonlinear dynamics of free waves produced a weak upshift of the spectral peak for k p h<1.36, instead of the downshift commonly observed in deep water.
Although transverse perturbation may be expected to lead to focusing of wave energy in finite water depths, and hence generate large amplitude waves, our simulations did not show any substantial variation of the statistical properties of directional sea states in water of arbitrary depth, k p h=O(1), where only weak deviation from Gaussian statistics were observed.Nevertheless, for k p h=0.8, a weak effect of wave directionality on the upper tail of the probability density function was observed.At this relative depth, however, it is rather difficult to reach a firm conclusion on the significance of the directional properties of shallow water waves, as nonlinearity higher than third order might be significant in the description of wave statistics.
It is interesting to mention that deviations from Gaussian statistics occurred in the upper tail of the probability density function (i.e.positive elevation), while the lower tail (negative elevation) fitted the Normal distribution reasonably well.Significant deviations in the lower tail of the probability density function (for p(η)<0.01)were only observed for k p h<1.36.

Fig. 3 .
Fig. 3. Example of temporal evolution of the mean directional spread: case N =90.