## Abstract

The temporal evolution of the energy spectrum of a field of random surface gravity waves in deep water is investigated by means of direct numerical simulations of the deterministic primitive equations. The detected rate of change of the spectrum is shown to be proportional to the cubic power of the energy density and to agree very well with the nonlinear energy transfer Snl as predicted by Hasselmann. Despite the fact that use of various asymptotic relations that are valid only for t → ∞ or integration with respect to time over a time scale much longer than O[period × (ak)−2] is necessary in the derivation of Hasselmann’s Snl, it is clearly demonstrated that the rate of change of the spectrum given by the numerical simulation agrees very well with Hasselmann’s Snl at every instance of ordinary time scale comparable to the period. The result implies that the four-wave resonant interactions control the evolution of the spectrum at every instant of time, whereas nonresonant interactions do not make any significant contribution even in a short-term evolution. It is also pointed out that the result may call for a reexamination of the process of derivation of the kinetic equation for the spectrum.

## 1. Introduction

In the standard treatment of ocean waves, the ocean wave field is considered as a superposition of an infinite number of component waves, which have various amplitudes, frequencies, and directions of propagation. These component waves interact with each other through the nonlinearity of the boundary conditions at the free surface. To describe this complicated nonlinear system with infinite degrees of freedom, we usually resort to a statistical description and introduce the notion of the energy spectrum. In terms of the wavenumber spectrum ɛ(k) or the directional spectrum Φ(ω, θ), the energy density E (i.e., the energy of the wave field per unit horizontal area) is expressed as

where k is the two-dimensional vector wavenumber in the horizontal plane, ω is the linear frequency given by ω(k) = (g|k|)1/2, and θ is the direction of propagation. Frequency spectrum Ψ(ω) is obtained by integrating Φ(ω, θ) with respect to θ:

If the wave field is homogeneous in space and there is no energy input from the wind or dissipation because of wave breaking, the only source of the change of the spectrum is the energy transfer Snl due to nonlinear interactions between component waves. In such a situation the evolution of the spectrum is described by the kinetic equation

For Snl, Hasselmann (1962) has derived a complicated yet explicit expression as follows:

where N(k) is the wave action density defined as N(k) = ɛ(k)/ω(k), ωi = ω(ki) (i = 1, 2, 3, 4), T1234 is a complicated function of k1, k2, k3, and k4, and δ is Dirac’s delta function. The existence of two delta functions for both k and ω implies that the energy transfer occurs only among those quartets of waves that satisfy the resonance condition

In most software packages for numerical wave prediction that are currently used, this form of Snl is employed as the theoretical basis for the modeling of nonlinear energy transfer.

Tanaka (2001a) performed a series of direct numerical simulations in order to investigate the validity of Hasselmann’s theory. He constructed random wave fields corresponding to the Joint North Sea Wave Project (JONSWAP) or Pierson–Moskowitz (P–M) spectra with cos2θ-type directionality and traced their temporal evolution deterministically according to the primitive equations for nonlinear surface gravity waves. Even though all the resonant and nonresonant interactions up to four-wave processes are included in his simulations, the rate of change of the energy spectrum that he detected agrees very well with Hasselmann’s Snl both qualitatively and quantitatively. But some doubt has been cast on his result, mainly because of the shortness of his simulations. In his simulations the wave fields evolve only for 20 periods, and this length of time appears definitely too short for any sensible discussion on the spectral change, whose characteristic time scale is generally believed to be O[T/(ak)4], where T and ak are the characteristic period and steepness of the wave field, respectively.

More recently Janssen (2003) reviewed the process of derivation of Snl and proposed a new form for Snl in which not only the contributions from resonant interactions but also that from nonresonant four-wave interactions are retained. Janssen’s new Snl reduces to Hasselmann’s Snl when t → ∞. Janssen claims that the time required to approach the asymptotic form may be so large that in the meantime the spectrum would experience a considerable change, and for this change the contributions from nonresonant interactions would not be negligible. However, this argument by Janssen apparently contradicts Tanaka’s result because Tanaka (2001a) observed the rate of change of the spectrum due to resonant interactions from the beginning of the evolution except for the transient phase, which lasts just the first several periods.

Then the aims of the paper are as follows:

1. By repeating similar numerical simulations as Tanaka (2001a) with improved accuracy and extended duration of time, we want to clarify the role of the resonant and nonresonant interactions in the short-term evolution of the spectrum of ocean waves, and thereby confirm the validity of the results of Tanaka (2001a).

2. Some researchers (e.g., Annenkov and Shrira 2006) have presented arguments that one needs to trace the evolution of the wave field for O[T/(ak)4] ≈ thousands of periods in order to cancel out the effects of nonresonant interactions and extract the change of the spectrum that is brought about solely by the four-wave resonant interactions. The typical time scale of the evolution of the spectrum is O[T/(ak)4], and there is no doubt about it. But this does not necessarily mean that we need to trace the evolution for such a long time to observe the effects of resonant interactions. We believe that one of the most important findings of Tanaka (2001a) is the fact that the change of the spectrum due to four-wave resonant interactions can be detected from a short-term evolution of the wave field for just 20 periods. We want to confirm this fact again because we think this has an important implication as follows.

The structure of the paper is as follows. In section 2, we first review the process of derivation of nonlinear transfer Snl starting from the Zakharov equation. We explain in section 3 the method that we employ to evaluate the rate of change of the spectrum by direct numerical simulations, and the results obtained are presented in section 4. Conclusions and discussion are given in section 5.

## 2. Kinetic equation

