## Abstract

A discretized equation for the evolution of random surface wave fields on deep water is derived from Zakharov's equation, allowing for a general treatment of the stability and long-time behaviour of broad-banded sea states. It is investigated for the simple case of degenerate four-wave interaction, and the instability of statistically homogeneous states to small inhomogeneous disturbances is demonstrated. Furthermore, the long-time evolution is studied for several cases and shown to lead to a complex spatio-temporal energy distribution. The possible impact of this evolution on the statistics of freak wave occurrence is explored.

This article is part of the theme issue ‘Nonlinear water waves’.

## 1. Introduction

Alber [1] has studied the stability of narrow-banded random waves to inhomogeneous disturbances by means of an equation derived from the Davey–Stewartson equation. The long-time behaviour of slightly unstable modes for this equation was subsequently investigated by Janssen [2]. In recent years, Alber's equation has also been studied extensively by Regev *et al.* [3], Ribal *et al.* [4] and Stiassnie *et al.* [5].

Alber [1, p. 544ff] addresses the question whether ocean-wave spectra are stable to inhomogeneous perturbations via a bandwidth criterion. While the fundamental equation used in his study is derived under the assumption of narrow-banded waves, a homogeneous spectrum is found to be stable to inhomogeneous disturbances if its bandwidth is sufficiently broad. This somewhat unsatisfying state of affairs invites a direct study of the instability of broad spectra, without the limiting processes for narrow bandwidth employed in the derivation of Alber's equation.

Some headway was made towards clarifying this situation by Stiassnie [6], who studied the stability of uni- and bimodal spectra, and found regions of instability even for well-separated modes on the basis of a random counterpart of an equation obtained by Rasmussen & Stiassnie [7]. In the present work, a rather general treatment of the stability and long-time behaviour of random seas is initiated, and applied herein to the most basic case of degenerate four-wave interaction. It should be emphasized that the intention of the present work is not to attempt to model ocean waves, but rather to study the leading-order effects arising in the statistically inhomogeneous description of random wave interaction.

## 2. A discretized equation for broad-banded random seas

The starting point of the present study is the reduced Zakharov equation introduced in [8] (see also the discussion by Craig & Wayne [9]) for the nonlinear interaction of deterministic waves. This is the foundation of many statistical model equations for random wave fields, including the kinetic equation (KE) [10], the modified kinetic equation (MKE) [11], the generalized kinetic equation (GKE) [12], as well as the Alber equation (AE) [13].

The discretized version of the Zakharov equation is obtained by substituting for discrete wave vectors **k**_{i}, where *δ* is the Dirac delta distribution, and integrating with respect to **k**, to yield
2.1
Here * denotes a complex conjugate, is a Kronecker delta function with the property that
and *T*_{npqr}=*T*(**k**_{n},**k**_{p},**k**_{q},**k**_{r}) is the symmetrized kernel introduced in [14] (see also [15, Ch. 14]). Here, the complex amplitudes *b*_{i}(*t*) are functions of time *t* only, related in the present discretized formulation to the free surface elevation *η*(**x**,*t*), to leading order, by
2.2
For deep water, which is used throughout this manuscript, the frequency *ω*_{i}=*ω*(**k**_{i}) is related to the wavevector **k**_{i} by the dispersion relation where *g*=9.81 *m* s^{−2} is the acceleration of gravity.

As the aim of the present study is an examination of the stability and long-time evolution of wave spectra, it is necessary to develop a suitable equation for the one-time, two-component spectral correlation functions
2.3
where 〈 〉 denotes an ensemble average. The stochastic model equations mentioned above (KE, GKE, MKE and AE), while requiring various underlying assumptions about the statistical homogeneity of the wave-field (or lack thereof), and the order of the approximation, are all derived by making use of a quasi-Gaussian hypothesis (to replace higher-order averages by products of lower-order averages) (see [16]), and the same is employed herein. Subsequently, the derivation proceeds in a manner analogous to [1] or [13]. The time evolution of the *r*_{nm} is given by
2.4
Define Substituting this into (2.4) yields
2.5
where *Δ*_{i,j,k,l}=*ω*_{i}+*ω*_{j}−*ω*_{k}−*ω*_{l}. Equation (2.5) is the discrete analogue of [13, Eq. (19)], and seems to be a rather convenient model equation for future studies of inhomogeneous water–wave fields, as was pointed out by Janssen [11].

