Single-spectrum prediction of kurtosis of water waves in a non-conservative model

TL;DR

This paper investigates the relationship between kurtosis, bandwidth, and other spectral properties of water waves under wind forcing, proposing a simple spectral measure to predict rogue wave likelihood in non-conservative models.

Contribution

It introduces a robust quadratic relation between kurtosis and bandwidth, enabling single-spectrum assessment of rogue wave potential in water waves with wind forcing.

Findings

Kurtosis and BFI are related quadratically but depend on forcing and damping.

A simple quadratic relation exists between kurtosis and spectral bandwidth.

The evolution of kurtosis can be predicted using bandwidth and damping coefficients.

Abstract

We study statistical properties after a sudden episode of wind for water waves propagating in one direction. A wave with random initial conditions is propagated using a forced-damped higher order Nonlinear Schr\"odinger equation (NLS). During the wind episode, the wave action increases, the spectrum broadens, the spectral mean shifts up and the Benjamin-Feir index (BFI) and the kurtosis increase. Conversely, after the wind episode, the opposite occurs for each quantity. The kurtosis of the wave height distribution is considered the main parameter that can indicate whether rogue waves are likely to occur in a sea state, and the BFI is often mentioned as a means to predict the kurtosis. However, we find that while there is indeed a quadratic relation between these two, this relationship is dependent on the details of the forcing and damping. Instead, a simple and robust quadratic relation does exist between the kurtosis and the bandwidth. This could allow for a single-spectrum assessment of the likelihood of rogue waves in a given sea state. In addition, as the kurtosis depends strongly on the damping and forcing coefficients, by combining the bandwidth measurement with the damping coefficient, the evolution of the kurtosis after the wind episode can be predicted.

Tables

Table 1. Table 1: Parameters in numerical experiments

$f_0=1.667$ Hz	${\sigma }_k\left(0\right)=0.2$	$\epsilon \left(0\right)=0.08$	$L=60{\lambda }_0$
$n_x=2^{10}$	$G=2$	$T_{\text{on}}=5T_0$	$N_{\text{sim}}$ = 250

Table 2. Table 2: Parameters of wind forcing r 𝑟 r and damping d 𝑑 d for each of the 18 numerical experimental conditions

#	$d$	$r$	#	$d$	$r$
1	0.01	0.04	2	0.01	0.2
3	0.01	1	4	0.01	3
5	0.01	5	6	0.02	1
7	0.025	0.2	8	0.025	1
9	0.025	3	10	0.025	5
11	0.05	0.2	12	0.05	1
13	0.05	3	14	0.05	5
15	0.1	0.2	16	0.1	1
17	0.1	3	18	0.1	5

Equations

$$B=2\sqrt{2}\frac{{\sigma }_k(t)}{\epsilon (t)}$$

$$\frac{\partial a}{\partial t}+\frac{{\omega }_0}{2k_0}\frac{\partial a}{\partial x}=i\frac{{\omega }_0}{8k_0^2}\frac{{\partial }^{2}a}{\partial x^2}+\frac{1}{2}i{k}_0^2{\omega }_{0}a|a{|}^2$$

$$\frac{1}{2}\Gamma a-2k_0^2\nu a+\frac{3i}{4k_0}\Gamma \frac{\partial a}{\partial x}-4i{k}_0\nu \frac{\partial a}{\partial x}$$

$$+i{k}_{0}a\frac{\partial \bar{\varphi }}{\partial x}-\frac{3}{2}{k}_0{\omega }_0|a{|}^2\frac{\partial a}{\partial x}-\frac{1}{4}{k}_0{\omega }_0{a}^2\frac{\partial a^{\ast }}{\partial x}$$

$$+\frac{{\omega }_0}{16k_0^3}\frac{{\partial }^{3}a}{\partial x^3}-i\frac{5{\omega }_0}{128k_0^4}\frac{{\partial }^{4}a}{\partial x^4}+2\nu \frac{{\partial }^{2}a}{\partial x^2}$$

