Rogue Waves and Large Deviations in Deep Sea

TL;DR

This paper develops a probabilistic and computational framework combining large deviations theory and Monte Carlo methods to predict rogue waves in deep sea, identifying early warning signs and specific wave patterns.

Contribution

It introduces an efficient optimization-based approach to estimate rogue wave probabilities and precursors using the modified nonlinear Schrödinger equation with random initial conditions.

Findings

Rogue waves are linked to unlikely wave configurations triggering large disturbances.

Specific wave shape patterns serve as early precursors for rogue waves.

The method effectively estimates tail probabilities of extreme sea surface elevations.

Abstract

The appearance of rogue waves in deep sea is investigated using the modified nonlinear Schr\"odinger (MNLS) equation in one spatial-dimension with random initial conditions that are assumed to be normally distributed, with a spectrum approximating realistic conditions of a uni-directional sea state. It is shown that one can use the incomplete information contained in this spectrum as prior and supplement this information with the MNLS dynamics to reliably estimate the probability distribution of the sea surface elevation far in the tail at later times. Our results indicate that rogue waves occur when the system hits unlikely pockets of wave configurations that trigger large disturbances of the surface height. The rogue wave precursors in these pockets are wave patterns of regular height but with a very specific shape that is identified explicitly, thereby allowing for early detection. The method proposed here combines Monte Carlo sampling with tools from large deviations theory that reduce the calculation of the most likely rogue wave precursors to an optimization problem that can be solved efficiently. This approach is transferable to other problems in which the system's governing equations contain random initial conditions and/or parameters.

Figures

Equations

$${\partial }_{t}u+\frac{1}{2}{\partial }_{x}u+\frac{i}{8}{\partial }_x^{2}u-\frac{1}{16}{\partial }_x^{3}u+\frac{i}{2}|u{|}^{2}u$$

$$+\frac{3}{2}|u{|}^2{\partial }_{x}u+\frac{1}{4}{u}^2{\partial }_x\bar{u}-\frac{i}{2}\left|{\partial }_x\right||u{|}^2=0,x\in [0,L],$$

$$u_0(x)=\sum _{n\in \mathbb{Z}}{e}^{i{k}_{n}x}(2{\hat{C}}_n{)}^{1/2}{\theta }_n,{\hat{C}}_n=A{e}^{-k_n^2/(2{\Delta }^2)},$$

$$P_{T}(z)\equiv\mathbb{P}\,\big{(}F(u(T))\geq z\big{)},$$

$$F(u(T))=\underset{x\in [0,L]}{\max }|u(T,x)|.$$

$$\Vert u_0{\Vert }_C^2=\sum _{n\in \mathbb{Z}}\frac{|{\hat{a}}_n{|}^2}{{\hat{C}}_n},{\hat{a}}_n=\frac{1}{L}{\int }_0^L{e}^{-i{k}_{n}x}{u}_0(x)dx.$$

$$\frac{1}{2}\underset{u_0\in \Omega (z)}{\min }\Vert u_0{\Vert }_C^2\equiv I_T(z),$$

$$P_T(z)\asymp \exp \left(-I_T(z)\right).$$

$$\underset{u_0}{\min }\left(\frac{1}{2}\Vert u_0{\Vert }_C^2-\lambda F(u(T,u_0))\right)\equiv S_T(\lambda ),$$

$$I_T(z(\lambda ))=\frac{1}{2}\Vert u_0^{\star }(\lambda ){\Vert }_C^2,z(\lambda )=F(u(T,u_0^{\star }(\lambda ))).$$

$$S_T(\lambda )=\underset{z\in \mathbb{R}}{\sup }(\lambda z-I_T(z))=\underset{z\in \mathbb{R}}{\sup }(\lambda z-\frac{1}{2}\underset{u_0\in \Omega (z)}{\inf }\Vert u_0{\Vert }_C^2),$$

$$P^{N}_{T}(z)=\mathbb{P}\Big{(}\max_{x\in[0,NL]}|u(T,x)|\geq z\Big{)},$$

$$P_T^N(z)=1-{\left(1-P_T(z)\right)}^N.$$