### a. The Zakharov equation

We consider the irrotational motion of water of infinite depth. The water is assumed to be inviscid and incompressible, and the flow is described by a velocity potential ϕ(x, z, t), which satisfies Laplace’s equation. Then the governing equations for the motion of nonlinear surface gravity waves are

where η(x, t) is the free surface displacement and ∇h ≡ (∂/∂x, ∂/∂y) is the gradient operator in the horizontal xy plane. The vertical coordinate z is pointing upward with its origin at the mean free surface. In terms of the velocity potential ψ(x, t) {=ϕ[x, η(x, t), t]} evaluated at the free surface, the boundary conditions in (7) and (8) can be rewritten as

where W(x, t) is the vertical velocity at the free surface.

Zakharov (1968) proved that the system in (10) and (11) can be described as a Hamiltonian system as follows:

where the Hamiltonian H is the total energy

Zakharov also rewrites the system to a simpler Hamiltonian system

by introducing a new canonical variable b(k, t) defined by

which will be called the complex amplitude function. Here η̂(k) and ψ̂(k) are the Fourier transform of η(x) and ψ(x), respectively. Then H can be expressed by a series expansion in b(k) as

where U(1)012 and V(1)0123 are some complicated functions. Krasitskii (1994) found an optimal canonical variable c(k), which is related to b(k) through a nonlinear canonical transform as

In terms of c(k), H is expressed as

and the resultant canonical equation reads

which is known as the Zakharov equation. For the detailed expressions for T0123 and other kernel functions, see Krasitskii (1994).

### b. Derivation of the kinetic equation

In this study Hasselmann’s Snl and Janssen’s modification to it are of central importance, and we briefly review the process of their derivations following Janssen (2003). Similar derivation can also be found, for example, in Yuen and Lake (1982) and Zakharov et al. (1992).

The second and the fourth moment of c(k) are introduced by

where the angle brackets denote an ensemble average and Djklm is the fourth cumulant. The delta function in Bij is a result of the assumption of spatial homogeneity. From (19) the following evolution equations for these moments can be derived:

If we apply the random phase approximation in (21) and require D = 0, we have ∂Ni/∂t = 0, and no spectral evolution results. So we need to apply the closure at a higher order. Requiring the sixth cumulant to vanish and integrating (22) with respect to t under the assumption that D(0) = 0, we obtain the following expression for D:

where Δω = ωi + ωjωkωl, and

where

Substitution of the expression (23) for D into (21) results in the following equation:

and this is the new kinetic equation derived by Janssen (2003). If we replace the resonance function Riω, t) by its asymptotic form πδω), which is valid as t → ∞, (26) reduces to the standard kinetic equation

derived by Hasselmann (1962).

Note here that, although (27) is a differential equation with respect to t, it has been derived by using some asymptotic relations that are valid only for t → ∞. This fact is apparent in its derivation via Janssen’s new kinetic equation in (26) as above, but also in the original derivation of (27) by Hasselmann (1962), the use of various asymptotic relations valid only for t → ∞ are essential. Then one may wonder what the “t” in (27) is. It should be some new time scale that corresponds to infinitely large time when measured in terms of an ordinary time scale, which is comparable to the period of the wave field. However, in the context of wave forecast, one treats (27) as if it holds at every instant of an ordinary time scale and integrates it numerically with using a step size of some 10 min. Furthermore, Tanaka (2001a) observes a rate of change of the spectrum consistent with (27) in his direct numerical simulations of the deterministic primitive equations at every instant of an ordinary time scale, but this cannot be regarded as a trivial consequence of (27) at all in view of the method of derivation of (27).

Very recently Annenkov and Shrira (2006) derived another modification to the standard kinetic equation and investigated the role of nonresonant interactions based on the equation.

## 3. Numerical method

Our method of detecting the rate of change of the spectrum by direct numerical simulations consists of the following three steps:

1. Construct an initial wave field corresponding to a prescribed directional spectrum.

2. Trace the evolution of the wave field deterministically by direct numerical simulation of the governing equations for nonlinear surface gravity waves.

3. Estimate the directional spectrum at later times from the deterministic variables b(k, t) or equivalently [η(x, t), ψ(x, t)] in order to detect the change of the spectrum.

In the following we explain the methods to perform these steps briefly. More detailed explanation can be found in Tanaka (2001a).

### a. Numerical setup

We trace the temporal evolution of the wave field according to (10) and (11). For solving the Dirichlet problem of Laplace’s equation for ϕ(x, z, t) to obtain W(x, t), we use the high-order spectral method (West et al. 1987; Dommermuth and Yue 1987; Tanaka 2001b), and the integration in t is performed by the fourth-order Runge–Kutta method. In the high-order spectral method, (10) and (11) are expanded in η and ψ and truncated at some order M, and we choose M = 3 here. This choice of M implies that we take account of all the three- and four-wave nonlinear interactions, irrespective of whether they are resonant or nonresonant. This is completely equivalent to truncating the Hamiltonian H of (16) at O(b4) and solving the corresponding canonical equation.

The basin we treat numerically is a square with length L on the physical xy plane, and the wave field is assumed to be periodic both in x and y with period L. It is understood that the x axis is taken along the principal direction of propagation. The basin is covered by a rectangular mesh with Nx × Ny nodes. In all of the simulations, Nx = 213 = 8192 and Ny = 212 = 4096.