$$\eta (x,t)=\mathrm{Re}\{a(x,t)\exp (i({\omega }_{0}t-k_{0}x)\}$$

$$T$$

$$X$$

$$A$$

$$\gamma$$

$$\frac{\partial \bar{\varphi }}{\partial x}=-\frac{{\omega }_0}{2}\mathcal{H}[|a_x{|}^2],\text{so that}\frac{\partial \bar{\Phi }}{\partial X}=-\frac{1}{\sqrt{2}}\epsilon \mathcal{H}[|A_X{|}^2]$$

$$\mathcal{F}[\mathcal{H}[u]]=-i\mathrm{sign}(\mathrm{k})\mathcal{F}[\mathrm{u}],$$

$$\begin{split}&\underbrace{i\frac{\partial A}{\partial T}+\frac{1}{2}\frac{\partial^{2}A}{\partial X^{2}}+|A|^{2}A=iA(r-d)}_{\text{Forced/Damped NLS}}\\ &-\varepsilon\underbrace{\frac{\partial A}{\partial X}\bigg{[}4d-3r\bigg{]}}_{\text{HOT}}\\ +&i\varepsilon\underbrace{\bigg{[}-6|A|^{2}\frac{\partial A}{\partial X}-A^{2}\frac{\partial A^{*}}{\partial X}+\frac{1}{2}\frac{\partial^{3}A}{\partial X^{3}}-2iA\mathcal{H}[|A|^{2}_{X}]\bigg{]}}_{\text{Dysthe}}\\ +&\varepsilon^{2}\underbrace{\left[\frac{5}{8}\frac{\partial^{4}A}{\partial X^{4}}+4id\,\frac{\partial^{2}A}{\partial X^{2}}\right]}_{\text{HOT dispersion correction}}\end{split}$$

$$r$$

$$d$$

$$\omega$$

$${\omega }_r$$

$$\beta$$

$$\frac{k_{\max }}{k_0}=\frac{3\Gamma }{16\nu k_0^2}-1$$

$${\beta }_{\max }=\Gamma \left[\frac{9\Gamma }{128k_0^2\nu }-\frac{1}{4}\right]$$

$$\frac{k_{\max }}{k_0}=\frac{3r}{4d}-1,{\beta }_{\max }{t}_0=\frac{r}{2}\left[\frac{9r}{8d}-1\right]$$

$${\beta }_{\text{P}}=\frac{\delta \xi u_{\ast }^2}{c_{\text{p}}{K}^2}k-4\nu k^2$$

$$\Gamma =\alpha \delta {\omega }_0{\left(\frac{u_{\ast }}{c_{\text{p}}}\right)}^2$$

$${\beta }_{\text{P}}=\frac{\xi }{\alpha K^2}\frac{\Gamma }{k_0}k-4\nu k^2$$

$${\beta }_{\text{P}}=\frac{\Gamma }{2}-2\nu k_0^2+\left[\frac{\Gamma }{2k_0}-4\nu k_0\right]k-2\nu k^2$$

$$\epsilon (t)=k_0\text{rms}[\eta (x,t)]=\frac{k_0}{\sqrt{2}}\text{rms}[|a(x,t)|]$$

$$\text{rms}[\eta (x,t){]}^2={\int }_0^{\infty }{S}_{\eta }(k,t)dk$$

$${\sigma }_k(t)=\left\{\frac{1}{k_0}\sqrt{\frac{{\int }_{-2k_0}^{2k_0}{\left(k-k_{\text{m}}\right)}^2{S}_{a}dk}{{\int }_{-2k_0}^{2k_0}{S}_{a}dk}}\right\}$$

$$N=\frac{1}{L}{\int }_0^L|A{|}^{2}dx$$

$$K(t)=\left\{\frac{⟨(a-\bar{a}{)}^4⟩}{⟨(a-\bar{a}{)}^2{⟩}^2}\right\}-3.24$$

$${\kappa }_{\text{m}}(t)=\frac{P}{N},P=\frac{i}{2L}{\int }_0^L\left(A{A}_X^{\ast }-A_X{A}^{\ast }\right)dx$$

$$K=\frac{\pi }{\sqrt{3}}{B}^2$$

$$i{A}_x=\beta \frac{{\partial }^{2}A}{\partial t^2}+\alpha |A{|}^{2}A$$

$$K(x)=K(x_0)+\frac{1}{⟨N⟩}\left(\Omega (x{)}^2-\Omega (x_0{)}^2\right)$$

Full text

††preprint: APS/123-QED

I Introduction

In gravity water waves, sensitivity to initial conditions is so strong that after a typical time of $\mathcal{O}(10^{3})$ wave periods, no information is left on initial conditions [1], even in the absence of irreversible processes such as wave breaking. Therefore, deterministic approaches to water waves cannot give a complete picture, and a complementary statistical approach is needed.

In general, when studying statistics of a sea state or laboratory experiment, homogeneous and stationary theory is assumed. This analysis needs to be extended in order to include non-conservative systems. In this study, we examine the effects of wind forcing. During an episode of wind, the system is out of equilibrium due to its growth in energy. After the wind, it can be considered closed again, provided that the dissipation timescale is much longer than the wave period ($T_{\text{diss}}\gg T_{0}$). In water waves, viscous dissipation is part of the physical system and always present. A steady state, apart from small dissipation, can be reached after non-resonant effects have saturated.

We consider unidirectional, i.e., long-crested waves, where only non-resonant interactions are allowed. The dominant nonlinear phenomenology stems from modulational (or Benjamin-Feir) instability (MI). This is a simplification, since realistic winds are not 1D, and seas are typically multi-directional in nature [2]. The directional spread of energy of the waves decreases the occurrence of rogue waves due to modulation instability [3, 4]. In addition, the effect of directional waves is to reduce nonlinear focusing related to MI [5, 6, 7, 8]. Moreover, short-crested seas have weakly non-Gaussian wave statistics [9].

The goal of this work is twofold. On the one hand it is to seek for characteristic signatures of wind forcing and dissipation on the wave statistics. From a deterministic point of view, it is known that the spectral mean, bandwidth and steepness are functions of dissipation and wind forcing [10], where the dissipation and forcing have opposite signatures. It is not known to date how this behavior translates to ensemble simulations. On the other hand, we aim at defining a proxy to predict kurtosis.

Kurtosis, the fourth moment of the wave height distribution, is seen as the main indicator for the presence of rogue waves, as it shows the ’fatness’ of the tails of the wave-height distribution. However, measuring the kurtosis with a reasonable precision requires measuring an ensemble of spectra. Thus, a first step to prediction is finding a parameter that would be measurable from one single spectrum and could be used to estimate the kurtosis, and therefore the likelihood of rogue waves in a given sea state, whether numerical or experimental. Subsequently, we will relate this parameter to the influence of damping on the spectrum in order to predict the evolution of kurtosis at a future time in the swell evolution.

The main candidate as an estimation parameter for kurtosis has been the Benjamin-Feir index (BFI), defined as [11]

[TABLE]

where $\sigma_{k}$ is the bandwidth, and $\varepsilon$ the characteristic wave steepness. Note there is a factor 2 difference with the definition of the BFI for a time-like NLS [12]. For narrow-band waves propagating in one direction kurtosis is deemed to depend quadratically on the BFI, as derived by  [12]. For narrow-band waves, the asymptotic solution of Nonlinear Schrödinger equation (NLS) is fully determined by the BFI value. For the Dysthe or Zakharov equation [13, 14] the BFI is still relevant if the spectrum is sufficiently narrow. Here we investigate whether the quadratic relation proposed by Janssen [12] is valid after an episode of wind.

In addition to the BFI, the bandwidth too is quadratically proportional to the kurtosis for a conservative NLS [15]. In addition, correlations between bandwidth and various statistical variables, including kurtosis, were observed in wave tank experiments without wind forcing [16]. Here, we show that a quadratic relation still holds in our non conservative and higher-order model. Finally, we show how the different terms in the evolution equation affect the evolution and statistics of the wave spectrum.

We extend the analysis conducted in the framework of the NLS in [11]. In contrast to this study, we include high-order nonlinear and dispersion terms of the modified NLS [13]. And, at the same order in steepness as the nonlinear terms, we include effects of wind forcing and dissipation [10]. This allows for the description of broader spectra. Indeed, after an episode of wind the spectrum is broadened and NLS is not sufficient to describe the evolution of the waves. The Dysthe modification of the NLS is better suited for describing broad-banded sea states, giving results comparable to models without bandwidth constraints such as high-order spectral methods applied to the truncated Euler equations [7].

For our investigation, we run a deterministic model, with as initial condition a Gaussian spectrum with random phases. Waves are propagated in one direction. The first part of the propagation serves to let the nonlinear interactions of the system reach equilibrium. In this stage no wind forcing is present, only dissipation. Then, the wind episode occurs, where wind and dissipation act simultaneously. Afterwards, the wave energy decreases again due to dissipation. To obtain statistics, an ensemble of such simulations is performed, for different wind and dissipation strengths.

II Model

II.1 Model equations

Our model to propagate the envelope $a(x,t)$ is the space-like version of the forced/damped-Modified NLS (MNLS) equation [10], with higher order terms added to suppress unphysical resonances and improve numerical stability. See Appendix A for a full derivation and notations. In dimensional coordinates, the envelope equation reads

[TABLE]

Here, $x$ is the propagation direction, $k_{0}$ the carrier wave-number, $\omega_{0}$ the carrier frequency, $\nu$ the viscosity, $\Gamma$ the wind input. The surface elevation (without bound modes) can be calculated as

[TABLE]

We define the following adimensional variables:

[TABLE]

The wave-induced mean current $\frac{\partial\bar{\phi}}{\partial x}$ can be written in terms of the Hilbert transform $\mathcal{H}$ of the wave envelope as

[TABLE]

where the Hilbert transform is defined as

[TABLE]

$\mathcal{F}$ being the Fourier transform. Eq. (II.1) then reduces to the dimensionless form:

[TABLE]

where we set the forcing and damping coefficients, $r$ and $d$, in analogy to Ref. [11]:

[TABLE]

We observe that terms proportional to $\varepsilon$ are of high-order with respect to the conventional NLS (on the left-hand side), that is derived from the Euler equations in the incompressible irrotational limit through the multiple-scale method at third-order in steepness $\varepsilon$. In the present model, damping and forcing terms, represented by factors proportional to $d$ and $r$, respectively, appear both the leading- and higher-order level (proportional to $\varepsilon$).

The two terms proportional to $\varepsilon^{2}$ represent higher-order corrections to the dispersion. The first term with 4th-order derivative $\frac{5}{8}\frac{\partial^{4}A}{\partial X^{4}}$ has been used in [17] to eliminate numerical instabilities that are due to the appearance of high harmonics. As discussed in the next section, the second term, $4id\,\frac{\partial^{2}A}{\partial X^{2}}$ allows the spectrum to be cut at high wave-numbers, as required by the presence of viscosity.

II.2 Linear stability analysis

We insert an eigenmode of the envelope $a$ of the form: $a~{}=~{}\hat{A}e^{i(\omega t-kx)}$, with $\omega\in\mathbb{C}$ and $k\in\mathbb{R}^{+}$ into the linearized forced/damped-MNLS equation (the linear part of Eq. (LABEL:eqn_ModelAdim)), yielding

[TABLE]

Note that this linear dispersion relation can also be obtained from Eq. A. The real $\omega_{r}$ and imaginary $\beta$ parts are, respectively:

[TABLE]

Note that only the non-conservative wind and viscosity terms have an influence on the growth rate $\beta$. In particular, we find that the most unstable mode is

[TABLE]

and the corresponding maximum growth rate, calculated in $k_{\mathrm{max}}$, reads

[TABLE]

Note that the higher order wind contribution is proportional to $k$ in Eq. (II.8), which implies an asymmetric growth of positive modes (implying the wind has a direction) until they reach $k_{\mathrm{max}}$. At this point, the $5^{\text{th}}$ order viscosity contribution proportional to $-k^{2}$ becomes the dominant term, strongly damping the very high $k$ modes. A similar reasoning where the input of the wind is naturally bounded is given in [18], using nonlinear damping.

In terms of our two dimensionless parameters, $r$ and $d$, these conditions read:

[TABLE]

This predicted most unstable mode is in agreement with the most unstable positive mode of the simulations for the full equation, Eq. (LABEL:eqn_ModelAdim).

In the limit of small wind forcing ($r\rightarrow 0$), $k_{\mathrm{max}}\rightarrow-k_{0}$ and $\beta_{\mathrm{max}}\rightarrow 0$. It means that without energy input, the most unstable mode of the envelope is $-k_{0}$ (thus, the corresponding surface elevation mode is $-k_{0}+k_{0}=0$) with null growth rate, all other modes are damped. That is, the wave is damped to a flat surface.

In the limit of small viscosity ($d\rightarrow 0$), both $k_{\mathrm{max}}$ and $\beta_{\mathrm{max}}$ go to infinity and the growth rate $\beta$ becomes an unbounded linear function of $k$. Notice that, in general $\beta(-k_{0})=-\frac{1}{4}\Gamma<0$. This is attributed to the Taylor expansion of the growth rate. It contributes to improve the numerical stability of our solver, as it prevents unstable growth for this mode and all modes below.

II.3 Comparison with Plant growth rate

Plant and Wright [19] derived the following linear growth rate for the evolution of the intensity:

[TABLE]

where $\delta$ is the density ratio between air and water, $c_{p}$ the phase velocity, $u_{*}$ the friction velocity and $K=0.41$ the Von Kármán constant. The parameter $\xi$ in Eq. (II.14) was empirically estimated in [19] to $\xi\approx 3.3$. This value turns out to be only slightly wind-speed dependent. By inserting the Miles growth rate:

[TABLE]

where the empirical constant $\alpha\approx 32.5$ [20]. The Plant growth rate becomes

[TABLE]

The factor in the first term is approximately $\xi/\alpha K^{2}\sim 0.6$. Shifting the wave-number by $k\rightarrow k+k_{0}$, dividing the overall equation by a factor two (to move from a growth rate for the wave intensity to one for the wave envelope), and setting $\xi/\alpha K^{2}\simeq 1$ gives:

[TABLE]

which is equal to the growth rate of Eq. (II.10), except for the coefficient of $\Gamma k/k_{0}$ ($1/2$ instead of $3/4$). Hereby we show that our model automatically generates the viscous correction to the Miles growth rate that is included in the Plant formula, albeit with a slightly different factor for the forcing term, attributed to the approximation inherent to our asymptotic expansion.

III Numerical simulations

Table 1 lists the parameters used in the simulations. Since our simulations extend the work of [11], we adopt the same parameters where possible to facilitate comparison. The frequency of the carrier wave, $f_{0}=1.667$ Hz, is the same for all simulations, and corresponds to $k_{0}\approx 11$ rad $\text{m}^{-1}$, where the dispersion relation $k_{0}=\left(\frac{2\pi f_{0}}{g}\right)^{2}$ has been used.

In [11] the calculation is based on the surface elevation consisting only of free waves. In the case of unidirectional waves, indeed only free waves are dominant [21]. When considering only free waves, the calculation of $\eta$ (Eq. (II.2)) simply shifts the envelope spectrum by $+k_{0}$. However, doing so imposes a boundary on the spectrum at $k=0$ for the surface elevation, or $k=-k_{0}$ for the envelope. Due to the wind forcing, our spectrum is slightly broader than the interval $[-k_{0},k_{0}]$, as can be seen from Figure 11 in Appendix C. Moving to the surface elevation would therefore introduce asymmetries. Therefore, we base our calculation on the envelope, which in addition is the variable we solve our model for (Eq. (LABEL:eqn_ModelAdim)).

Secondly, [11] limits the calculation of the bandwidth to a range $k\in[0,2k_{0}]$ for the surface elevation, or to the interval $[-k_{0},k_{0}]$ for the envelope. In order to be consistent with all other quantities that are subsequently calculated we not only limit the calculation of the bandwidth, but all other quantities too. Due to an increase in bandwidth due to wind forcing, and because our model is one order higher in steepness, we truncate the spectrum of the envelope to a wider interval: $[-2k_{0},2k_{0}]$. See Appendix C for further discussion on how results depend on various options of truncation.

The system [Eq. (LABEL:eqn_ModelAdim)] is given as initial condition a Gaussian shaped power spectral density, with random phases with uniform distribution between 0 and $2\pi$ independently and identically in each spectral bin. This initial condition is defined by the wave steepness and bandwidth. We define the characteristic wave steepness

[TABLE]

where rms refers to root mean square. The initial steepness$\varepsilon(0)=\frac{a_{0}k_{0}}{\sqrt{2}}=0.08$ in all simulations. The steepness can be linked to the initial power spectrum $S_{\eta}(k,0)$ through the relation

[TABLE]

evaluated at $t=0$. From the generated initial surface elevation $\eta(0)$ we now calculate the envelope $a(0)$, and its power spectrum $S_{a}$.

Following [11], the bandwidth $\sigma_{k}(t)$ corresponds to the variance of a Gaussian function, that is, the second moment of the distribution. Unlike [11], however, our model is not symmetric in the spectral domain, therefore we calculate the bandwidth with respect to the spectral mean $k_{m}$ instead of $k_{0}$. Indeed, the bandwidth, or, spectral width, is equal to the standard deviation of the distribution, which is always computed with respect to its mean:

[TABLE]

where $k_{m}(t)$ is defined for the envelope, as such $k_{m}(0)=0$. The initial bandwidth $\sigma_{k}(0)=0.2$. The computational length $L$ is set to $60$ times the carrier-wave wavelength, corresponding to a physical tank length of $L\approx 33$ m, and is discretized over $n_{x}=2^{10}$ equispaced grid points. The numerical scheme is based on the interaction picture with an adaptive time-step: the linear terms are solved in the Fourier space, and the nonlinear terms by means of an Embedded Runge-Kutta 4(3) scheme [22].

The Gaussian initial condition, or homogeneous wave field [12], is integrated for a long enough time that in the absence of dissipation or forcing, an equilibrium condition is reached. The nonlinear terms act on the Gaussian spectrum. In the absence of damping or forcing, the system equilibriates after a time $T=5\,T_{\text{nl}}\simeq 125\,T_{0}$. Here, $T_{0}=2\pi/\omega_{0}$ is the wave period, $T_{\text{nl}}=1/(\varepsilon^{2}\omega_{0})$ [11]. In our system, we give the modes the same time to re-organize, however, dissipation is active in this first part of the propagation. The wind is switched on at $T_{\text{on}}=5T_{\text{nl}}$, and the wind episode begins. The wind is turned off, when the wave amplitude has increased by a factor $G=2$ from point $T_{\text{on}}$. That is, $N(T_{\text{off}})=G^{2}N(T_{\text{on}})=4N(T_{\text{on}})$, where the norm or wave action

[TABLE]

is based on the envelope in Eq. (LABEL:eqn_ModelAdim). The input energy of the wind, determined by $r$, remains constant during the time when the wind is active.

A total of 18 different combinations of $r$ and $d$ were considered (see Appendix B). For each of these, $N_{\mathrm{sim}}=250$ realizations with random initial conditions were performed. The results of 6 cases are displayed and discussed in detail. We consider two different damping values: weak ($d=0.01$) and strong ($d=0.05$), and three different wind strengths: $r=0.2,1,3$. The other data-sets ($d=0.02,0.025,0.1$ and $r=5$) are used for calculation of fits when indicated.

IV Results

IV.1 Temporal evolution

Figure 1a shows the time evolution of the norm, or wave-action. From $T=0$ to $T_{\text{on}}$ the wave energy decreases, the rate depending on the damping coefficient $d$. During the wind episode ($T_{\text{on}}$ to $T_{\text{off}}$) the wave action increases exponentially. The higher the value of $r$, the faster the increase. The duration of the wind action for $d=0.05$, $r=0.2$ (solid blue line) is indicated by by red tick marks. From $T_{\text{off}}$ to the end of the simulation, the rate of energy decrease is again determined by $d$. The norm asymptotically reaches the same value for a given $d$, irrespective of $r$.

A similar pattern of increase and decrease that follows the norm is observed for the kurtosis, Figure 1b. The kurtosis is calculated as

[TABLE]

where $\bar{a}=\langle a\rangle$. And $3.24$ is the fourth standardized moment of the Rayleigh distribution. The kurtosis reaches high values during the wind forcing, and tends back to 0 at long times $T$ due to dissipation.

In [11] a distinction is made between fast ($r=5$) and adiabatic ($r=0.2$) pumping, and a higher final kurtosis is found for the latter. In contrast, we observe no quantitative difference between fast and slow forcing. The final kurtosis value is largely determined by the viscosity, where a stronger damping (higher $d$) gives a stronger decay rate of the kurtosis. In addition, when comparing the correlations between statistical quantities in Section IV.2, where values of $r=5$ are taken into account, a distinction between fast and slow pumping is not revealed.

The adimensional spectral mean for the envelope, $\kappa_{\text{m}}=k_{\text{m}}x_{0}$, is calculated as

[TABLE]

Figure LABEL:sub@fig:MeanTime shows the mean up-shifts during forcing and downshifts during damping as discussed in [10]. A similar trend, but in the sense of widening and narrowing is seen for the bandwidth (Figure 1d).

Since our propagation equation does not have a natural bound to the energy input, which in reality is of course provided by wave breaking, care must be taken not to go outside of the validity of the model. For all simulations, the characteristic steepness $\varepsilon$ (Figure LABEL:sub@fig:steepTime), never exceeds 0.25. We can therefore safely assume there are no wave breaking events [23]. More recent models, such as the compact Zakharov equation [24], do provide limitations to the wave growth and allows for the study of the behaviour of sharply peaked wave-forms and the inclusion of wave-breaking.

The BFI (Figure LABEL:sub@fig:BFITime), increases when there is forcing, and decreases when there is damping. However, for strong forcing, there seems to be an overshoot, as observed in [11], followed by a rapid decrease. The BFI follows a different pattern in time than the kurtosis, mean and width. The next section (IV.2) will show that while the bandwidth correlates well with the kurtosis, the steepness does not (Figure  LABEL:sub@fig:steepTime). Since the BFI is the ratio of these two quantities, this causes a deviation from the trend of the temporal evolution of the kurtosis.

Figure 2 shows the development of the surface elevation at different points in the evolution for $d=0.01$, $r=0.2$ (weak damping, weak forcing) and for $d=0.05$, $r=3$ (strong damping, strong forcing). In agreement with the evolution of the spectral mean (Fig. LABEL:sub@fig:MeanTime), a stronger downshift is visible at $T_{\text{end}}$ in the case for $d=0.05$ than for $d=0.01$. The wind input is stopped at a fixed norm, therefore, the surface elevation at $T_{\text{off}}$ has a similar average amplitude for both cases.

Figure 3 shows the spectrum at different points in the evolution for the simulation with parameters $d=0.05$, $r=3$. At $T_{\text{off}}$ the wind has caused a broadening, upshift and increase of spectral energy, and the spectrum tends more towards the JONSWAP spectrum of the ocean, see the discussion on this point in Section V. In the last section of the simulation the spectrum is damped and shifted downward further.

In summary, during the wind episode, the norm increases, the spectrum broadens, the spectral mean shifts up and BFI and kurtosis increase. Conversely, after the wind episode, when only dissipation is present, the opposite happens; where the norm decreases, the bandwidth decreases, the spectral mean shifts down and kurtosis and BFI decrease. Naturally, properties of the spectrum such as norm, spectral mean and width are correlated, as they are influenced by $r$ and $d$ in the same direction. Since wind and dissipation have opposite effects on the wave amplitude, when they are balanced their effects cancel out and nonlinear interactions become dominant [25]. In Section V we will discuss the role of each term in Eq. (LABEL:eqn_ModelAdim) in more detail.

IV.2 Statistical quantities

One defining criterion for rogue waves is that the the wave height $H$ exceeds the significant wave height $H_{s}=4\sigma$ by a factor 2 ($H/H_{s}>2$). Since $|a|\sim 2H$, this corresponds to $|a|/H_{s}>1$. In general, the probability of finding a wave that exceeds $x$ times $H_{s}$, the exceedance probability $P(|a|>H_{s})$, is a measure of how dangerous the sea state is. The Rayleigh probability distribution corresponds to the envelope height distribution of a Gaussian, linear, process. The Tayfun ([26]) and Fedele-Tayfun ([27]) distributions, take into account second- and third-order nonlinearities respectively, that in the limit of deep water and narrow-banded waves, depend on the characteristic wave steepness [3].

As the initial condition is a Gaussian spectrum, the envelope follows the Rayleigh distribution at $T=0$. Due to the nonlinear processes, the wave-height distribution has moved away from the Rayleigh distribution at $T_{\text{on}}$ (before the wind input), and is more closely described by the Tayfun distribution (Figure LABEL:sub@fig:exceedBefore), due to the nonlinearities in our system. When measured at the end of the wind episode ($T_{\text{off}}$), the tails are much larger (Figure LABEL:sub@fig:exceedAfter), determined by $r$. As the wave-steepness has increased, the Tayfun distribution continues to give a good description of the tail. The Tayfun-Fedele distribution over-estimates the exceedance probability, as we do not include bound modes in our statistical analysis. As the kurtosis is maximal right after the wind input, the distribution of the wave height will be more peaked compared to before the wind input. The tails get a disproportionately higher weight as compared to the mean of the distribution, causing an inflection point in the exceedance probability plot.

We now turn our attention to the indicator for rogue waves, the kurtosis, and its possible predictors: the BFI and the bandwidth. Because of the presence of viscosity, the simulations cannot be compared at the same time after the wind episode as in [11], since different values of viscosity give rise to a different long-term behavior, as shown in Figure 1a. The results (kurtosis, BFI, spectral mean, steepness and band width) should instead be compared for the same value of wave action, i.e. the energy content within the wave field. Recall that we only consider free waves and $|k|<2k_{0}$. Times $T_{\text{on}}$ and $T_{\text{off}}$ (red lines) and the range where the norm is sampled after the wind input (grey area) are indicated for the case $d=0.05$, $r=0.2$ (solid blue line) in Figure LABEL:sub@fig:NormTime.

IV.3 Kurtosis versus BFI

Starting from the analysis of the NLS, a quadratic relation between the kurtosis and the BFI was derived in [12, 28]:

[TABLE]

In our results, we observe this quadratic relation between BFI and kurtosis after the wind episode (Figure 5). However, the large spread of data points in the figure shows the coefficients of Eq. (IV.3), instead of being constant as predicted in Eq. (IV.3), depend on viscosity, forcing strength and duration. In addition, for the truncated spectrum $\mid k\mid/k_{0}\leq 2$, the sign of the quadratic coefficient is inverted for lower damping values $d=0.01$ (circles), giving an opposite curvature to that of $d=0.05$ (diamonds). See Figure 12(a) in Appendix C for results for the full spectrum. Note that the bandwidth $\sigma_{k}$, and consequently the BFI, is calculated with respect to the spectral mean. If the BFI is calculated with respect to $k_{0}$, the quadratic relation is not recovered.

IV.4 Kurtosis versus bandwidth

A quadratic relationship between the bandwidth and the kurtosis is derived in [15] for a conservative system: the NLS propagating in space, without forcing/damping:

[TABLE]

By calculating the expected value of the Hamiltonian $\langle H\rangle$ and the norm $\langle N\rangle$, the kurtosis is found to depend quadratically on the bandwidth $\Omega$:

[TABLE]

Since our system propagates in time, and since we compare the quantities with respect to the energy content within the wave field, i.e. the norm $N$, we can write the relation in the general form

[TABLE]

This relationship is quite closely followed by all our data-sets (Figure 6), irrespective of the values of $r$ and $d$ (see list with parameters in Appendix B). A least-squares fit yields an $R^{2}=0.92$, for coefficients $c_{0}~{}=~{}-1.17~{}\pm~{}0.07$, $c_{1}~{}=~{}26.95~{}\pm~{}0.35$.

Although our system does not reach a steady-state statistical distribution, because it is inherently dissipative, it can be considered to evolve over a sequence of quasi-stable states sharing the same properties as those obtained in the literature by a kinetic approach. This hints at the quasi-homogeneity of the distribution, which is the main hypothesis used in [15] to derive the quadratic dependence of the kurtosis on the bandwidth. The quadratic dependence occurs in our parametric range, even for norm values sampled just after the wind episode.

Appendix C shows that for the non-truncated spectrum, for the spectral mean too, a quadratic relationship with the kurtosis exists. However, when the spectrum is truncated to a much narrower region ($[-k_{0},k_{0}]$in this case), there is a larger spread of data-points around this quadratic relation. The spectral mean is highly correlated with the bandwidth. As both increase (decrease) during forcing (damping). However, for the truncated spectrum, only the bandwidth remains as a reliable predictor.

In other words, in our mathematical model, the kurtosis is influenced by truncating the spectrum at the same rate as the bandwidth is. Since the spectral mean is mostly influenced by the balance of the modes close to $k_{0}$, it is less influenced by truncating the spectrum than the bandwidth is. Therefore, while there is a clear quadratic relation between the spectral mean and the kurtosis in the full model, this relation is lost when the spectrum is truncated (see Appendix C, Figure 12).

IV.5 Kurtosis prediction

In the previous section we have established that the bandwidth is strongly correlated to the kurtosis in a given one-dimensional sea-state. In addition, Figure 1b shows that the kurtosis exponentially decays as a function of time after the wind episode, with the decay rate depending on the damping coefficient $d$. Combining these two findings, the single-spectrum measurement of the bandwidth can be used to predict the kurtosis at a future time, given a certain damping coefficient $d$. As in the previous section, we only analyze data after the wind episode, i.e. the swell. Therefore, we assume that the kurtosis of the swell at any time $T_{s}=T-T_{\text{off}}$ after the start of the swell $T_{\text{s,0}}=T_{\text{off}}$, can be expressed as:

[TABLE]

The fitting parameters $a$ and $b$ are extracted from all investigated 18 combinations of (r,d), where we perform this fit on the times between $T_{\text{off}}$ and $T_{\text{end}}$. Figure 7 indeed shows that parameter $a$ is a linear function of the initial bandwidth (at $T_{\text{off}}$), and that the decay rate of the kurtosis ($b$) linearly depends on the damping coefficient $d$. In our study we have varied the value of $d$ for illustrational purposes, but in reality this is a measurable property [29, 30]. Thus we conclude that, assuming an evolution in time of the form of Eq. (IV.7), the kurtosis can be predicted by the initial bandwidth.

The relationship between the kurtosis and the initial bandwidth stems from the fact that the kurtosis scales linearly with the norm, as can be seen from the exponential decays of both in Figure 1. The norm experiences an exponential decay with decay rate $d$. Therefore, predicting the kurtosis at a future time $T$, for a given damping coefficient $d$, corresponds to a certain decrease in the norm and can also be seen as moving down the curve in Figure 6.

V Contribution of terms

To gain insight into the contribution of the different terms in the model, we make a comparison between the following variations: $NLS$, $NLS+Dysthe$, $NLS+HOT$ and the full model $NLS+Dysthe+HOT$, see labels in Eq. (LABEL:eqn_ModelAdim). Here $HOT$ denotes the higher order terms in $\varepsilon$, and $Dysthe$ the terms preceded by $i\varepsilon$. In all cases dissipation $d$ and forcing $r$ at leading order are included, with values $r=1,d=0.1$. The higher order dispersion terms $\left[\frac{5}{8}\frac{\partial^{4}A}{\partial X^{4}}+4id\,\frac{\partial^{2}A}{\partial X^{2}}\right]$ are included in the higher order terms, but do not influence results and only provide numerical stability. Comparison for the spectral mean, spectral peak, and kurtosis as a function of time are displayed in Figure 8.

The $NLS$ is symmetric and therefore the spectral mean and spectral peak stay equal to $k_{0}$. In the case of $NLS+Dysthe$, we observe upshift of the spectral mean, downshift of the spectral peak, as also predicted by a simple three-wave system  [31, 32]. Interestingly, the downshift of the peak is permanent, recurrence is not observed. Comparing this behavior to that of $NLS+HOT$ and the full model ($NLS+Dysthe+HOT$), makes it clear that the higher order terms are responsible for the permanent downward trend of the spectral peak and mean.

For the kurtosis, in the case where these effects are not included, we verify that results are in agreement with those obtained in [11] without damping ($d=0$). The maximal kurtosis is strongly influenced by which terms are included. It is increased by the $HOT$, and decreased by the $Dysthe$ terms. The addition of these effects gives a maximal kurtosis for the full model, equal to that of the damped/forced NLS. After the wind episode, the NLS has a linear decrease of kurtosis, while the full model relaxes in a strong exponential way.

The influence of the different terms on the spectrum is compared in Figure 9, at the end of the wind input $T=T_{\text{off}}$. As we let the complete model ($NLS+HOT+Dysthe$) act, the Gaussian initial spectrum develops into a JONSWAP spectrum. The Dysthe terms affect the distribution for $k/k_{0}<1$ (Figure 9b), while $HOT$ terms are important to reproduce the JONSWAP spectrum at $k/k_{0}>1$ (Figure 9c).

It is interesting to see how the various terms in the model affect the quadratic relation proposed by Janssen  [12]. Both the Dysthe and NLS equations without forcing ($d=r=0$) give values of kurtosis in agreement with Eq. (IV.3), indicated by the dark yellow and red squares in Figure 10. The value for the kurtosis and BFI were obtained in the steady state, for the given steepness value of $\varepsilon=0.08$. In non-conservative conditions, ($d>0,r>0$), a quadratic relation cannot always be found, or is inverted, for the damped/forced NLS simulations (dark blue and green diamonds). Instead, the nonlinear Dysthe terms need to be included (red diamonds). This explains why NLS models without higher order terms [11] do not observe the quadratic relation, while it is observed in experimental settings [33, 28]. Interestingly, similar behavior in optics and water waves is found for the kurtosis, in terms of its behavior as a function of propagation [34]. However, as the results included only two only different values of the BFI, a conclusion on the quadratic behavior cannot be made.

VI Conclusion

We investigate the predicted quadratic relation between the BFI and kurtosis of the wave height distribution [12]. As shown in [11], a forced NLS model could not reproduce this prediction. While we do find a quadratic relation between the kurtosis and the BFI with our higher order model, it has to be parameterized depending on the details of the damping and forcing (rate and duration). Instead, we show that the bandwidth is a good candidate to predict the kurtosis, as it provides a general relation for the whole range of damping and forcing investigated in our work.

For a conservative system, the steepness remains roughly constant. As the BFI is the ratio of the bandwidth and the steepness, it is roughly proportional to the bandwidth. Hence, it is no surprise that for the conservative NLS a quadratic relation was found between both BFI and bandwidth with respect to kurtosis.

The bandwidth is calculated as the variance of the wavenumber distribution (Eq. (III.3)). To our knowledge there is no fundamental relation between the kurtosis of the distribution of the wave-height $a$, and the variance of the distribution of the wave number $k$. Therefore this relation must come from the model equations, as is derived in [15] for the conservative NLS. Our observation is that the bandwidth is a robust predictor of the kurtosis behavior also for our non-conservative higher order NLS model, for example after a wind episode.

In addition, we demonstrate that the evolution of the kurtosis is strongly influenced by the damping and forcing rates (Figure 1b). Therefore, once the kurtosis is estimated at the end of the wind episode, i.e. the start of the swell, its subsequent evolution can by predicted by the damping coefficient. In this way, the evolution of the kurtosis can be predicted based on the single spectrum measurement of the bandwidth and on the damping coefficient.

Acknowledgements.

The authors gratefully acknowledge financial support of the Swiss National Science Foundation (Projects Nos. 200021- 155970 and 200020-175697).

Appendix A Model

By denoting the free-surface elevation $\eta(x,t)$ and the velocity potential $\phi(x,z,t)$, where $z$ is the depth coordinate and $x$ the longitudinal propagation direction, the dispersive part of our model equation can be obtained from the following linearized Euler system with viscosity and wind forcing [35]

[TABLE]

where $\Gamma$ is the Miles growth rate due to wind forcing and $\nu$ is the kinematic viscosity, while $\omega_{0}$ and $k_{0}$ are radial frequency and wavenumber of the carrier wave, respectively. These are the Laplace equation within the fluid column, the rigid condition at the bottom, the kinematic and dynamic boundary conditions at the free surface. Positive values for $k$ are selected by imposing the boundary condition at the bottom to the solutions of the Laplace equation. Expanding the above system according to the normal modes

[TABLE]

the eigenvalue problem reduces to

[TABLE]

whose only non-trivial solutions correspond to the case where the determinant of the matrix is zero, leading to the following dispersion relation:

[TABLE]

or, equivalently:

[TABLE]

We then apply the method suggested in Ref. [36] that allows to find dispersive terms at all orders in steepness for an evolution equation written in Fourier space as:

[TABLE]

where $\hat{a}(\ell,t)$ is the Fourier transform of the envelope field $a(x,t)$ and $\epsilon=a_{0}k_{0}$ is the steepness, $a_{0}$ being the initial wave amplitude. Thus, by Taylor expanding $\omega(k_{0}+\epsilon\ell)$ about $\epsilon=0$, and setting $\nu=\epsilon^{2}\nu^{\prime}$, $\Gamma=\epsilon^{2}\Gamma^{\prime}$, one obtains:

[TABLE]

Moving to the real space by using $\frac{\partial}{\partial x}=-i\epsilon l$, omitting terms beyond fifth-order in steepness, and multiplying by $-i$, the following evolution equation can be obtained:

[TABLE]

Note that while the viscosity series is finite, i.e. bound to fourth order in $\epsilon$, the wind series is not, since the wind growth-rate parameter $\Gamma$ occurs under the square root in Eq. (A.5). Since a natural cutoff in the spectrum appears in physical fluids for high wave-numbers, meaning that viscosity dominates at very small scales, we neglect the wind terms at $O(\epsilon^{4})$ in order to mimic such natural behavior. We will see in section II.3 that this choice allows us to reproduce the empirical formula obtained by Plant in Ref. [19] with a good agreement.

Including the nonlinear Dysthe terms [13] in the previous evolution equation (A.8), we finally obtain the forced/damped-modified NLS equation:

[TABLE]

Appendix B Simulation parameters

Appendix C Truncating of the spectrum

Figure 11 shows the power spectrum for the surface elevation and the envelope. After the wind input, $T_{\text{off}}$ (when the energy is maximal) the spectrum is quite broad. On a numerical level, the fact that the spectrum for the surface elevation $S_{\eta}$ does not tend to zero can lead to artifacts in calculating quantities that rely on the balance between the upper and lower side of the spectrum. It can introduce spurious asymmetries, such as in the calculation for the bandwidth and the spectral mean.

As mentioned in [11], the calculation for the BFI is quite different if it is based on the full spectrum (figure 8 of [11]) versus the spectrum cut at $k=2k_{0}$ for the surface elevation, corresponding to the interval $[-k_{0},k_{0}]$ for the envelope (as used in the analysis of [11]).To keep our calculations consistent, we chose to bound not only the calculation of the BFI, but of all calculated quantities based on the spectrum. Figure 12 demonstrates the effect of different truncation bandwidths for the spectrum.

In addition, as our model is of higher order in steepness, we chose an upper bound of $2k_{0}$ instead of $1k_{0}$. In order to retain the symmetric behavior of the envelope, and allowing a more direct comparison with the forced NLS, which is a symmetric equation, we limit the interval for the envelope to $[-2k_{0},2k_{0}]$. The main results are repeated in Figure 12b for comparison. While the modes $k<-k_{0}$ are technically nonphysical, Figure 11 shows these modes are several orders of magnitude smaller than the modes $k>k_{0}$. In addition, as shown in Eq II.10, our approximation gives a negative growth rate for $k=-k_{0}$, such that this mode will not grow unbounded. Therefore, as expected, an asymmetric interval cutting the modes $k<-k_{0}$, that is $[-k_{0},2k_{0}]$ (Figure 12(c)) yields very similar results to those presented in Section IV of this paper, but loses the aforementioned symmetry.