$$p\equiv\mathbb{P}\Big{(}\max_{t\in[0,T_{I}]}\max_{x\in[0,L]}|u(t,x)|\geq z\Big{)}\sim 1-(1-P(z))^{T_{I}/\tau_{c}}\,,$$

$$P(z)=\mathbb{P}(\varphi (\theta )>z),$$

$$P(z)=(2\pi {)}^{-D/2}{\int }_{\Omega (z)}{e}^{-\frac{1}{2}|\theta {|}^2}d\theta ,$$

$${\theta }^{\star }(z)=\underset{\theta \in \partial \Omega (z)}{\mathrm{argmin}}|\theta {|}^2,$$

$$\frac{1}{2}|{\theta }^{\star }(z){|}^2$$

$$\underset{z\to \infty }{\lim }\frac{1}{2}|{\theta }^{\star }(z){|}^2=\infty .$$

$${\theta }^{\star }(z)=\lambda \nabla \varphi ({\theta }^{\star }(z))$$

$$\frac{{\theta }^{\star }(z)}{|{\theta }^{\star }(z)|}=\frac{\nabla \varphi ({\theta }^{\star }(z))}{|\nabla \varphi ({\theta }^{\star }(z))|}=\hat{n}(z).$$

$$|\theta {|}^2$$

$$=|{\theta }^{\star }(z){|}^2+2|{\theta }^{\star }(z)|z+z^2,$$

$$(2\pi {)}^{-D/2}{\int }_{P(z)}{e}^{-\frac{1}{2}|\theta {|}^2}d\theta$$

$$=(2\pi {)}^{-D/2}{\int }_0^{\infty }{e}^{-\frac{1}{2}|{\theta }^{\star }(z){|}^2-|{\theta }^{\star }(z)|u-\frac{1}{2}{u}^2}du{\int }_{{\mathbb{R}}^{N-1}}{e}^{-\frac{1}{2}|\eta {|}^2}d\eta$$

$$=(2\pi {)}^{-1/2}{e}^{-\frac{1}{2}|{\theta }^{\star }(z){|}^2}{\int }_0^{\infty }{e}^{-|{\theta }^{\star }(z)|u-\frac{1}{2}{u}^2}du$$

$$=(2\pi {)}^{-1/2}|{\theta }^{\star }(z){|}^{-1}{e}^{-\frac{1}{2}|{\theta }^{\star }(z){|}^2}{\int }_0^{\infty }{e}^{-v-\frac{1}{2}|{\theta }^{\star }(z){|}^{-2}{v}^2}dv$$

$$\sim (2\pi {)}^{-1/2}|{\theta }^{\star }(z){|}^{-1}{e}^{-\frac{1}{2}|{\theta }^{\star }(z){|}^2}\text{as z→∞.}$$

$$\underset{z\to \infty }{\lim }|{\theta }^{\star }(z){|}^{-2}\log \left((2\pi {)}^{-D/2}{\int }_{P(z)}{e}^{-\frac{1}{2}|\theta {|}^2}d\theta \right)=-\frac{1}{2}.$$

$$\begin{gathered} {\int }_{P(z)}{e}^{-\frac{1}{2}|\theta {|}^2}d\theta \\ \asymp e^{-\frac{1}{2}|{\theta }^{\star }(z){|}^2}\text{as z→∞.} \end{gathered}$$

$$E(u_0,\lambda )\equiv \frac{1}{2}\Vert u_0{\Vert }_C^2-\lambda F(u(T,u_0)).$$

$$\frac{\delta E}{\delta u_0}=C^{-1}{u}_0-\lambda J^T(T,u_0)\frac{\delta F}{\delta u}$$

$${\partial }_{t}u=b(u),u(t=0)=u_0,$$

$${\partial }_{t}J=\frac{\delta b}{\delta u}J,J(t=0)=\text{Id}.$$

Full text