## 3. The degenerate four-wave example

Rather than analyse (2.5) directly, it is convenient to make the simplifying assumption that only three distinct waves are present, i.e. the sums over *p*,*q*,*r* are allowed to take on only values in {*a*,*b*,*c*},^{1} corresponding to the wave vectors **k**_{a},**k**_{b} and **k**_{c}, so that 2**k**_{a}=**k**_{b}+**k**_{c}. This leads to the following system of six ordinary differential equations:
3.1
3.2
3.3
3.4
3.5
3.6
It may be shown that the system (3.1)–(3.6) has the following conserved quantities:
3.7
In the absence of inhomogeneous terms *R*_{ij}(*i*≠*j*) this system undergoes no evolution at this order *O*(*R*^{2}). This leads to the question whether the homogeneous state (with *R*_{aa},*R*_{bb} and *R*_{cc} terms only) is stable to inhomogeneous perturbations.

### (a) Linear stability analysis

Linearizing in the inhomogeneous terms decouples (3.4) and (3.5) from the remaining equations:
3.8
and
3.9
where the coefficients are given by
In what follows, the notation *Δ*=*Δ*_{bcaa} will be employed where there is no risk of confusion. The linear system (3.8)–(3.9) has a solution in the form
3.10
which is stable if and unstable if *σ* has a non-vanishing imaginary part. For given values of the homogeneous terms *R*_{aa},*R*_{bb} and *R*_{cc}, this may be determined by computing a discriminant *B*−*A*^{2}/4 with
3.11
To further investigate the presence of possibly unstable modes requires the specification of initial conditions for the homogeneous amplitudes (or spectral action densities) *R*_{aa},*R*_{bb} and *R*_{cc}. The specification of the three waves is performed by writing, with no loss of generality, **k**_{a}=(1,0)^{2} and **k**_{b}=(1+*p*,*q*) for some *p*,*q*. Consequently, to satisfy 2**k**_{a}=**k**_{b}+**k**_{c} requires **k**_{c}=(1−*p*,−*q*). The steepness *ε*_{i} of the three waves **k**_{i},*i*∈{*a*,*b*,*c*} is the second defining parameter, which is related to the amplitude *a*_{i} by *a*_{i}|**k**_{i}|=*ε*_{i}. In fact, for simplicity, *ε*_{c} is taken to be zero throughout. The two initially present ‘mother waves’ **k**_{a} and **k**_{b} will be shown to be unstable to the effects of small inhomogeneities, and thereby give rise to the growth of the daughter wave **k**_{c}. A deterministic counterpart to this situation was recently explored experimentally by Bonnefoy *et al.* [17].

The initial complex amplitudes are assumed Gaussian, and written as
3.12
and
3.13
where the magnitudes |*b*_{a}|,|*b*_{b}|,|*μ*_{a}| and |*μ*_{b}| are Rayleigh distributed, the phases *ϕ*_{a},*ϕ*_{b} and *ϕ* are uniformly distributed over [0,2*π*), and |*μ*_{i}|≪|*b*_{i}|.

Then, neglecting the products of small terms,
3.14
Specifying a value for *R*_{ii}(0) then means specifying an average value for |*b*_{i}|, which may be related in turn to the value of the amplitude *a*_{i} (see [15, p. 854]) by
3.15

Vanishing *ε*_{c} then implies that *R*_{cc}(*t*=0)=0. In the calculations to follow, four degenerate quartets with three distinct combinations of the initial steepness of **k**_{a} and **k**_{b} are considered: case (a) *ε*_{a}=*ε*_{b}=0.15, cases (b) and (c) *ε*_{a}=0.15,*ε*_{b}=0.05 and case (d) *ε*_{a}=0.05,*ε*_{b}=0.15.