To remove the aliasing errors that are generated in the high-order spectral method, the maximum mode numbers kmax and lmax in the x and y directions are determined as kmax = Nx/(M + 1) = 2048 and lmax = Ny/(M + 1) = 1024, respectively. This method of removal of aliasing error reduces to the well-known “3/2 rule” when M = 2, which is often used in computational fluid dynamics (e.g., Canuto et al. 1988).

In most computations, we require that the maximum mode number kmax correspond to the 16th harmonic of the peak of the spectrum. This in turn fixes the mode number kp at the peak as kp = kmax/16 = 128, implying that the total basin on the physical x plane is a square with area 128λp × 128λp, with λp(=2π) being the wavelength corresponding to the peak of the spectrum. On the other hand, the corresponding region in the k space is a rectangle −16 ≤ kx ≤ 16, −8 ≤ ky ≤ 8, and it is covered by a uniform square mesh with a mesh size Δk = 1/128. Note that k, l, kmax, lmax, and kp denote mode numbers and take integral numbers, while kx and ky are components of a vector wavenumber k and take real values. Here we have employed the normalization of space and time such that the gravity g and the angular frequency ωp at the peak of the spectrum are both unity. This normalization is used throughout this study. As the mode numbers k and l run from −kmax to kmax, and from −lmax to lmax, respectively, the total number of component waves whose temporal evolutions are followed amounts to more than 8 million (4097 × 2049, more precisely).

### b. Construction of initial wave field

Because the wave field is assumed to be periodic both in x and y with period L, η(x) and ψ(x) are represented by discrete Fourier transforms as

where the wavenumber k is discretized as

Then the complex amplitude function b(k) is defined as

only for discrete values of k. Each mesh point in the k space represents a component wave with corresponding wavenumber, and bk determines its amplitude and phase.

Equating the two different expressions for E (=) as

gives an approximate relation between |bk| and Φ(ω, θ) as

and this relation determines the distribution of {bk} corresponding to a prescribed Φ(ω, θ). The initial phases of {bk} are given by a homogeneous random number in [0, 2π].

We assume that Φ(ω, θ) is expressed as Φ(ω, θ) = Ψ(ω)G(θ). For the frequency spectrum Ψ(ω), we employ the JONSWAP spectrum

with

and the bimodal spectrum employed by Masuda (1980):

where

For the directional dependence G(θ), we employ cos2θ, cos4θ, and the “Donelan type” (Donelan et al. 1985) as given by

The energy density E is varied in the range 0.001 ≤ E ≤ 0.005. This range corresponds to 2.0 m ≤ H1/3 ≤ 4.5 m if the peak period is 8 s and the corresponding wavelength is 100 m.

### c. Estimate of spectrum

By connecting two different expressions for E as

we obtain an approximate expression for the frequency spectrum Ψ(ω) at discrete values of ω with an interval δω as follows:

where Σ′k,l denotes summation over those modes that satisfy the condition

We set δω = 0.05 (remember that ωp = 1). An interval [ωδω/2, ω +δω/2] corresponds to a circular ring on the k plane with radius ω2 and width 2ωδω. Because the area of the ring is 4πω3δω and the area of one mesh on the k plane is (Δk)2 = 1/1282, the number of component waves included in a circular ring is approximately 4πω3δω × 1282. Then the single marker at ω = 1 on a graph of Ψ(ω) like that shown in Fig. 1 below represents more than 10 000 component waves, while the marker at ω = 2 represents more than 80 000 component waves.

Fig. 1.

Frequency spectrum Ψ(ω) at t = 0, 25Tp, and 100Tp (JONSWAP, cos2θ, E = 0.003, case 0).

Fig. 1.

Frequency spectrum Ψ(ω) at t = 0, 25Tp, and 100Tp (JONSWAP, cos2θ, E = 0.003, case 0).

The directional spectrum Φ(ω, θ) with resolutions δω and δθ is estimated similarly by

where Σ″k,l denotes summation over those modes that satisfy

In this study we set δω = 0.05 and δθ = 10°.

## 4. Results

### a. Statistics of η(x, t)

Figure 1 shows the frequency spectrum Ψ(ω, t) obtained from the simulation in which the initial wave field has the JONSWAP spectrum with E = 0.003 and G(θ) ∝ cos2θ. In the figure Ψ(ω, t) is plotted at t = 0, 25Tp, and 100Tp, where Tp is the period corresponding to the peak of the spectrum and is equal to 2π according to our normalization. Note that the spectrum does not change significantly and remains practically the same as that given initially throughout the evolution up to 100Tp. Note also that the spectrum remains smooth even though it is not an ensemble average over many realizations but is obtained from just a single realization. This smoothness of the spectrum has been brought about by the large number of component waves that are represented by a single point on the plot of the spectrum as discussed in section 3c.

Figure 2 shows the same spectrum as shown in Fig. 1 in a log–log scale. Although the spectrum remains approximately the same up to t = 100Tp as shown in Fig. 1 as far as the absolute amount of change is concerned, the relative change is seen to be significant in the higher-frequency region. The accumulation of energy near the cutoff frequency has been brought about by the nonlinear energy transfer. As we mentioned in section 1, we have not included the effect of wave breaking (as well as energy input from the wind) in our numerical simulations. For the actual ocean wave spectra, the dissipation due to wave breaking is known to affect a wide range of frequency, including the energy-containing region around the peak of the spectrum. Therefore if we introduce a model of wave breaking in our simulation, it would render the task of extracting the change of spectrum due to nonlinear energy transfer impossible.

Fig. 2.

Same as in Fig. 1, but in log–log scale.

Fig. 2.

Same as in Fig. 1, but in log–log scale.