Rogue waves, long considered a figment of sailor’s imagination, are now recognized to be a real, and serious, threat for boats and naval structures muller2005rogue ; white1998chance . Oceanographers define them as deep water waves whose crest-to-trough height $H$ exceeds twice the significant wave height $H_{s}$, which itself is four times the standard deviation of the ocean surface elevation. Rogue waves appear suddenly and unpredictably, and can lead to water walls with vertical size on the order of $20$–$30$ m haver2004possible ; nikolkina2011rogue , with enormous destructive power. Although rare, they tend to occur more frequently than predicted by linear Gaussian theory onorato-residori-bortolozzo-etal:2013 ; nazarenko:2016 . While the mechanisms underlying their appearance remain under debate akhmediev2009extreme ; akhmediev2010editorial ; onorato2016origin , one plausible scenario has emerged over the years: it involves the phenomenon of modulational instability benjamin-feir:1967 ; zakharov:1968 , a nonlinear amplification mechanism by which many weakly interacting waves of regular size can create a much larger one. Such an instability arises in the context of the focusing nonlinear Schrödinger (NLS) equation zakharov:1968 ; kuznetsov:1977 ; peregrine:1983 ; akhmediev1987exact ; osborne2000nonlinear ; zakharov-ostrovsky:2009 ; onorato:2009 or its higher order variants dysthe:1979 ; stiassnie:1984 ; trulsen1996modified ; craig2010hamiltonian ; gramstad2011hamiltonian , which are known to be good models for the evolution of a unidirectional, narrow-banded surface wave field in a deep sea. Support for the description of rogue waves through such envelope equations recently came from experiments in water tanks onorato2004observation ; chabchoub2011rogue ; chabchoub2012super ; goullet2011numerical , where Dysthe’s MNLS equation in one spatial dimension dysthe:1979 ; stiassnie:1984 was shown to accurately describe the mechanism creating coherent structures which soak up energy from its surroundings. While these experiments and other theoretical works lo1985numerical ; cousins-sapsis:2015A give grounds for the use of MNLS to describe rogue waves, they have not addressed the question of their likelihood of appearance. Some progress in this direction has been recently made in cousins2016reduced , where a reduced model based on MNLS was used to estimate the probability of a given amplitude within a certain time, and thereby compute the tail of the surface height distribution. These calculations were done using an ansatz for the solutions of MNLS, effectively making the problem two-dimensional. The purpose of this paper is to remove this approximation, and study the problem in its full generality. Specifically, we consider the MNLS with random initial data drawn from a Gaussian distribution nazarenko2011wave . The spectrum of this field is chosen to have a width comparable to that of the JONSWAP spectrum hasselmann:1973 ; onorato:2001 obtained from observations in the North Sea. We calculate the probability of occurrence of a large amplitude solution of MNLS out of these random initial data, and thereby also estimate the tail of the surface height distribution. These calculations are performed within the framework of large deviations theory (LDT), which predicts the most likely way by which large disturbances arise and therefore also explains the mechanism of rogue wave creation. Our results are validated by comparison with brute-force Monte-Carlo simulations, which indicate that rogue waves in MNLS are indeed within the realm of LDT. Our approach therefore gives an efficient way to assess the probability of large waves and their mechanism of creation.

I Problem setup

Our starting point will be the MNLS equation for the evolution of the complex envelope $u(t,x)$ of the sea surface in deep water dysthe:1979 , in terms of which the surface elevation reads \eta(t,x)=\Re\big{(}u(t,x)e^{i(k_{0}x-\omega_{0}t)}\big{)} (here $k_{0}$ denotes the carrier wave number, $\omega_{0}=\sqrt{gk_{0}}$, and $g$ is the gravitational acceleration). Measuring $u$ and $x$ in units of $k^{-1}_{0}$ and $t$ in $\omega_{0}^{-1}$ we can write MNLS in non-dimensional form as

[TABLE]

where the bar denotes complex conjugation. We will consider Eq. 1 with random initial condition $u_{0}(x)\equiv u(0,x)$, constructed via their Fourier representation,

[TABLE]

where $k_{n}=2\pi n/L$, $\theta_{n}$ are complex Gaussian variables with mean zero and covariance $\mathbb{E}\theta_{n}\bar{\theta}_{m}=\delta_{m,n}$, $\mathbb{E}\theta_{n}\theta_{m}=\mathbb{E}\bar{\theta}_{n}\bar{\theta}_{m}=0$. This guarantees that $u_{0}(x)$ is a Gaussian field with mean zero and $\mathbb{E}(u_{0}(x)\bar{u}_{0}(x^{\prime}))=2\sum_{n\in\mathbb{Z}}e^{ik_{n}(x-x^{\prime})}\hat{C}_{n}$.To make contact with the observational data, the amplitude $A$ and the width $\Delta$ in Eq. 2 are picked so that $\hat{C}_{n}$ has the same height and area as the JONSWAP spectrum hasselmann:1973 ; onorato:2001 – see Supporting Information for details.