Figure 1 depicts the unstable modes in the first quadrant of the (*p*,*q*)-plane. The shaded region contains the (*p*,*q*) values such that a wave-field initially containing two waves **k**_{a}=(1,0) and **k**_{b}=(1+*p*,*q*) with fixed steepness *ε*_{a} and *ε*_{b}, and small inhomogeneities will exhibit an initial growth of the inhomogeneous contributions with time, and hence give rise to further interaction among the waves in an averaged sense. Case (d) with *ε*_{a}<*ε*_{b} has negative discriminant, and so exhibits no instability; it is not depicted. The modes with the largest growth rate are found to be parallel for both (a) (*p*,*q*)=(0.285,0) and (b) (*p*,*q*)=(0.32,0). In moving from the point of largest growth along the region of instability, the growth rate decreases, e.g. in case (b) (figure 1*b*), the growth rate at (*p*,*q*)=(0.32,0) (red dot) is nearly four times larger than that for case (c) where (*p*,*q*)=(0.8,0.31) (blue dot).

These results confirm those obtained by Stiassnie [6] using a modified discretization of the Zakharov equation, and consequently may be seen as a preliminary extension of Alber's [1] results to widely separated modes. It may also be observed that the instability range and growth rate coefficients tend to be smaller for larger bandwidth. It should be emphasized that the instability here observed is that for the wave–action densities of two ‘mother waves’ **k**_{a} and **k**_{b} to small inhomogeneities. This should not be confused with cosmetically similar deterministic results [18,19] for the instability of a single Stokes’ wave to modulational perturbations.

### (b) Long-time evolution

Having established that there are degenerate quartets—even those with widely separated modes—for which initially small inhomogeneous perturbations grow with time, the behaviour of these unstable cases should be investigated. An analytical treatment of the subsequent evolution is complicated by the high dimension of the nonlinear system (3.1)–(3.6). Consequently, numerical integration of the system was carried out to analyse the long-time evolution of the four cases presented in §3a, and summarized in table 1.

In addition to a choice of steepness *ε*_{a},*ε*_{b} and wavenumber (*p*,*q*), initial conditions for the inhomogeneous disturbances must then be chosen. As delineated previously in §3a, one possible choice consists in assuming that the complex amplitudes are given as in (3.12)–(3.13):
3.16
and
3.17
where |*μ*_{a}|,|*μ*_{b}| are small compared with |*b*_{a}|,|*b*_{b}| and *ϕ*,*ϕ*_{a},*ϕ*_{b} are uniformly distributed random phases over [0,2*π*). Hence, upon averaging,
The values 〈|*b*_{a}|〉 and 〈|*b*_{b}|〉 are fixed by specifying a wavevector and steepness; see (3.15). It remains to specify the magnitude of the inhomogeneous contributions |*μ*_{a}|,|*μ*_{b}| which are chosen such that
3.18
for small *δ*=*o*(1). Note that if either 〈|*μ*_{a}|〉=0 or 〈|*μ*_{b}|〉=0, the initial inhomogeneous disturbance *R*_{ab}(0) vanishes identically and there is no subsequent evolution. In what follows, the terms 〈|*μ*_{i}|〉^{2}, smaller by a factor *δ*^{2} than 〈|*b*_{i}|〉^{2}, appearing in the initial conditions for the homogeneous terms *R*_{aa} and *R*_{bb} are neglected.

The evolution of the four cases (a)–(d) is depicted in figure 2. For cases (a) and (b) the evolution is for the most unstable modes (having the largest growth rate), while for case (c) very widely separated modes with a smaller growth rate are selected (figure 1). For case (d), the analysis of §3a showed that all modes are stable when the prescribed steepnesses are (*ε*_{a},*ε*_{b},*ε*_{c})=(0.05,0.15,0), and predicted no subsequent interaction. This is also demonstrated by numerical integration, as depicted in figure 2*d* for the parallel quartet (*p*,*q*)=(0.3,0). Note that all plots are 1000 s long, i.e. 498 periods of the wave **k**_{a}, although calculations for up to 10 000 s were performed.