The piling up of energy at the high-frequency end of the spectrum is of course a nuisance and would sometimes terminate the simulation because of a numerical overflow. To prevent this from happening, an artificial dissipation is sometimes introduced that is tuned so that it would be effective only within the region of highest frequencies. The necessity of such an artificial dissipation varies according to the length of time for which the evolution is to be traced as well as the magnitude of the nonlinear energy transfer toward higher frequencies. For example, if we were to trace the evolution of the wave field for thousands of periods in order to reproduce the equilibrium power-law spectrum as done by Yokoyama (2004), introduction of such an artificial dissipation would be indispensable. In the case of our simulations, however, we have not included any artificial dissipation from the following reasons: 1) for the range of E and the length of time that we are concerned here, the accumulation of energy at the highest frequencies does not cause any serious problems, and the simulations finish successfully; 2) even if the accumulation of energy at the highest frequencies appears significant in terms of relative magnitude, its absolute magnitude is still insignificant, as shown in Fig. 1, and does not affect the rate of change of the spectrum in the energy-containing region around the peak of the spectrum that we are most concerned with.

Figure 3 shows the evolution of the skewness S and the kurtosis K of η(x, t) obtained from the same simulation as that shown in Fig. 1. Since the initial wave field is constructed by superposing a huge number of independent linear waves as explained in section 3b, it is approximately a Gaussian field and hence S ≈ 0 and K ≈ 3 at t = 0. Although S and K show large fluctuations in the initial transient phase during which bound waves are generated, this transition finishes in several periods and S and K remain nearly constant afterward. This behavior suggests that the initial wave field has evolved to a nonlinear wave field that includes bound wave components in a way consistent to the spectrum with a finite amount of energy during the first several periods.

Fig. 3.

Evolution of skewness S and kurtosis K of η (the same case as Fig. 1).

Fig. 3.

Evolution of skewness S and kurtosis K of η (the same case as Fig. 1).

Figures 4 and 5 show the behavior of S and K for the first 25 periods for various values of E. Each line is obtained by taking the average of 10 simulations, all of which have the same initial spectrum but differ only in the initial phases of the component waves. The time required for the initial transition to finish does not depend on E, and the process finishes by t ≈ 10Tp for all values of E. The values that S and K settle down to after the initial transition are larger for larger E. Figure 6 is a log–log plot of S and = K − 3 as a function of E. The values plotted in the figure are their mean values evaluated during 9Tpt ≤ 25Tp. It is clearly observed that S and are scaled as SE1/2 and E, and this scaling is consistent with the analytical expression for derived by Janssen (2003).

Fig. 4.

Evolution of skewness S of η for various E (JONSWAP, cos2θ).

Fig. 4.

Evolution of skewness S of η for various E (JONSWAP, cos2θ).

Fig. 5.

Evolution of kurtosis K of η for various E (JONSWAP, cos2θ).

Fig. 5.

Evolution of kurtosis K of η for various E (JONSWAP, cos2θ).

Fig. 6.

Skewness S and kurtosis K vs E (JONSWAP, cos2θ).

Fig. 6.

Skewness S and kurtosis K vs E (JONSWAP, cos2θ).

Note here that there is a great quantitative difference between one- and two-dimensional propagations. Figure 7 shows the values to which settles down after the initial transition as a function of E. The lower curve is the same as that shown in Fig. 6, while the upper curve is obtained from the numerical simulations in which the wave field has the Pierson–Moskowitz spectrum at t = 0 and the propagation is restricted to be uni-directional. Although K is an increasing function of E(=) for both situations, there is a remarkable quantitative difference between them, and of unidirectional propagation is almost one order of magnitude larger than that of the case of two-dimensional propagation when E is the same. If we trace the evolution of the unidirectional wave field corresponding to the JONSWAP spectrum as shown in Fig. 1, we observe a much larger change in the spectrum during the same period of time.

Fig. 7.

Kurtosis K vs E for 1D and 2D propagations (1D: P–M; 2D: JONSWAP; cos2θ).

Fig. 7.

Kurtosis K vs E for 1D and 2D propagations (1D: P–M; 2D: JONSWAP; cos2θ).

This kind of enhancement or exaggeration of the effects of nonlinearity appears in various aspects in the behavior of one-dimensional wave fields. Recently we found by direct numerical simulations that the occurrence probability of freak waves increases monotonically with E and K in the case of one-dimensional propagation, but it does not show such a clear dependence on E or K when the propagation is extended to two-dimensional with cos2 θ-type directional spread. The result will be reported elsewhere in more detail in the near future. Thus the results obtained by those studies that restrict the propagation in one dimension can sometimes be misleading, and we should exercise appropriate caution in interpreting the results of such studies.

Figure 8 shows the probability density function (PDF) of η(x, t) obtained from the case with the JONSWAP spectrum with E = 0.003 and G(θ) ∝ cos2θ. In the figure we use a normalized variable ξ, which is defined by ξ = η/E. As expected from the behavior of S and K as shown above, the PDF of η shifts from the initial Gaussian shape to a shape with positive S and positive by t = 10Tp, and remains practically the same for the rest of the evolution. This fact supports again our view that the transition from the initial Gaussian field, which is constructed as a superposition of a huge number of linear free waves, to a consistent and quasi-steady nonlinear wave field completes by t = 10Tp, and after this transition is over, various statistics of η change in conjunction with the secular change of the spectrum with a time scale of O(T/E2).

Fig. 8.

Probability density function of η at t = 0, 10Tp, and 100Tp (JONSWAP, cos2 θ, e = 0.003, case 0).

Fig. 8.