Because the initial data for Eq. 1 are random, so is the solution at time $t>0$, and our aim is to compute

[TABLE]

where $\mathbb{P}$ denotes probability over the initial data and $F$ is a scalar functional depending on $u$ at time $T>0$. Even though our method is applicable to more general observables, here we will focus on

[TABLE]

II Large deviations theory approach

A brute force approach to calculate Eq. 3 is Monte-Carlo sampling: Generate random initial conditions $u_{0}(x)$ by picking random $\theta_{n}$’s in Eq. 2, evolve each of these $u_{0}(x)$ deterministically via Eq. 1 up to time $t=T$ to get $u(T,x)$, and count the proportion that fulfill $F(u(T))\geq z$. While this method is simple, and will be used below as benchmark, it looses efficiency when $z$ is large, which is precisely the regime of interest for the tails of the distribution of $F(u(T))$. In that regime, a more efficient approach is to rely on results from LDT which assert that Eq. 3 can be estimated by identifying the most likely initial condition that is consistent with $F(u(T))\geq z$. To see how this result comes about, recall that the probability density of $u_{0}$ is formally proportional to $\exp(-\frac{1}{2}\|u_{0}\|^{2}_{C})$, where $\|u_{0}\|^{2}_{C}$ is given by

[TABLE]

To calculate Eq. 3 we should integrate this density over the set $\Omega(z)=\{u_{0}:F(u(T,u_{0}))\geq z\}$, which is hard to do in practice. Instead we can estimate the integral by Laplace’s method. As shown in Material and Methods, this is justified for large $z$, when the probability of the set $\Omega(z)$ is dominated by a single $u_{0}(x)$ that contributes most to the integral and can be identified via the constrained minimization problem

[TABLE]

which then yields the following LDT estimate for Eq. 3

[TABLE]

Here $\asymp$ means that the ratio of the logarithms of both sides tends to 1 as $z\to\infty$. As discussed in Material and Methods, a multiplication prefactor can be added to (7) but it does not affect significantly the tail of $P_{T}(z)$ on a logarithmic scale.

In practice, the constraint $F(u(T,u_{0}))\geq z$ can be imposed by adding a Lagrange multiplier term to Eq. 6, and it is easier to use this multiplier as control parameter and simply see a posteriori what value of $z$ it implies. That is to say, perform for various values of $\lambda$ the minimization

[TABLE]

over all $u_{0}$ of the form in Eq. 2 (no constraint), then observe that this implies the parametric representation

[TABLE]

where $u^{\star}_{0}(\lambda)$ denotes the minimizer obtained in Eq. 8. It is easy to see from Eqs. 6 and 8 that $S_{T}(\lambda)$ is the Legendre transform of $I_{T}(z)$ since:

[TABLE]

III Results

We considered two sets of parameters. In Set 1 we took $A=5.4\cdot 10^{-5}k_{0}^{-2}$ and $\Delta=0.19k_{0}$. Converting back into dimensional units using $k^{-1}_{0}=36$ m consistent with the JONSWAP spectrum hasselmann:1973 ; onorato:2001 , this implies a significant wave height $H_{s}=4\sqrt{C(0)}=3.3$ m classified as a rough sea WMO:2016 . It also yields a Benjamin-Feir index BFI$\,=2\sqrt{2C(0)}/\Delta=0.34$, onorato:2001 ; janssen:2003 , meaning that the modulational instability of a typical initial condition is of medium intensity. In Set 2 we took $A=3.4\cdot 10^{-4}k_{0}^{-2}$ and $\Delta=0.19k_{0}$, for which $H_{s}=8.2$ m is that of a high sea and the BFI is 0.85, meaning that the modulational instability of a typical initial condition is stronger.