Cases (a) and (b) exhibit strikingly different behaviour. For both, an initial inhomogeneous disturbance is only present in *R*_{ab}, and is one order of magnitude smaller than any of the homogeneous initial values. A period of visually indiscernible interaction for times *t*<100 is observed, during which the terms |*R*_{ab}|,|*R*_{ac}|,|*R*_{bc}| grow only slowly.

Thereafter, the initially small inhomogeneities for the unstable quartet (**k**_{a},**k**_{a},**k**_{b},**k**_{c}) grow to magnitudes comparable to the homogeneous modes. Likewise, the ‘daughter wave’ *R*_{cc} grows and interacts with the other modes. Case (a) exhibits slow, nearly-periodic recurrence in time. Case (b), on the other hand, undergoes a faster, non-recurring evolution with time. A similar division of the dynamics into simple recurrence and more complex behaviour has been observed for the evolution of narrow-banded random seas under inhomogeneous perturbations via the AE by Stiassnie *et al.* [5]. Case (c) presents the evolution of the modes (*p*,*q*)=(0.8,0.31) with steepnesses (*ε*_{a},*ε*_{b})=(0.15,0.05), which corresponds to a very large separation distance between the modes (∥**k**_{b}∥=1.83∥**k**_{a}∥). Despite this, significant, though slower, interaction is still observed. Throughout the numerical computations, the invariants (3.7) exhibit a relative error of 3.1×10^{−14},2.7×10^{−14} and 9.3×10^{−11} for *I*_{1},*I*_{2} and *I*_{3}, respectively.

While the homogeneous states associated with initial values of *R*_{aa},*R*_{bb},*R*_{cc} may be naturally assumed to be given, the source of the small inhomogeneous disturbances is less clear. However, it is shown in figure 3 that the dependence of the subsequent evolution on the size of the initial inhomogeneities has a particularly simple form. For clarity, only *R*_{aa} of case (a) is shown as a function of time, for different sizes of the initial inhomogeneous disturbances. The solid lines correspond to *R*_{ab}(0)=2.5×10^{−3}|*b*_{a}||*b*_{b}|, the dashed lines to *R*_{ab}(0)=4×10^{−4}|*b*_{a}||*b*_{b}| and the dotted lines to *R*_{ab}(0)=1×10^{−4}|*b*_{a}||*b*_{b}|, or values of *δ*=0.05,0.03 and 0.02, respectively. Provided the inhomogeneous disturbances are initially small, the subsequent evolution is unaltered by changes in their size, except for a shift in the emergence of interaction, i.e. in the duration of the ‘warm-up time.’

## 4. Spatiotemporal distribution of energy

In terms of the complex amplitudes *b*_{i}, the free surface of the fluid is written as
4.1
accurate to leading order, and where c.c. stands for the complex conjugate terms (see (2.2)). The variance of the surface elevation (proportional to the energy) associated with this wave field may be written as
4.2
where *r*_{ij} is defined as in (2.3).

The contours of the variance for cases (a) and (b) are presented in figure 4. It may be observed that case (a) has a larger homogeneous background energy than case (b), due to the larger steepness of wave **k**_{b} (*ε*_{b}=0.15 versus *ε*_{b}=0.05). The variance is seen to be periodic in *x* with period 2*π*/*p*, and quite naturally follows the long-time behaviour: for case (a) the variance is recurrent in time, while for case (b) the temporal evolution is more complex and non-recurrent. From these distributions of free-surface variance in time and space, it is possible to obtain statistical information for the wave-fields.

In each case, the resulting variations in the spatio-temporal distribution of energy may be interpreted differently: for the recurrent case (a) it suffices to select one cycle of interaction, while for case (b) the data are truncated up to the point of initial energy transfer (which depends on *δ* (see (3.18)), and here occurs at time *t*≈100). It is thus possible to present empirical (PDF) and cumulative distribution functions (CDFs) for the variance, as given in figure 5.

In figure 5*a*, the interaction from time 100 to 480 is used as the basis for the computation of the PDF and CDF. In figure 5*b*, the PDF and CDF are computed from the interactions taking place after the initial warm-up period (from time 110 to 1000). In each case, the variance (4.2) has been normalized by the homogeneous variance, *E*_{n}=*E*/*E*_{h}, with
cf. the first invariant (3.7). For case (a) *E*_{h}=0.018 m^{2}, for case (b) *E*_{h}=0.012 m^{2}.