Probability density function of η at t = 0, 10Tp, and 100Tp (JONSWAP, cos2 θ, e = 0.003, case 0).

### b. Nonlinear energy transfer

As shown in Fig. 1, the spectrum experiences only a very small change, and we need to estimate the rate of change of the spectrum based on this very small change. We pick up the spectral density at ω = 0.95 and 1.10 in the evolution of the spectrum as shown in Fig. 1 and plot them as functions of t/Tp in Fig. 9. It can be observed that, for both values of ω, Ψ(ω, t) changes almost linearly in t, implying that the rate of change of Ψ(ω, t), which will be denoted by T1(ω), remains almost constant during the course of the evolution. If we fit a straight line by the least squares method to the curve for ω = 0.95 over the range 9Tpt ≤ 100Tp, the slope of the straight line is 4.41 × 10−6 and the correlation coefficient R = 0.997. If we apply the least squares fitting only to the range 9Tpt ≤ 25Tp, the slope of the best-fit straight line is 4.43 × 10−6 and R = 0.969. This closeness to a straight line of Ψ(ω, t) as a function of t, which is observed for t ≥ 9Tp, strongly suggests that T1(ω) is constant so long as the spectrum remains constant and does not depend on t explicitly as suggested by the new form (26) of Snl derived by Janssen (2003). The above result also implies that we can estimate T1(ω) very accurately from the change of Ψ(ω, t) during 9Tpt ≤ 25Tp instead of 9Tpt ≤ 100sTp. It should be noted that Fig. 9 is a result of just a single simulation and any ensemble averaging over realizations is not applied yet.

Fig. 9.

Evolution of Ψ(ω) at ω = 0.95 and 1.10 (JONSWAP, cos2θ, E = 0.003, case 0).

Fig. 9.

Evolution of Ψ(ω) at ω = 0.95 and 1.10 (JONSWAP, cos2θ, E = 0.003, case 0).

Figure 10 shows T1(ω) thus obtained from 10 simulations for the same initial spectrum, that is, JONSWAP, E = 0.003, and cos2 θ directionality. They differ only in the set of random numbers specifying the initial phase of each component wave. As expected from the nice linearity of Ψ(ω) with respect to t as depicted in Fig. 9, the dispersion among T1(ω) obtained from different simulations is small, and a reliable ensemble mean can be obtained even if the number of realizations is modest. The problem of the relation between the number of realizations and the density of modes in k space, which are necessary to secure a prescribed reliability of T1(ω), is investigated by Tanaka and Yokoyama (2004), and their result suggests that the ensemble averaging is not indispensable and even a single realization would do, provided that the distribution of modes in the k space is sufficiently dense.

Fig. 10.

T1(ω) obtained from 10 realizations with the same spectrum The thick dashed line denotes the mean of 10 realizations (JONSWAP, E = 0.003, cos2θ).

Fig. 10.

T1(ω) obtained from 10 realizations with the same spectrum The thick dashed line denotes the mean of 10 realizations (JONSWAP, E = 0.003, cos2θ).

Figure 11 shows T1(ω)/Ψ3p for various values of E, where Ψp = Ψ(ω = 1, t = 0). The coincidence of T1(ω)/Ψ3p for different values of E implies that T1(ω) is actually in proportion to E3 as Hasselmann’s theory predicts. It should be noted that in our numerical simulations all of the three-wave and four-wave nonlinear interactions are taken into account, irrespective of whether they are resonant or nonresonant, and that it is not guaranteed at all by the framework of the study that T1(ω) is scaled by E3. The weird oscillation observed in T1(ω) for E = 0.001 would probably be due to the smallness of the change in Ψ(ω). When E = 0.001 the change in Ψ(ω) between t = 9Tp and t = 25Tp is almost two orders of magnitude smaller than the change that we see in Fig. 1.

Fig. 11.

Normalized T1(ω) for various E (JONSWAP, cos2θ).

Fig. 11.

Normalized T1(ω) for various E (JONSWAP, cos2θ).

Figure 12 shows Int = ∫0 |T1(ω)| as a function of E. This figure also shows clearly the fact that T1(ω) scales like E3 for the whole range of E treated here. As expected from the shape of T1(ω), the peak frequency ωp of the JONSWAP spectrum downshifts in time. Figure 13 shows the change of ωp obtained from the direct numerical simulations. The peak frequency ωp is estimated by approximating Ψ(ω) around the peak by a quadratic curve that interpolates the nearby three points around the peak. It can be observed that the rate of downshift of ωp is proportional to E2, which is again consistent with the fact that T1(ω) is scaled as E3.

Fig. 12.

Int = ∫0 |T1(ω)| vs E (JONSWAP, cos2θ).

Fig. 12.

Int = ∫0 |T1(ω)| vs E (JONSWAP, cos2θ).

Fig. 13.

Downshift of peak frequency ωp in time (JONSWAP, cos2θ).

Fig. 13.

Downshift of peak frequency ωp in time (JONSWAP, cos2θ).

Since we do not have a numerical code for evaluating the nonlinear energy transfer Snl predicted by Hasselmann (1962), we compare T1(ω), which we obtained by direct numerical simulations, with Hasselmann’s Snl as evaluated by Masuda (1980) and Resio and Perrie (1991). Their corresponding figures are electronically scanned and saved as PostScript files, and then the curves included in the file are digitized by using the reading of the pointer in Ghostscript software. The result is shown in Fig. 14. In the figure T1(ω) is the average of T1(ω)/Ψ3p obtained for various values of E, and Snl are also normalized such that Ψp = 1. The figure clearly shows that the T1(ω) obtained from the direct simulations is not just scaled as E3 but is nothing but the nonlinear energy transfer due to four-wave resonant interactions as that predicted by Hasselmann (1962). Among the three results shown in the figure only Masuda’s result differs quantitatively by a factor of about 2. In the figure Masuda’s result has been multiplied by 2. We do not know the origin of this discrepancy at the moment, but judging from the nice agreement of the functional form as a whole, there would probably be some misapprehension in the way of normalization of the result.