Fig. 1 (top) shows the time evolution of $|u(t,x)|$ starting from an initial condition from Set 1 optimized so that $\max_{x}|u(T,x)|=8$ m at $T=20$ min. For comparison, Fig. 1 (bottom) shows $|u(t,x)|$ for a typical initial condition drawn from its Gaussian distribution. To illustrate what is special about the initial conditions identified by our optimization procedure, in Fig. 2 we show snapshots of the surface elevation $\eta(t,x)$ at three different times, $t=0,10,20$ min (black lines), using the constraint that $\max_{x}|u(T,x)|\geq 4.8$ m at $T=20$ min. Additionally, we averaged all Monte-Carlo samples achieving $\max_{x}|u(t,x)|\geq 4.8$ m, translated to the origin. Snapshots of this mean configuration are shown in Fig. 2 (blue lines). They agree well with those of the optimized solution (black lines). The one standard deviation spread around the mean Monte-Carlo realization (light blue) is reasonably small, especially around the rogue wave at final time. This indicates that the event $\max_{x}|u(T,x)|\geq 4.8$ m is indeed realized with probability close to 1 by starting from the most likely initial condition consistent with this event, as predicted by LDT. The usefulness of LDT is confirmed in Figs. 3 (a,b) depicting the probabilities of $\max_{x}|u(T,x)|$ for both Sets 1 and 2 calculated via LDT optimization (lines), compared to Monte-Carlo sampling (dots). As can be seen, the agreement is remarkable, especially in the tail corresponding to the rogue wave regime. As expected, the Monte-Carlo sampling becomes inaccurate in the tail, since there the probabilities are dominated by unlikely events. The LDT calculation, in contrast, remains efficient and accurate far in the tail.

The probabilities plotted in Fig. 3 (a,b) show several remarkable features. First, they indicate that, as $T$ gets larger, their tails fatten significantly. For example, in Set 1 $P_{T=20\,\text{min}}(6\,\text{m})\approx 10^{-5}$, which is $5$ orders of magnitude larger than initially, $P_{t=0\,\text{min}}(6\,\text{m})\approx 10^{-10}$. Secondly, the probabilities converge to a limiting density for large $T$. This occurs after some decorrelation time  $\tau_{c}\approx 10$ min in Set 1 and $\tau_{c}\approx 3$ min in Set 2. Similarly, the LDT results converge. In fact, this convergence can be observed at the level of the trajectories generated from the optimal $u_{0}^{\star}$. As Fig. 4 shows, reading these trajectories backward from $t=T$, their end portions coincide, regardless on whether $T=20$ min, $T=15$ min, or $T=10$ min. The implications of these observations, in particular on the mechanism of creation of rogue waves and their probability of appearance within a time window, will be discussed in Interpretation below.

Before doing so, let us discuss the scalability of our results to larger domain sizes, referring the reader to the Supporting Information for more details. As shown above, the optimization procedure based on large deviation theory predicts that the most likely way a rogue wave will occur in the domain is via the apparition of a single large peak in $|u(t,x)|$. In the set-up considered before, this prediction is confirmed by the brute-force simulations using Monte-Carlo sampling. It is clear, however, that for increased domain size, e.g. by taking a domain size of $NL$ with $N\gg 1$, it will become increasingly likely to observe multiple peaks, for the simple reason that large waves can occur independently at multiple sufficiently separated locations. In these large domains, the large deviation predictions remain valid if we look at the maximum of $|u(t,x)|$ in observation windows that are not too large (that is, about the size of the domain $L$ considered above). However, they deteriorate if we consider this maximum in the entire domain of size $NL$, in the sense that the value $\mathbb{P}\left(\max_{x\in[0,NL]}|u(t,x)|\geq z\right)$ predicted by large deviation theory matches that from Monte-Carlo sampling at values of $z$ that are pushed further away in the tails. This is an entropic effect, which is easy to correct for: events in different subwindows must be considered independent, and their probabilities superposed. That is, if we denote by

[TABLE]

it can be related to $P_{T}(z)=\mathbb{P}(\max_{x\in[0,L]}|u(T,x)|\geq z)$ via

[TABLE]

This formula is derived in the Supporting Information and shown to accurately explain the numerical results. For efficiency $L$ is chosen to be the smallest domain size for which boundary effects can be neglected, in the sense that the shape of the optimal trajectories does no longer change if $L$ is increased further. In effect, this provides us with a method to scale up our results to arbitrary large observation windows.

IV Interpretation