In the recurrent case (a), spatio-temporal variations of the normalized variance (and thus the energy) range between 0.32 and 1.82, while case (b) exhibits a larger range of 0.015–2.79, indicating greater departures from the energy of the homogeneous state. In each case, the dynamics associated with inhomogeneous disturbances have the potential to significantly alter the wave-field. While the time domain in case (b) reaches only *t*=1000 (or approx. 500 periods; see above), even considerably longer data do not materially alter the distribution, or the main conclusions.

### (a) Freak wave statistics

It has been noted that the absence of inhomogeneous disturbances leads to an absence of evolution for the present equations. Consequently, if the complex amplitudes are assumed initially Gaussian, their statistics will remain unchanged at this order, and for the timescales considered herein. For a purely homogeneous wave field, the probability of the wave height *H* non-dimensionalized by the RMS wave height of the homogeneous case *H*_{rms0} exceeding a value *H*_{0} is calculated directly from the Rayleigh distribution
4.3
Taking into account the fluctuations in free-surface variance described above, the probability of exceeding waves of height *H*_{0} throughout the spatial and temporal evolution of the inhomogeneous wave field may be calculated akin to [4]:
4.4

where *f*_{E}(*E*_{n}) is the PDF obtained from the data of figure 5. These probabilities for cases (a) and (b) are plotted in figure 6. It is observed that the probabilities for encountering freak waves with heights equal to twice the significant wave height *H*_{s} (equivalent to *H*_{0}=2.85) are increased significantly when the wave field evolves as a consequence of inhomogeneous disturbances. For case (a) the probability increases by a factor of 4, and for case (b) by a factor of 13. For extremely high freak waves (exceeding 3*H*_{s}), the dynamics of case (a) yield an encounter probability some two orders of magnitude higher than the homogeneous value of 1×10^{−8}, while case (b) predicts an increase by three orders of magnitude when compared to a homogeneous ocean.

## 5. Conclusion

Hitherto the stability of random seas has been studied primarily from the standpoint of narrow-banded spectra. By employing a discretized, stochastic analogon of the reduced Zakharov equation, which contains no narrow-band restrictions, it seems possible to investigate directly some simple examples concerning the instability of broad spectra to inhomogeneous disturbances. This may be understood as a preliminary attempt to extend the analysis of Alber [1] and others. In addition, it presents a possible point of departure for a study of wave turbulence beyond the KE paradigm, which is restricted to homogeneous sea states.

While the equation presented is very general, its potential to produce new insight is best demonstrated by treating a degenerate quartet of waves. Even in this simple case, considering the interaction of three modes in a model of the nonlinear evolution of a random sea subject to inhomogeneous perturbations, instability is observed for suitable choices of the two mother waves **k**_{a} and **k**_{b}. For the unstable cases discussed, two waves **k**_{a} and **k**_{b} suffice to generate the third wave **k**_{c} satisfying 2**k**_{a}=**k**_{b}+**k**_{c}, through mutual interaction. In the cases investigated, the unstable region is of finite extent, and the maximum initial growth rate of the inhomogeneous terms occurs for modes with no transverse components (i.e. *q*=0). Two cases (a) and (b) evolving from the points of maximal growth and one case of modes with smaller initial growth (c) were investigated. Significant interaction was observed for all cases, even with the very widely separated mother waves **k**_{a}=(1,0), **k**_{b}=(1.8,0.31) of case (c).

While the cases considered herein, with *ε*_{c}=0, have unstable modes lying below the curve of exact resonance 2*ω*_{a}=*ω*_{b}+*ω*_{c}, other cases have been obtained from computations of the discriminant, including those lying wholly above the exact resonance curve and having fastest growing modes with *p*=0. Exact classifications of the various regimes of instability, and of the subsequent long-time evolution of the unstable modes remain interesting challenges for future work. Likewise, while higher-order inhomogeneous equations along the lines of [13, eqs. (28),(29)] quickly become cumbersome, future work may also be aimed at introducing *O*(*R*^{3}) terms to the equations for the homogeneous terms (see [20]).