Fig. 14.

Comparison of normalized T1(ω) with Hasselmann’s theory (JONSWAP, cos2θ).

Fig. 14.

Comparison of normalized T1(ω) with Hasselmann’s theory (JONSWAP, cos2θ).

The two-dimensional nonlinear energy transfer T2(ω, θ) can also be evaluated similarly by treating the directional spectrum Φ(ω, θ; t) as a function of t, and fitting a straight line by the least squares method. Figure 15 is the surface plot of T2(ω, θ), which is obtained from our simulations for the case with JONSWAP spectrum with E = 0.005 and G(θ) ∝ cos2 θ. This shows perfect correspondence with the contour plot of Hasselmann’s Snl for the same spectrum evaluated by Masuda (1980).

Fig. 15.

T2(ω, θ) for JONSWAP with cos2θ.

Fig. 15.

T2(ω, θ) for JONSWAP with cos2θ.

So far we have shown the results only for JONSWAP with cos2 θ-type directionality. But we have applied the method to other spectral shape Ψ(ω) and directionality G(θ) as mentioned in section 3b, and confirmed the agreement of T1(ω) and T2(ω, θ) with Hasselmann’s Snl. As an example of these, we show the result for the bimodal spectrum (35) with E = 0.003 and G(θ) ∝ cos2θ. The Hasselmann’s Snl for this spectrum is calculated by Masuda (1980). The evolution of Ψ(ω) during 25Tp is shown in Fig. 16, and Fig. 17 shows T1(ω) and Hasselmann’s Snl as evaluated by Masuda (1980) multiplied by 2. The correspondence between T1(ω) obtained by the direct numerical simulations and Hasselmann’s Snl is obvious again.

Fig. 16.

Bimodal spectrum and its evolution up to 25TP [G(θ) ∝ cos2θ].

Fig. 16.

Bimodal spectrum and its evolution up to 25TP [G(θ) ∝ cos2θ].

Fig. 17.

Normalized T1(ω) and Hasselmann’s Snl for bimodal spectrum as shown in Fig. 16.

Fig. 17.

Normalized T1(ω) and Hasselmann’s Snl for bimodal spectrum as shown in Fig. 16.

## 5. Conclusions and discussion

We have repeated the direct numerical simulations of Tanaka (2001a) to detect the rate of change of the spectrum of the field of random gravity waves and compared it with the theoretical predictions by Hasselmann (1962) and Janssen (2003). Our simulations have been improved from those of Tanaka (2001a) in various aspects as follows:

• The number of mesh points in both x and y directions are doubled, hence the density of modes in the k space is also doubled in both directions. This improves the reliability of the spectrum, which is estimated by smoothing over nearby modes in the k space.

• For some of the simulations, the temporal evolution is traced up to t = 100Tp, which is 5 times as long as the simulations in Tanaka (2001a).

• The step size Δt of the Runge–Kutta method is almost halved, hence the conservation of the integrals of motions, such as the total energy, is greatly improved. Take the case with the JONSWAP spectrum with E = 0.003 and cos2 θ-type directional spread, for example; the total energy remained constant within 0.48% up to t = 20Tp in the simulation of Tanaka (2001a), while it remains constant within 0.034% up to t = 20Tp and 0.21% up to 100Tp in the simulation of the present study.

• A bimodal frequency spectrum is treated in addition to the JONSWAP spectrum. As to the directional dependence G(θ), cos4θ type and Donelan type are newly included as well as the cos2θ type, which was treated previously.

### a. Conclusions

The conclusions of the study can be summarized as follows:

1. Starting from a Gaussian wave field with skewness S = 0 and kurtosis K = 3, which is constructed by superposing a huge number of free waves with random initial phases, the wave field evolves to a nonlinear wave field that includes bound wave components in a way consistent to the energy spectrum with a finite amount of energy. The time required for this transition process does not depend on the energy density E and is less than 10 periods for the whole range of E treated. After this transition finishes, various statistics of η, including the skewness, kurtosis, and the probability density function, change very slowly in conjunction with the secular change of the spectrum with a time scale of O(T/E2) and remain practically constant throughout the evolution up to 100Tp.

2. The spectral density Ψ(ω) at a fixed ω evolves in time t almost linearly except for the first several periods. This linearity implies the existence of nonlinear energy transfer, which is completely determined by the spectrum autonomously and does not depend on t explicitly as suggested by Janssen’s new form (26) of Snl.

3. The nonlinear energy transfer that is evaluated based on the evolution of the spectrum given by the direct numerical simulations during 9Tpt ≤ 25Tp agrees very well with Hasselmann’s Snl for all the frequency spectra Ψ(ω) and directional dependence G(θ) treated in the study. This fact implies that the evolution of the spectrum is completely governed by the four-wave resonant interactions except for the initial transition phase of several periods. Janssen (2003) argues that the time required for his resonance function Riω, t) to approach its asymptotic form πδω) would be so large that in the meantime there would occur a considerable change of the spectrum for which the nonresonant interactions would make some significant contributions. But our results suggest that the approach of Riω, t) to πδω) finishes in just several periods, and there is practically no change of the spectrum in such a short time of transition. Therefore we would expect that the change of the spectrum due to nonresonant interactions would be of some importance only for the cases of one-dimensional propagation where resonant interactions are not allowed.