The convergence of $P_{T}(z)$ towards a limiting function $P(z)$ has important consequences for the significance and interpretation of our method and its results. Notice first that this convergence can be explained if we assume that the probability distribution of the solutions to Eq. 1 with Gaussian initial data converges to an invariant measure. In this case, for large $T$, the Monte-Carlo simulations will sample the value of $\max_{x}|u|$ on this invariant measure, and the optimization procedure based on LDT will do the same. The timescale $\tau_{c}$ over which convergence occurs depends on how far this invariant measure is from the initial Gaussian measure of $u_{0}(x)$. Interestingly the values we observe for $\tau_{c}$ are in rough agreement with the time scales predicted by the semi-classical limit of NLS that describes high-power pulse propagation bertola2013universality ; tikan2017universality . As recalled in the Supporting Information, this approach predicts that the timescale of apparition of a focusing solution starting from a large initial pulse of maximal amplitude $U_{i}$ and length-scale $L_{i}$ is $\tau_{c}=\sqrt{T_{\text{nl}}T_{\text{lin}}}$, where $T_{\text{nl}}=\left(\tfrac{1}{2}\omega_{0}k_{0}^{2}U_{i}^{2}\right)^{-1}$ is the nonlinear timescale for modulational instability and $T_{\text{lin}}=8\omega_{0}^{-1}k_{0}^{2}L_{i}^{2}$ is the linear timescale associated to group dispersion. Setting $U_{i}=H_{s}$ (the size at the onset of rogue waves) and $L_{i}=\sqrt{2\pi}\Delta^{-1}$ (the correlation length of the initial field) gives $\tau_{c}\simeq 18$ min for Set 1 and $\tau_{c}\simeq 8$ min for Set 2, consistent with the convergence times of $P_{T}(z)$. This observation has implications in terms of the mechanism of apparition of rogue waves, in particular their connection to the so-called Peregrine soliton, that has been invoked as prototype mechanism for rogue waves creation peregrine:1983 ; akhmediev2009waves ; shrira2010makes ; akhmediev2013recent ; onorato-residori-bortolozzo-etal:2013 ; toenger2015emergent , in particular for water waves chabchoub2011rogue ; chabchoub2012super ; chabchoub2016tracking , plasmas bailung2011observation and fiber optics kibler2010peregrine ; suret2016single ; tikan2017universality . This connection is discussed in the Supporting Information.

Our findings also indicate that, even though the assumption that $u_{0}(x)$ is Gaussian is incorrect in the tail (that is, $P_{T=0}(z)$ is not equal to the limiting $P(z)$ in the tail), it contains the right seeds to estimate $P(z)$ via $P_{T}(z)$ if $T\gtrsim\tau_{c}$111This convergence occurs on the timescale $\tau_{c}$ which is much smaller than the mixing time for the solutions of Eq. 1, i.e.  the time it would take from a given initial condition, rather than an ensemble thereof, to sample the invariant measure.. Altogether this is consistent with the scenario put forward by Sapsis and collaborators in mohamad2016probabilistic ; farazmand-sapsis:2017 to explain how extreme events arise in intermittent dynamical systems and calculate their probability: they occur when the system hits small instability pockets which trigger a large transient excursion. In this scenario, as long as the initial probability distribution in these pockets is accurate, the dynamics will permit precise estimation of the distribution tail. In some sense, the distribution of the initial condition plays a role of the prior distribution in Bayesian inference222Note in particular that the Gaussian field in Eq. 2 is the random field that maximizes entropy given the constraint on its covariance $C(x)$., and the posterior can be effectively sampled by adding the additional information from the dynamics over short periods of time during which instabilities can occur. In mohamad2016probabilistic , this picture was made predictive by using a two-dimensional ansatz for the initial condition $u_{0}(x)$ to avoid having to perform sampling in high-dimension over the original $u_{0}(x)$. What our results show is that this approximation can be avoided altogether by using LDT to perform the calculations directly with the full Gaussian initial condition in Eq. 2.

Interestingly, we can use the results above to calculate the probability of occurrence of rogue waves in a given time window. More precisely, the probability $p(z,T_{I})$ that a rogue wave of amplitude larger than $z$ be observed in the domain $[0,L]$ during $[0,T_{I}]$ (i.e. that $\max_{t\in[0,T_{I}]}\max_{x\in[0,L]}|u(t,x)|\geq z$) can be estimated in terms of $P(z)$ and $\tau_{c}$ as

