Evolution and Statistical Analysis of Internal Random Wave Fields within the Benjamin–Ono Equation

: This study investigates the numerical evolution of an initially internal random wave ﬁeld characterized by a Gaussian spectrum shape using the Benjamin–Ono (BO) equation. The research focuses on analyzing various properties associated with the random wave ﬁeld, including the transition to a steady state of the spectra, statistical moments, and the distribution functions of wave amplitudes. Numerical simulations are conducted across different Ursell parameters, revealing intriguing ﬁndings. Notably, it is observed that the spectra of the wave ﬁeld converge to a stationary state in a statistical sense, while exhibiting statistical characteristics that deviate from a Gaussian distribution. Moreover, as the Ursell parameter increases, the positive skewness of the wave ﬁeld intensiﬁes, and the kurtosis increases. The investigation also involves the computation of the probability of rogue wave formation, revealing deviations from the Rayleigh distribution. Notably, the study uncovers distinct types of internal rogue waves, speciﬁcally referred to as the “two sisters” and “three sisters” phenomena.


Introduction
Ocean internal waves are fascinating phenomena characterized by underwater oscillations with amplitudes ranging from 50 to 100 m, and sometimes even more [1][2][3].These waves are generated when there is a disturbance at the interface between water layers of varying densities.This disturbance can propagate over vast distances, causing the layers to oscillate in relation to their initial equilibrium state.One of the main causes of internal waves is the deterministic mechanism triggered by tidal flows interacting with bathymetric features like seamounts or continental shelves; other mechanisms may include instability of ocean currents in zones with strong shear flow (such as an ocean gulf), or directly by wind stress [4].This phenomenon is commonly referred to as an internal tide.As the tides flow over these underwater obstacles, they create disturbances that propagate as internal waves throughout the ocean.In addition to deterministic mechanisms, internal waves can also be excited by random events such as tsunamis or severe ocean storms.These powerful and unpredictable phenomena have the ability to generate significant disturbances in the ocean, leading to the formation of large-scale internal waves.
Extensive research has been conducted to investigate weakly nonlinear models that depict the behavior of internal waves.Notably, the Korteweg-de Vries (KdV) equation and the Gardner equation have received significant attention in this regard [4][5][6].The KdV and Gardner equations are applicable to shallow water scenarios, while the Intermediate Long Wave (ILW) equation is suitable for fluids of finite depth [7][8][9].Additionally, the Benjamin-Ono (BO) equation pertains to deep water dynamics [7,[10][11][12][13][14][15][16][17].These well-established models exhibit intriguing features, such as the existence of periodic and solitary wave solutions that endure over time.However, it is crucial to acknowledge that these model equations have inherent limitations that constrain their applicability to more generalized problems.
The widely recognized Benjamin-Ono (BO) equation [7,11] is frequently employed to examine the behavior of a perturbed interface between two inviscid fluids.This interface is characterized by a flat rigid lid and infinite depth.In this equation, various parameters and functions play important roles.The thickness of the upper fluid layer with density ρ 1 is represented by h 1 , while ρ 2 denotes the density of the lower fluid.The ratio of densities between the lighter upper layer and the heavier lower layer is denoted as ρ r = ρ 1 /ρ 2 < 1.The linear speed c 0 is [7,11] where g is the acceleration of gravity.More details of the geometry of the problem are depicted in Figure 1.The elevation of the interface in the position y and time τ is denoted by η(y, τ), and H denotes the Hilbert transform defined as [7,11] H[η(y,     < l a t e x i t s h a 1 _ b a s e 6 4 = " L B T T t s w S D L S q S 7 n 5 The BO equation is used to model several problems such as ship wakes, flow of water over rocks, the formation of storms in the ocean, and even atmospheric problems such as atmospheric flows encountering obstacles [18,19].These problems are typically approached as deterministic, where the initial state of the perturbed interface is precisely known, and the evolution of the phenomenon over time can be computed using the BO equations.However, in many cases, the scarcity of data for initializing deterministic models presents a challenge.In such instances, assuming random initial states becomes preferable, allowing for the study and formulation of statistically reliable predictions [20][21][22][23][24]. Analogously, the medium in which a wave propagates can also be considered random [25][26][27][28].This approach is commonly employed to investigate rogue waves, which are characterized as large waves that abruptly emerge within the wave field.Rogue waves are characterized by their exceptionally high wave heights and their abrupt appearance.The physical mechanisms behind their formation can be attributed to two primary factors: dispersion compression and geometric focusing.Dispersion compression arises due to the varying propagation speeds of waves with different wavelengths.In essence, waves of different lengths move at different speeds, causing them to temporarily converge and create a rogue wave.Geometric focusing, on the other hand, occurs when multiple wave systems converge from different directions, amplifying their combined effect and giving rise to rogue waves.Additionally, changes in environmental conditions, such as variable water depth and strong, non-uniform currents, can contribute to the concentration of energy within surface waves [29][30][31].Rogue waves have been extensively catalogued, and they have been shown to pose a significant threat to coastal areas, often causing substantial damage.To ensure the safety of both small boats navigating in shallow waters and people on the shoreline, in-depth investigations into rogue wave phenomena, including the study of their lifetimes (as indicated by Didenkulova and Pelinovsky [32]), are essential.Rogue waves have been part of maritime folklore for centuries, but their systematic measurement only began in the mid-20th century.This measurement process involves recording the displacement of the water-air interface at specific measurement points, utilizing various equipment such as floating moored buoys, laser altimeters mounted on platforms, or overpressure sensors located on the seabed.The first digitally documented instance of a rogue wave, famously known as the "Draupner wave," was recorded on 1 January 1995, at the Draupner platform in the North Sea.This rogue wave reached an astounding maximum wave height of 25.6 m and caused minor damage to the platform, which was situated high above sea level.Understanding the dynamics of rogue waves is crucial for enhancing maritime safety and improving the design of marine platforms.Although rogue waves are typically associated with large surface waves, they can also manifest within internal waves [5].However, detecting these internal rogue waves with the naked eye can be challenging, and they occasionally induce distinct surface responses.In general, internal waves exhibit higher amplitudes compared to surface waves.A rogue internal wave can unexpectedly displace less dense water beneath a neutrally buoyant submarine, causing it to descend to depths beyond the hull's pressure capacity.These waves maintain well-defined patterns as they propagate along the interface, in contrast to the more chaotic nature of turbulence.
In the context of surface waves, the study conducted by Pelinovsky and Sergeeva [22] focused on numerically examining the evolution of an initially random wave field with a Gaussian spectrum shape using the KdV equation.Their findings revealed that the irregular wave field, resulting from the wave evolution within the KdV model, does not adhere to Gaussian statistics.Instead, its statistical properties are influenced by the Ursell parameter, which represents the ratio of nonlinear effects to dispersion.Consequently, the wave field becomes asymmetric, exhibiting sharper crests and leading to a positive third moment.Importantly, the study demonstrates the existence of a steady state for statistical characteristics, including skewness, kurtosis, distribution functions, and spectral density.Through computations, it was observed that both statistical moments and distribution functions evolve until they reach a certain bound level.A similar effect is observed in the evolution of a random wave spectrum.In a related study, Didenkulova et al. [20] conducted direct numerical simulations of nonlinear wave evolution within the framework of the KdV equation for cases involving bimodal wave spectra models.They investigated the coexistence effect of an additional wave system on the evolution of wave statistical characteristics, spectral shapes, and the resulting equilibrium state.Furthermore, the study demonstrated that the presence of a low-frequency spectral component leads to more asymmetric waves with more extreme statistical properties.
The objective of this study is to investigate the behaviour of the evolution of a random internal wave field through the BO equation.As our model is deterministic, we conduct numerous simulations with randomly generated initial conditions with identical statistical characteristics.More precisely, our approach involves assuming that the initial conditions possess a power spectrum that adheres to a Gaussian distribution.This assumption serves as our initial step in unraveling the intricate nature of nonlinear dynamics, laying the foundation for the analysis of internal wave turbulence.Notably, in cases of more complex stratification, the Garrett-Munk spectrum typically constitutes the fundamental spectrum for internal waves.However, in scenarios featuring a two-layer stratification, our selected spectrum becomes a more suitable choice due to its relevance and appropriateness.For the mechanism of turbulence on internal waves, the readers are referred to [33][34][35][36].Our approach is related with a case of integrable turbulence of internal waves because the BO equation is an integrable model.Through a series of numerical simulations, we investigate spectrum evolution and stabilization, the third and the fourth statistical moments of the random wave field, and the distribution function for the crest amplitudes.In addition, we pay special attention to the formation of "rogue internal waves", finding the so-called "three sister" and other extensions.
For reference, this article is organized as follows: The BO equation is presented in Section 2. In Section 3, we describe the numerical methods, and the results are presented in Section 4. The final considerations are presented in Section 5.

The Benjamin-Ono Equation
Our study aims to investigate the underlying cause of the dominant dynamic, focusing on the nonlinearity of dispersion.To achieve this, it is convenient to utilize dimensionless variables by introducing the BO equation.Equation ( 1) is typically initialized with [22] η(y, 0) = A s F(Kk 0 y) sin(k 0 y).
Here, A s represents a standard wave amplitude, F denotes the wave envelope with spectrum width K, and k 0 corresponds to the carrier wave number.Thus, we introduce dimensionless variables Notice that with these choices, the crests of the wave field represented by η align with the troughs of the wave field denoted by u and vice versa.Substituting the dimensionless variables (5) into Equation ( 1), we derive the dimensionless BO equation where U r represents the Ursell parameter, given by U r = A s ρ r /h 1 k 0 .This parameter governs the nonlinearity of the BO equation, with higher values of U r indicating predominantly nonlinear dynamics, while lower values of U r suggest predominantly linear dynamics.Additionally, the Ursell parameter depends on the thickness of upper layer: hence, variations on the thickness of the upper layer may change the wave regime significantly.
To incorporate randomness into the problem, we introduce a zero-mean random state modeled as a Fourier series with M harmonics.The initial condition is given by Here, S(k) represents the initial power spectrum, k i = i∆k, with ∆k being the sampling wave number, and ϕ i denotes a random variable uniformly distributed in the interval (0, 2π).The length of the initial realization is L = 2π/∆k.We assume that the initial power spectrum follows a Gaussian distribution: The wave characteristics are determined by several parameters: the dimensionless peak wave number is k 0 = 1, the spectral width is denoted as K = 0.18, and the relative energy is indicated by the values of Q.All scenarios share the same total energy of the waves, which is defined by the variance σ 2 = 0.25 Within the framework of the integrable Benjamin-Ono (BO) equation, an intriguing characteristic emerges: the variance remains constant throughout the motion, offering a stable reference point for analyzing its behavior.In our spectral domain, we opt for a dimension featuring 256 harmonics.This choice accommodates a gradual spectrum decay in the higher wave number (k) region, allowing us to capture a more realistic representation of the wave dynamics.It is worth noting that these specific parameter selections closely align with those employed by Pelinovsky and Sergeeva [22] as well as Didenkulova et al. [20] in their respective studies.By adopting these consistent parameters, we facilitate a robust basis for comparing our research findings with the results presented in the subsequent sections, thus enhancing the overall comprehensibility and validity of our analysis.

Numerical Methods
In this section, we present the numerical method used to compute solutions of the BO Equation ( 6).Equation ( 6) is solved numerically by employing a Fourier pseudospectral method that incorporates an integrating factor, a technique akin to the one presented by Trefethen in [37].This integrating factor effectively resolves the linear component of Equation ( 6), thereby circumventing numerical instabilities stemming from the dispersive term.It is worth noting that for a given function f (ξ), we can compute the nonlocal operator H using the formula where F −1 denotes the Fourier inverse transform.
Rewriting Equation ( 6) in the Fourier frequency space, it takes the form The integrating factor is defined as Defining Û(k, t) = E(k, t) û(k, t) and replacing it in Equation (10), we obtain the equation where F denotes the Fourier transform, and F −1 is its inverse.It is important to note that once we have computed Û(k, t), we can recover the original function u as The solution of Equation ( 11) is performed within a periodic computational domain [0, L] using a uniform grid consisting of an even number of points denoted as N.The spatial points are discretized as x j = (j − 1)∆x, j = 1, 2, . . ., N, where ∆x = L/N, (12) and the frequencies are discretized as Fourier transforms and the Hilbert operator are computed using the Fast Fourier Transform (FFT), while spatial derivatives are evaluated spectrally following the methods outlined by Trefethen in [37].The time advance is computed using the Runge-Kutta fourth-order method (RK4) with time step ∆t.The parameters used in the following simulations are L = 2π/∆k, where ∆k = 0.023, and there is a uniform grid containing N = 2 12 points.For the time evolution of the equation, we utilize the classical fourth-order Runge-Kutta method with discrete time steps of size ∆t = 0.005.For a more detailed resolution study and comprehensive understanding of a similar numerical method, readers are referred to the work of Flamarion et al. [38].Furthermore, it is noteworthy that the reliability of the identical numerical method has been confirmed in a recent study by Flamarion and Pelinovsky [39].In this context, we use the fact that the BO Equation ( 6) conserves both the total mass (m(t)) and the momentum (p(t)) to verify the accuracy of the chosen numerical method.The total mass is defined as and the momentum is defined as Numerical simulations are controlled by retaining of the first and second moments with precision of machine and 10 −9 , respectively.Figure 2 illustrates the conservation of mass and momentum over time for different simulations.As observed, both quantities remain conserved throughout the simulations.Notice that for the mass conservation, we have machine precision.

Wave Field
The evolution of the wave record is depicted in Figure 3 at various time instances.Initially, the wave field exhibits a relatively narrow spectrum, resulting in smooth group structures of the waves.However, as time progresses, the wave profile becomes asymmetric, transitioning from smooth crests to sharp ones.This signifies an increase in the skewness of the wave field beyond its initial value, indicating a departure from a symmetric distribution, which will be discussed later.Moreover, the chosen parameter U r = 6.7 corresponds to a significant nonlinearity, leading to the emergence of large-amplitude waves during the interaction.To comprehend the role of the Ursell parameter in the dynamics, we analyze the trajectory patterns displayed in Figure 4. On the time-space plane, we observe that for higher values of U r , the wave propagation results in the formation of prominent peak amplitudes, as indicated in Figure 4.The distribution of crest amplitudes presented in Figure 5 shows that nonlinearity is essential for the existence of waves with large amplitudes: in particular, in the formation of rogue waves, which will be discussed below.In recent years, there has been a growing fascination with rogue waves, an intriguing type of nonlinear wave also known as freak waves.These waves have captured the attention of researchers from various scientific disciplines due to their unique characteristics.Initially observed in the deep ocean, rogue waves exhibit abnormal behavior, with amplitudes two to three times higher than the surrounding waves.What makes them particularly captivating is their sudden formation, seemingly appearing out of nowhere.
Rogue waves, characterized by their sudden and exceptionally high amplitudes, have emerged as a compelling subject of extensive investigation across diverse scientific disciplines.They have piqued the interest of researchers in fields such as oceanography [29,30], where understanding their behavior is critical for maritime safety and coastal protection.Furthermore, rogue waves have made their presence felt in areas as diverse as optical fibers [40][41][42], impacting the reliability of communication networks and signal transmission.In the realm of quantum physics, rogue waves have found relevance in Bose-Einstein condensates [43][44][45], offering insights into the behavior of ultracold atomic systems.Even financial markets have not remained untouched, with rogue-wave-like events manifesting themselves in market fluctuations and risk assessment [46].These broad-ranging applications underscore the pervasive influence of rogue wave phenomena in both natural and man-made systems.Mathematically, the criterion for the occurrence of freak waves can be described by the equation [30] where A f r represents the height of the freak wave, and H s represents the height of the "significant" wave field.The value of H s is determined by averaging one-third of the largest waves observed in the given context, such as in the field of oceanology.Figure 6 (left) displays different simulated wave fields at different time records, where prominent peaks significantly surpass the surrounding waves, representing various types of freak waves.Zooming in on the left panels, Figure 6 (right) provides a closer look at their distinctive appearance.For instance, the top-right image showcases a freak wave of substantial amplitude known as the "two sisters" and "three sisters".The available data, along with theoretical studies, provide ample evidence for the existence of rogue waves in diverse shapes.For example, these types of waves were also documented and their lifetimes estimated in a recent study by Didenkulova and Pelinovsky [32], where they were examined within the context of wind waves utilizing the KdV equation.It is worth mentioning that series of big waves with more than four peaks were also observed.The occurrence of freak waves is strongly influenced by the parameter U r , with a higher frequency expected in wave fields characterized by stronger nonlinearity and smaller values of U r .
By studying and understanding the characteristics and behavior of these rogue waves, researchers aim to shed light on their formation mechanisms and develop methods to predict and mitigate their potentially hazardous effects.

Spectra
The evolution of the spectrum is studied across various U r values, ranging from 0.1 (representing nearly linear progression) to 6.7 (indicating a highly nonlinear wave behavior).To ensure reliable statistical outcomes, 300 realizations are averaged over different time periods.As anticipated, the presence of nonlinearity triggers a transformation in the spectrum, leading to its widening and eventual convergence into a stable state (refer to Figure 7).The nature of this stable state, determined by the Ursell parameter, displays an asymmetrical shape with a noticeable shift of energy towards lower frequencies, commonly known as the spectrum downshift effect.In cases of larger Ursell parameter values (as illustrated in (refer to Figure 7) (top)), the spectral density becomes more evenly distributed across smaller wave numbers (k-values).The prominence of spectrum flatness intensifies under strong nonlinearity (U r = 6.7), indicating a higher energy level in the wave field and an increased significance of nonlinear effects.Consequently, a broader range of frequencies is necessary to accurately represent the wave spectrum.This inclination towards spectrum flatness aligns with the concept of statistical equilibrium in the absence of external inputs or outputs.These findings corroborate and strengthen the outcomes reported by Pelinovsky and Sergeeva [22] within the framework of the KdV equation.

Statistical Moments
To gain a deeper understanding of the wave fields interactions described in the previous subsection, we focus on examining four specific integrals, corresponding to four statistical moments More precisely, we focus on two statistical quantities that characterize the wave spectrum.The kurtosis excess (κ) and the skewness (ς) are defined as The kurtosis is a measure that indicates the tail heaviness of the spectrum.In simpler terms, it quantifies the degree of peakedness in the distribution and characterizes the influence of large waves on the overall distribution.A positive kurtosis implies a substantial contribution from large waves.Skewness, on the other hand, measures the asymmetry of the spectrum relative to the mean.Specifically, it represents the statistical measure of vertical asymmetry in the wave field, with its sign indicating the ratio of crests to troughs.A positive skewness signifies that crests are larger than troughs.
The behavior of statistical moments provides insights into the presence of a stationary state and its evolution over time.During a transition period of approximately 10 to 20 nonlinear time units, both moments of the wave field converge towards nearly constant values (see Figure 8).Regardless of the conditions, the skewness of the wave field remains positive, indicating that positive waves (crests) have larger amplitudes compared to negative waves (troughs).Furthermore, the asymptotic value of skewness increases as the Ursell parameter rises.
When the Ursell number is 2, the kurtosis shows oscillatory patterns centered around zero.This scenario resembles the findings reported by Onorato et al. [47] and Tanaka [48] in the context of deep water and the findings of Dutykh [21] for different depths using the Euler equations.In their studies, it was demonstrated that in deep water, the calculated kurtosis values exhibit oscillations around zero.For highly nonlinear random wave processes, some values surpass zero, indicating an increased likelihood of encountering large-amplitude waves, including freak waves.On the other hand, for smaller values, the kurtosis becomes negative, suggesting a lower probability of encountering such extreme events compared to what is expected for Gaussian processes.Conversely, under strong nonlinearity, the asymptotic value of kurtosis exceeds zero, indicating a higher probability of encountering large waves, with significant contributions from smaller waves within the wave field.

Crest Distribution
In this section, we compare the exceedance probability distributions of wave crests for different U r values with the Rayleigh distribution of amplitudes for a narrow-band Gaussian process The instant distributions of wave crests, denoted as u max (representing local maxima of u between consecutive zero-crossing points), are computed for a dataset comprising approxi-mately 150,000 waves.Initially, the wave phases are randomly distributed, and the crest distributions, denoted as P(u max ), exhibit a reasonable agreement with the Rayleigh law for narrow-banded waves (refer to Figure 9 (left)).The probability distributions evolve over time.In the quasi-equilibrium stage (Figure 9 (right)), the wave crest distributions are compared with the theoretical predictions provided by Equation (19).Across all five cases, the distributions exhibit qualitatively similar behavior in relation to the reference Rayleigh curve for small amplitude waves.The asymptotic distribution surpasses the Rayleigh distribution, indicating an increased probability of encountering higher wave crests.Qualitatively, the shape of the amplitude distribution function aligns with the behavior of skewness and kurtosis depicted in Figure 8.The positive waves exhibit larger amplitudes than the negative waves, as indicated by the skewness, while the kurtosis highlights the significant contribution of small waves to the overall distribution.

Conclusions
In this work, we have explored the statistical properties of a random internal wave field with initially identical statistical properties.Based on the numerical simulations, our findings indicate that under strong nonlinearity, different types of freak waves might exist, including the intriguing three sisters phenomenon.We have observed that the spectral evolution of the wave field reaches a steady state across all Ursell parameters, but with increased nonlinearity, the spectral density becomes more uniformly distributed.This outcome aligns with previous research conducted by Pelinovsky and Sergeeva [22].Furthermore, we have computed the kurtosis and skewness of the wave field.The positive skewness values indicate the asymmetry of the wave field, with sharper crests.As for the kurtosis, strong nonlinearity (large Ursell parameter values) demonstrates a significant contribution from larger waves in the wave dynamics.In contrast, strong dispersion (small Ursell parameter values) leads to kurtosis oscillating around zero or becoming negative.Additionally, we have examined the wave crest distribution for various Ursell parameter values, comparing it with the Rayleigh distribution.Our analysis reveals that as the Ursell number increases, the probability distribution function deviates slightly from the theoretical Rayleigh distribution, indicating differences in the distributions.

< l a t e x i t s h a 1 _ b a s e 6 4 =
" g 0 7 N g 3 K v I 9 F W M O X z w L W R 8 G A a o r k = " > A A A B + X i c b V D L S s N A F J 3 U V 6 2 v q E s 3 g 6 1 Q Q U p S 8 L E R C m 5 c V r A P a E K Z T C f t 0 M m D m Z tC D P 0 T N y 4 U c e u f u P N v n L Z Z a O u B C 4 d z 7 u X e e 7 x Y c A W W 9 W 0 U 1 t Y 3 N r e K 2 6 W d 3 b 3 9 A / P w q K 2 i R F L W o p G I Z N c j i g k e s h Z w E K w b S 0 Y C T 7 C O N 7 6 b + Z 0 J k 4 p H 4 S O k M X M D M g y 5 z y k B

h 1 <
l a t e x i t s h a 1 _ b a s e 6 4 = " 4 7 H a 1 n Y n i / h n T T Z d d + 9 0 H I S S x v Q = " > A A A B 7 3 i c b V D L S g N B E O z 1 G e M r 6 t H L Y C J 4 C r s B H 8 e A F 4 8 R z A O S J c x O Z p M h s 7 P r T K + w h P y E F w + K e P V 3 v P k 3 T p I 9 a

1 < l a t e x i t s h a 1 _ b a s e 6 4 =
r l e y e M o w C m c w Q V 4 c A 1 1 u I M G N I G B h G d 4 h T f n 0 X l x 3 p 2 P R e u a k 8 + c w B 8 4 n z 9 2 W I + C < / l a t e x i t > " f

⇢ 1 < l a t e x i t s h a 1 _ b a s e 6 4 =
" v e 6 u 7 O B N w l I l d B n Z / z g k Z p O 1 l Z 4 = " > A A A B 7 3 i c b V D J S g N B E K 1 x j X G L e v T S m A i e w k z A 5 R j w 4 j G C W S A Z Q k + n J 2 n S y 9 j d I 4 Q h P + H F g y J e / R 1 v / o 2 d Z A 6 a + K D g 8 r l e y e M o w C m c w Q U E c A 1 1 u I M G N I E A h 2 d 4 h T f v 0 X v x 3 r 2 P R e u a l 8 + c w B 9 4 n z / 8 y I 8 y < / l a t e x i t > ⇢ 2 r / 4 i b / 4 b t 2 0 O 2 v p g 4 P H e D D P z / F h w b V z 3 2 y m s r W 9 s b h W 3 S z u 7 e / s H 5 c O j t o 4 S x b D

Figure 2 .
Figure 2. The conserved quantities of BO equation for different simulations.On the left is the mass conservation and on the right is the momentum.

Figure 3 .
Figure 3.The evolution of the wave field at different times with U r = 6.7.

Figure 4 .
Figure 4.The evolution of the wave field for different values of the parameter U r .From (top) to (bottom) and (left) to (right) U r = 6.7,U r = 4.0, U r = 2.0 and U r = 1.0.

Figure 5 .
Figure 5. (Left): Maximum of wave amplitudes over time.(Right): Distribution of maximum crest amplitudes over 300 realizations for different values of U r .

Figure 7 .
Figure 7. Averaged evolution of spectra of the wave fields at different times.From (top) to (bottom) and (left) to (right) U r = 6.7,U r = 4.0, U r = 2.0, U r = 1.0,U r = 0.5 and U r = 0.1.

Figure 8 .
Figure 8. Temporal evolution of the skewness and kurtosis for different values of the parameter U r averaged over 300 realizations.On the right, the zooming of the figures on the left.

Figure 9 .
Figure 9.The exceedance probability distributions for wave crests at the initial moment (t = 0) and the asymtotic distribution (t = 100) (right) for different Ursell numbers over 300 realizations.The black solid line with circles corresponds to the Rayleigh distribution of the narrow-band Gaussian process.