4. For each spectrum and the value of E, we have taken an ensemble average over 10 runs that differ from each other only in the initial phases of component waves in order to improve the statistical reliability of the results. However, as can be expected from the linearity of the evolution of spectral density at fixed ω‘s as shown in Fig. 9 or the small dispersion among T1(ω) obtained from different runs with the same spectrum as shown in Fig. 10, the ensemble average with respect to the sets of initial phases of component waves is not essential provided that the density of modes in the k space is sufficiently large. The equivalence between the ensemble averaging with respect to the sets of initial phases and the smoothing or the coarse graining in the k space in the evaluation of the spectrum, which is suggested by Tanaka and Yokoyama (2004), has been reconfirmed.

### b. Discussion

Some researchers (e.g., Annenkov and Shrira 2006) have presented arguments that one needs to trace the evolution of the wave field for O[T/(ak)4] ≈ thousands of periods in order to extract the change of the spectrum that is brought about by the four-wave resonant interactions. On the other hand, we have succeeded in detecting the nonlinear energy transfer due to four-wave resonant interactions as predicted by Hasselmann from short-term direct numerical simulations that trace the evolution of the wave field for just a few tens of periods, during which the spectrum practically remains unchanged. Owing to the great progress of the numerical wave forecast during the last few decades, we are now predicting the ocean wave spectra globally with ever-increasing accuracy. The result of the present study implies that it is possible to investigate, for a prescribed spectrum, various important characteristics of a nonlinear wave field that require information on the phases of component waves by performing short-term direct numerical simulations, and thus to construct a correspondence between the stochastic description of the wave field in terms of the spectrum and the deterministic description in terms of (η, ψ). Along the same approach, we are now investigating the correspondence between the spectrum and the occurrence probability of freak waves, and the result will be reported elsewhere in the near future.

As discussed in section 2b, some asymptotic relations valid only for t → ∞ are used in the original derivation of (27) by Hasselmann (1962). The use of these asymptotic relations is indispensable in order to delete the effects of all the nonresonant interactions and extract the contributions from resonant interactions. Yuen and Lake (1982) derived (27) again via Zakharov’s equation more systematically by using the method of multiple time-scale analysis. They introduced a new “slow time” τ = ɛ2t, where ɛ is a small parameter of order ak. Then the contributions from nonresonant interactions have been removed by integrating with respect to τ and making use of the “delta calculus,” something like

In any case Hasselmann’s kinetic (27) has been derived by making use of some asymptotic relation that is valid when t → ∞ or integration with respect to t over a very long interval of time (much longer than hundreds of periods if ak ≈ 0.1). But our direct numerical simulations based on the primitive governing equations for nonlinear surface gravity waves have revealed that the spectrum actually evolves in time at the rate predicted by (27) at every instance of the ordinary time scale. This result implies that the change of the spectrum is completely governed by the four-wave resonant interactions, not in the sense of long-term average but at every second.

The spectrum that is directly treated is the k spectrum both in the theoretical study and our numerical simulations. (The ω spectrum is derived from the k spectrum by using the linear dispersion relation.) However there is one important difference between the theory and our numerical study in the definition of the k spectrum. In the theory the k spectrum is defined by the ensemble average, while in our numerical study it is estimated by the so-called “direct method”, that is, the method of estimating the spectrum by combining the fast Fourier transform and the averaging or the smoothing over modes in the k space. This is one of the most standard methods used frequently in data analysis in various fields and would probably be the only possible method that we can rely on in practical applications because it is utterly impossible to prepare a large number of seas that are statistically equivalent to take the ensemble average.

If we seriously consider the assumptions used in the derivation of (27), we would not be allowed to predict the spectrum of the immediate future, say a minute later, based on (27). However, it would be irrational to expect that we can predict the spectrum thousands of periods later even though we cannot predict it a few periods later. Actually one treats (27) as if it holds at every instant of ordinary time scale and integrates it numerically with using a step size of about 10 min. Thus there seems to be some kind of inconsistency between the process of derivation of (27) and the way it is used. However our numerical result strongly suggests that (27) is valid at every instant of time for the spectrum that is defined by the direct method. In this sense our result seems to suggest that it might be possible to derive (27) without relying on asymptotic relations for t → ∞ or ensemble average but relying only on the smoothing in the k space. We do not know the solution to this problem at the moment, and further investigations should be needed.

In this paper we have shown for weakly nonlinear wave fields in deep water that the spectrum defined by the direct method evolves in time solely because of the effect of four-wave resonant interactions at every instance of ordinary time scale. Although we have not performed the simulations for the case of finite depth yet, the same conclusion would probably hold so long as the depth is moderately deep. However the situation might be essentially different when the water is shallow. Zakharov (1999) argues that the applicability of the weak turbulence theory itself might be strictly limited when the water becomes shallow from the following reason. According to the weak-turbulence theory developed by Zakharov and his colleagues, the energy spectrum evolves in time because of four-wave resonant interactions, and this fact is described mathematically by a kinetic equation. The spectrum that appears in the kinetic equation is defined as the correlation function of a canonical variable that is related to the original observable variables such as η by some very complicated nonlinear canonical transformation. According to the analysis by Zakharov (1999), the relative difference [n(k) − N(k)]/n(k) between the spectrum N(k) defined in terms of the new canonical variable and the spectrum n(k) defined in terms of observable variables is O[(ak)2] when the water is deep, and we do not need to make distinction between the two spectra in ordinary conditions. When the water is shallow, on the other hand, Zakharov’s analysis shows that [n(k) − N(k)]/n(k) is O[(ak)2/(kh)5]. This implies that the condition of applicability for a weakly nonlinear statistical theory of waves on shallow water is (ak)2 ≪ (kh)5, which can practically never be satisfied when kh is as small as, say, 0.2.