[TABLE]

where we used the fact that rogue waves can be considered independent on timescales larger than $\tau_{c}$ and assumed $T_{I}\gg\tau_{c}$. The function $p$ is plotted in Fig. 5 as a function of $z$ and $T_{I}$. For example for Set 1, Eq. 13 indicates a 50% chance to observe a rogue wave of height $z=4$ m (that is, about 8 m from crest-to-trough) after 11 hours (using $\tau_{c}=10$ min and $P(z=4~{}\text{m})=1.1\cdot 10^{-2}$); if we wait 30 hours, the chance goes up to 85%. Similarly, for Set 2 the chance to observe a wave of 11 m height is about 50% after 3 hours and about 85% after 8 hours ($\tau_{c}=3$ min and $P(z=11~{}\text{m})=1.2\cdot 10^{-2}$).

V Concluding remarks

We have shown how an optimization problem building on LDT can be used to predict the pathway and likelihood of appearance of rogue waves in the solutions of MNLS fed by random initial data consistent with observations. This setup guarantees accuracy of the core of the initial distribution, which in turn permits the precise estimation of its tail via the dynamics. Our results give quantitative estimate for the probabilities of observing high amplitude waves within a given time window. These results also show that rogue waves have very specific precursors, a feature that was already noted in farazmand2017reduced in the context of a reduced model and could potentially be used for their early detection.

VI Materials and Methods

VI.1 Laplace method and large deviations

For the reader’s convenience, here we recall some standard large deviations results that rely on the evaluation of Gaussian integrals by Laplace’s method and are at the core of the method we propose. It is convenient to rephrase the problem abstractly and consider the estimation of

[TABLE]

where $\theta\in\mathbb{R}^{D}$ are Gaussian random variables with mean zero and covariance Id, and $\phi:\mathbb{R}^{D}\to\mathbb{R}$ is some real valued function – as long as we truncate the sum in Eq. 2 to a finite number of modes, $|n|\leq M$, the problem treated in this paper can be cast in this way, with $\theta$ playing the role of $C^{-1/2}u_{0}$ and $\phi(\cdot)$ that of $F(u(T,C^{-1/2}\cdot))$. The probability $P(z)$ in Eq. 14 is given by

[TABLE]

where $\Omega(z)=\{\theta:\phi(\theta)>z\}$. The interesting case is when this set does not contain the origin, $0\not\in\Omega(z)$, which we will assume is true when $z>0$. We also make two additional assumptions:

The point on the boundary $\partial\Omega(z)$ that is closest to the origin is isolated: Denoting this point as

[TABLE]

we assume that

[TABLE] 2. 2.

The connected piece of $\partial\Omega(z)$ that contains $\theta^{\star}(z)$ is smooth with a curvature that is bounded by a constant independent of $z$.

The point $\theta^{\star}(z)$ satisfies the Euler-Lagrange equation for Eq. 16, with the constraint incorporated via a Lagrange multiplier term:

[TABLE]

for some Lagrange multiplier $\lambda$. This implies that

[TABLE]

where $\hat{n}(z)$ denotes the inward pointing unit vector normal to $\partial\Omega(z)$ at $\theta^{\star}(z)$. If we move inside the set $\Omega(z)$ from $\theta^{\star}(z)$ in the direction of $\hat{n}(z)$, the norm $|\theta|^{2}$ increases under the assumptions in Eq. 17. Indeed, setting $\theta=\theta^{\star}(z)+\hat{n}(z)u$ with $u\geq 0$, we have

[TABLE]

where we used Eq. 19. In fact, if we were to perform the integral in that direction, the natural variable of integration would be to rescale $u\to u/|\theta^{\star}(z)|$. In particular, if we were to replace $\Omega(z)$ by the half space $P(z)=\{\theta\,|\,\hat{n}(z)\cdot(\theta-\theta^{*}(z))>0\}$, it would be easy to estimate the integral in Eq. 15 by introducing a local coordinate system around $\theta^{*}(z)$, whose first coordinate is in the direction of $\hat{n}(z)$. Indeed this would give:

[TABLE]

The last approximation goes beyond a large deviations estimate (i.e. it includes the prefactor), and it implies

[TABLE]

This log-asymptotic estimate is often written as

[TABLE]

Interestingly, while the asymptotic estimate in Eq. 21 does not necessarily apply to the original integral in Eq. 15 (that is, the prefactor may take different forms depending on the shape of $\partial\Omega(z)$ near $\theta^{\star}(z)$), the rougher log-asymptotic estimate in Eq. 23 does as long as the the boundary $\partial\Omega(z)$ is smooth, with a curvature that is bounded by a constant independent of $z$. This is because because the contribution (positive or negative) to the integral over the region between $\Omega(z)$ and $P(z)$ is subdominant in that case, in the sense that the log of its amplitude is dominated by $|\theta^{\star}(z)|$. This is the essence of the large deviations result that we apply in this paper.

VI.2 Numerical aspects

To perform the calculations, we solved Eq. 1 with $L=40\pi$ and periodic boundary conditions, and checked that this domain is large enough to make the effect of these boundary conditions negligible (see Supporting Information). The spatial domain was discretized using $2^{12}$ equidistant gridpoints, which is enough to resolve the solution of Eq. 1. To evolve the field $u(t,x)$ in time we used a pseudo-spectral second order exponential time-differencing (ETD2RK) method cox:2002 ; kassam:2005 .

When performing the Monte-Carlo simulations, we used $10^{6}$ realizations of the random initial data constructed by truncating the sum in Eq. 2 over the $M=23$ modes with $-11\leq n\leq 11$, i.e $-3\Delta\leq k_{n}\leq 3\Delta$: these modes carry most of the variance, and we checked that adding more modes to the initial condition did not affect the results in any significant way (see Supporting Information).

VI.3 Optimization procedure

As explained above, the large deviation rate function $I_{T}(z)$ in Eq. 6 can be evaluated by solving the dual optimization problem in Eq. 8, which we rewrite as $S_{T}(\lambda)=\inf_{u_{0}}E(u_{0},\lambda)$, where we defined the cost function

[TABLE]

We performed this minimization using steepest descent with adaptive step (line-search) and preconditioning of the gradient borzi2011computational . This involves evaluating the (functional) gradient of $E_{T}(u_{0},\lambda)$ with respect to $u_{0}$. Using the chain rule, this gradient can be expressed as (using compact vectorial notation)

[TABLE]

where $J(t,u_{0})=\delta u(t,u_{0})/\delta u_{0}$ is the Jacobian of the transformation $u_{0}\to u(t,u_{0})$. Collecting all terms on the right-hand-side of the MNLS Eq. 1 into $b(u)$, this equation can be written as

[TABLE]

and it is easy to see that in this notation $J(t,u_{0})$ satisfies

[TABLE]

Consistent with what was done in the Monte-Carlo sampling, to get the results presented above we truncated the initial data $u_{0}$ over $M=23$ modes using the representation

[TABLE]

This means that minimization of Eq. 24 was performed in the $2M-1=45$ dimensional space spanned by the modes $\hat{a}_{n}$, accounting for invariance by an overall phase shift – to check convergence we also repeated this calculation using larger values of $M$ and found no noticeable difference in the results (see Supporting Information).

In practice, the evaluation of the gradient in Eq. 25 was performed by integrating both $u(t)$ and $J(t)$ up to time $t=T$. Eq. 27 was integrated using the same pseudo-spectral method as for Eq. 1 on the same grid. To perform the steepest descent step, we then preconditioned the gradient through scalar multiplication by the step-independent, diagonal metric with the components of the spectrum $\hat{C}_{n}$ as diagonal elements.

Acknowledgment

We thank W. Craig and M. Onorato for helpful discussions, and O. Bühler, M. Mohamad, and T. Sapsis for interesting comments. We also thank the anonymous reviewer for drawing our attention to the semi-classical theory for the Nonlinear Schrödinger Equation. GD is supported by the joint Math PhD program of Politecnico and Università di Torino. EVE is supported in part by the Materials Research Science and Engineering Center (MRSEC) program of the National Science Foundation (NSF) under award number DMR-1420073 and by NSF under award number DMS-1522767.