The long-time behaviour of the unstable cases is complex and seems in most instances to be non-recurring. Insofar as the initial inhomogeneous disturbances are sufficiently small with respect to the homogeneous terms, the subsequent dynamics may be shown to be largely independent of their actual size; a shift in the magnitude of the initial inhomogeneities gives rise only to a shift in the time of initial interaction, a type of ‘warm-up’. The physical mechanism behind these inhomogeneities, possibly related to forcing by the wind, remains a major issue for future research.

Another fundamental issue which remains to be addressed is the relationship between the present stochastic model and simulations of the deterministic Zakharov equation. In formulating the stochastic model as detailed in §2, the relevant statistical properties are assumed to hold throughout the entire subsequent evolution (see [1, §3.1]). On the other hand, for Monte Carlo simulations of a deterministic equation with random initial phases, such assumptions are introduced only at the initial time. Various authors have reported different degrees of agreement between such Monte Carlo simulations and numerical solutions of the derived stochastic equations: for the (homogeneous) KE, Stiassnie & Shemer [21] noted the lack of even qualitative agreement for a toy model of four modes; Annenkov & Shrira [12] subsequently modified the Monte Carlo simulations by considering clusters to achieve agreement for four modes, and also reported good agreement for many-mode simulations close to equilibrium; for the (homogeneous) GKE, Gramstad & Stiassnie [20] found good agreement using one- and two-dimensional JONSWAP spectra; for the (inhomogeneous) AE, Onorato *et al.* [22] found reasonable qualitative agreement for the onset of instability when comparing with Monte Carlo simulations of the nonlinear Schrödinger equation using JONSWAP spectra. For our part, the agreement between equations (3.1)–(3.6) and Monte Carlo simulations of the Zakharov equation for a degenerate quartet remains an open problem. The inclusion of higher-order effects, or simply more modes may be necessary. Nevertheless, we believe the dynamics of the simplest occurring interaction—the degenerate quartet here considered—provide a blueprint for more complex, physically relevant cases to be considered in future work.

The fact that the presence of inhomogeneous disturbances in an otherwise homogeneous wave field sets up fluctuations in wave energy in (*x*,*t*)-space indicates that the corresponding statistics of wave height will also be altered. For recurrent dynamics, the resulting energy distribution may be extended to arbitrary times, while for non-recurrent dynamics the computations have been performed for an interaction time of approximately 500 wave periods. Assuming that the wave heights of the homogeneous sea are Rayleigh distributed, the presence of inhomogeneities yields up to an 11-fold increase in encounter probability for freak waves twice the significant wave height. It is left to future work to investigate the behaviour of the fundamental equation (2.5) for a larger number of modes, thus more closely approximating a real sea spectrum. The present investigation of the most basic case of degenerate quartets already gives a clear indication of the potential instability of broad sea spectra to inhomogeneous disturbances, and the richness of the dynamics appearing for the fundamental stochastic equation (2.5). Considerable recent work has been devoted to resonant interactions and Hamiltonian formulations for problems related to water waves, including effects of capillarity, vorticity or two fluid layers [23–27], and the considerations pursued herein for the Zakharov equation for deep water waves may well be extended to other primitive equations describing various physical or mathematical problems.

## Data accessibility

This article has no additional data.

## Authors' contributions

The present study was conceived and executed collaboratively by both M.S. and R.S. Both the authors read and approved the manuscript.

## Competing interests

We declare we have no competing interests.

## Funding

This research was supported by the Israel Science Foundation (grant nos. 261/17 and 464/13).

## Footnotes

One contribution of 19 to a theme issue ‘Nonlinear water waves’.

↵1 These three waves are usually called a ‘degenerate quartet’ because the wave

**k**_{a}is counted twice.↵2 Hence |

**k**_{a}|=1 m^{−1}serves to determine, together with the gravitational acceleration*g*=9.81 m s^{−2}, the scales of the example cases.

- Accepted May 30, 2017.

- © 2017 The Author(s)

Published by the Royal Society. All rights reserved.