As the water depth decreases and the dispersion gets weaker, it is expected that the non- (but quasi) resonant three-wave interactions begin to play more and more important roles in the evolution of the spectrum than in the deep-water case. Harmonic generation or generation of spectral humps at second and higher harmonics of the peak of the spectrum may be typical examples of those phenomena that occur only when the water is shallow. Young and Eldeberky (1998) points out that there can be some room for argument even in the way to define the spectrum itself. Thus there seems to be many fundamental aspects related to the spectral representation of nonlinear wave field that should be clarified before we actually plunge into the investigation of interesting but difficult problem of the relative importance of the nonresonant three-wave processes and the resonant four-wave processes in shallow water.

## Acknowledgments

This work was supported by a Grant-in-Aid for Scientific Research from the Japan Society for the Promotion of Science (JSPS). Numerical calculations were performed on the computer facilities of the Information Technology Center of Nagoya University.

## REFERENCES

REFERENCES
Annenkov
,
A. Y.
, and
V. I.
Shrira
,
2006
:
Role of non-resonant interactions in evolution of nonlinear random water wave fields.
J. Fluid Mech.
,
561
,
181
207
.
Canuto
,
C.
,
M. Y.
Hussaini
,
A. Q.
Quarteroni
, and
T. A.
Zang
,
1988
:
Spectral Methods in Fluid Dynamics.
Springer-Verlag, 557 pp
.
Dommermuth
,
D. G.
, and
D. K. P.
Yue
,
1987
:
A high-order spectral method for the study of nonlinear gravity waves.
J. Fluid Mech.
,
184
,
267
288
.
Donelan
,
M. A.
,
J.
Hamilton
, and
W. H.
Hui
,
1985
:
Directional spectra of wind-generated waves.
Philos. Trans. Roy. Soc. London
,
315A
,
509
562
.
Hasselmann
,
K.
,
1962
:
On the non-linear energy transfer in a gravity-wave spectrum. Part 1. General theory.
J. Fluid Mech.
,
12
,
481
500
.
Janssen
,
P. A. E. M.
,
2003
:
Nonlinear four-wave interactions and freak waves.
J. Phys. Oceanogr.
,
33
,
863
884
.
Krasitskii
,
V. P.
,
1994
:
On reduced equations in the Hamiltonian theory of weakly nonlinear surface waves.
J. Fluid Mech.
,
272
,
1
20
.
Masuda
,
M.
,
1980
:
Nonlinear energy transfer between wind waves.
J. Phys. Oceanogr.
,
10
,
2082
2093
.
Resio
,
D.
, and
W.
Perrie
,
1991
:
A numerical study of nonlinear energy fluxes due to wave-wave interactions. Part 1. Methodology and basic results.
J. Fluid Mech.
,
223
,
603
629
.
Tanaka
,
M.
,
2001a
:
Verification of Hasselmans’s energy transfer among surface gravity waves by direct numerical simulations of primitive equations.
J. Fluid Mech.
,
444
,
199
221
.
Tanaka
,
M.
,
2001b
:
A method of studying nonlinear random field of surface gravity waves by direct numerical simulation.
Fluid Dyn. Res.
,
28
,
41
60
.
Tanaka
,
M.
, and
N.
Yokoyama
,
2004
:
On the effects of discretization of the spectrum in numerical study of water wave turbulence.
Fluid Dyn. Res.
,
34
,
199
216
.
West
,
B. J.
,
K. A.
Brueckner
,
R. S.
Janda
,
M.
Milder
, and
R. L.
Milton
,
1987
:
A new numerical method for surface hydrodynamics.
J. Geophys. Res.
,
92
,
11803
11824
.
Yokoyama
,
N.
,
2004
:
Statistics of gravity waves obtained by direct numerical simulation.
J. Fluid Mech.
,
501
,
169
178
.
Young
,
I. R.
, and
Y.
Eldeberky
,
1998
:
Observations of triad coupling of finite depth wind waves.
Coastal Eng.
,
33
,
137
154
.
Yuen
,
H. C.
, and
B. M.
Lake
,
1982
:
Nonlinear dynamics of deep-water gravity waves.
,
22
,
67
229
.
Zakharov
,
V. E.
,
1968
:
Stability of periodic waves of finite amplitude on the surface of a deep fluid.
J. Appl. Mech. Tech. Phys.
,
2
,
190
194
.
Zakharov
,
V. E.
,
1999
:
Statistical theory of gravity and capillary waves on the surface of a finite-depth fluid.
Eur. J. Mech. B/Fluids
,
18
,
326
344
.
Zakharov
,
V. E.
,
V. S.
L’vov
, and
G.
Falkovich
,
1992
:
Kolmogorov Spectra of Turbulence I—Wave Turbulence.
Springer, 277 pp
.

## Footnotes

Corresponding author address: Mitsuhiro Tanaka, Department of Mathematical and Design Engineering, Faculty of Engineering, Gifu University, 1-1 Yanagido, Gifu 501-1193, Japan. Email: tanaka@cc.gifu-u.ac.jp